! Last change: DOS 24 Jul 2000 8:52 pm ! *** copyright 2000 *** ! *** filename chlsoi.f95 *** John F. Monahan ** ! ********************** PROGRAM pchlsoi ! COMPUTE CHOLESKY FACTORS AND INVERSE OF POSITIVE DEFINITE MATRIX ! ONLY USES LOWER TRIANGULAR PART OF MATRIX ! ! DECLARATIONS IMPLICIT NONE REAL, DIMENSION(10,10) :: A,ACOPY REAL, DIMENSION(10) :: B REAL DET INTEGER I,J,K,N,IDET ! FORMAT STATEMENTS 20 FORMAT(2X,4F6.2) 21 FORMAT(1X,10F13.8) 22 FORMAT(2X,'DETERMINANT IS ',F9.6,' * ( 2**',I6,' )') 25 FORMAT(' INTERMEDIATE SOLUTION Y THAT SOLVES L*Y = B ') 26 FORMAT(' INPUT MATRIX ') 27 FORMAT(' FACTORED MATRIX -- LOWER TRIANGULAR L ') 28 FORMAT(' INVERSE MATRIX ') 29 FORMAT(' IS IT THE IDENTITY? ') ! INTERFACE BLOCK INTERFACE SUBROUTINE CHLSKY(A,N,DET,IDET) REAL, DIMENSION(:,:), INTENT(IN OUT) :: A REAL, INTENT(OUT) :: DET INTEGER, INTENT(IN) :: N INTEGER, INTENT(OUT) :: IDET END SUBROUTINE CHLSKY SUBROUTINE CHLSOI(A,N) REAL, DIMENSION(:,:), INTENT(IN OUT) :: A INTEGER, INTENT(IN) :: N END SUBROUTINE CHLSOI END INTERFACE ! OPEN( UNIT=5, FILE='chlsky.dat') OPEN( UNIT=6, FILE='chlsoi.out') ! READ THE MATRIX A N = 4 DO I = 1,N READ(5,20) (A(I,J),J=1,N) ENDDO ! LOOP ON I ! WRITE THE MATRIX BACK OUT WRITE(6,26) DO I = 1,N WRITE(6,21) (A(I,J),J=1,N) ENDDO ! LOOP ON I ! SAVE A COPY OF MATRIX ACOPY = A ! COMPUTE CHOLESKY FACTORIZATION CALL CHLSKY(A,N,DET,IDET) ! WRITE OUT DETMINANT IN DET*(2**IDET) FORM WRITE(6,22) DET,IDET ! WRITE(6,27) DO I = 1,N WRITE(6,21) (A(I,J),J=1,N) ENDDO ! LOOP ON I ! ! COMPUTE THE INVERSE AND OVERWRITE CALL CHLSOI(A,N) WRITE(6,28) DO I = 1,N WRITE(6,21) (A(I,J),J=1,N) ENDDO ! LOOP ON I ! CHECK IT OUT -- GET A-INV * A WRITE(6,29) DO I = 1,N DO J = 1,N DET = 0. DO K = 1,N DET = DET + A(I,K)*ACOPY(K,J) ENDDO ! LOOP ON K B(J) = DET ENDDO ! LOOP ON J WRITE(6,21) (B(J),J=1,N) ENDDO ! LOOP ON I STOP END program pchlsoi SUBROUTINE ADJUST(D,I) ! NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I) INTEGER, PARAMETER :: IBIG = -2147483644 REAL, INTENT (IN OUT) :: D INTEGER, INTENT(IN OUT) :: I IF ( D .GT. 0.0 ) THEN DO WHILE ( D .LT. 1.0 ) D = D*16. I = I-4 ENDDO DO WHILE ( D .GT. 16.0 ) D = D/16. I = I+4 ENDDO ELSE I = IBIG ! IF D < 0 THEN I = - BIG ENDIF RETURN END SUBROUTINE ADJUST SUBROUTINE CHLSKY(A,N,DET,IDET) ! CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT) FACTORIZATION ! A = L * ( L - TRANSPOSE ) WHERE L IS LOWER TRIANGULAR ! A IS OVERWRITTEN BY L IN ITS LOWER TRIANGULAR PART ! SUBROUTINE ADJUST KEEPS DETERMINANT FROM EXPLODING USING ! 1 LE DET LE 16 AND DETERMINANT OF A IS DET*2**IDET ! ! J F MONAHAN (DEC,1983) DEPT OF STAT, N C S U, RALEIGH, N C 27650 USA ! REVISED JANUARY, 2000 FOR FORTRAN 95 IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A REAL, INTENT(OUT) :: DET INTEGER, INTENT(IN) :: N INTEGER, INTENT(OUT) :: IDET REAL T,S INTEGER I,J,K,IM1,JM1 ! INITIALIZE FOR DETERMINANT DET = 1. IDET = 0 ! DO I-TH ROW DO I = 1,N T = A(I,I) ! FIRST ROW IS TRIVIAL IF(I .GT. 1) THEN IM1 = I-1 DO J=1,IM1 A(J,I) = 0. ! WORK ON (I,J)-TH ELEMENT S = A(I,J) IF(J .GT. 1) THEN JM1 = J-1 DO K = 1,JM1 S = S - A(J,K) * A(I,K) ENDDO ! LOOP ON K ENDIF A(I,J)= S / A(J,J) ! WORK ON DIAGONAL ELEMENT T = T - A(I,J) * A(I,J) ENDDO ! LOOP ON J ! FINISHED WITH LOOP ON J ! UPDATE DETERMINANT WITH DIAGONAL ENDIF DET = DET * T ! KEEP DET FROM EXPLODING CALL ADJUST(DET,IDET) ! ABORT IF DET IS NEGATIVE OR ZERO IF(IDET.LT.-2000000000) RETURN ! FINISH A POSITIVE DIAGONAL A(I,I) = SQRT(T) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLSKY SUBROUTINE CHLSOI(A,N) ! ! OVERWRITES L WITH INVERSE(L')*INVERSE(L) WHERE L IS A LOWER ! TRIAGULAR (N BY N) MATRIX STORED IN A ! ! TO COMPUTE THE INVERSE OF A POSITIVE DEFINITE MATRIX A ! FIRST: FACTOR A = L * L' WHERE L IS LOWER TRIANGULAR ! USING CHLSKY(A,N,D,IDET) ! SECOND: COMPUTE THE INVERSE(A) = INVERSE(L') * INVERSE(L) ! USING CHLSOI(A,N) ! ! J F MONAHAN (JULY,1984) DEPT OF STATISTICS, N C S U, RALEIGH, N C 2769 ! REVISED JANUARY, 2000 FOR FORTRAN 95 REAL, DIMENSION(:,:), INTENT(IN OUT) :: A INTEGER, INTENT(IN) :: N REAL S INTEGER I,J,K,IM1,KP1 ! FIRST INVERT L IN PLACE DO K = 1,N ! SOLVE L * X = E(K) A(K,K) = 1./A(K,K) ! BRANCH AROUND IF K EQUALS N IF( K .LT. N ) THEN KP1 = K + 1 DO I = KP1,N S = 0. IM1 = I - 1 DO J = K,IM1 S = S - A(I,J)*A(J,K) ENDDO ! LOOP ON J A(I,K) = S/A(I,I) ENDDO ! LOOP ON I ENDIF ENDDO ! LOOP ON K ! DONE WITH INVERSION, NOW COMPUTE INNER PRODUCT DO J = 1,N DO I = J,N ! GET (I,J) ELEMENT OF MATRIX IS THE INNER PRODUCT ! OF COLUMN I AND COLUMN J S = 0. DO K = I,N S = S + A(K,I)*A(K,J) ENDDO ! LOOP ON K ! STORE THE INNER PRODUCT A(I,J) = S A(J,I) = S ENDDO ! LOOP ON I ENDDO ! LOOP ON J RETURN END SUBROUTINE CHLSOI ! *** end of filename chlsoi.f95 *********************