Page 110 - Compact Numerical Methods For Computers

The symmetric positive definite matrix again 99
A difficulty in this is that the quantities /p have to be stored or they will
be overwritten by A , during the work. This is overcome by using a working
i-1 j-1
vector to store the needed quantities temporarily.
Because of the re-numeration we also have the matrix of coefficients in the
form

(8.24)

T T
This permits, via Z = Z and Y = – X , the computation involved in (8.23), where
A j1 for j < n-k
A =
1j (8.25)
- A for j > n-k
j l
without recourse to a full matrix. By using row-wise storage of the lower triangle
of A, the algorithm below is obtained. Note that after n re-numerations the array
elements are in the correct order. Also by counting backwards (step 1) from n to
1 the counter will contain (n - k + 1).

Algorithm 9. Bauer-Reinsch inversion of a positive definite symmetric matrix

procedure brspdmi(n : integer; {order of problem}
var avector : smatvec; {matrix in vector form}
var singmat : boolean); {singular matrix flag}
{alg09.pas ==
Bauer - Reinsch inversion of a symmetric positive definite matrix in
situ, Lower triangle of matrix is stored as a vector using row
ordering of the original matrix.

Copyright 1988 J.C.Nash
}
var
i,j,k,m,q : integer;
s,t : real;
X : vector;
begin
writeln(‘alg09.pas -- Bauer Reinsch inversion’);
singmat := false; {initial setting for singularity flag}
for k := n downto 1 do {STEP 1}
begin
if (not singmat) then
begin
s := avector[1]; {STEP 2}
if s>0.0 then {STEP 3}
begin
m := 1; {STEP 4}
for i := 2 to n do {STEP 5}
begin {STEP 6}
q := m; m := m+i; t := avector[q+1]; X[i] := -t/s;
{Note the use of the temporary vector X[]}
if i>k then X[i] := -X[i];{take account of Eqn. 8.22 -- STEP 7}