module Dynamic_HEVI_LinArg_lapack
  ! 
  ! LAPACK を利用
  !

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private
  
  !変数の定義
  integer               :: N = 128           !係数行列/改行列の次数, 整合寸法
  integer               :: M = 256           !方程式の組数
  integer               :: NUD = 1           !係数行列の上三角部分の帯幅
  integer               :: NLD = 1           !係数行列の下三角部分の帯幅
  integer               :: NAL = 1           !LU 分解の結果 L の整合寸法
  integer               :: NA = 3            !NUD + NLD + 1
  real(8), allocatable  :: D(:)              !係数行列の対角成分
  real(8), allocatable  :: DL(:)             !係数行列の下三角部分
  real(8), allocatable  :: DU(:)             !係数行列の上三角部分
  real(8), allocatable  :: DU2(:)            !LU 分解の結果 L (1 次元配列)
  integer, allocatable  :: IPIV(:)           !部分ピボット交換の情報を格納

  !値を保管する.
  save N, M, NUD, NLD, NAL, NA
  save D, DL, DU, DU2, IPIV

  !総称名を定義して, 総称名を public にする.
  interface LinArg_LUfact
     module procedure Dynamic_HEVI_LinArg_lapack_LUfact
  end interface
  interface LinArg_LUsolv
     module procedure Dynamic_HEVI_LinArg_lapack_LUsolv
  end interface
  public LinArg_LUfact
  public LinArg_LUsolv

contains

!!!--------------------------------------------------------------------!!!
  subroutine Dynamic_HEVI_LinArg_lapack_LUfact(RegXMin, RegXMax, RegZMin, RegZMax, A, B, C)
    !
    ! 実 3 項行列の LU 分解(倍精度). 
    ! LAPACK DGTTRF ルーチンを利用. 
    !
    !  DGTTRF computes an LU factorization of a real tridiagonal matrix A
    !  using elimination with partial pivoting and row interchanges.
    !
    !  The factorization has the form
    !     A = L * U
    !  where L is a product of permutation and unit lower bidiagonal
    !  matrices and U is upper triangular with nonzeros in only the main
    !  diagonal and first two superdiagonals.
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer, intent(in) :: RegXMin                 ! X 方向の計算領域の下限
    integer, intent(in) :: RegXMax                 ! X 方向の計算領域の上限
    integer, intent(in) :: RegZMin                 ! Z 方向の計算領域の下限
    integer, intent(in) :: RegZMax                 ! Z 方向の計算領域の上限
    real(8), intent(in) :: A(RegZMin+1:RegZMax)    !係数行列の対角成分 (作業配列)
    real(8), intent(in) :: B(RegZMin+2:RegZMax)    !係数行列の上三角部分 (作業配列)
    real(8), intent(in) :: C(RegZMin+1:RegZMax-1)  !係数行列の下三角部分 (作業配列)
    integer             :: INFO                    !解のコンディションチェック
    
    !変数の初期化
    N   = RegZMax - RegZMin !係数行列/改行列の次数, 整合寸法
    M   = RegXMax - RegXMin !方程式の組数
    NUD = 1                 !係数行列の上三角部分の帯幅
    NLD = 1                 !係数行列の下三角部分の帯幅
    NAL = NLD               !LU 分解の結果 L の整合寸法
    NA  = NUD + NLD + 1
    INFO = 0

    !配列の割り当て
    allocate( D(N), DL(N-1), DU(N-1), DU2(N-2), IPIV(N) )

    !A, B, C の定義
    D  = A
    DU = B
    DL = C
        
    !解行列の計算. LAPACK を使用. 
    call DGTTRF(N, DL, D, DU, DU2, IPIV, INFO)
    
!    !解のコンディションをチェック. 
!    if (INFO /= 0) then
!      call MessageNotify("Error", "lapack_linear", "INFO is not 0")
!      stop
!    end if
    
  end subroutine Dynamic_HEVI_LinArg_lapack_LUfact
  

!!!--------------------------------------------------------------------!!!
  subroutine Dynamic_HEVI_LinArg_lapack_LUsolv( X )
    !
    ! LU 分解された実 3 項行列の連立 1 次方程式(倍精度). 
    ! LAPACK DGTTRS ルーチンを利用. 
    !
    ! DGTTRS solves one of the systems of equations
    !    A*X = B  or  A'*X = B,
    ! with a tridiagonal matrix A using the LU factorization computed
    ! by DGTTRF.
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(inout) :: X(M, N)              !定数/解行列
    real(8)                :: TX(N, M)             !定数/解行列の転置行列
    integer                :: NRHS        
    integer                :: INFO
    integer                :: LDB
    character(1),parameter :: TRANS = 'N'

    !変数の初期化
    NRHS = M
    INFO = 0
    LDB  = N
    TX = transpose( X )
    
    !解行列の計算. LAPACK を使用. 
    call DGTTRS(TRANS, N, NRHS, DL, D, DU, DU2, IPIV, TX, LDB, INFO)

    !解の出力
    X = transpose( TX )
    
    !解のコンディションをチェック. 
!    if (INFO /= 0) then
!      call MessageNotify("Error", "lapack_linear", "INFO is not 0")
!      stop
!    end if
     
  end subroutine Dynamic_HEVI_LinArg_lapack_LUsolv
  
end module Dynamic_HEVI_LinArg_lapack
