/***************************************************************************/
/* */
/* Algorithm 254. EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX */
/* BY THE QR METHOD. */
/* */
/* P. A. Businger, University of Texas, Austin, Texas. */
/* */
/***************************************************************************/
/* Translated to PL/I by R. A. Vowels, 11 February 2006. */
/* The results are the same as the sample results given with the algorithm
(to 5 figures for single precision), except that
the eigenvectors for x(3) and x(5) have their signs reversed.
/* Businger gives the following test data:
/* 5
/* 4 6
/* 3 0 7
/* 2 4 6 8
/* 1 3 5 7 9
/* and the following results:
/* eigenvalues:
/* 1 = 22.406875306
/* 2 = 7.5137241547
/* 3 = 4.8489501203
/* 4 = -1.0965951820
/* 5 = 1.3270455995
/* eigenvectors:
/* x(1) = ( 0.24587793851, 0.30239603954, 0.45321452335, 0.57717715229, 0.55638458400) */
/* x(2) = ( 0.55096195546, 0.70944033954, -0.34017913315, -0.083410953290, -0.26543567685) */
/* x(3) = ( 0.54717279573, -0.31256992008, 0.61811207635, -0.11560659356, -0.45549374666) */
/* In his Certification, J. H. Welsch states that in his results the signs
were reversed for x(3) and x(4).
I don't know whether Welsch's article has a typo, and that he
really means x(3) and x(5). However, it is clear that
three versions of the algorithm produce different results.
I will try to find another algorithm to determine who is correct.
In Algorithm 297, Boothroyd does not mention any problems
with this algorithm. His algorithm 297 calls this Algorithm.
I ran Algorithm 85, which uses the Jacobi method, and
obtained the same eigenvalues (though in a different order),
and the same eigenvectors (again in a different order),
and some of these also have the signs reversed (different equations
from either of the above).
*/
/* Uses Householder's method and the QR algorithm to */
/* find all n eigenvalues and eigenvectors of the real symmetric */
/* matrix whose lower triangular part is given in the array g. The */
/* computed eigenvalues are stored as the diagonal elements */
/* g(i, i) and the eigenvectors as the corresponding columns of the */
/* array x. The original contents of the lower triangular part of g */
/* are lost during the computation whereas the strictly upper */
/* trianguar part of g is left untouched. */
symmetric_QR2:
procedure (n, g, x);
declare (g(5,5), x(5,5)) float, n fixed binary;
sum: procedure (i, m, n, a, p, q, b, r, s) returns (float);
declare (i, m, n, p, q, r, s) fixed binary; declare (a(5,5), b(5,5)) float;
declare su float;
su = 0;
do i = m to n;
su = su + a(p,q)*b(r,s);
end;
return (su);
end sum;
sum2: procedure (i, m, n, a, p, q, b, r) returns (float);
declare (i, m, n, p, q, r) fixed binary; declare (a(5,5), b(2:5)) float;
declare s float;
s = 0;
do i = m to n;
s = s + a(p,q)*b(r);
end;
return (s);
end sum2;
/* Reduces the given real symmetric n by n matrix G to tri-diagonal form */
/* using n-2 elementary orthogonal transformations */
/* (i-2ww') = (I - gamma ww'). Only the lower triangular part of G need */
/* be given. The computed diagonal and sub-diagonal elements of the */
/* reduced matrix are stored in a(1:n) and b(1:n-1) respectively. The */
/* transformations on the right are also applied to the n by n matrix X. */
/* The columns of the strict lower triangular part of G are replaced by */
/* the non-zero portion of the vectors U, NORM is set equal to the */
/* infinity norm of the reduced matrix. */
Householder_tridiagonalization2:
procedure (n, g, a, b, x, norm);
declare n fixed binary, (g(5,5), a(5), b(0:5), x(5,5), norm) float;
declare (i, j, k) fixed binary,
(t, sigma, alpha, beta, gamma, absb, p(2:5)) float;
norm = 0; absb = 0;
do k = 1 to n-2;
a(k) = g(k,k);
sigma = sum(i, k+1, n, g, i, k, g, i, k);
t = absb+abs(a(k)); absb = sqrt(sigma);
norm = max(norm, t+absb); alpha = g(k+1,k);
if alpha < 0 then beta = absb; else beta = -absb;
b(k) = beta; /* multiple assignment here */
if sigma ^= 0 then
do;
gamma = 1/(sigma-alpha*beta);
g(k+1, k) = alpha - beta;
do i = k+1 to n;
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));
end;
t = 0.5 * gamma*sum2(i, k+1, n, g, i, k, p, i);
do i = k+1 to n; p(i) = p(i) - t*g(i,k); end;
do i = k+1 to n;
do j = k+1 to i;
g(i,j) = g(i,j)-g(i,k)*p(j)-p(i)*g(j,k);
end;
end;
do i = 2 to n;
p(i) = gamma*sum(j, k+1, n, x, i, j, g, j, k);
end;
do i = 2 to n;
do j = k+1 to n;
x(i,j) = x(i,j) - p(i)*g(j,k);
end;
end;
end;
end;
a(n-1) = g(n-1, n-1); a(n) = g(n,n); b(n-1) = g(n, n-1);
t = abs(b(n-1));
norm = max(norm, absb+abs(a(n-1))+t);
norm = max(norm, t+abs(a(n)));
end Householder_tridiagonalization2;
declare (i, j, k, m, m1, its, its2, II) fixed binary,
(t, norm, eps, sine, cosine, lambda, mu, a0, a1, b0, beta, x0, x1) float;
declare (a(5), b(0:5), c(0:4), cs(4), sn(4)) float;
/* set x equal to the identity matrix. */
do i = 1 to n;
x(i,i) = 1;
do j = i+1 to n;
x(i,j) = 0; x(j,i) = 0;
end;
end;
call Householder_tridiagonalization2 (n, g, a, b, x, norm);
eps = norm * 1.5e-11;
eps = norm * 1.5e-6; /* single precision */
/* The tolerance used in the QR iteration is set equal to the */
/* product of the infinity norm of the reduced matrix and the */
/* relative machine precision (here assumed to be 1.5e-11 which is */
/* appropriate for a machine with a 36-bit mantissa) */
b(0) = 0; mu = 0; m = n;
its = 0;
inspect:
if m = 0 then return;
its = its + 1;
put skip list ('ITERATION', its);
i = m-1; k = m-1; m1 = m-1;
if abs(b(k)) <= eps then
do;
g(m,m) = a(m); mu = 0; m = k; go to inspect;
end;
PUT SKIP;
DO II = 0 TO N; PUT LIST (B(II)); END;
PUT SKIP LIST ('EPS=', EPS);
its2 = 0;
i = i-1; /* k = i; */
do while (abs(b(i)) > eps);
k = i;
PUT SKIP LIST (I, abs(b(i)));
its2 = its2 + 1;
if its2 > 10 then stop;
i = i - 1;
end;
if abs(a(m)-mu) < 0.5*abs(a(m)) | (m1 = k) then
lambda = a(m) + 0.5*b(m1);
else
lambda = 0;
mu = a(m); a(k) = a(k)-lambda; beta = b(k);
do j = k to m1; /* Transformation on the left. */
a0 = a(j); a1 = a(j+1)-lambda; b0 = b(j);
t = sqrt(a0*a0+beta*beta);
cosine = a0/t; cs(j) = a0/t; sine = beta/t; sn(j) = beta/t;
a(j) = cosine*a0 + sine*beta; a(j+1) = -sine*b0+cosine*a1;
b(j) = cosine*b0 + sine*a1; beta = b(j+1);
b(j+1) = cosine*beta; c(j) = sine*beta;
end;
b(k-1) = 0; c(k-1) = 0;
do j = k to m1; /* Transformation on the right */
sine = sn(j); cosine = cs(j);
a0 = a(j); b0 = b(j);
b(j-1) = b(j-1)*cosine + c(j-1)*sine;
a(j) = a0*cosine + b0*sine + lambda;
b(j) = -a0*sine + b0*cosine; a(j+1) = a(j+1)*cosine;
do i = 1 to n;
x0 = x(i,j); x1 = x(i,j+1);
x(i,j) = x0*cosine + x1*sine; x(i,j+1) = -x0*sine + x1*cosine;
end;
end;
a(m) = a(m)+lambda; go to inspect;
return:
end symmetric_QR2;