Page 101 - Compact Numerical Methods For Computers

90 Compact numerical methods for computers

Algorithm 8. Choleski back-substitution (cont.)
}
var
i,j,q : integer;
begin {Forward substitution phase -- STEP 1}
if a[1]=0.0 then x[1]:=0.0 {to take care of singular case}
else x[1]:=x[1]/a[1];
{STEP 2}
if n>1 then {do Steps 3 to 8; otherwise problem is trivial}
begin
{STEP 3}
q:=1; {to initialize the index of matrix elements}
for i:=2 to n do {STEP 4}
begin {STEP 5}
for j:=1 to (i-1) do
begin
q:=q+1;x[i]:=x[i]-a[q]*x[j];
end; {loop on j}
q:=q+1; {STEP 6 -- to give index of diagonal element of row i}
if a[q]=0.0 then x[i]:=0.0 {to handle singular case}
else x[i]:=x[i]/a[q]; {STEP 7}
end; {loop on i -- STEP 8}
end; {non-trivial case. This completes the forward substitution}
{STEP 9 -- Back substitution phase}
if a[n*(n+1) div 2]=0.0 then x[n]:=0.0 {for singular case}
else x[n]:=x[n]/a[n*(n+1) div 2];
if n>1 then {STEP 10} {test for trivial case; otherwise do steps 11 to 15}
begin {STEP 11}
for i:=n down to 2 do
begin {STEP 12}
q:=i*(i-1) div 2; {to give base index for row i}
for j:=1 to (i-1) do x[j]:=x[j]-x[i]*a[q+j]; {STEP 13}
if a[q]=0.0 then x[i-1]:=0.0 {for singular case}
else x[i-1]:=x[i-1]/a[q]; {STEP 14}
end; {loop on i -- STEP 15}
end; {non-trivial case -- STEP 16}
end; {alg08.pas == Choleski back-substitution cholback}

Note that this algorithm will solve consistent sets of equations whose coefficient matrices are
symmetric and non-negative definite. It will not detect the cases where the equations are not
consistent or the matrix of coefficients is indefinite.

7.3. SOME ORGANISATIONAL DETAILS
The algorithm given for the Choleski decomposition uses the most direct form for
computing the array indices. In practice, however, this may not be the most
efficient way to program the algorithm. Indeed, by running a FORTRAN version of
algorithm 7 against the subroutine of Healy (1968) and that of Kevin Price
(Nash 1984a, pp 97-101) on an IBM 370/168 it became quite clear that the latter two
codes were markedly faster in execution (using the FORTRAN G compiler) than my
own. On the Data General NOVA and Hewlett-Packard 9830, however, this