/***************************************************************************/
/* */
/* ACM ALGORITHM 150 Invert Positive Definite Symmetric Matrix. */
/* */
/* H. Rutishauser, Eigd. Technische Hochschule, Zurich, Switzerland. */
/* */
/***************************************************************************/
/* Translated from Algol to PL/I by R. A. Vowels, 1 December 2008. */
/* SYMINV2 obtains the inverse of a positive definite symmetric matrix A */
/* of order N by a method that is similar to that given by Busing and Levy */
/* (Comm. ACM 5 (1962), 446) but requires no interchange of rows and */
/* columns nor storage space for an additional matrix Q, yet is */
/* numerically equivalent. The procedure required the upper triangular */
/* part of A to be given and overwrites it by the upper triangular part of */
/* the inverse which is again denoted by A. All pivots are chosen on the */
/* diagonal, and if all further diagonal elements that are elegible as */
/* pivots vanish (this is impossible for a positive definite matrix A) */
/* then the condition SINGULAR is raised. */
/* This procedure is approximately twice as fast as ACM Algorithm 120 */
/* which inverts a regular matrix. */
SYMINV2: PROCEDURE (A);
declare a(*,*) float (18);
declare bigajj float (18);
declare (i, j, k, n) fixed binary;
declare (p, q) (hbound(a,1)) float (18);
declare r(hbound(a,1)) bit (1) aligned;
n = hbound(a,1);
r = '1'b;
grand_loop:
do i = 1 to n;
/* Search for pivot: */
bigajj = 0;
do j = 1 to n;
if r(j) & abs(a(j,j)) > bigajj then
do;
bigajj = abs(a(j,j));
k = j;
end;
end;
if bigajj = 0 then signal condition(singular);
/* Preparation of elimination step i: */
r(k) = '0'b;
q(k) = 1/a(k,k);
p(k) = 1;
a(k,k) = 0;
do j = 1 to k-1;
p(j) = a(j,k);
if r(j) then q(j) = -a(j,k) * q(k) ; else q(j) = a(j,k) * q(k);
a(j,k) = 0;
end;
do j = k+1 to n;
if r(j) then p(j) = a(k,j); else p(j) = -a(k,j);
q(j) = -a(k,j) * q(k);
a(k,j) = 0;
end;
/* Elimination proper: */
do j = 1 to n;
do k = j to n;
a(j,k) = a(j,k) + p(j) * q(k);
end;
end;
end grand_loop;
END SYMINV2;