module Dynamic_HEVI_LinArg_matrix
  ! 
  ! MATRIX を利用
  !

  !暗黙の型宣言禁止
  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  :: AL2(:,:)      !LU 分解の結果 L (2 次元配列)  
  real(8), allocatable  :: AU2(:,:)      !LU 分解の結果 U (2 次元配列)
  integer, allocatable  :: IP(:)         !部分ピボット交換の情報を格納

  !値を保管する.
  save N, M, NUD, NLD, NAL, NA
  save AL2, AU2, IP

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

!!!--------------------------------------------------------------------!!!
  subroutine Dynamic_HEVI_LinArg_matrix_LUfact(RegXMin, RegXMax, RegZMin, RegZMax, A, B, C)
    !
    !実 3 項行列の LU 分解(倍精度). 
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer, intent(in) :: RegXMin                ! X 方向の計算領域の下限
    integer, intent(in) :: RegXMax                ! X 方向の計算領域の上限
    integer, intent(in) :: RegZMin                ! Z 方向の計算領域の下限
    integer, intent(in) :: RegZMax                ! Z 方向の計算領域の上限
    integer, intent(in) :: A(RegZMin+1:RegZMax)   !係数行列の対角成分
    integer, intent(in) :: B(RegZMin+2:RegZMax)   !係数行列の上三角部分
    integer, intent(in) :: C(RegZMin+1:RegZMax-1) !係数行列の下三角部分
    real(8)             :: TX(N, M)               !定数/解行列 (ダミー)
    real(8), parameter  :: EPS = 1.0d-16          !特異性判定基準値
    integer             :: IOPT(2)                !制御コード
    real(8)             :: WK1(NLD+NUD+1,2)       !作業配列
    real(8)             :: WK2(N,2)               !作業配列
    integer             :: IER
    integer             :: NM

    !変数の初期化
    N   = RegZMax - RegZMin !係数行列/改行列の次数, 整合寸法
    M   = RegXMax - RegXMin !方程式の組数
    NUD = 1                 !係数行列の上三角部分の帯幅
    NLD = 1                 !係数行列の下三角部分の帯幅
    NAL = NLD               !LU 分解の結果 L の整合寸法
    NA  = NUD + NLD + 1
    NM  = N

    !配列の割り当て
    allocate( AL2(NAL, N), AU2(NA, N), IP(N) )

    !AU2 の定義
    AU2           = 0.0d0
    AU2(1, 2:N  ) = C
    AU2(2, 1:N  ) = A
    AU2(3, 1:N-1) = B

    !IOPT の定義
    IOPT(1)     = 2       !LU 分解だけを行う
    IOPT(2)     = 1       !行スケーリングを行わない. 計算の安定性を優先
    
    !解行列の計算. 
    call HDLGBM(AU2, N, NLD, NUD, NA, TX, M, NM, EPS, &
      &         IOPT, AL2, NAL, IP, WK1, WK2, IER )
    
    !コンティションチェック
!    if (IER /= 0) then 
!      write(*,*) "MATRIX code is eregular !!"
!      stop
!    end if
    
  end subroutine Dynamic_HEVI_LinArg_matrix_LUfact
  

!!!--------------------------------------------------------------------!!!
  subroutine Dynamic_HEVI_LinArg_matrix_LUsolv( X )
    !
    !  LU 分解された実 3 項行列の連立 1 次方程式(倍精度). MATRIX 利用
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(inout) :: X(M, N)              !定数/解行列
    real(8)                :: TX(N, M)             !定数/解行列の転置行列
    real(8), parameter     :: EPS = 1.0d-16        !特異性判定基準値
    integer                :: IOPT(2)              !制御コード
    real(8)                :: WK1(NLD+NUD+1,2)     !作業配列
    real(8)                :: WK2(N,2)             !作業配列
    integer                :: IER
    integer                :: NM

    !変数の初期化
    IOPT(1)     = 3       ! LU 分解の結果を利用し解を求める
    IOPT(2)     = 1       !行スケーリングを行わない. 計算の安定性を優先
    TX = transpose( X )
    NM = N

    !解行列の計算. 
    call HDLGBM(AU2, N, NLD, NUD, NA, TX, M, NM, EPS, &
      &         IOPT, AL2, NAL, IP, WK1, WK2, IER )

    !解の出力
    X = transpose( TX )

    !コンティションチェック
!    if (IER /= 0) then 
!      write(*,*) "MATRIX code is eregular !!"
!      stop
!    end if
   
  end subroutine Dynamic_HEVI_LinArg_matrix_LUsolv
  
end module Dynamic_HEVI_LinArg_matrix
