begin;
comment eigenvalues of a real symmetric matrix
by the QR method. Algorithm 253, P.A. Businger, CACM 8 (1965) 217
;
real procedure EXTERN(I, M, N, A); code 31234;
integer N;
N := READ;
begin;
integer I, J;
real array A[1 : N, 1 : N];
procedure SYMMETRICQR1(N, G);
value N;
integer N;
array G;
comment uses Housholders' method and the QR algorithm to
find all n eigenvalues of the real symmetric matrix whose lower
triangular part is given in the array g[1:n,1:n]. The computed
eigenvalues are stored as the diagonal elements g[i,i]. The
original contents of the lower triangular part of g are lost during
the computation whereas the strictly upper triagular part of g
is left untouched;
begin;
real procedure SUM(I, M, N, A);
value M, N;
integer I, M, N;
real A;
begin;
real S;
S := 0;
for I := M step 1 until N do S := S + A;
SUM := S;
end SUM;
real procedure MAX(A, B);
value A, B;
real A, B;
MAX := if A > B then A else B;
procedure HOUSHOLDERTRIDIAGONALIZATION1(N, G, A, BQ, NORM);
value N;
integer N;
array G, A, BQ;
real NORM;
comment nonlocal real procedure sum, max;
comment reduces the given real symmetric n by n matrix g
to tridiagonal form using n - 2 elementary orthogonal trans-
formations (I-2ww') = (I-gamma uu'). Only the lower tri
angular part of g need be given. The diagonal elements and
the squares of the subdiagonal elements of the reduced matrix
are stored in a[1:n] and bq[1:n-1] respectively. norm is set
equal to the infinity norm of the reduced matrix. The columns
of the strictly lower triagular part of g are replaced by the
nonzero portions of the vectors u
;
begin;
integer I, J, K;
real T, ABSB, ALPHA, BETA, GAMMA, SIGMA;
array P[2 : N];
NORM := ABSB := 0;
for K := 1 step 1 until N - 2 do
begin;
A[K] := G[K, K];
SIGMA := BQ[K] := SUM(I, K + 1, N, G[I, K] 
2);
T := ABSB + ABS(A[K]);
ABSB := SQRT(SIGMA);
NORM := MAX(NORM, T + ABSB);
if SIGMA 
0 then begin;
ALPHA := G[K + 1, K];
BETA := if ALPHA < 0 then ABSB else -ABSB;
GAMMA := 1 ÷ (SIGMA - ALPHA 
BETA);
G[K + 1, K] := ALPHA - BETA;
for I := K + 1 step 1 until N do
P[I] := GAMMA (SUM(J, K + 1, I, G[I, J] G[J, K]) + SUM(J, I + 1, N, G[J, I] G[J, K]));
T := 0.5 GAMMA SUM(I, K + 1, N, G[I, K] P[I]);
for I := K + 1 step 1 until N do
P[I] := P[I] - T G[I, K];
for I := K + 1 step 1 until N do
for J := K + 1 step 1 until I do
G[I, J] := G[I, J] - G[I, K] P[J] - P[I] G[J, K];
end ;
end K;
A[N - 1] := G[N - 1, N - 1];
BQ[N - 1] := G[N, N - 1] 2;
A[N] := G[N, N];
T := ABS(G[N, N - 1]);
NORM := MAX(NORM, ABSB + ABS(A[N - 1]) + T);
NORM := MAX(NORM, T + ABS(A[N]));
end HOUSHOLDERTRIDIAGONALIZATION1;
integer I, K, M, M1;
real NORM, EPSQ, LAMBDA, MU, SQ1, SQ2, U, PQ, GAMMA, T;
array A[1 : N], BQ[0 : N - 1];
HOUSHOLDERTRIDIAGONALIZATION1(N, G, A, BQ, NORM);
EPSQ := 2.25 mod -22 NORM 2;
comment The tollerance used in the QR iteration depends
on the square of the relative machine precision. Here 2.25%-22
is used which is appropriate for a machine with a 36-bit
mantissa
;
MU := 0;
M := N;
INSPECT: if M = 0 then goto RETURN else I := K := M1 := M - 1;
BQ[0] := 0;
if BQ[K] < EPSQ then begin;
G[M, M] := A[M];
MU := 0;
M := K;
goto INSPECT;
end ;
for I := I - 1 while BQ[I] > EPSQ do K := I;
if K = M1 then begin;
comment treat 2 * 2 block separately
;
MU := A[M1] A[M] - BQ[M1];
SQ1 := A[M1] + A[M];
SQ2 := SQRT((A[M1] - A[M]) 2 + 4 BQ[M1]);
LAMBDA := 0.5 (if SQ1 > 0 then SQ1 + SQ2 else SQ1 - SQ2);
G[M1, M1] := LAMBDA;
G[M, M] := MU ÷ LAMBDA;
MU := 0;
M := M - 2;
goto INSPECT;
end ;
LAMBDA := if ABS(A[M] - MU) < 0.5 ABS(A[M]) then A[M] + 0.5 SQRT(BQ[M1]) else 0.0;
MU := A[M];
SQ1 := SQ2 := U := 0;
for I := K step 1 until M1 do
begin;
comment shortcut single QR iteration
;
GAMMA := A[I] - LAMBDA - U;
PQ := if SQ1 1 then GAMMA 2 ÷ (1 - SQ1) else (1 - SQ2) BQ[I - 1];
T := PQ + BQ[I];
BQ[I - 1] := SQ1 T;
SQ2 := SQ1;
SQ1 := BQ[I] ÷ T;
U := SQ1 (GAMMA + A[I + 1] - LAMBDA);
A[I] := GAMMA + U + LAMBDA;
end I;
GAMMA := A[M] - LAMBDA - U;
BQ[M1] := SQ1 (if SQ1 1 then GAMMA 2 ÷ (1 - SQ1) else (1 - SQ2) BQ[M1]);
A[M] := GAMMA + LAMBDA;
goto INSPECT;
RETURN: ;
end SYMMETRICQR1;
for I := 1 step 1 until N do
for J := 1 step 1 until I do
A[I, J] := READ;
SYMMETRICQR1(N, A);
for I := 1 step 1 until N do
begin;
NLCR;
ABSFIXT(2, 0, I);
PRINT(A[I, I]);
end ;
end ;
end;