Page 140 - Compact Numerical Methods For Computers

Real symmetric matrices 129
Algorithm 14. A Jacobi algorithm for eigensolutions of a real symmetric matrix (cont.)
}
{STEP 0 -- via the calling sequence of the procedure, we supply the matrix
and its dimensions to the program.}
var
count, i, j, k, limit, skipped : integer;
c, p, q, s, t : real;
ch : char;
oki, okj, rotn : boolean;
begin
writeln(‘alg14.pas -- eigensolutions of a real symmetric’);
writeln(‘matrix via a Jacobi method’);
if initev then {Do we initialize the eigenvectors to columns of
the identity?}
begin
for i := l to n do
begin
for j := 1 to n do V[i,j] := 0.0;
V[i,i] := 1.0; {to set V to a unit matrix -- rotated to become
the eigenvectors}
end; {loop on i;}
end; {initialize eigenvectors}
count := 0;
limit := 30; {an arbitrary choice following lead of Eberlein}
skipped := 0; {so far no rotations have been skipped. We need to set
skipped here because the while loop below tests this variable.}
{main loop}
while (count<=limit) and (skipped<((n*(n-1)) div 2) ) do
{This is a safety check to avoid indefinite execution of the algorithm.
The figure used for limit here is arbitrary. If the program terminates
by exceeding the sweep limit, the eigenvalues and eigenvectors computed
may still be good to the limitations of the machine in use, though
residuals should be calculated and other tests made. Such tests are
always useful, though they are time- and space-consuming.}
begin
count := count+1; {to count sweeps -- STEP 1}
write(‘sweep’,count,‘ ’);
skipped := 0; {to count rotations skipped during the sweep.}
for i := 1 to (n-1) do {STEP 2}
begin {STEP 3}
for j := (i+1) to n do {STEP 4}
begin
rotn := true; {to indicate that we carry out a rotation unless
calculations show it unnecessary}
p := 0.5*(A[i,j]+A[j,i]); {An average of the off-diagonal elements
is used because the right and left multiplications by the
rotation matrices are non-commutative in most cases due to
differences between floating-point and exact arithmetic.}
q := A[i,i]-A[j,j]; {Note: this may suffer from digit cancellation
when nearly equal eigenvalues exist. This cancellation is not
normally a problem, but may cause the algorithm to perform more
work than necessary when the off-diagonal elements are very
small.}