! crout.f90  solve simultaneous equation  A X = Y  by Crout method

subroutine crout(n, sz, A, Y, X)
  implicit none
  integer, intent(in) :: n   
  integer, intent(in) :: sz   
  double precision, dimension(sz,sz), intent(inout) :: A
  double precision, dimension(sz), intent(inout) :: Y
  double precision, dimension(sz), intent(inout) :: X
  double precision :: sum
  integer :: i, j, k, m 

  do m=1,n
    ! below diagonal computations
    do i=m,n
      sum = 0.0
      do k=1,m-1
        sum = sum + A(i,k)*A(k,m) ! m was j
      end do ! k
      A(i,m) = A(i,m) - sum
    end do ! i

    ! above diagonal computations
    do j=m+1,n
      sum = 0.0
      if(m>1) then
        do k=1,m-1
          sum = sum + A(m,k)*A(k,j)
        end do ! k
      end if
      A(m,j) = (A(m,j)-sum)/A(m,m)
    end do ! j
    sum = 0.0
    if(m>1) then
      do k=1,m-1
        sum = sum + A(m,k)*Y(k) ! m was i
      end do ! k
    end if
    Y(m) = (Y(m)-sum)/A(m,m)
  end do ! m

  ! back substitute
  do i=n,1,-1
    sum = 0.0
    if(i<n) then
      do k=i+1,n
        sum = sum + A(i,k)*X(k)
      end do ! k
    end if
    X(i) = Y(i) - sum
  end do ! i
end subroutine crout