Page 100 - Compact Numerical Methods For Computers

!The Choleski decomposition 89
Algorithm 7. Choleski decomposition in compact storage (cont.)
}
var
i,j,k,m,q: integer;
s : real; {accumulator}
begin
singmat := false; {singmat will be set true if matrix found
computationally singular}
for j := 1 to n do {STEP 1}
begin {STEP 2}
q := j*(i+1) div 2; {index of the diagonal element of row j}
if j>1 then {STEP 3}
begin {prepare for the subtraction in Eqn. (7.8). This is not needed
for the first column of the matrix.}
for i := j to n do {STEP 4}
begin
m := (i*(i-1) div 2)+j; s := a[m];
for k := 1 to (j-1) do s := s-a[m-k]*a[q-k];
a[m] := s;
end; {loop on i}
end; {of STEP 4}
if a[q]<=0.0 then {STEP 5}
begin {matrix singular}
singmat := true;
a[q] := 0.0; {since we shall assume matrix is non-negative definite)}
end;
s := sqrt(a[q]); {STEP 7}
for i := j to n do {STEP 8}
begin
m := (i*(i-1) div 2)+j;
if s=0.0 then a[m] := 0 {to zero column elements in singular case}
else a[m] := a[m]/s; {to perform the scaling}
end; (loop on i)
end; {loop on j -- end-loop is STEP 9}
end; {alg07.pass == Choleski decomposition choldcmp}

This completes the decomposition. The lower-triangular factor L is left in the vector a in
row-wise storage mode.

Algorithm 8. Choleski back-substitution
procedure cholback(n: integer; {order of problem}
a: smatvec; {the decomposed matrix}
var x: rvector); {the right hand side}
{alg08.pass ==
Choleski back substitution for the solution of consistent sets of
linear equations with symmetric coefficient matrices.
Copyright 1988 J.C.Nash