! -*- mode: f90; coding: utf-8 -*-
!------------------------------------------------------------------------
! Copyright(c) 2002-2019 SPMODEL Development Group. All rights reserved.
!
! 履歴  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 版変更
!       2013/12/14  佐々木洋平 interface 部を修正
!       2019/09/20  佐々木洋平 コメントを doxygen 用に修正, lapackと統合
!------------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2002-2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief 行列の LU 分解による線形連立方程式の解法: LAPACK 版
!> @details
!> spml/lumatrix モジュールは, LU 分解法により連立 1 次方程式を解くための
!> Fortran90 関数を提供する.
!>
!> 他のスペクトル計算用モジュールの中で登場する境界値問題を解くために
!> 用いられている.
!>
!> ベクトル計算機を意識して, 同じ次数の複数個の連立 1 次方程式
!> \f[
!>     A[ij] X [j] = B[i]
!> \f]
!> の解を同時に複数個の右辺ベクトル \f$B[i]\f$ に対して求めることが
!> できるようになっている.
module lumatrix
  use dc_types, only: DP
  use dc_message, only: MessageNotify
  private
  public LUDecomp
  public LUSolve
  !------------------------------------------------------------------------
  !> 与えられた行列の LU 分解を行い, ピボットを格納する.
  !>
  !> * LU 分解された結果は `Alu` に上書きされる.
  !>   そのピボット情報は `kp` に格納される.
  !>
  !> * LUSolve を用いる前にこのサブルーチンを呼んで `Alu` と `kp` を
  !>   計算しておく.
  interface LUDecomp
     module procedure ludecomp21
     module procedure ludecomp32
  end interface LUDecomp
  !-----------------------------------------------------------------------
  !> 連立 1 次方程式の解を求める
  !>
  !> @warning
  !> LUSolve を用いる前に LUDecompを呼んで Alu を LU 分解し,
  !> ピボット情報 kp を計算しておかねばならない.
  interface LUSolve
     module procedure lusolve211
     module procedure lusolve212
     module procedure lusolve322
     module procedure lusolve323
  end interface LUSolve
contains
  !-----------------------------------------------------------------------
  !> NDIM x NDIM の行列 ALU を LU 分解する.
  !> @note
  !> 結果は 入力行列に上書き
  subroutine ludecomp21(alu,kp)
    !> 入出力行列
    real(DP), intent(inout) :: alu(:,:)
    !> ピボット
    integer, intent(out)   :: kp(size(alu,1))
    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
    call dgetrf( nn, nn, alu, nn, kp, icon )
    if ( icon /= 0 ) &
      call MessageNotify('E','ludecomp21','LU decompostion not succeeded.')
  end subroutine ludecomp21
  !-----------------------------------------------------------------------
  !> JDIM 個 の NDIM x NDIM の行列を一度に LU 分解.
  !> @note
  !> 結果は入力行列に上書き
  subroutine ludecomp32(alu,kp)
    !> 入力行列．JDIM が行列の個数, 後ろ二つの NDIMxNDIM が行列の形状
    real(DP), intent(inout) :: alu(:,:,:)
    !> ピボット
    integer, intent(out)    :: kp(size(alu,1),size(alu,2))
    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 decompostion（クラウト法）
    !$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
  !-----------------------------------------------------------------------
  !> ALU(NDIM,NDIM), KP(NDIM), B(NDIM)
  !> NDIM x NDIM 型行列の連立方程式
  !> A X = B を 1 個の B に対して計算する.
  function lusolve211(alu,kp,b)
    !> 入力／LU 行列
    real(DP), intent(in)  :: alu(:,:)
    !> ピボット
    integer, intent(in)  :: kp(:)
    !> 右辺ベクトル
    real(DP), intent(in)  :: b(:)
    !> 解
    real(DP) :: 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 ) then
      call MessageNotify('E','lusplve211',&
        & 'LU decompostion not succeeded.')
    endif
  end function lusolve211
  !-----------------------------------------------------------------------
  !> ALU(NDIM,NDIM), KP(NDIM), B(JDIM,NDIM)
  !> NDIM x NDIM 型行列の連立方程式
  !> A X = B を JDIM 個の B に対して計算する.
  function lusolve212(alu,kp,b)
    !> 入力／LU 行列
    real(DP), intent(in)  :: alu(:,:)
    !> ピボット
    integer,  intent(in)  :: kp(:)
    !> 右辺ベクトル
    real(DP), intent(in)  :: b(:,:)
    !> 解
    real(DP) :: 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
  !-----------------------------------------------------------------------
  !> ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(JDIM,NDIM)
  !> NDIM x NDIM 型行列を JDIM 個並べた連立方程式
  !> A X = B をひとつの B の並びに対して計算する.
  function lusolve322(alu,kp,b)
    !> 入力／LU 行列
    real(DP), intent(in)  :: alu(:,:,:)
    !> ピボット
    integer,  intent(in)  :: kp(:,:)
    !> 右辺ベクトル
    real(DP), intent(in)  :: b(:,:)
    !> 解
    real(DP) :: 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
  !-----------------------------------------------------------------------
  !> ALU(JDIM,NDIM,NDIM), KP(JDIM,NDIM), B(IDIM,JDIM,NDIM)
  !> NDIM x NDIM 型行列を JDIM 個並べた連立方程式
  !> A X = B を IDIM 個の B に対して計算する.
  function lusolve323(alu,kp,b)
    !> 入力／LU 行列
    real(DP), intent(in)  :: alu(:,:,:)
    !> ピボット
    integer, intent(in)  :: kp(:,:)
    !> 右辺ベクトル
    real(DP), intent(in)  :: b(:,:,:)
    !> 解
    real(DP) :: 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
