program cholesky
  implicit none
  integer, parameter    :: n = 4  ! number of rows/cols in matrix
  integer               :: i,j    ! loop counter
  real, dimension(4,4)  :: A      ! input matrix
  real, dimension(4)    :: p      ! diagonal elements of Cholesky factor L 
  real choldet, det

  data A / &
       +5.000000,  -1.000000,  -1.000000,  -1.000000, &
       -1.000000,   5.000000,  -1.000000,  -1.000000, &
       -1.000000,  -1.000000,   5.000000,  -1.000000, &
       -1.000000,  -1.000000,  -1.000000,   5.000000/
  
  ! Det = 432.000000
  
  ! wykonaj rozkład cholesky'ego
  call choldcmp(A,n,p)
  !call choldc(A,n,p)

  ! wypisz wynik
  print *,'Wynik:'
  print *, p(1)
  do i = 2,n
     print "(4f10.4)", A(i,1:i-1),p(i)
  end do

  ! Wyznacznik a
  det = product(p)
  det = det*det
  print *, 'Det = ', det, choldet(a,n)

end program cholesky


! ======================================================================
subroutine choldc (a,n,p)
  implicit none
  integer                    :: i,j,k,n
  real, dimension(n,n)       :: a
  real, dimension(n)         :: p
  real                       :: asum
  !
  ! given a positive definite symmetric matrix a(1:n,1:n), with
  ! physical dimensions np, this routine constructs its Cholesky
  ! decomposition, A=L L^T. On input, only the upper triangle of
  ! a need to be given; it is not modified. The Cholesky factor L
  ! is returned in the lower triangle of a, except for its diagonal
  ! elements which are returned in p(1:n).
  ! (c) Numerical recipes
  !
  do i = 1,n
     do j = i,n
        asum = a(i,j)
        do k = i-1,1,-1
           asum = asum-a(i,k)*a(j,k)
        end do
        if(i.eq.j) then
           if (asum.le.0.) stop 'choldc failed'
           p(i) = sqrt(asum)
        else
           a(j,i) = asum / p(i)
        end if
     end do
  end do
  return
end subroutine choldc


! ======================================================================
subroutine choldcmp(A,n,p)
  implicit none
  integer :: n      ! liczba kolumn/wierszy
  real    :: A(n,n) ! macierz
  real    :: p(n)   ! diagonalne elementy faktora Cholesky'ego L
  integer :: j
  !
  ! Given a positive-definite symmetric matrix a(1:n,1:n),
  ! this routine constructs its Cholesky decomposition, A = L·L'.
  ! The Cholesky factor L is returned in the lower triangle of
  ! a, except for its diagonal elements which are returned in p(1:n).
  !
  do j = 1,n
     ! składowe diagonalne
     p(j) = sqrt(A(j,j) - dot_product(A(j,1:j-1),A(j,1:j-1)))
     if (p(j) .eq. 0) stop 'choldcmp failed.'
     ! skladowe pozadiagonalne
     if (j < n) then
        A(j+1:n,j) = (A(j+1:n,j) - &
             matmul(A(j+1:n,1:j-1),A(j,1:j-1))) / p(j)
     end if
  end do
  return
end subroutine choldcmp


! ======================================================================
function choldet(a,n)
  integer n
  real a(n,n)
  real, pointer :: p(:)
  real, pointer :: b(:,:)
  real d
  real choldet
  !
  ! Determinant of the matrix from Cholesky decomposition
  ! Uses: choldc(mat,int)
  !
  allocate( p(n), b(n,n) )
  d=1.0
  b = a
  call choldc(b,n,p)
  do i = 1, n
     d = d * p(i)
  end do
  choldet = d * d
  deallocate(p, b)
  return
end Function choldet

