module m_cholesky
contains
  subroutine cholesky(Aa,nn) 
    integer, intent(in) :: nn 
    real, intent(inout), dimension(1:nn,1:nn) :: Aa 
    integer :: i, j, k 
    real :: sum=0 
!    print*, '<><><><><><><> Cholesky is working for you<><><><><><><>'
    do i=1, nn 
       do j=1, i-1 
          sum=0 
          do k=1, j-1 
             sum=sum + Aa(i,k)*Aa(j,k) 
          end do
          Aa(i,j)=(1/Aa(j,j))*(Aa(i,j)-sum)  
       end do
       sum=0 
       do k=1, i-1 
          sum=sum+(Aa(i,k))**2 
       end do
       Aa(i,i)=sqrt(Aa(i,i)-sum)
       do j=i+1, nn ! Lower triangle sqrt matrix
          Aa(i,j)=0 
       end do
    end do
!    print*, '<><><><><><><> Cholesky made good job!! <><><><><><><>'
    return  
  end subroutine cholesky

!         Exponential covariance function
  function exponential(h)
    real, intent(in) :: h
    real exponential

    if(h<0) stop 'm_cholesky exponential covariance: negative distance'
    exponential=exp(-h)
  end function exponential

!         Gaussian covariance function
  function gaussian(h)
    real, intent(in) :: h
    real gaussian

    if(h<0) stop 'm_cholesky gaussian covariance: negative distance'
    gaussian=0.999*exp(-h*h)
    if (h.eq.0.0) gaussian=1.0
  end function gaussian

!         spherical covariance function
  function spherical(h)
    real, intent(in) :: h
    real spherical

    if(h<0) stop 'm_cholesky spherical covariance: negative distance'
    if(h>1.) then
       spherical=0.
    else
       spherical=1.-1.5*h+0.5*h*h*h
    endif
  end function spherical

end module m_cholesky