!--
!------------------------------------------------------------------------
! Copyright(c) 2002-2013 SPMODEL Development Group. All rights reserved.
!------------------------------------------------------------------------
!表題  lumatrix : 行列の LU 分解による線形連立方程式の解
!
!      spml/lumatrix モジュールは, LU 分解法により連立 1 次方程式を解くための
!      Fortran90 関数を提供する.
!
!      他のスペクトル計算用モジュールの中で登場する境界値問題を解くために
!      用いられている.
!
!      並列計算機を意識して, 同じ次数の複数個の連立1 次方程式
!
!          A[ij]^(n) X [j]^(n) = B[i]^(n)
!
!      の解を同時に複数個の右辺ベクトル B[i]^(n)b に対して求めることが
!      できるようになっている.
!
!
!履歴  2002/01/20  竹広真一
!      2002/06/10  竹広真一  ベクトル長問題回避のため lusol2 を使用
!      2005/01/10  竹広真一  msgdmp -> MessageNotify に変更
!      2006/03/04  竹広真一  コメントを RDoc 用に修正
!      2009/01/29  佐々木洋平  コメントを RDoc 用に修正
!      2009/08/06  竹広真一    ludecomp21 用ルーチン OMP 版変更
!      2010/03/17  竹広真一    Lapack 版作成
!      2013/12/14  佐々木洋平 interface 部を修正
!
!++
module lumatrix
  !
  != lumatrix
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/lumatrix モジュールは, LU 分解法により連立 1 次方程式を解くため
  ! のFortran90 関数を提供する.
  !
  ! 他のスペクトル計算用モジュールの中で登場する境界値問題を解くために
  ! 用いられている.
  !
  ! ベクトル計算機を意識して, 同じ次数の複数個の連立1 次方程式
  !
  !     A[ij]^(n) X [j]^(n) = B[i]^(n)
  !
  ! の解を同時に複数個の右辺ベクトル B[i]^(n)b に対して求めることが
  ! できるようになっている.
  !
  !== 変数・手続き群の要約
  !
  ! LUDecomp    :: 行列の LU 分解を行う
  ! LUSolve     :: 連立 1 次方程式の解を求める
  !
  private
  public LUDecomp, LUSolve

  interface LUDecomp
    !
    !=== 与えられた行列の LU 分解を行い, ピボットを格納する.
    !
    ! * LU 分解された結果は alu に上書きされる.
    !   そのピボット情報は kp に格納される.
    !
    ! * LUSolve を用いる前にこのサブルーチンを呼んで alu と kp を
    !   計算しておく.
    !
    ! * 入力する行列とピボット用配列の次元によって内部でサブルーチンを
    !   使い分けている. ユーザーインターフェースは共通であるので
    !   下部ルーチンである ludecomp21, ludecomp32 を呼ぶ必要はない.
    !
    !=== 引数と結果の型
    !
    ! * Alu が 2 次元配列(与える係数行列が 1 つ)の場合
    !
    !     ! ALU(NDIM,NDIM), KP(NDIM)
    !     ! NDIM x NDIM の行列を LU 分解.
    !     ! LU 行列は 入力行列に上書きされる.
    !
    !     real(8), intent(inout) :: alu(:,: )         ! 入力／LU 行列
    !     integer, intent(out)   :: kp(size(alu,1))   ! ピボット
    !
    !
    ! * Alu が 3 次元配列(与える係数行列が複数)の場合
    !
    !     ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM)
    !     ! NDIM x NDIM の行列 JDIM 個を一度に LU 分解.
    !     ! LU 行列は 入力行列に上書きされる.
    !
    !     real(8), intent(inout) :: alu(:,:,:)      ! 入力／LU 行
    !     integer, intent(out)   :: kp(size(alu,1),size(alu,2)) ! ピボット
    !
    module procedure ludecomp21, ludecomp32
  end interface LUDecomp

  interface LUSolve
    !
    ! 連立 1 次方程式の解を求める
    !
    !  * LUSolve を用いる前に LUDecompを呼んで Alu を LU 分解し,
    !    ピボット情報 kp を計算しておかねばならない.
    !
    !  * 入力する行列とピボット用配列の次元によって内部でサブルーチンを
    !    使い分けている. ユーザーインターフェースは共通であるので
    !    下部ルーチンである lusolve??? を呼ぶ必要はない.
    !
    ! 引数と結果の型
    !
    !  *  Alu が 2 次元配列(与える係数行列が 1 つ),
    !     b が 1 次元配列(与える右辺ベクトルが 1 つ)の場合
    !
    !     ! ALU(NDIM,NDIM), KP(NDIM), B(NDIM)
    !     ! NDIM x NDIM 型行列の連立方程式
    !     ! A X = B を 1 個の B に対して計算する.
    !
    !     real(8), intent(in)  :: alu(:,:)              ! 入力／LU 行列
    !     integer, intent(in)  :: kp(:)                 ! ピボット
    !     real(8), intent(in)  :: b(:)                  ! 右辺ベクトル
    !
    !     real(8) :: lusolve(size(b))                   ! 解
    !
    !  * Alu が 2 次元配列(与える係数行列が 1 つ),
    !    b が 2 次元配列(与える右辺ベクトルが複数)の場合
    !
    !     ! ALU(NDIM,NDIM), KP(NDIM), B(JDIM,NDIM)
    !     ! NDIM x NDIM 型行列の連立方程式
    !     ! A X = B を JDIM 個の B に対して計算する.
    !
    !     real(8), intent(in)  :: alu(:,:)              ! 入力／LU 行列
    !     integer, intent(in)  :: kp(:)                 ! ピボット
    !     real(8), intent(in)  :: b(:,:)                ! 右辺ベクトル
    !
    !     real(8) :: lusolve(size(b,1),size(b,2))       ! 解
    !
    !
    !  * Alu が 3 次元配列(与える係数行列が複数),
    !    b が 2 次元配列(与える右辺ベクトルが 1 つ)の場合
    !
    !     ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(JDIM,NDIM)
    !     ! NDIM x NDIM 型行列を JDIM 個並べた連立方程式
    !     ! A X = B をひとつの B の並びに対して計算する.
    !
    !     real(8), intent(in)  :: alu(:,:,:)            ! 入力／LU 行列
    !     integer, intent(in)  :: kp(:,:)               ! ピボット
    !     real(8), intent(in)  :: b(:,:)                ! 右辺ベクトル
    !
    !     real(8) :: lusolve(size(b,1),size(b,2))             ! 解
    !
    !
    !  * Alu が 3 次元配列(与える係数行列が複数),
    !    b が 3 次元配列(与える右辺ベクトルが複数)の場合
    !
    !     ! 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(:,:,:)                ! 入力／LU 行列
    !     integer, intent(in)  :: kp(:,:)                   ! ピボット
    !     real(8), intent(in)  :: b(:,:,:)                  ! 右辺ベクトル
    !
    !     real(8) :: lusolve(size(b,1),size(b,2),size(b,3)) ! 解
    !
    !
    module procedure lusolve211, lusolve212, lusolve322, lusolve323
  end interface LUSolve

contains
  subroutine ludecomp21(alu,kp)
    !
    ! ALU(NDIM,NDIM), KP(NDIM)
    ! NDIM x NDIM の行列を LU 分解.
    ! LU 行列は 入力行列に上書きされる.
    !
    use dc_message

    real(8), intent(inout) :: alu(:,:)                  ! 入力／LU 行列
    integer, intent(out)   :: kp(size(alu,2))           ! ピボット
    integer  :: icon
    integer  :: kk, nn

    kk = size(alu,1) ;  nn = size(alu,2)

    if ( kk /= nn ) then
      call MessageNotify('E','ludecomp21',&
        'The first dimension differs from the second')
    endif

    !" 行列のLU 分解
    call dgetrf( nn, nn, alu, nn, kp, icon )
    if ( icon /= 0 ) &
      call MessageNotify('E','ludecomp21','LU decompostion not succeeded.')

  end subroutine ludecomp21

  subroutine ludecomp32(alu,kp)
    !
    ! ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM)
    ! NDIM x NDIM の行列 JDIM 個を一度に LU 分解.
    ! LU 行列は 入力行列に上書きされる.
    !
    use dc_message

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

    integer :: icon
    integer :: jj, kk, nn
    integer :: j

    jj = size(alu,1) ; kk = size(alu,2) ;  nn = size(alu,3)

    if ( kk /= nn ) then
      call MessageNotify('E','ludecomp32',&
        'The second dimension differs from the third')
    endif

    !" 行列のLU 分解（クラウト法）
    !$omp parallel do private(icon)
    do j=1,jj
      call dgetrf( nn, nn, alu(j,:,:), nn, kp(j,:), icon )
      if ( icon /= 0 ) &
        call MessageNotify('E','ludecomp32','LU decompostion not succeeded.')
    end do
    !$omp end parallel do

  end subroutine ludecomp32

  function lusolve211(alu,kp,b)
    !
    ! ALU(NDIM,NDIM), KP(NDIM), B(NDIM)
    ! NDIM x NDIM 型行列の連立方程式
    ! A X = B を 1 個の B に対して計算する.
    !
    use dc_message

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

    real(8) :: lusolve211(size(b))                ! 解
    integer :: kk, nn
    integer :: icon

    kk = size(alu,1) ;  nn = size(alu,2)

    if ( kk /= nn ) then
      call MessageNotify('E','ludsolve211',&
        'The first dimension is greater than the second')
    elseif ( nn /= size(b)) then
      call MessageNotify('E','lusolve211', &
        'The input vector length differs from the second dimension of the matrix')
    endif

    lusolve211 = b
    call dgetrs( 'N', nn, 1, alu, nn, kp, lusolve211, nn, icon )
    if ( icon /= 0 ) &
      call MessageNotify('E','lusplve211','LU decompostion not succeeded.')

  end function lusolve211

  function lusolve212(alu,kp,b)
    !
    ! ALU(NDIM,NDIM), KP(NDIM), B(JDIM,NDIM)
    ! NDIM x NDIM 型行列の連立方程式
    ! A X = B を JDIM 個の B に対して計算する.
    !
    use dc_message

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

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

    integer :: kk, nn
    integer :: j
    integer :: icon

    kk = size(alu,1) ;  nn = size(alu,2)

    if ( kk /= nn ) then
      call MessageNotify('E','ludsolve212',&
        'The first dimension differs from the second')
    elseif ( nn /= size(b,2)) then
      call MessageNotify('E','lusolvep212', &
        'The vector length differs from the second dimension of the matrix')
    endif

    lusolve212 = b

    !$omp parallel do private(icon)
    do j=1,size(b,1)
      call dgetrs( 'N', nn, 1, alu, nn, kp, lusolve212(j,:), nn, icon )
      if ( icon /= 0 ) &
        call MessageNotify('E','lusplve212','LU decompostion not succeeded.')
    enddo
    !$omp end parallel do

  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 の並びに対して計算する.
    !
    use dc_message

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

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

    integer  :: jj, kk, nn
    integer  :: j
    integer :: icon

    jj = size(alu,1) ; kk = size(alu,2) ;  nn = size(alu,3)

    if ( kk /= nn ) then
      call MessageNotify('E','ludsolve322',&
        'The second dimension differs from the third')
    elseif ( nn /= size(b,2)) then
      call MessageNotify('E','lusolvep322', &
        'The vector length differs from the second dimension of the matrix')
    endif

    lusolve322 = b

    !$omp parallel do private(icon)
    do j=1,jj
      call dgetrs( 'N', nn, 1, alu(j,:,:), nn, kp(j,:), lusolve322(j,:), nn, icon )
      if ( icon /= 0 ) &
        call MessageNotify('E','lusplve322','LU decompostion not succeeded.')
    end do
    !$omp end parallel do

  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 に対して計算する.
    !
    use dc_message

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

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

    integer  :: jj, kk, nn
    integer  :: i, j
    integer :: icon

    jj = size(alu,1) ; kk = size(alu,2) ;  nn = size(alu,3)


    if ( kk /= nn ) then
      call MessageNotify('E','ludsolve323',&
        'The first dimension differs from the second')
    elseif ( nn /= size(b,3)) then
      call MessageNotify('E','lusolvep323', &
        'The vector length differs from the second dimension of the matrix')
    endif

    lusolve323 = b

    !$omp parallel do private(icon)
    do i=1,size(b,1)
      do j=1,jj
        call dgetrs( 'N', nn, 1, alu(j,:,:), nn, kp(j,:), lusolve323(i,j,:), nn, icon )
        if ( icon /= 0 ) &
          call MessageNotify('E','lusplve323','LU decompostion not succeeded.')
      enddo
    enddo
    !$omp end parallel do

  end function lusolve323

end module lumatrix