Page 63 - Compact Numerical Methods For Computers

52 Compact numerical methods for computers
Algorithm 3. Givens’ reduction of a real rectangular matrix (cont.)

b := abs(c); if abs(s)>b then b := abs(s);
if b>0 then
begin {rotation is needed}
c := c/b; s := s/b; {normalise elements to avoid over- or under-flow}
p := sqrt(c*c+s*s); {STEP 4}
s := s/p;
if abs(s)>=eps then {STEP 5}
begin {need to carry out rotations} {STEP 6}
c := c/p;
for i := 1 to nCo1 do {STEP 7 -- rotation of A}
begin
p := A[j,i]; A[j,i] := c*p+s*A[k,i]; A[k,i] := -s*p+c*A[k,i];
end; {loop on i for rotation of A}
for i := 1 to nRow do {STEP 8 -- rotation of Q. Note: nRow not nCo1.}
begin
p := Q[i,j]; Q[i,jl := c*p+s*Q[i,k]; Q[i,k] := -s*p+c*Q[i,k];
end; {loop on i for rotation of Q}
end; {if abs(s)>=eps}
end; {if b>0}
end; {loop on k -- the end-loop is STEP 9}
end; {loop on j -- the end-loop is STEP 10}
end; {alg03.pas == Givens’ reduction}

After the Givens’ procedure is complete, the array A contains the triangular factor R in rows 1
to mn. If nRow is less than nCol, then the right-hand (nCo1 – nRow) columns of array A contain the
transformed columns of the original matrix so that the product Q*R = A, in which R is now
trapezoidal. The decomposition can be used together with a back-substitution algorithm such as
algorithm 6 to solve systems of linear equations.
The order in which non-zero elements of the working array are transformed to zero is not
unique. In particular, it may be important in some applications to zero elements row by row
instead of column by column. The ,file alg03a.pas on the software disk presents such a row-wise
variant. Appendix 4 documents driver programs DR03.PAS and DR03A.PAS which illustrate
how the two Givens’ reduction procedures may be used.
Example 4.1. The operation of Givens’ reduction
The following output of a Data General ECLIPSE operating in six hexadecimal
digit arithmetic shows the effect of Givens’ reduction on a rectangular matrix. At
each stage of the loops of steps 1 and 2 of algorithm 3 the entire Q and A matrices
are printed so that the changes are easily seen. The loop parameters j and k as
well as the matrix elements c = A[j,j] and s = A[k,j] are printed also. In this
example, the normalisation at step 3 of the reduction is not necessary, and the
sine and cosine of the angle of rotation could have been determined directly from
the formulae (4.5), (4.6) and (4.7).
The matrix chosen for this example has only rank 2. Thus the last row of the
FINAL A MATRIX is essentially null. In fact, small diagonal elements of the
triangular matrix R will imply that the matrix A is ‘nearly’ rank-deficient.
However, the absence of small diagonal elements in R, that is, in the final array A,
do not indicate that the original A is of full rank. Note that the recombination of