!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module LAPACK_Linear
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: lapack_linear.f90,v 1.2.2.2 2005/11/14 01:55:12 kitamo Exp $ 
!   * Tag Name: $Name: arare3m-20051114 $
!   * Change History: 
!
!== Overview 
!
!LAPACK 線形計算のためのパッケージ型モジュール
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!== Future Plans
!
!  * 三重対角行列の解法しかコーディングされていない
!
                                                                 !=end

module LAPACK_linear
  implicit none
  
  private
                                                                 !=begin
  !=== Public Interface
  public lapack_dgtsv_init  !初期化
  public lapack_dgtsv       !実 3 項行列の連立 1 次方程式(倍精度)
  public lapack_dgttrf      !実 3 項行列の LU 分解
  public lapack_dgttrs      !LU 分解済の実 3 項行列の連立 1 次方程式(倍精度)
                                                                 !=end
  integer                   :: N = 10       !値はダミー

  save N

contains
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Initialize module
  !
  !LAPACK の DGTSV ルーチンの初期化
  !
  subroutine lapack_dgtsv_init(M)

    !=== Input
    integer, intent(in) :: M  !配列の大きさ
                                                                 !=end
    N = M
  end subroutine lapack_dgtsv_init
  
  
                                                                 !=begin
  !== Procedure Interface 
  !
  !=== Execute Subroutine for the module
  !
  !実 3 項行列の連立 1 次方程式(倍精度)
  !LU 分解を行ってから方程式を解く
  ! 
  subroutine lapack_dgtsv(A, B, C, D)
    !=== Dependency
    use dc_message, only: MessageNotify
                                                                 !=end
    !=== 暗黙の型宣言禁止
    implicit none
                                                                 !=begin
    !=== In/Out
    real(8), intent(inout)         :: A(N)      !係数行列
    real(8), intent(inout)         :: B(N-1)    !係数行列
    real(8), intent(inout)         :: C(N-1)    !係数行列
    real(8), intent(inout)         :: D(1,N)    !定数/解行列
                                                                 !=end
    !=== Work
    integer                        :: NRHS
    integer                        :: INFO
    integer                        :: LDB

    !=== 変数の初期化
    NRHS = 0; INFO = 0; LDB = N
    
    !=== 解行列の計算. LAPACK を使用. 
    call DGTSV(N, NRHS, C, A, B, D, LDB, INFO)

    !=== 解のコンディションをチェック. 
    if (INFO /= 0) then
       call MessageNotify("Error", "lapack_linear", &
&                        "INFO is not 0")
       stop
    end if
    
  end subroutine lapack_dgtsv


  !== Procedure Interface 
  !
  !=== Execute Subroutine for the module
  !
  !実 3 項行列の連立 1 次方程式(倍精度)
  !
  subroutine lapack_dgttrf(A, B, C, B2, IPIV)
    !=== Dependency
    use dc_message, only: MessageNotify
                                                                 !=end
    !=== 暗黙の型宣言禁止
    implicit none
                                                                 !=begin
    !=== In/Out
    real(8), intent(inout)         :: A(N)      !係数行列
    real(8), intent(inout)         :: B(N-1)    !係数行列
    real(8), intent(inout)         :: C(N-1)    !係数行列
    real(8), intent(out)           :: B2(N-2)   !係数行列
    real(8), intent(out)           :: IPIV(N)   !部分ピボット交換の情報を格納
                                                                 !=end
    !=== Work
    integer                        :: INFO

    !=== 変数の初期化
    INFO = 0
    
    !=== 解行列の計算. LAPACK を使用. 
    call DGTTRF(N, C, A, B, B2, IPIV, INFO)

    !=== 解のコンディションをチェック. 
    if (INFO /= 0) then
       call MessageNotify("Error", "lapack_linear", &
&                        "INFO is not 0")
       stop
    end if

  end subroutine lapack_dgttrf
  
  
  !== Procedure Interface 
  !
  !=== Execute Subroutine for the module
  !
  !LU 分解されている実 3 項行列の連立 1 次方程式(倍精度)
  !
  !
  subroutine lapack_dgttrs(A, B, C, D, B2, IPIV)
    !=== Dependency
    use dc_message, only: MessageNotify
                                                                 !=end
    !=== 暗黙の型宣言禁止
    implicit none
                                                                 !=begin
    !=== In/Out
    real(8), intent(inout)         :: A(N)      !係数行列
    real(8), intent(inout)         :: B(N-1)    !係数行列
    real(8), intent(inout)         :: C(N-1)    !係数行列
    real(8), intent(inout)         :: D(1,N)    !定数/解行列
    real(8), intent(inout)         :: B2(N-2)    !係数行列
    real(8), intent(in)            :: IPIV(N)   !部分ピボット交換の情報を格納
                                                                 !=end
    !=== Work
    integer                        :: NRHS
    integer                        :: INFO
    integer                        :: LDB
    character(1),parameter         :: TRANS = 'N'

    !=== 変数の初期化
    NRHS = 0; INFO = 0; LDB = N
    
    !=== 解行列の計算. LAPACK を使用. 
    call DGTTRS(TRANS, N, NRHS, C, A, B, B2, IPIV, D, LDB, INFO)

    !=== 解のコンディションをチェック. 
    if (INFO /= 0) then
       call MessageNotify("Error", "lapack_linear", &
&                        "INFO is not 0")
       stop
    end if
    
  end subroutine lapack_dgttrs

end module LAPACK_linear