Page 71 - Compact Numerical Methods For Computers

60 Compact numerical methods for computers
Algorithm 4. Givens’ reductions, singular-value decomposition and least-squares sol-
ution (cont.)

{Now come important convergence test considerations.
First we will decide if rotation will exchange order of rows.}
{STEP 16 If q<r then goto step 19}
{Check if the rows are ordered.}
if q>= r then
begin {STEP 17 Rows are ordered, so try convergence test.}
if not ((q<=e2*svs[1]) or (abs(p)c=tol*q)) then
{First condition checks for very small row norms in BOTH rows,
for which no rotation makes sense. The second condition
determines if the inner product is small with respect to the
larger of the rows, which implies a very small rotation angle.}
begin {revised STEP 18}
{columns are in order, but not converged to smallinner product.
Calculate angle and rotate.}
p := p/q; r := 1-r/q; vt := sqrt(4*p*p + r*r);
c := sqrt(0.5*(1+r/vt)); s := p/(vt*c);
rotnsub; {STEP 19 in original algorithm}
count := count+1;
end;
end
else {q<r, columns out of order -- must rotate}
{revised STEP 16. Note: r > q, and cannot be zero since both are
sums of squares for the svd. In the case of a real symmetric
matrix, this assumption must be questioned.}
begin
p := p/r; q := q/r-1; vt := sqrt(4*p*p + q*q);
s := sqrt(0.5*(1-q/vt));
if p<0 then s := -s;
c := p/(vt*s);
rotnsub; {STEP 19 in original algorithm}
count := count+1;
end;
{Both angle calculations have been set up so that large numbers
do not occur in intermediate quantities. This is easy in the svd
case, since quantities q and r cannot be negative.}
{STEP 20 has been removed, since we now put the number of
rotations in count, and do not count down to zero.}
end; {loop on k -- end-loop is STEP 21}
end; (loop on j)
sweep := sweep +1;
writeln( ‘Sweep’,sweep,‘ ’,count,‘rotations performed’);
{Set EstColRank to largest column index for which
svs[column index] > (svs[1]*tol + tol*tol)
Note how Pascal expresses this more precisely.}
while (EstRowRank >= 3) and (svs[EstRowRank] <= svs[1]*tol+tol*tol)
do EstRowRank := EstRowRank-1;
until (sweep>slimit) or (count=0); {STEP 22}
{Singular value decomposition now ready for extraction of information
and formation of least squares solution.}
writeln(‘Singular values and principal components’);
for j := 1 to n do {STEP 23}
begin
s := svs[j];