Page 69 - Compact Numerical Methods For Computers

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

begin
for i := m to tcol do {Note: starts at column m, not column 1.}
begin
r := W[j,i];
W[j,i] := r*c+s*W[k,i];
W[k,i] := -r*s+c*W[k,i];
end;
end; {rotnsub}
begin {Givsvd}
writeln(‘alg04.pas -- Givens’,chr(39),
‘reduction, svd, and least squares solution’);
Write(‘Order of 1s problem and no. of right hand sides:’);
readln(infile,n,nRHS); {STEP 0}
if infname<>‘con’ then writeln(n,‘ ’,nRHS);
write(‘Enter a number to indicate end of data’);
readln(infile,endflag);
if infname<>‘con’ then writeln(endflag);
tcol := n+nRHS; {total columns in the work matrix}
k := n+1; {current row of interest in the work array during Givens’ phase}
for i := l to n do
for j := 1 to tcol do
W[i,j] := 0.0; {initialize the work array}
for i := 1 to nRHS do rss[i] := 0.0; {initialize the residual sums of squares}
{Note that other quantities, in particular the means and the total
sums of squares will have to be calculated separately if other
statistics are desired.}
eps := calceps; {the machine precision}
tol := n*n*eps*eps; {a tolerance for zero}
nobs := 0; {initially there are no observations]
{STEP 1 -- start of Givens’ reduction}
enddata := false; {set TRUE when there is no more data. Initially FALSE.}
while (not enddata) do
begin {main loop for data acquisition and Givens’ reduction}
getobsn(n, nRHS, W, k, endflag, enddata); {STEP 2}
if (not enddata) then
begin {We have data, so can proceed.} {STEP 3}
nobs := nobs+1; {to count the number of observations}
write(‘Obsn’,nobs,‘ ’);
for j := 1 to (n+nRHS) do
begin
write(W[k,j]:10:5,‘ ’);
if (7 * (j div 7) = j) and (j<n+nRHS) then writeln;
end;
writeln;
for j := 1 to (n+nRHS) do
begin {write to console file}
end;
for j := 1 to n do {loop over the rows of the work array to
move information into the triangular part of the
Givens’ reduction} {STEP 4}
begin
m := j; s := W[k,j]; c := W[j,j]; {select elements in rotation}