cholesky.f90

SUBROUTINE CHOLESKY(A_IN,A_OUT,N) 

! TRANSLATED AND SIMPLIFIED FROM EMBEDDED MATLAB FUNCTION. THIS SIMPLIFIED
! VERSION ONLY FUNCTIONS FOR REAL VALUES OF THE INPUT MATRIX A, FOR A POSITIVE
! DEFINITE MATRIX A, AND FOR SQUARE MATRIX A.

IMPLICIT NONE

INTEGER,INTENT(IN) :: N
REAL,INTENT(IN) :: A_IN(N,N)
REAL,INTENT(OUT) :: A_OUT(N,N)

!LOCALS
INTEGER :: R,C,J,K,W
REAL :: T,AKJ,AJJ,TEMP

!SET A_OUT EQUAL TO A_IN
A_OUT=A_IN
!SET ALL ENTRIES BELOW DIAGONAL TO 0. (TESTED WITH A 3X3 SYMMETRIC MATRIX)
DO C=1,N
  DO R=C+1,N
    A_OUT(R,C)=0. 
  END DO
END DO

!I TRANSLATED THIS 'BLIND' WITHOUT FIGURING OUT WHAT WAS GOING ON
DO J=1,N
  T=0.
  IF(J.GT.1)THEN
    DO K=1,J-1
      AKJ=A_OUT(K,J)
      T=T+AKJ*AKJ
    END DO
  END IF
  AJJ=A_OUT(J,J)-T
  IF(AJJ.LT.0.)THEN
    PRINT *,'Cholesky: Input matrix must be positive definite! Aborting.' 
    STOP
  ELSEIF(AJJ.EQ.0)THEN
    PRINT *, 'Cholesky: Ajj=0! Covariance matrix has 0.0 diagonal entries, ',&
       'has all identical entries or some other illegal condition. Aborting.'
    STOP
  END IF
  AJJ=SQRT(AJJ)
  A_OUT(J,J)=AJJ
  IF(J.LT.N)THEN
    DO K=J+1,N
      TEMP=0.
      IF(J.GT.1)THEN
        DO W=1,J-1
          TEMP=TEMP+A_OUT(W,J)*A_OUT(W,K)
        END DO
      END IF
      A_OUT(J,K)=(A_OUT(J,K)-TEMP)/AJJ
    END DO
  END IF
END DO

END SUBROUTINE CHOLESKY