Page 62 - Compact Numerical Methods For Computers

Handling larger problems 51

Thus the zeros below R multiply the last (m – n) elements of
(4.12)

where d 1 is of order n and d 2 of order (m – n). Thus
T
T
A Ax = R Rx
T T
= R d + 0d = A b. (4.13)
l
2
These equations are satisfied regardless of the values in the vector d 2 by
solutions x to the triangular system
Rx = d . (4.14)
1
This system is trivial to solve if there are no zero elements on the diagonal of R.
Such zero elements imply that the columns of the original matrix are not linearly
independent. Even if no zero or ‘small’ element appears on the diagonal, the
original data may be linearly dependent and the solution x to (4.14) in some way
‘unstable’.
Algorithm 3. Givens’ reduction of a real rectangular matrix

procedure givens(nRow,nCol : integer; (size of matrix to be decomposed)
var A, Q: rmatrix); (the matrix to be decomposed
which with Q holds the resulting decomposition)
{alg03.pas ==
Givens’ reduction of a real rectangular matrix to Q * R form.
This is a conventional version, which works one column at a time.
Copyright 1988 J. C. Nash
}
var
i, j, k, mn: integer; {loop counters and array indices}
b, c, eps, p, s : real; {angle parameters in rotations}
begin
writeln(‘alg03.pas -- Givens’,chr(39),’ reduction -- column-wise’);
{STEP 0 -- partly in the procedure call}
mn := nRow; if nRow>nCol then mn := nCo1; {mn is the minimum of nRow
and nCo1 and gives the maximum size of the triangular matrix
resulting from the reduction. Note that the decomposition is
still valid when nRow<nCol, but the R matrix is then
trapezoidal, i.e. an upper trianglular matrix with added
columns on the right.}
for i := 1 to nRow do
begin {set Q to a unit matrix of size nRow}
for j := 1 to nRow do Q[i,j] := 0.0;
Q[i,i] := 1.0;
end; {loop on i -- Q now a unit matrix}
eps := calceps; {the machine precision}
for j := 1 to (mn-1) do {main loop on diagonals of triangle} {STEP 1}
begin {STEP 2}
for k := (j+1) to nRow do {loop on column elements}
begin {STEP 3}
c := Au[j,j]; s := A[k,j]; {the two elements in the rotation}