Page 68 - Compact Numerical Methods For Computers

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

One could save the rotations and carefully combine them to produce the U
matrix. However, this algorithm uses plane rotations not only to zero
elements in the data in the Givens’ reduction and to orthogonalize rows of
the work array in the svd portion of the code, but also to move the data
into place from the (n+1)st row of the working array into which the data
is read. These movements i.e. of the observation number nobs,
would normally move the data to row number nobs of the original
matrix A to be decomposed. However, it is possible, as in the array given
by data file ex04.cnm
3 1 <--- there are 3 columns in the matrix A
and 1 right hand side
-999 <--- end of data flag
1 2 3 1 <--- the last column is the RHS vector
2 4 7 1
2 2 2 1
5 3 1 1
-999 0 0 0 <--- end of data row
that this movement does not take place. This is because we use a complete
cycle of Givens’ rotations using the diagonal elements W[i,j], j := 1 to
n, of the work array to zero the first n elements of row nobs of the
(implicit) matrix A. In the example, row 1 is rotated from row 4 to row 1
since W is originally null. Observation 2 is loaded into row 4 of W, but
the first Givens’ rotation for this observation will zero the first TWO
elements because they are the same scalar multiple of the corresponding
elements of observation 1. Since W[2,2] is zero, as is W[4,2], the second
Givens’ rotation for observation 2 is omitted, whereas we should move the
data to row 2 of W. Instead, the third and last Givens’ rotation for
observation 2 zeros element W[4,3] and moves the data to row 3 of W.
In the least squares problem such permutations are irrelevant to the final
solution or sum of squared residuals. However, we do not want the
rotations which are only used to move data to be incorporated into U.
Unfortunately, as shown in the example above, the exact form in which such
rotations arise is not easy to predict. Therefore, we do not recommend
that this algorithm be used via the rotations to compute the svd unless
the process is restructured as in Algorithm 3. Note that in any event a
large data array is needed.
The main working matrix W must be n+l by n+nRHS in size.
Copyright 1988 J. C. Nash
}
var
count, EstRowRank, i, j, k, m, slimit, sweep, tcol : integer;
bb, c, e2, eps, p, q, r, s, tol, trss, vt : real;
enddata : boolean;
endflag : real;
procedure rotnsub; {to allow for rotations using variables local
to Givsvd. c and s are cosine and sine of
angle of rotation.}
var
i: integer;
r: real;