/* WR-17 */
/* SYMMETRIC DECOMPOSITION OF A POSITIVE DEFINITE MATRIX - CHOLESKY. */
/* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation",
Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 17-21. */
/* Input:
n = order of the matrix a.
a = elements of the positive definite matrix A, in lower triangular form,
stored as a linear array of n(n+1)/2 elements, as follows:
array A stored as
------- ---------
a11 a1
a21 a22 a2 a3
a31 a32 a33 a4 a5 a6
a41 a42 a43 a44 a7 a8 a9 a10
Output:
a = elements of the lower triangle L of the Cholesky decomposition of A,
stored as for the input.
d1 and d2 = two numbers which together give the determinant of A. */
/* The lower triangle of a positive definite symmetric matrix, A, is
stored row by row in a 1 x n(n+1)/2 array a[i],i = 1(1)n(n+1)/2.
The Cholesky decomposition A = LU, where U is the transpose of L,
is performed and L is overwritten on A. The reciprocals of the diagonal
elments are stored instead of the elements themselves. The determinant,
d1 x 2^d2, of A is also computed. The procedure will fail
if A, modified by the rounding erors, is not positive definite. */
CHOLDET2:PROC(N,A,D1,D2,FAIL) REORDER; /*PAGE 17*/
DCL (N,D2,I,J,K,P,Q,R) FIXED BIN(31,0),
A(*) BIN FLOAT(21), /*LOWER TRIANGLE*/
FAIL LABEL, X BIN FLOAT (53),
D1 FLOAT BIN(21),
(LB,UB) BIN FIXED(31);
LB=LBOUND(A,1); UB=LB+N-1;
D1=1E0;D2=0;P=LB-1;
DO I=LB TO UB;
Q=P+1;R=LB-1;
DO J=LB TO I;
X=A(P+1);
DO K=Q TO P;
R=R+1;X=X-A(K)*A(R);
END;
R=R+1;P=P+1;
IF I=J THEN DO;
D1=D1*X;
IF X=0E0 THEN DO;
D2=0;GO TO FAIL;
END;
DO WHILE(ABS(D1)>=1E0);
D1=D1*.0625E0;D2=D2+4;
END;
DO WHILE(ABS(D1)<.0625E0);
D1=D1*16E0;D2=D2-4;
END;
IF X<0E0 THEN GO TO FAIL;
A(P)=1E0/SQRT(X);
END;
ELSE A(P)=X*A(R);
END;
END;
END CHOLDET2;
/* p. 18. */
/* Input:
n = order of the matrix a.
r = number of right-hand sides for which Ax = b is to be solved.
a = elements of the lower triangle L of the Cholesky decomposition of a
positive definite matrix A as produced by choldet2.
b = the n x r matrix consisting of the r right-hand sides.
Output:
b = the n x r matrix consisting of the r solution vectors. */
/* Solves Ax = b, where A is a positive definite symmetric matrix and
b is an n x r matrix of r right-hand sides. The procedure cholsol2
must be preceded by choldet2 in which L is produced in a[i], from A.
Ax = b is solved in two steps, Ly = b and UX = y. The matrices y
and then x are overwritten on b. */
CHOLSOL2:PROC(N,R,A,B) REORDER; /* PAGE 18 */
DCL (N,R,I,J,K,P,Q,S) BIN FIXED(31,0),
A(*) BIN FLOAT(21),
B(*,*) FLOAT BIN (21),
X FLOAT BIN (53),
(LB,UB,LBR,UBR) BIN FIXED(31);
LB=LBOUND(A,1); UB=LB+N-1;
LBR=LBOUND(B,2); UBR=LBR +R-1;
DO J=LBR TO UBR;
/* SOLUTION OF L y = b */
P=LB;
DO I=LB TO UB;
Q=I-1;X=B(I,J);
DO K=LB TO Q;
X=X-A(P)*B(K,J);P=P+1;
END;
B(I,J)=X*A(P);P=P+1;
END;
/* SOLUTION OF U x = y */
DO I=UB TO LB BY -1;
S,P=P-1;Q=I+1;X=B(I,J);
DO K=UB BY -1 TO Q;
X=X-A(S)*B(K,J);S=S-K+1;
END;
B(I,J)=X*A(S);
END;
END;
END CHOLSOL2;
/* p. 19 */
/* Input:
n = order of the matrix a.
a = elements of the positive definite matrix A, stored as a linear array
of n(n+1)/2 elements.
Output:
a = elements of inv(A), stored as linear array as described for the input. */
/* The lower triangle of a positive definite matrix, A, is stored row by
row in a 1 x n(n+1)/2 array a[i], i = 1(1)n(n+1)/2. The Cholesky
decomposition A = LU, where U is the transpose of L, is performed
and L is overwritten on A. The reciprocals of the diagonal elements
are stored instead of the elements themselved. L is then replaced by
its inverse and this in turn is replaced by the lower triangle of the
inverse of A. The procedure will fail is A, modified by the rounding
errors, is not positive definite. */
CHOLINVERSION2:
PROC(N,A,FAIL) REORDER; /*PAGE 19*/
DCL (N,I,J,K,P,Q,R,S,T,U) BIN FIXED (31,0),
Y BIN FLOAT(21), X BIN FLOAT(53),
A(*) BIN FLOAT (21),
(LB,UB) BIN FIXED (31),
FAIL LABEL;
/*FORMATION OF L */
LB=LBOUND(A,1); UB=LB+N-1;
P=LB-1;
DO I=LB TO UB;
Q=P+1;R=LB-1;
DO J=LB TO I;
X=A(P+1);
DO K=Q TO P;
R=R+1;X=X-A(K)*A(R);
END;
R=R+1;P=P+1;
IF I=J THEN DO;
IF X<=0E0 THEN GO TO FAIL;
A(P)=1E0/SQRT(X);
END;
ELSE A(P)=X*A(R);
END;
END;
/* INVERSION OF L */
P,R,T=LB-1;
DO I=LB+1 TO UB;
P=P+1;R=R+I-LB+1;Y=A(R+1);
DO J=LB+1 TO I;
P=P+1;S,T=T+1;U=I-LB-1;X=0E0;
DO K=R BY -1 TO P;
X=X-A(K)*A(S);S=S-U;U=U-1;
END;
A(P)=X*Y;
END;
END;
/*CALCULATION OF THE INVERSE OF A */
P=LB-1;
DO I=LB TO UB;
R=P+I-LB+1;
DO J=LB TO I;
S=R;T,P=P+1;X=0E0;
DO K=I-LB+1 TO N;
X=X+A(S)*A(T);S=S+K;T=T+K;
END;
A(P)=X;
END;
END;
END CHOLINVERSION2;