! Last change: DOS 24 Jul 2000 9:01 pm ! *** copyright 2000 *** ! *** filename chlzoi.f95 *** John F. Monahan ** ! ********************** PROGRAM PCHLZOI ! COMPUTE CHOLESKY FACTORS AND INVERSE OF POSITIVE DEFINITE MATRIX ! ONLY USES LOWER TRIANGULAR PART OF MATRIX ! ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! DECLARATIONS IMPLICIT NONE INTEGER, PARAMETER :: LENGTH = 10 REAL, DIMENSION( (LENGTH*(LENGTH+1))/2 ) :: A REAL DET INTEGER I,J,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 ') ! OPEN( UNIT=5, FILE='chlsky.dat') OPEN( UNIT=6, FILE='chlzoi.out') ! READ THE MATRIX A N = 4 DO I = 1,N READ(5,20) (A( (I*(I-1))/2 + J ),J=1,I) ENDDO ! LOOP ON I ! WRITE THE MATRIX BACK OUT WRITE(6,26) DO I = 1,N WRITE(6,21) (A( (I*(I-1))/2 + J ),J=1,I) ENDDO ! LOOP ON I ! COMPUTE CHOLESKY FACTORIZATION CALL CHLZKY(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*(I-1))/2 + J ),J=1,I) ENDDO ! LOOP ON I ! ! COMPUTE THE INVERSE AND OVERWRITE CALL CHLZOI(A,N) WRITE(6,28) DO I = 1,N WRITE(6,21) (A( (I*(I-1))/2 + J ),J=1,I) ENDDO ! STOP END PROGRAM PCHLZOI 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 CHLZOI(A,N) ! ! OVERWRITES L WITH INVERSE(L')*INVERSE(L) WHERE L IS A LOWER ! TRIANGULAR (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 JULY 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A 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))/2 ) = 1./A( (K*(K+1))/2 ) ! 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*(I-1))/2 + J )*A( (J*(J-1))/2 + K ) ENDDO ! LOOP ON J A( (I*(I-1))/2 + K ) = S/A( (I*(I+1))/2 ) 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*(K-1))/2 + I )*A( (K*(K-1))/2 + J ) ENDDO ! LOOP ON K ! STORE THE INNER PRODUCT A( (I*(I-1))/2 + J ) = S ENDDO ! LOOP ON I ENDDO ! LOOP ON J RETURN END SUBROUTINE CHLZOI SUBROUTINE CHLZKY(A,N,DET,IDET) ! CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT) FACTORIZATION ! A = L * ( L - TRANSPOSE ) WHERE L IS LOWER TRIANGULAR ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! 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 JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL, INTENT(OUT) :: DET 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+1))/2 ) ! FIRST ROW IS TRIVIAL IF(I .GT. 1) THEN IM1 = I-1 DO J=1,IM1 ! WORK ON (I,J)-TH ELEMENT S = A( (I*(I-1))/2 + J ) IF(J .GT. 1) THEN JM1 = J-1 DO K = 1,JM1 S = S - A( (J*(J-1))/2 + K ) * A( (I*(I-1))/2 + K ) ENDDO ! LOOP ON K ENDIF A( (I*(I-1))/2 + J )= S / A( (J*(J+1))/2 ) ! WORK ON DIAGONAL ELEMENT T = T - A( (I*(I-1))/2 + J ) * A( (I*(I-1))/2 + 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+1))/2 ) = SQRT(T) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZKY ! *** end of filename chlzoi.f95 *********************