Page 141 - Compact Numerical Methods For Computers

130 Compact numerical methods for computers
Algorithm 14. A Jacobi algorithm for eigensolutions of a real symmetric matrix (cont.)

t := sqrt(4.0*p*p+q*q);
if t=0.0 then {STEP 5}
begin {STEP 11 -- If t is zero, no rotation is needed.}
rotn := false; {to indicate no rotation is to be performed.}
end
else
begin {t>0.0}
if q>=0.0 then {STEP 6}
begin {rotation for eigenvalue approximations already in order}
{STEP 7 -- test for small rotation}
oki := (abs(A[i,i])=abs(A[i,i])+l00.0*abs(p));
okj := (abs(A[j,j])=abs(A[j,j])+l00.0*abs(p));
if oki and okj then rotn := false
else rotn := true;
{This test for a small rotation uses an arbitrary factor of
100 for scaling the off-diagonal elements. It is chosen to
ensure “small but not very small” rotations are performed.}
if rotn then
begin {STEP 8}
c := sqrt((t+q)/(2.0*t)); s := p/(t*c);
end;
end {if q>=0.0}
else
begin {q<0.0 -- always rotate to bring eigenvalues into order}
rotn := true; {STEP 9}
s := sqrt((t-q)/(2,0*t));
if p<0.0 then s := -s;
c := p/(t*s);
end; {STEP 10}
if 1.0=(1.0+abs(s)) then rotn := false; {test for small angle}
end; {if t=0.0}
if rotn then {STEP 11 -- rotate if necessary}
begin {STEP 12}
for k := 1 to n do
begin
q := A[i,k]; A[i,k] := c*q+s*A[j,k]; A[j,k] := -s*q+c*A[j,k];
end; {left multiplication of matrix A}
{STEP 13}
for k := l to n do
begin {right multiplication of A and V}
q := A[k,i]; A[k,i] := c*q+s*A[k,j]; A[k,j] := -s*q+c*A[k,j];
{STEP 14 -- can be omitted if eigenvectors not needed}
q := V[k,i]; V[k,i] := c*q+s*V[k,j]; V[k,j] := -s*q+c*V[k,j];
end; {loop on k for right multiplication of matrices A and V}
end {rotation carried out}
else
{STEP 11 -- count skipped rotations}
skipped := skipped+1; {to count the skipped rotations}
end; {loop on j} {STEP 15}
end; {loop on i. This is also the end of the sweep. -- STEP 16}
writeln(‘ ’,skipped,‘ / ’,n*(n-1) div 2,‘ rotations skipped’);
end; {while -- main loop}
end; {alg14.pas = evJacobi -- STEP 17}