Page 70 - Compact Numerical Methods For Computers

Handling larger problems 59
Algorithm 4. Givens’ reductions, singular-value decomposition and least-squares sol-
ution (cont.)

bb := abs(c); if abs(s)>bb then bb := abs(s);
if bb>0.0 then
begin {can proceed with rotation as at least one non-zero element}
c := c/bb; s := s/bb; p := sqrt(c*c+s*s); {STEP 7}
s := s/p; {sin of angle of rotation}
if abs(s)>=tol then
begin {not a very small angle} {STEP8}
c := c/p; {cosine of angle of rotation}
rotnsub; {to perform the rotation}
end; {if abs(s)>=tol}
end; {if bb>0.0}
end; {main loop on j for Givens’ reduction of one observation} {STEP 9}
{STEP 10 -- accumulate the residual sums of squares}
write(‘Uncorrelated residual(s):’);
for j := 1 to nRHS do
begin
rss[j] := rss[j]+sqr(W[k,n+j]); write(W[k,n+j]:10,‘ ’);
if (7 * (j div 7) = j) and (j < nRHS) then
begin
writeln;
end;
end;
writeln;
{NOTE: use of sqr function which is NOT sqrt.}
end; {if (not enddata)}
end; {while (not enddata)}
{This is the end of the Givens’ reduction part of the program.
The residual sums of squares are now in place. We could find the
least squares solution by back-substitution if the problem is of full
rank. However, to determine the approximate rank, we will continue
with a row-orthogonalisation.}
{STEP 11} {Beginning of svd portion of program.}
m := 1; {Starting column for the rotation subprogram}
slimit := n div 4; if slimit< then slimit := 6; {STEP 12}
{This sets slimit, a limit on the number of sweeps allowed.
A suggested limit is max([n/4], 6).}
sweep := 0; {initialize sweep counter}
e2 := l0.0*n*eps*eps; {a tolerance for very small numbers}
tol := eps*0.1; {a convergence tolerance}
EstRowRank := n; {current estimate of rank};
repeat
count := 0; {to initialize the count of rotations performed}
for j := 1 to (EstRowRank-1) do {STEP 13}
begin {STEP 14}
for k := (j+1) to EstRowRank do
begin {STEP 15}
p := 0.0; q := 0.0; r := 0.0;
for i := 1 to n do
begin
p := p+W[j,i]*W[k,i]; q := q+sqr(W[j,i]); r := r+sqr(W[k,i]);
end; {accumulation loop}
svs[j] := q; svs[k] := r;