!--
!----------------------------------------------------------------------
!     Copyright (c) 2011 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  複素固有値問題サブルーチン (SSL2)
!
!履歴  2011/11/05  竹広真一
!
!++
module ssl2_ceigen
  !
  != ssl2_ceigen
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id: 
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  ! 
  !== 概要
  !
  ! spml/ssl2_ceigen は eigmatrix モジュールの公開サブルーチン ceigen
  ! 用の固有値計算用共通インターフェースを与える.
  !
  !   * 行列 CMAT の i 番目固有値を eigen(i) に格納
  !   * 対応する固有ベクトルを eigvec(:,i) に格納
  !   * 格納する固有値の数は引数 eigen の大きさで決まる
  !   * 固有値の順番は sort と order で定められる. 
  !   * sort によって順番を定めるために用いる量を指定する. 
  !     実部(R), 実部の絶対値(RA), 虚部(I), 虚部の絶対値(IA)
  !   * reverse によって小さい順(.false.), 大きい順(.true.)を指定できる.
  !   * デフォルトは sort='R', reverse=.false.
  !   * 行列 CMAT は保存されない.
  !
  ! 内部では DCEIG2/SSL2 ルーチンによる複素行列の固有値/固有ベクトル計算を
  ! 行っている. が, ユーザーは用いているライブラリとサブルーチンを意識
  ! することなく使うことができる. 
  use dc_message, only : MessageNotify

  implicit none
  private
  public dceig2_ssl2

contains
  subroutine dceig2_ssl2(cmat,eigen,eigvec,info,sort,reverse)
    !
    ! このサブルーチンは固有値計算用共通インターフェースを与える
    ! eigmatrix モジュールの公開サブルーチン ceigen として用いられる. 
    !
    !   * 行列 CMAT の i 番目固有値を eigen(i) に格納
    !   * 対応する固有ベクトルを eigvec(:,i) に格納
    !   * 格納する固有値の数は引数 eigen の大きさで決まる
    !   * 固有値の順番は sort と order で定められる. 
    !   * sort によって順番を定めるために用いる量を指定する. 
    !     実部(R), 実部の絶対値(RA), 虚部(I), 虚部の絶対値(IA)
    !   * reverse によって小さい順(.false.), 大きい順(.true.)を指定できる.
    !   * デフォルトは sort='R', reverse=.false.
    !   * 行列 CMAT は保存されない.
    !
    ! 内部では DCEIG2/SSL2 ルーチンによる複素行列の固有値/固有ベクトル計算を
    ! 行っている. が, ユーザーは用いているライブラリとサブルーチンを意識
    ! することなく使うことができる. 
    !
    interface 
       function indexx(arrin)
         implicit none
         real(8), dimension(:), intent(in)  :: arrin
         integer, dimension(size(arrin))    :: indexx
       end function indexx
    end interface

   !------------ 引数 ------------
    complex(8), dimension(:,:)                :: cmat      ! 入力正方行列
    complex(8), intent(out), dimension(:)     :: eigen     ! 固有値実数部
    complex(8), 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(8), dimension(size(cmat,1))              :: ceig  ! 固有値
    complex(8), dimension(size(cmat,1),size(cmat,1)) :: cvec  ! 固有ベクトル
    real(8), dimension(size(cmat,1))                 :: vw    ! 作業変数
    integer, dimension(size(cmat,1))                 :: ivw   ! 作業変数
    integer, dimension(size(cmat,1))                 :: index ! 並び変え用
    integer, parameter :: mode=0   ! DCEIG1へ渡すスイッチ(平衡化省略)

    integer :: nm, i, j

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

    !------- DCEIG2/SSL2 による計算 ------
    call dceig2(cmat,nm,nm,mode,ceig,cvec,vw,ivw,info)

    !------- サブルーチンエラー処理 -------
    if ( info /= 0 ) then
       call MessageNotify('W','DCEIG2_SSL2',&
            '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) = cvec(:,j)
    enddo

  end subroutine dceig2_ssl2

end module ssl2_ceigen
