/* WR-52 */ /* SYMMETRIC DECOMPOSITION OF POSITIVE DEFINITE BAND MATRICES */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 52-54. */ /* The lower half of a positive definite symmetric band matrix, A, with m lines on either side of the diagonal is stored as an n x (m+1) array, a[i,k], i = 1(1)n, k = 0(1)m, a[i,m] being the diagonal elements. The Cholesky decomposition A = L U, where U is the transpose of L, is performed and L is stored in l[i,k] in the same form as A. The reciprocals of the diagonal elements are stored instead of the elements themselves. A is retained so that the solution obtained can subsequently be improved. However, L and A can be identified in the call of the procedure. 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. */ /* Input: n = order of the matrix A. m = number of lines of A on either side of the diagonal. a = elements of the positive definite band matrix A. They must be given as described in Section 5. [see text] Output: l = the elements of the lower-triangle of the Cholesky decomposition of A stored as described in section 5. [see text] d1 and d2 = two numbers which together give the determinant of A. fail = exit used when A, possibly as a result of the rounding errors is not positive definite. */ CHOBANDDET: PROC(N,M,A,L,D1,D2,FAIL) OPTIONS (REORDER); /*PAGE 52*/ /*00000010*/ DCL (M,N,D2,I,J,K,P,Q,R,S,LBN,UBN,LBM,UBM) BIN FIXED(31,0), /*00000020*/ D1 BIN FLOAT (21), /*00000030*/ (A,L)(*,*) BIN FLOAT(21), /*00000040*/ FAIL LABEL, /*00000050*/ Y BIN FLOAT (53); /*00000060*/ LBN=LBOUND(A,1); UBN=LBN+N-1; /*00000070*/ LBM=LBOUND(A,2); UBM=LBM+M; /*00000080*/ D1=1E0;D2=0; /*00000090*/ DO I=LBN TO UBN; /*00000100*/ IF I-LBN+1>M THEN P=LBM; ELSE P=UBM-(I-LBN+1)+1; /*00000110*/ IF P=LBM THEN R=I-M; ELSE R=LBN; /*00000120*/ DO J=P TO UBM; /*00000130*/ S=J-1; Q=UBM -J+P; Y=A(I,J); /*00000140*/ DO K=P TO S; /*00000150*/ Y=Y-L(I,K)*L(R,Q);Q=Q+1; /*00000160*/ END; /*00000170*/ IF J=UBM THEN DO; /*00000180*/ D1=D1*Y; /*00000190*/ IF Y=0E0 THEN DO; /*00000200*/ D2=0;GO TO FAIL; /*00000210*/ END; /*00000220*/ DO WHILE (ABS(D1)>=1E0); /*00000230*/ D1=D1*.0625E0;D2=D2+4; /*00000240*/ END; /*00000250*/ DO WHILE (ABS(D1)<.0625E0); /*00000260*/ D1=D1*16E0;D2=D2-4; /*00000270*/ END; /*00000280*/ IF Y<0E0 THEN GO TO FAIL; /*00000290*/ L(I,J)=1E0/SQRT(Y); /*00000300*/ END; /*00000310*/ ELSE DO; /*00000320*/ L(I,J)=Y*L(R,UBM); R=R+1; /*00000330*/ END; /*00000340*/ END; /*00000350*/ END; /*00000360*/ END CHOBANDDET; /*00000370*/ /* page 53 */ /* Solves A x = b, where A is a positive definite symmetric band matrix with m lines on either side of the diagonal and b is an n x r matrix of r right-hand sides. The procedure CHOBANDSOL must be preceded by CHOBANDDET in which L is produced in l[i,k], from A. A x = b is solved in two steps, L y = b and U x = 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. */ /* Input: n = order of the matrix A. m = number of lines of A on either side of the diagonal. r = the number of right-hand sides for which A x = b is to be solved. l = the elements of the lower-triangle of the Cholesky decomposition of the positive definite matrix A as produced by procedure CHOBANDDET. b = the m x r matrix formed by the matrix of right-hand sides. Output: x = the n x r matrix fornmed by the solution vectors. */ CHOBANDSOL: PROC(N,M,R,L,B,X) OPTIONS (REORDER); /*PAGE 53*/ /*00000010*/ DCL (M,N,R,I,J,K,P,Q,S,LBN,UBN,LBM,UBM,LBR,UBR) BIN FIXED (31,0), /*00000020*/ Y BIN FLOAT (53), /*00000030*/ (A(LBN:UBN,LBM:UBM) CTL, /*00000040*/ (L,B,X)(*,*)) BIN FLOAT(21); /*00000050*/ LBN=LBOUND(L,1); UBN=LBN+N-1; /*00000060*/ LBM=LBOUND(L,2); UBM=LBM+M; /*00000070*/ LBR=LBOUND(B,2); UBR=LBR+R-1; /*00000080*/ S=UBM-1; /*00000090*/ DO J=LBR TO UBR; /*00000100*/ /*SOLUTION OF LY=B*/ /*00000110*/ DO I=LBN TO UBN; /*00000120*/ IF I-LBN+1>M THEN P=LBM; ELSE P=UBM-(I-LBN+1)+1; /*00000130*/ Q=I; Y=B(I,J); /*00000140*/ DO K=S BY -1 TO P; /*00000150*/ Q=Q-1; Y=Y-L(I,K)*X(Q,J); /*00000160*/ END; /*00000170*/ X(I,J)=Y*L(I,M); /*00000180*/ END; /*00000190*/ /*SOLUTION OF UX=Y*/ /*00000200*/ DO I=UBN TO LBN BY -1; /*00000210*/ IF UBN-I>M THEN P=LBM; ELSE P=UBM-UBN+I; /*00000220*/ Y=X(I,J);Q=I; /*00000230*/ DO K=S BY -1 TO P; /*00000240*/ Q=Q+1; Y=Y-L(Q,K)*X(Q,J); /*00000250*/ END; /*00000260*/ X(I,J)=Y*L(I,M); /*00000270*/ END; /*00000280*/ END; /*00000290*/ END CHOBANDSOL; /*00000300*/