begin comment eigenvalues of a real symmetric matrix
by the QR method. Algorithm 253, P.A. Businger, CACM 8 (1965) 217;
integer n;
n:= read;
begin integer i,j;
real array a[1:n,1:n];
procedure symmetric QR1(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 Housholder tridiagonalization 1(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 Housholder tridiagonalization 1;
integer i,k,m,m1;
real norm,epsq,lambda,mu,sq1,sq2,u,pq,gamma,t;
array a[1:n],bq[0:n-1];
Housholder tridiagonalization 1(n,g,a,bq,norm);
epsq:= 2.25%-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 symmetric QR 1;
for i:= 1 step 1 until n do
for j:= 1 step 1 until i do
a[i,j]:= read;
symmetric QR 1(n,a);
for i:= 1 step 1 until n do
begin NLCR; ABSFIXT(2,0,i); PRINT(a[i,i]) end
end
end