Page 135 - Compact Numerical Methods For Computers

124 Compact numerical methods for computers
Algorithm 13. Eigensolutions of a real symmetric matrix via the singular-value
decomposition (cont.)

W : wmatrix; {working array}
var Z: rvector); {to return the eigenvalues}
{alg13.pas ==
eigensolutions of a real symmetric matrix via the singular value
decomposition by shifting eigenvalues to form a positive definite
matrix.
This algorithm replaces Algorithm 13 in the first edition of Compact
Numerical Methods.
Copyright 1988 J.C.Nash
}
var
count, i, j, k, limit, skipped : integer;
c, p, q, s, shift, t : real ; {rotation angle quantities}
oki, okj, rotn : boolean;
ch : char;
begin
writeln(‘alg13.pas-- symmetric matrix eigensolutions via svd’);
{Use Gerschgorin disks to approximate the shift. This version
calculates only a positive shift.}
shift:=0.0;
for i:=1 to n do
begin
t:=A[i,i];
for j:=1 to n do
if i<>j then t:=t-abs(A[ij]);
if t<shift then shift:=t; {looking for lowest bound to eigenvalue}
end; {loop over rows}
shift:=-shift; {change sign, since current value < 0 if useful}
if shift0.0 then shift:=0.0;
writeln(‘Adding a shift of ’,shift,‘ to diagonal of matrix.’);
for i:=1 to n do
begin
for j:=1 to n do
begin
W[i,j]:=A[i,j]; {copy matrix to working array}
if i=j then W[i,i]:=A[i,i]+shift; {adding shift in process}
if initev then {initialize eigenvector matrix}
begin
if i=j then W[i+n,i]:=0.0
else
begin
W[i+n,j]:=0.0;
end,
end; {eigenvector initialization}
end; {loop on j}
end; {loop on i}
NashSVD(n, n, W, Z); {call alg01 to do the work}
for i:=1 to n do
begin
Z[i]:=sqrt(Z[i])-shift; {to adjust eigenvalues}
for j:=1 to n do
V[i,j]:=W[n+i,j]; {extract eivenvectors}
end; {loop on i}
end; {alg13.pas == evsvd}