/* WR-15 */ /* 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. 15-17. */ /* Input: n = order of the matrix a. a = elements of the positive definite matrix A, in upper triangular form. Only the diagonal and above-diagonal elements are required. Output: a = elements of the lower triangle L of the Cholesky decomposition of A (excepting the diagonal elements). These are stored below the diagonal. p = a vector consisting of the reciprocals of the diagonal elements of L. d1 and d2 = two numbers which together give the determinant of A. */ /* The upper triangle of a positive definite symmetric matrix, A, is stored in the upper triangle of an m x n array a[i,j], i = 1(1)n, j = 1(1)n. The Cholesky decomposition A = LU, where U is the transpose of L, is performed and stored in the remainder of the array a except for the reciprocals of the diagonal ekements which are stored in p[i], i = 1(1)n, instead of the elements themselves. A is retained so that the solution obtained dcan be sunsequently improved. The determinant, d1 x 2^d2, of A is also computed. The procedure will fail if A, modified by the rounding errors, is not positive definite. */ CHOLDET1:PROC(N,A,P,D1,D2,FAIL); /*PAGE 15*/ /*00000010*/ DCL (N,D2) BIN FIXED (31,0),(I,J,K) BIN FIXED(31,0), /*00000020*/ D1 BIN FLOAT(21), X BIN FLOAT(53), /*00000030*/ A(*,*) BIN FLOAT(21), /*00000040*/ P(*) BIN FLOAT(21), /*00000050*/ FAIL LABEL, (LB,UB) BIN FIXED(31); /*00000060*/ LB=LBOUND(A,1); UB=LB+N-1; /*00000070*/ D1=1E0;D2=0; /*00000080*/ DO I=LB TO UB; /*00000090*/ DO J=I TO UB; /*00000100*/ X=A(I,J); /*00000110*/ DO K=I-1 BY -1 TO LB; /*00000120*/ X=X-A(J,K)*A(I,K); /*00000130*/ END; /*00000140*/ IF J=I THEN DO; /*00000150*/ D1=D1*X; /*00000160*/ IF X=0E0 THEN DO; /*00000170*/ D2=0;GO TO FAIL; /*00000180*/ END; /*00000190*/ DO WHILE(ABS(D1)>=1E0); /*00000200*/ D1=D1*.0625E0;D2=D2+4; /*00000210*/ END; /*00000220*/ DO WHILE(ABS(D1)<.0625E0); /*00000230*/ D1=D1*16E0;D2=D2-4; /*00000240*/ END; /*00000250*/ IF X<0E0 THEN GO TO FAIL; /*00000260*/ P(I)=1E0/SQRT(X); /*00000270*/ END; /*00000280*/ ELSE A(J,I)=X*P(I); /*00000290*/ END; /*00000300*/ END; /*00000310*/ END CHOLDET1; /*00000320*/ /* p. 16 */ /* 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 choldet1. p = a vector consisting of the reciprocals of the diagonal elements of L as produced by choldet1. b = the n x r matrix consisting of the r right-hand sides. Output: x = 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 cholsol1 must be preceded by choldet1 in which L is produced in a[i,j] and p[i], from A. Ax = b is solved in two steps, Ly = b and Ux = y. The matrix b is retained in order to facilitate the refinement of x, but x is overwritten on y. However, x and b can be identified in the call of the procedure. */ CHOLSOL1:PROCEDURE(N,R,A,P,B,X) REORDER; /* PAGE 15 */ DCL (R,N) BIN FIXED(31,0), (A,B,X) (*,*) BIN FLOAT(21), P(*) BIN FLOAT(21), (I,J,K) BIN FIXED(31,0), Z BIN FLOAT(53), (LB,UB,LBC,UBC) BIN FIXED(31); LB=LBOUND(A,1); UB=LB+N-1; LBC=LBOUND(A,2); UBC=LBC+R-1; DO J=LBC TO UBC; /* SOLUTION OF L y = b */ DO I=LB TO UB; Z=B(I,J); DO K=I-1 BY -1 TO LB; Z=Z-A(I,K)*X(K,J); END; X(I,J)=Z*P(I); END; /* SOLUTION OF U x = y */ DO I=UB TO LB BY -1; Z=X(I,J); DO K=I+1 TO UB; Z=Z-A(K,I)*X(K,J); END; X(I,J)=Z*P(I); END; END; END CHOLSOL1; /* p. 16. */ /* Input: n = order of the matrix a. a = elements of the positive definite matrix A, in upper triangular form. Only the diagonal and above-diagonal elements are required. Output: a = elements of inv(A). */ /* 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 = 1(1)n + 1, j = 1(1)n. The Cholesky decomposition A = LU, where U is the transpose of L, is performed and L is stored in the remainder of the array a. The reciprocals of the diagonal elements are stored instead of the elements themselves. L is then replaced by its inverse and this in turn is replaced by the lower triangle of the inverse of A. A is retained so that the inverse can be subsequently improved. The procedure will fail if A, modified by the rounding errors, is not positive definite. */ CHOLINVERSION1: PROC(N,A,FAIL); /*PAGE 16*/ /*00000010*/ DCL(N,I,J,K,I1,J1) BIN FIXED (31,0), /*00000020*/ A(*,*) BIN FLOAT(21), /*00000030*/ FAIL LABEL, Z BIN FLOAT (53), /*00000040*/ (X, Y) BIN FLOAT(21), /*00000050*/ (LB,UB) BIN FIXED(31); /*00000060*/ LB=LBOUND(A,1); UB=LB + N-1; /*00000070*/ /* FORMATION OF L */ /*00000080*/ DO I=LB TO UB; /*00000090*/ I1=I+1; /*00000100*/ DO J=I TO UB; /*00000110*/ J1=J+1;X=A(I,J); /*00000120*/ DO K=I-1 BY -1 TO LB; /*00000130*/ X=X-A(J1,K)*A(I1,K); /*00000140*/ END; /*00000150*/ IF J=I THEN DO; /*00000160*/ IF X<=0E0 THEN GO TO FAIL; /*00000170*/ A(I1,I),Y=1E0/PRECISION(SQRT(X),53); /*00000180*/ END; /*00000190*/ ELSE A(J1,I)=X*Y; /*00000200*/ END; /*00000210*/ END; /*00000220*/ /* INVERSION OF L */ /*00000230*/ DO I=LB TO UB; /*00000240*/ DO J=I+1 TO UB; /*00000250*/ Z=0E0;J1=J+1; /*00000260*/ DO K=J-1 BY -1 TO I; /*00000270*/ Z=Z-A(J1,K)*A(K+1,I); /*00000280*/ END; /*00000290*/ A(J1,I)=Z*A(J1,J); /*00000300*/ END; /*00000310*/ END; /*00000320*/ /* CALCULATION OF INVERSE OF A */ /*00000330*/ DO I=LB TO UB; /*00000340*/ DO J=I TO UB; /*00000350*/ Z=0;J1=N+1; /*00000360*/ DO K=J+1 TO J1; /*00000370*/ Z=Z+A(K,J)*A(K,I); /*00000380*/ END; /*00000390*/ A(J+1,I)=Z; /*00000400*/ END; /*00000410*/ END; /*00000420*/ END CHOLINVERSION1; /*00000430*/