!******************************************************************************
! subroutine ludecomp
!******************************************************************************
!     Arrangement of array elements
!----------------------------------------------------------------------
!  A x = b
!          a(1,1)  a(1,2)  a(1,3)  a(1,4) ...  a(1,n)  |   x(1)  |   b(1)
!          a(2,1)  a(2,2)  a(2,3)  a(2,4) ...  a(2,n)  |   x(2)  |   b(2)
!          a(3,1)  a(3,2)  a(3,3)  a(3,4) ...  a(3,n)  |   x(3)  |   b(3)
!          a(4,1)  a(4,2)  a(4,3)  a(4,4) ...  a(4,n)  |   x(4)  |   b(4)
!          a(5,1)  a(5,2)  a(5,3)  a(5,4) ...  a(5,n)  |   x(5)  |   b(5)
!          ...     ...     ...     ...    ...  ...     |   ...   |   ...
!          a(n,1)  a(n,2)  a(n,3)  a(n,4) ...  a(n,n)  |   x(n)  |   b(n)
!----------------------------------------------------------------------
!  A = L * U
!          a(1,1)  a(1,2)  a(1,3)  a(1,4) ...  a(1,n)
!          a(2,1)  a(2,2)  a(2,3)  a(2,4) ...  a(2,n)
!          a(3,1)  a(3,2)  a(3,3)  a(3,4) ...  a(3,n)
!          a(4,1)  a(4,2)  a(4,3)  a(4,4) ...  a(4,n)
!          a(5,1)  a(5,2)  a(5,3)  a(5,4) ...  a(5,n)
!          ...     ...     ...     ...    ...  ...
!          a(n,1)  a(n,2)  a(n,3)  a(n,4) ...  a(n,n)
!----------------------------------------------------------------------
!  matrix L
!          1       0       0       0      ...  0
!          a(2,1)  1       0       0      ...  0
!          a(3,1)  a(3,2)  1       0      ...  0
!          a(4,1)  a(4,2)  a(4,3)  1      ...  0
!          a(5,1)  a(5,2)  a(5,3)  a(5,4) ...  0
!          ...     ...     ...     ...    ...  ...
!          a(n,1)  a(n,2)  a(n,3)  a(n,4) ...  a(n,n)
!----------------------------------------------------------------------
!  matrix L
!          a(1,1)  a(1,2)  a(1,3)  a(1,4) ...  a(1,n)
!          0       a(2,2)  a(2,3)  a(2,4) ...  a(2,n)
!          0       0       a(3,3)  a(3,4) ...  a(3,n)
!          0       0       0       a(4,4) ...  a(4,n)
!          0       0       0       0      ...  a(5,n)
!          ...     ...     ...     ...    ...  ...
!          0       0       0       0      ...  a(n,n)
!******************************************************************************

module ludecomp_module

  ! モジュール引用 ; USE statements
  !
  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP,  & ! 倍精度実数型. Double precision.
    &                 TOKEN  ! キーワード.   Keywords. 

  implicit none

  private

  public :: ludecomp_prep_simple_many
  public :: ludecomp_solve_simple_many

contains

  !**************************************************************************************

  subroutine ludecomp_prep_simple_many( mm, nn, mat, ms, me )

    integer , intent(in   ) :: mm, nn
    real(DP), intent(inout) :: mat( mm, nn, nn )
    integer , intent(in   ) :: ms, me


    !
    ! local variables
    !
    real(DP) :: lmat( ms:me, nn, nn )
    real(DP) :: ratio( ms:me )
    integer  :: i, j, k, m


    do j = 1, nn
       do i = 1, nn
          do m = ms, me
             lmat( m, i, j ) = 0.0d0
          end do
       end do
    end do
    do i = 1, nn
       do m = ms, me
          lmat( m, i, i ) = 1.0d0
       end do
    end do

    do k = 1, nn-1
       do i = k+1, nn
          do m = ms, me
             ratio( m ) = mat( m, i, k ) / mat( m, k, k )
          end do

          do j = 1, nn
             do m = ms, me
                mat( m, i, j ) = mat( m, i, j ) - mat( m, k, j ) * ratio( m )
             end do
          end do

          do m = ms, me
             lmat( m, i, k ) = ratio( m )
          end do
       end do

    end do

    !
    ! assemble into 1 matrix
    !
    do j = 1, nn
       do i = j+1, nn
          do m = ms, me
             mat( m, i, j ) = lmat( m, i, j )
          end do
       end do
    end do


  end subroutine ludecomp_prep_simple_many

  !**************************************************************************************

  subroutine ludecomp_solve_simple_many( mm, nn, mat, vec, ms, me )

    integer , intent(in   ) :: mm, nn
    real(DP), intent(in   ) :: mat( mm, nn, nn )
    real(DP), intent(inout) :: vec( mm, nn )
    integer , intent(in   ) :: ms, me


    !
    ! local variables
    !
    real(DP) :: tmp( ms:me )
    integer  :: i, j, m


    ! solve matrix

    do i = 1+1, nn
       do m = ms, me
          tmp( m ) = 0.0d0
       end do
       do j = 1, i-1
          do m = ms, me
             tmp( m ) = tmp( m ) + mat( m, i, j ) * vec( m, j )
          end do
       end do
       do m = ms, me
          vec( m, i ) = vec( m, i ) - tmp( m )
       end do
    end do


    do i = nn, 1, -1
       do m = ms, me
          tmp( m ) = 0.0d0
       end do
       do j = i+1, nn
          do m = ms, me
             tmp( m ) = tmp( m ) + mat( m, i, j ) * vec( m, j )
          end do
       end do
       do m = ms, me
          vec( m, i ) = ( vec( m, i ) - tmp( m ) ) / mat( m, i, i )
       end do
    end do



  end subroutine ludecomp_solve_simple_many

  !**************************************************************************************

end module ludecomp_module