Program cholesky
  Implicit None

  Integer, parameter :: pr=4
  Integer, parameter :: dim=3
  Real(kind=pr), dimension(:,:), allocatable :: A

  Allocate(A(dim,dim))

  A(1,1)=4;A(1,2)=-2;A(1,3)=4
  A(2,1)=-2;A(2,2)=10;A(2,3)=-5
  A(3,1)=4;A(3,2)=-5;A(3,3)=6

  Print *,A
  print *,decomp(A)

  Deallocate(A)

  Contains

  Function decomp(X)
    Real(kind=pr), Dimension(:,:), Intent(in) :: X
    Real(kind=pr), Dimension(dim,dim)         :: decomp
    Integer                                   :: i,j,k,n
    Real(kind=pr)                             :: temp

    n=dim
    decomp=0

    Do j=1,n

       temp=0._pr
       if (j>=2) Then
          Do k=1,j-1
             temp=temp+decomp(j,k)**2
          End Do
       End if

       decomp(j,j)=sqrt(X(j,j)-temp)

       Do i=j+1,n

          temp=0._pr
          Do k=1,j-1
             temp=temp+decomp(i,k)*decomp(j,k)
          End Do

          decomp(i,j)=(X(i,j)-temp)/decomp(j,j)

       End Do

    End Do
  End Function decomp



End Program cholesky