内部では ZGEEV/LAPACK ルーチンによる実行列の固有値/固有ベクトル計算を
行っている. が, ユーザーは用いているライブラリとサブルーチンを意識
することなく使うことができる.
subroutine dceigen_lapack(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 は保存されない.
!
! 内部では ZGEEV/LAPACK ルーチンによる実行列の固有値/固有ベクトル計算を
! 行っている. が, ユーザーは用いているライブラリとサブルーチンを意識
! することなく使うことができる.
!
!------------ 引数 ------------
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)) :: vl ! 左固有ベクトル
complex(8), dimension(size(cmat,1),size(cmat,1)) :: vr ! 右固有ベクトル
complex(8), dimension(size(cmat,1)*2) :: work ! 作業変数
real(8), dimension(size(cmat,1)*2) :: rwork ! 作業変数
integer, dimension(size(cmat,1)) :: index ! 並び変え用
character(len=1), parameter :: jobvl='N', jobvr='V' ! ZEGGV へ渡すスイッチ
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(dimag(ceig))
elseif ( sort == 'IA' ) then ! 固有値虚部の絶対値
index=indexx(abs(dimag(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 dceigen_lapack