Page 72 - Compact Numerical Methods For Computers

Handling larger problems 61
Algorithm 4. Givens’ reductions, singular-value decomposition and least-squares sol-
ution (cont.)
s := sqrt(s); svs[j] := s; {to save the singular value}
writeln(‘Singular value[‘,j,’]= ’,s);
if s>=tol then
begin
for i := 1 to n do W[j,i] := W[j,i]/s;
for i := 1 to n do
begin
if (8 * (i div 8) = i) and (i<n) then
begin
writeln;
end;
end; {for i=1...}
{principal component is a column of V or a row of V-transpose. W
stores V-transpose at the moment.}
writeln;
end; {if s>=tol}
{Principal components are not defined for very small singular values.}
end; {loop on j over the singular values}
{STEP 24 -- start least squares solution}
q := 0.0; {to ensure one pass of least squares solution}
while q>=0.0 do
begin
write(‘Enter a tolerance for zero (<0 to exit)’);
readln(infile,q);
if infname<>‘con’ then writeln(q);
if q>=0.0 then
begin
{For each value of the tolerance for zero entered as q we must
calculate the least squares solution and residual sum of squares
for each right hand side vector. The current elements in columns
n+1 to n+nRHS of the work array W give the contribution of each
principal coordinate to the least squares solution. However, we do
not here compute the actual principal coordinates for reasons
outlined earlier.}
for i := 1 to nRHS do {STEP 25}
begin
trss := rss[i]; {get current sum of squared residuals}
for j := 1 to n do {loop over the singular values -- STEP 26}
begin {STEP 27}
p := 0.0;
for k := l to n do
begin
if svs[k]>q then p := p+W[k,j]*W[k,n+i]/svs[k];
end; {loop over singular values}
B[j,i] := p; {to save current solution -- STEP 28}
writeln(‘Solution component [‘,j,’]= ’,p);
if svs[j]<=q then trss := trss+sqr(W[j,n+i]); {to adjust the
residual sum of squares in the rankdeficient case}
end; {loop on j -- end-loop is STEP 29}
writeln(‘Residual sum of squares = ’,trss);
end; {loop on i -- end-loop is STEP 30}
end; {if q>=0.0}
end; {while q>=0.0}
end; {alg04.pas -- Givens’ reduction, svd, and least squares solution}