! -*- mode: f90; coding: utf-8 -*-
!-----------------------------------------------------------------------
!> 汎用ユーティリティモジュール
!>
!> これまで複数のモジュールに分離していた諸関数をひとまとまりにした.
!>
module spml_utils
  use dc_message, only : MessageNotify
  use dc_types, only: DP
  implicit none
  private
  public pi
  public gammaln
  public factrl
#ifdef EIG_LAPACK
  public eigen
  public ceigen
#endif

  !> 円周率 40 桁: 3.14159265358979323846264338327950288
  real(DP), parameter :: pi = 3.141592653589793238d0

contains
  !----------------------------------------------------------
  !> Returns the vaule \f$ ln[Γ(x)] for x >= 0 \f$
  function gammaln(xx)
    real(DP), intent(IN) :: xx
    real(DP)             :: gammaln
    real(DP)             :: dlgamma

    gammaln = dlgamma(xx)

  end function gammaln

  !----------------------------------------------------------
  !> Returns the value \f$n!\f$ as a floating point number
  function factrl(n)
    integer, intent(IN)  :: n
    real(DP)             :: factrl

    factrl = exp(gammaln(n+1.0d0))

  end function factrl

#ifdef EIG_LAPACK
  !---------------------------------------------------------------------
  !> 実行列の固有値計算: DGEEV/LAPACK の wrapper
  !>
  !>   * 行列 AMAT の i 番目固有値を eigen_r(i), eigen_i(i) に格納
  !>   * 対応する固有ベクトルを eigvec_r(:,i), eigvec_i(:,i) に格納
  !>   * 格納する固有値の数は引数 eigen_r の大きさで決まる
  !>   * 固有値の順番は sort と order で定められる.
  !>   * sort によって順番を定めるために用いる量を指定する.
  !>     実部(R), 実部の絶対値(RA), 虚部(I), 虚部の絶対値(IA)
  !>   * reverse によって小さい順(.false.), 大きい順(.true.)を指定できる.
  !>   * デフォルトは sort='R', reverse=.false.
  !>
  subroutine eigen(                             &
    & amat,eigen_r, eigen_i, eigvec_r,eigvec_i, &
    & info, sort, reverse )


    !> 入力正方行列
    real(DP), dimension(:,:)                  :: amat
    !> 固有値実数部
    real(DP), intent(out), dimension(:)       :: eigen_r
    !> 固有値虚数部
    real(DP), intent(out), dimension(:)       :: eigen_i
    !> 固有ベクトル実部
    real(DP), intent(out), &
      & dimension(size(amat,1),size(eigen_r))   :: eigvec_r
    !> 固有ベクトル虚部
    real(DP), intent(out), &
      & dimension(size(amat,1),size(eigen_i))   :: eigvec_i
    !> ステータス
    integer, intent(out)                      :: info
    !> 並び変えの量
    character(len=2), intent(in), optional    :: sort
    !> 並び変えスイッチ
    logical, intent(in), optional             :: reverse

    real(DP), dimension(size(amat,1))              :: wr
    real(DP), dimension(size(amat,1))              :: wi
    real(DP), dimension(size(amat,1),size(amat,1)) :: vl
    real(DP), dimension(size(amat,1),size(amat,1)) :: vr
    real(DP), dimension(size(amat,1)*4)            :: work
    integer, dimension(size(amat,1))               :: index
    character(len=1), parameter                    :: jobvl='N'
    character(len=1), parameter                    :: jobvr='V'
    integer :: nm, i, j

    !------- 形状チェック ------
    if (size(amat,1) /= size(amat,2))then
      call MessageNotify('E','DEIGEN_LAPACK','Input matrix not square')
    else
      nm = size(amat,1)
    endif

    !------- DGEEV/LAPACK による計算 ------
    call DGEEV( jobvl, jobvr, nm, amat, nm, wr, wi, &
      & vl, nm, vr, nm, work, nm*4, info )

    !------- サブルーチンエラー処理 -------
    if ( info /= 0 ) then
      call MessageNotify('W','DEIGEN_LAPACK',&
        'Error in calculating eigenvalues/vectors...',i=(/info/) )
      return
    endif

    !------- 固有ベクトル入れ換え -------
    if ( present(sort) ) then
      if ( sort == 'RA' ) then          ! 固有値実部の絶対値
        index=indexx(abs(wr))
      elseif ( trim(sort) == 'I ' ) then ! 固有値虚部
        index=indexx(wi)
      elseif ( sort == 'IA' ) then      ! 固有値虚部の絶対値
        index=indexx(abs(wi))
      else
        index=indexx(wr)               ! defaultは固有値実部
      endif
    else
      index=indexx(wr)                  ! defaultは固有値実部
    endif

    if ( present(reverse) )then
      if ( reverse ) then               ! 大きい順
        index=index(size(index):1:-1)
      endif
    endif

    do i=1,size(eigen_r)
      j = index(i)
      eigen_r(i) = wr(j)
      eigen_i(i) = wi(j)

      if ( abs(wi(j)) <= epsilon(0.0D0) ) then
        eigvec_r(:,i) = vr(:,j)
        eigvec_i(:,i) = 0.0D0

      elseif ( wi(j) > 0.0D0 ) then
        eigvec_r(:,i) = vr(:,j)
        eigvec_i(:,i) = vr(:,j+1)
      else
        eigvec_r(:,i) = vr(:,j-1)
        eigvec_i(:,i) = -vr(:,j)
      endif
    enddo

  end subroutine eigen

  !---------------------------------------------------------------------
  !> 複素行列の固有値計算: ZGEEV/LAPACK の wrapper
  !>
  !>  * 行列 CMAT の i 番目固有値を eigen(i) に格納
  !>  * 対応する固有ベクトルを eigvec(:,i) に格納
  !>  * 格納する固有値の数は引数 eigen の大きさで決まる
  !>  * 固有値の順番は sort と order で定められる.
  !>  * sort によって順番を定めるために用いる量を指定する.
  !>    実部(R), 実部の絶対値(RA), 虚部(I), 虚部の絶対値(IA)
  !>  * reverse によって小さい順(.false.), 大きい順(.true.)を指定できる.
  !>  * デフォルトは sort='R', reverse=.false.
  !>  * 行列 CMAT は保存されない.
  !>
  subroutine ceigen(cmat,eigen,eigvec,info,sort,reverse)
    !> 入力正方行列
    complex(DP), dimension(:,:)                :: cmat
    !> 固有値実数部
    complex(DP), intent(out), dimension(:)     :: eigen
    !> 固有ベクトル実部
    complex(DP), intent(out), &
      dimension(size(cmat,1),size(eigen))     :: eigvec
    !> ステータス
    integer, intent(out)                      :: info
    !> 並び変えの量
    character(len=2), intent(in), optional    :: sort
    !> 並び変えスイッチ
    logical, intent(in), optional             :: reverse

   !------------ 作業変数 ------------
    complex(DP), dimension(size(cmat,1))              :: ceig
    complex(DP), dimension(size(cmat,1),size(cmat,1)) :: vl
    complex(DP), dimension(size(cmat,1),size(cmat,1)) :: vr
    complex(DP), dimension(size(cmat,1)*2)            :: work
    real(DP), dimension(size(cmat,1)*2)               :: rwork
    integer, dimension(size(cmat,1))                  :: index

    character(len=1), parameter :: jobvl='N'
    character(len=1), parameter :: jobvr='V'
    integer                     :: nm, i, j

    !------- 形状チェック ------
    if (size(cmat,1) /= size(cmat,2))then
      call MessageNotify('E','DCEIGEN_LAPACK','Input matrix not square')
    else
      nm = size(cmat,1)
    endif

    !------- DGEEV/LAPACK による計算 ------
    call ZGEEV( jobvl, jobvr, nm, cmat, nm, ceig,   &
      &         vl, nm, vr, nm, work, nm*2, rwork, info )

    !------- サブルーチンエラー処理 -------
    if ( info /= 0 ) then
      call MessageNotify('W','DCEIGEN_LAPACK',&
        'Error in calculating eigenvalues/vectors...',i=(/info/) )
      return
    endif

    !------- 固有ベクトル入れ換え -------
    if ( present(sort) ) then
      if ( sort == 'RA' ) then          ! 固有値実部の絶対値
        index=indexx(abs(dble(ceig)))
      elseif ( trim(sort) == 'I' ) then ! 固有値虚部
        index=indexx(aimag(ceig))
      elseif ( sort == 'IA' ) then      ! 固有値虚部の絶対値
        index=indexx(abs(aimag(ceig)))
      else
        index=indexx(dble(ceig))       ! defaultは固有値実部
      endif
    else
      index=indexx(dble(ceig))          ! defaultは固有値実部
    endif

    if ( present(reverse) )then
      if ( reverse ) then               ! 大きい順
        index=index(size(index):1:-1)
      endif
    endif

    do i=1,size(eigen)
      j = index(i)
      eigen(i)    = ceig(j)
      eigvec(:,i) = vr(:,j)
    enddo

  end subroutine ceigen
#endif

  !---------------------------------------------------------------------
  !> 与えられた配列を値の小さい順に sort する
  ! *need rewrite!!*
  function indexx(arrin)
    implicit none
    real(DP), dimension(:), intent(in)  :: arrin
    integer, dimension(size(arrin))    :: indexx

    integer :: ir, l, i, j, indxt
    real(DP) :: Q

    integer :: n

    n = size(arrin)

    indexx = (/(j, j=1,n)/)

    l  = n/2 + 1
    ir = n

10  continue

    if( l .gt. 1 )then
      l     = l - 1
      indxt = indexx(l)
      q     = arrin(indxt)
    else
      indxt = indexx(ir)
      q     = arrin(indxt)
      indexx(ir) = indexx(1)
      ir = ir - 1
      if( ir .eq. 1 )then
        indexx(1) = indxt
        return
      endif
    endif
    i = l
    j = l + l

20  if( j .le. ir )then
      if( j .lt. ir )then
        if( arrin( indexx(j) ) .lt. arrin( indexx(j+1) )  ) j = j + 1
      endif
      if(  q .lt. arrin( indexx(j) )  )then
        indexx( i ) = indexx( j )
        i = J
        j = j + j
      else
        j = ir + 1
      endif
      go to 20
    endif
    indexx( i ) = indxt
    go to 10
  end function indexx

end module spml_utils
