Page 86 - Compact Numerical Methods For Computers

Linear equations—a direct approach 75
If is zero, we cannot proceed, and ‘small’ are quite likely to occur during
subtractions involving digit cancellations, so that multipliers m ij that are large and
inaccurate are possible. However, we can ensure that multipliers m are all less
ij
than one in magnitude by permuting the rows of A' (and hence A) so that the
largest of for i = j, (j + 1), . . . , n, is in the diagonal or pivot position. This
modified procedure, called Gauss elimination with partial pivoting, has a large
literature (see Wilkinson (1965) for a discussion with error analysis). Since the
rows of A have been exchanged, this procedure gives the triangular decomposition
of a transformed matrix
PA = LR (6.27)

where P is the permutation matrix, simply a re-ordering of the columns of a unit
matrix appropriate to the re-ordering of the rows of A.
Particular methods exist for further minimising error propagation in the Gauss
elimination procedure by using double-precision accumulation of vector inner
products. These go by the names of Crout and Doolittle and are discussed, for
instance, by Dahlquist and Björck (1974) as well as by Wilkinson (1965). Since the
accumulation of inner products in several computer programming languages on a
variety of machines is a non-trivial operation (though on a few machines it is
simpler to accumulate in double than in single precision), these will not be
discussed here. Nor will the Gauss elimination with complete pivoting, which
chooses as pivot the largest element in the current matrix, thereby requiring both
column and row permutations. The consensus of opinion in the literature appears
to be that this involves too much work for the very few occasions when complete
pivoting is distinctly more accurate than partial pivoting.
The Gauss elimination with partial pivoting and back-substitution are now
stated explicitly. All work is done in an array A of dimension n by n + p, where p
is the number of right-hand sides b to be solved.
Algorithm 5. Gauss elimination with partial pivoting
Procedure gelim ( n : integer; {order of equations}
p : integer; {number of right hand sides}
var A : rmatrix; {equation coefficients in row order with right
hand sides built into the matrix columns n+1, . . ,n+p}
to1 : real); {pivot tolerance}

{alg05.pas == Gauss elimination with partial pivoting.
This form does not save the pivot ordering, nor does it keep
the matrix decomposition which results from the calculations.
Copyright 1988 J. C. Nash
}
var
det, s : real;
h,i,j,k: integer;
begin {STEP 0}
det := 1.0; {to initialise determinant value}
writeln(‘alg05.pas -- Gauss elimination with partial pivoting’);
for j := 1 to (n-1) do {STEP 1}