!c Description:
!c  SSL2 線形計算のためのパッケージ型モジュール
!c
!c Current Code Owner:
!c   sugiyama@gfd-dennou.org
!c
!c Histry:
!c   Version    Date          Comment
!c   -------    ----------    --------
!c   1.0        2004-04-28    作成
!c
!c Copyright (C) SUGIYAMA Ko-ichiro, 2003, All rights reserved

module ssl2_linear
  implicit none

  private
  public ssl2_lax_init       !初期化
  public ssl2_lax            !実行列の連立 1 次方程式(倍精度)
  public ssl2_lsx_init       !初期化
  public ssl2_lsx            !実行列の連立 1 次方程式(倍精度)
  public ssl2_lsix_init      !初期化
  public ssl2_lsix           !実行列の連立 1 次方程式(倍精度)

  real(8), allocatable      :: VW(:)
  real(8)                   :: EPSZ = 0.0d0 ! 標準値
  integer, allocatable      :: IP(:)
  integer, allocatable      :: IVW(:)
  integer                   :: ISW = 1
  integer                   :: K = 10       !値はダミー
  integer                   :: N = 10       !値はダミー

  save K, N, ISW, IP, EPSZ, VW, IVW

contains

  subroutine ssl2_lax_init(M)
    integer, intent(in) :: M

    if (allocated(IP)) deallocate(IP, VW)
    allocate(IP(M), VW(M))
    N = M
    K = M
  end subroutine ssl2_lax_init

  subroutine ssl2_lax(A, B, sw)
    real(8), intent(inout)         :: A(:,:)         !係数行列
    real(8), intent(inout)         :: B(:)           !定数/解行列
    logical, intent(out), optional :: sw
    integer                        :: ICON
    integer                        :: IS

    IS = 0; VW = 0.0d0; IP = 0; ICON = 0
    if (present(sw)) sw = .false.

    !配列の大きさチェック
    if (N /= size(B,1) .OR. N /= size(A,1) .OR. N /= size(A,2)) then
       write(*,*) "ssl2_lax; matrix size is wrong"
       stop
    end if

    !解行列の計算. SSL II を使用.
    call DLAX(A, K, N, B, EPSZ, ISW, IS, VW, IP, ICON)

    !解のコンディションをチェック.
    if (ICON == 30000) then
       write(*,*) 'ssl2_lax: ICON is 30000'
       if (present(sw)) then
          sw = .true.
       else
          stop
       end if
    elseif (ICON == 20000) then
       write(*,*) 'ssl2_lax: ICON is 20000'
       write(*,*) "matrix is singular!"
       if (present(sw)) then
          sw = .true.
       else
          stop
       end if
    end if
  end subroutine ssl2_lax

  subroutine ssl2_lsx_init(M)
    integer, intent(in) :: M

    N = M
  end subroutine ssl2_lsx_init

  subroutine ssl2_lsx(A, B, sw)
    real(8), intent(inout)         :: A(:)          !係数行列
    real(8), intent(inout)         :: B(:)          !定数/解行列
    logical, intent(out), optional :: sw
    integer                        :: ICON

    ICON = 0
    if (present(sw)) sw = .false.

    !配列の大きさチェック
    if (N /= size(B,1) .OR. N*(N+1)/2 /= size(A,1)) then
       write(*,*) "ssl2_lsx; matrix size is wrong"
    end if

    !解行列の計算. SSL II を使用.
    call DLSX(A, N, B, EPSZ, ISW, ICON)

    !解のコンディションをチェック.
    if (ICON == 30000) then
       write(*,*) 'ssl2_lsx: ICON is 30000'
       if (present(sw)) then
          sw = .true.
       else
          stop
       end if
    elseif (ICON == 20000) then
       write(*,*) 'ssl2_lsx: ICON is 20000'
       write(*,*) "matrix is singular!"
       if (present(sw)) then
          sw = .true.
       else
          stop
       end if
    end if
  end subroutine ssl2_lsx

  subroutine ssl2_lsix_init(M)
    integer, intent(in) :: M

    if (allocated(IP)) deallocate(IP, VW, IVW)
    allocate(IP(M), VW(2*M), IVW(M))
    N = M
    K = M
  end subroutine ssl2_lsix_init

  subroutine ssl2_lsix(A, B, sw)
    real(8), intent(inout)         :: A(:)           !係数行列
    real(8), intent(inout)         :: B(:)           !定数/解行列
    logical, intent(out), optional :: sw
    integer                        :: ICON

    ICON = 0; VW = 0.0d0; IP = 0; IVW = 0
    if (present(sw)) sw = .false.

    !配列の大きさチェック
    if (N /= size(B,1) .OR. N*(N+1)/2 /= size(A,1)) then
       write(*,*) "ssl2_lsx; matrix size is wrong"
    end if

    !解行列の計算. SSL II を使用.
    call DLSIX(A, N, B, EPSZ, ISW, VW, IP, IVW, ICON)

    !解のコンディションをチェック.
    if (ICON == 30000) then
       write(*,*) 'ssl2_lsix: ICON is 30000'
       if (present(sw)) then
          sw = .true.
       else
          stop
       end if
    elseif (ICON == 20000) then
       write(*,*) 'ssl2_lsix: ICON is 20000'
       write(*,*) "matrix is singular!"
       if (present(sw)) then
          sw = .true.
       else
          stop
       end if
    end if

  end subroutine ssl2_lsix

end module ssl2_linear