  subroutine ludecomp21(alu,kp)

    ! ALU(NDIM,NDIM), KP(NDIM)
    ! NDIM x NDIM の行列 JDIM 個を一度に LU 分解.
    ! ＬＵ行列は 入力行列に上書きされる.

    real(8), intent(inout) :: alu(:,: )                 ! 入力／ＬＵ行列
    integer, intent(out)   :: kp(size(alu,1))           ! ピボット

    if ( size(alu,1) > size(alu,2) ) then
       call msgdmp('E','ludecomp',&
            'The third dimension is less than the second')
    elseif( size(alu,1) < size(alu,2) ) then
       call msgdmp('W','ludecomp',&
            'The third dimension is grater than the second')
    endif

    !" 行列のＬＵ分解（部分ピボット選択）
    call LUMAKE( alu, kp, 1, size(alu,1) )

  end subroutine ludecomp21

  subroutine ludecomp32(alu,kp)

    ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM)
    ! NDIM x NDIM の行列 JDIM 個を一度に LU 分解.
    ! ＬＵ行列は 入力行列に上書きされる.

    real(8), intent(inout) :: alu(:,:,:)                  ! 入力／ＬＵ行列
    integer, intent(out)   :: kp(size(alu,1),size(alu,2)) ! ピボット

    if ( size(alu,2) > size(alu,3) ) then
       call msgdmp('E','ludecomp',&
            'The third dimension is less than the second')
    elseif( size(alu,2) < size(alu,3) ) then
       call msgdmp('W','ludecomp',&
            'The third dimension is grater than the second')
    endif

    !" 行列のＬＵ分解（部分ピボット選択）
    call LUMAKE( alu, kp, size(alu,1), size(alu,2) )

  end subroutine ludecomp32

  function lusolve211(alu,kp,b)

    ! ALU(NDIM,NDIM), KP(NDIM), B(NDIM)
    ! NDIM x NDIM 型行列の連立方程式
    ! A X = B を IDIM 個の B に対して計算する. 
    ! 解は右辺の入力ベクトルに上書きされる

    real(8), intent(in)  :: alu(:,:)              ! 入力／ＬＵ行列
    integer, intent(in)  :: kp(:)                 ! ピボット
    real(8), intent(in)  :: b(:)                  ! 右辺ベクトル

    real(8) :: lusolve211(size(b))                   ! 解

    lusolve211 = b
    call LUSOLV( lusolve211, alu , kp, &             !" ＬＵ分解による解の計算
                 1, 1, size(b) )

  end function lusolve211

  function lusolve212(alu,kp,b)

    ! ALU(NDIM,NDIM), KP(NDIM), B(JDIM,NDIM)
    ! NDIM x NDIM 型行列の連立方程式
    ! A X = B を IDIM 個の B に対して計算する. 
    ! 解は右辺の入力ベクトルに上書きされる

    real(8), intent(in)  :: alu(:,:)              ! 入力／ＬＵ行列
    integer, intent(in)  :: kp(:)                 ! ピボット
    real(8), intent(in)  :: b(:,:)                ! 右辺ベクトル

    real(8) :: lusolve212(size(b,1),size(b,2))       ! 解

    lusolve212 = b
    call LUSOLV( lusolve212, alu , kp, &             !" ＬＵ分解による解の計算
                 size(b,1), 1, size(b,2) )

  end function lusolve212

  function lusolve322(alu,kp,b)

    ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(JDIM,NDIM)
    ! NDIM x NDIM 型行列を JDIM 個並べた連立方程式
    ! A X = B を B に対して計算する. 
    ! 解は右辺の入力ベクトルに上書きされる

    real(8), intent(in)  :: alu(:,:,:)                   ! 入力／ＬＵ行列
    integer, intent(in)  :: kp(:,:)                      ! ピボット
    real(8), intent(in)  :: b(:,:)                       ! 右辺ベクトル

    real(8) :: lusolve322(size(b,1),size(b,2))             ! 解

    lusolve322 = b
    call LUSOLV( lusolve322, alu , kp, &           !" ＬＵ分解による解の計算
                 1, size(b,1), size(b,2) )

  end function lusolve322

  function lusolve323(alu,kp,b)

    ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(IDIM,JDIM,NDIM)
    ! NDIM x NDIM 型行列を JDIM 個並べた連立方程式
    ! A X = B を IDIM 個の B に対して計算する. 
    ! 解は右辺の入力ベクトルに上書きされる

    real(8), intent(in)  :: alu(:,:,:)                   ! 入力／ＬＵ行列
    integer, intent(in)  :: kp(:,:)                      ! ピボット
    real(8), intent(in)  :: b(:,:,:)                     ! 右辺ベクトル

    real(8) :: lusolve323(size(b,1),size(b,2),size(b,3)) ! 解

    lusolve323 = b
    call LUSOLV( lusolve323, alu , kp, &           !" ＬＵ分解による解の計算
                 size(b,1), size(b,2), size(b,3) )

  end function lusolve323

module lumatrix
  private 
  public ludecomp, lusolve

  interface ludecomp
     subroutine ludecomp21(alu,kp)
       real(8), intent(inout) :: alu(:,: )                 ! 入力／ＬＵ行列
       integer, intent(out)   :: kp(size(alu,1))           ! ピボット
     end subroutine ludecomp21

     subroutine ludecomp32(alu,kp)
       real(8), intent(inout) :: alu(:,:,:)                  ! 入力／ＬＵ行列
       integer, intent(out)   :: kp(size(alu,1),size(alu,2)) ! ピボット
     end subroutine ludecomp32
  end interface

  interface lusolve
     function lusolve211(alu,kp,b)

      ! ALU(NDIM,NDIM), KP(NDIM), B(NDIM)
      ! NDIM x NDIM 型行列の連立方程式
      ! A X = B を IDIM 個の B に対して計算する. 
      ! 解は右辺の入力ベクトルに上書きされる

       real(8), intent(in)  :: alu(:,:)              ! 入力／ＬＵ行列
       integer, intent(in)  :: kp(:)                 ! ピボット
       real(8), intent(in)  :: b(:)                  ! 右辺ベクトル
       real(8) :: lusolve211(size(b))                 ! 解
     end function lusolve211

     function lusolve212(alu,kp,b)

       ! ALU(NDIM,NDIM), KP(NDIM), B(JDIM,NDIM)
       ! NDIM x NDIM 型行列の連立方程式
       ! A X = B を IDIM 個の B に対して計算する. 
       ! 解は右辺の入力ベクトルに上書きされる

       real(8), intent(in)  :: alu(:,:)              ! 入力／ＬＵ行列
       integer, intent(in)  :: kp(:)                 ! ピボット
       real(8), intent(in)  :: b(:,:)                ! 右辺ベクトル

       real(8) :: lusolve212(size(b,1),size(b,2))       ! 解

     end function lusolve212

     function lusolve322(alu,kp,b)

       ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(JDIM,NDIM)
       ! NDIM x NDIM 型行列を JDIM 個並べた連立方程式
       ! A X = B を B に対して計算する. 
       ! 解は右辺の入力ベクトルに上書きされる

       real(8), intent(in)  :: alu(:,:,:)                   ! 入力／ＬＵ行列
       integer, intent(in)  :: kp(:,:)                      ! ピボット
       real(8), intent(in)  :: b(:,:)                       ! 右辺ベクトル

       real(8) :: lusolve322(size(b,1),size(b,2))             ! 解
     end function lusolve322

     function lusolve323(alu,kp,b)

       ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(IDIM,JDIM,NDIM)
       ! NDIM x NDIM 型行列を JDIM 個並べた連立方程式
       ! A X = B を IDIM 個の B に対して計算する. 
       ! 解は右辺の入力ベクトルに上書きされる

       real(8), intent(in)  :: alu(:,:,:)                   ! 入力／ＬＵ行列
       integer, intent(in)  :: kp(:,:)                      ! ピボット
       real(8), intent(in)  :: b(:,:,:)                     ! 右辺ベクトル

       real(8) :: lusolve323(size(b,1),size(b,2),size(b,3)) ! 解

     end function lusolve323
  end interface

end module lumatrix
