%PROCESS MACRO, MDECK;
/* ITERATIVE REFINEMENT OF THE SOLUTION OF A */
/* POSITIVE DEFINITE SYSTEM OF EQUATIONS. */
/* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation",
Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 34-39. */
/* Translated to PL/I from Algol 60 by R. A. Vowels, 5 December 2007. */
/* NOTE: the parameter EPS has been removed from the parameter list. */
/* Status: Compiled, but not executed. */
/* P. 34 - 39 */
%innerprod: procedure (initval, inc, finalval, c, a, b, cv, result);
declare (initval, inc, finalval, cv) character;
declare (c, a, b, result) character;
answer ('Temp = ' || c || ';' ) skip noscan col (16);
answer ('do ' || cv || ' = ' || initval || ' to ' || finalval || ' by ' ||
inc || ';' ) skip noscan col (16);
answer ('Temp = Temp + ' || a || '*' || b || ';' ) skip noscan col (19);
answer ('end;') skip col(16);
answer (result || ' = Temp' ) skip noscan col (16);
%end innerprod;
%activate innerprod;
accsolve: procedure (r, a, p, b, x, bb, l, ill) ;
declare (r, l) fixed binary (31),
(a, p, x, b, bb)(*,*) float (18),
Temp float (18),
ill label;
/* Solves Ax = b where A is an n x n positive definite symmetric
matrix and b is an n x r matrix of right hand sides, using the
procedure cholsol1. The procedure must be preceded by choldet1
in which L is produced in a[i,j] and p[i]. The residuals
bb = b - Ax are calculated and Ad=bb is solved, overwriting
d on bb. The refinement is repeated, as long as the maximum
correction at any stage is less than half that at the previous
stage, until the maximum correction is less than 2 eps times
the maximum abs (x). Exits to label ill if the solution fails
to improve. Uses the procedure innerprod which forms accurate
innerproducts. l is the number of iterations. */
declare (i, j, k, d2, n) fixed binary (31),
(d0, d1, c, cc, xmax, bbmax) float (18);
n = hbound(a,1);
do j = 1 to r;
do i = 1 to n;
x(i,j) = 0;
bb(i,j) = b(i,j);
end;
end;
l = 0;
d0 = 3;
L3: call cholsol1 (n, r, a, p, bb, bb) ;
l = l + 1;
d2, d1 = 0;
do j = 1 to r;
do i = 1 to n;
x(i,j) = x(i,j) + bb(i,j);
end;
end;
do j = 1 to r;
xmax, bbmax = 0;
do i = 1 to n;
if abs( x(i,j)) > xmax then xmax = abs(x(i,j));
if abs(bb(i,j)) > bbmax then bbmax = abs(bb(i,j));
innerprod (1, 1, i-1, -b(i,j), a(k,i), x(k,j), k, c);
innerprod (i, 1, n, c, a(i,k), x(k,j), k, c);
bb(i,j) = -c;
end;
if bbmax > d1*xmax then d1 = bbmax / xmax;
if bbmax > 2 * epsilon(xmax) * xmax then d2 = 1;
end;
if d1 > d0/2 then go to ill;
d0 = d1;
if d2 = 1 then goto L3;
end accsolve;
accinverse: procedure (n, a, l, fail, ill);
declare (n, l) fixed binary (31),
a(*,*) float (18),
Temp float (18),
(fail, ill) label;
/* The upper triangle of a positive definite symmetric matrix, A,
is stored in the upper triangle of an (n+1) x n array a[i,j],
i=l(l)n+l, j=1(1)n. X, the inverse of A, is formed in the
remainder of the array a[i,j] by the procedure cholinversionl.
The inverse is improved by calculating X=X+Z until the correction,
Z, is such that maximum abs(Z[i,j]) is less than 2eps times
maximum abs(X[i,j]), where Z=XB and B = I - AX. b is an n x n
array and z is a 1 x n array, X being overwritten a row at
a time. Exits to the label fail if A is not positive definite
and to the label ill if the maximum correction at any stage
is not less than half that at the previous stage. l is the
number of corrections applied. Uses the procedure innerprod. */
declare (i,j,k,j1) fixed binary (31);
declare (c, d, xmax, zmax, e) float (18);
declare ( b(n,n), z(n)) float (18);
e = 1; l = 0;
call cholinversion1 (n, a, fail);
L1: do i = 1 to n;
do j = 1 to n;
j1 = j + 1;
if j >= i then
do;
innerprod (1, 1, i, -(i=j), a(k,i), a(j1,k), k, c);
innerprod (i+1, 1, j, c, a(i,k), a(j1,k), k, c);
innerprod (j1, 1, n, c, a(i,k), a(k+1,j), k, c);
end;
else
do;
innerprod ( 1, 1, j, 0, a(k,i), a(j1,k), k, c);
innerprod (j1, 1, i, c, a(k,i), a(k+1,j), k, c);
innerprod (i+1,1, n, c, a(i,k), a(k+1,j), k, c);
end;
b(i,j) = -c;
end;
end;
xmax, zmax = 0;
do i = 1 to n;
do j = 1 to i;
innerprod(1, 1, i, 0, a(i+1,k), b(k,j), k, c);
innerprod(i+1, 1, n, c, a(k+1,i), b(k,j), k, z(j));
end;
do j = 1 to i;
c = abs(a(i+1,j));
d = abs(z(j));
if c > xmax then xmax = c;
if d > zmax then zmax = d;
a(i+1,j) = a(i+1,j)+z(j);
end;
end;
l = l + 1;
d = zmax/xmax;
if d > e/2 then go to ill;
e = d;
if d > 2 * epsilon(d) then go to L1;
end accinverse;