module igmcore_SimLinearAlgebraEq

  use dc_types, only: &
    & DP

  implicit none
  private

  public :: LUdecomp_TriDiagMat, sol_SimLinearAlEq_lutri
contains
subroutine LUdecomp_TriDiagMat( &
  & al, ad, au &
  & )

  real(DP), intent(inout) :: al(:)
  real(DP), intent(inout) :: ad(:)
  real(DP), intent(inout) :: au(:)

  integer :: i
  integer :: N

  N = size(ad)

  do i=1, N-1
    al(i+1) = al(i+1) / ad(i)
    ad(i+1) = ad(i+1) - au(i) * al(i+1)
  end do

end subroutine LUdecomp_TriDiagMat

function sol_SimLinearAlEq_lutri( &
  & al, ad, au, b &
  & ) result(answerVec)

  real(DP), intent(in) :: al(:)
  real(DP), intent(in) :: ad(:)
  real(DP), intent(in) :: au(:)
  real(DP), intent(in) :: b(:)
  real(DP) :: answerVec(size(ad,1))

  integer :: N
  integer :: i

  N = size(ad,1)

  answerVec(1) = b(1)

  do i=1, N-1
    answerVec(i+1) =  b(i+1) - answerVec(i) * al(i+1)
  end do

  answerVec(N) = answerVec(N) / ad(N)

  do i=N-1, 1, -1
    answerVec(i) = ( answerVec(i) - au(i) * answerVec(i+1) ) / ad(i)
  end do

end function sol_SimLinearAlEq_lutri

end module igmcore_SimLinearAlgebraEq
