subroutine chol1(n, a, na)
! Cholesky decomposition of matrix A
! Ulrich Schmitt, 2002-12-31
implicit none
integer :: i, j, k, n, na
real, dimension(na, n) :: a
real :: sum
a(1,1) = sqrt(a(1,1))
do j = 2, n
  a(1,j) = a(1,j)/a(1,1)
  a(j,1) = 0.
end do
do i = 2, n-1
  sum = a(i,i)
  do k = 1, i-1
    sum = sum - a(k,i)**2
  end do
  a(i,i) = sqrt(sum)
  do j = i+1, n
    sum = a(i,j)
    do k = 1, i-1
      sum = sum - a(k,i)*a(k,j)
    end do
    a(i,j) = sum/a(i,i)
!    a(j,i) = a(i,j) ! *** error: inconsistent, a(j, 1) missing
    a(j,i) = 0.
  end do
end do
sum = a(n,n)
do k = 1, n-1
  sum = sum - a(k,n)**2
end do
a(n,n) = sqrt(sum)
end

subroutine chol2(n, a, na, b, x)
! Solution of A * x = b
! (uses result of chol1)
! Ulrich Schmitt, 2002-12-31
implicit none
integer :: i, j, n, na
real, dimension(na, n) :: a
real, dimension(n) :: b, x
real :: sum
x(1) = b(1)/a(1,1)
do i = 2, n
  sum = 0.
  do j = 1, i-1
    sum = sum + a(j,i)*x(j)
  end do
  x(i) = (b(i) - sum)/a(i,i)
end do
x(n) = x(n)/a(n,n)
do i = n-1, 1, -1
  sum = 0.
  do j = i+1, n
    sum = sum + a(i,j)*x(j)
  end do
  x(i) = (x(i) - sum)/a(i,i)
end do
end
