内部では DM_*/SSL2TP ルーチンによる複素行列の固有値/固有ベクトル計算を
行っている. が, ユーザーは用いているライブラリとサブルーチンを意識
することなく使うことができる.
subroutine dm_ceigen_ssl2tp(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 は保存されない.
!
! 内部では DM_*/SSL2TP ルーチンによる複素行列の固有値/固有ベクトル計算を
! 行っている. が, ユーザーは用いているライブラリとサブルーチンを意識
! することなく使うことができる.
!
!------------ 引数 ------------
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),size(cmat,1)) :: cmatrix ! 行列保存
complex(8), dimension(size(cmat,1)) :: ceig ! 固有値
complex(8), dimension(size(cmat,1),size(cmat,1)) :: cvec ! 固有ベクトル
real(8), dimension(size(cmat,1)) :: dv ! 作業変数
integer, dimension(size(cmat,1)) :: ip ! 作業変数
integer, dimension(size(cmat,1)) :: ind ! 作業変数
complex(8), dimension(size(cmat,1),size(cmat,1)+1) :: zaw ! 作業変数
integer, dimension(size(cmat,1)) :: index ! 並び変え用
integer, parameter :: mode=0 ! DCEIG1へ渡すスイッチ(平衡化省略)
integer :: nm, i, j, neig
!------- 形状チェック ------
if (size(cmat,1) /= size(cmat,2))then
call MessageNotify('E','DM_CEIGEN_SSL2','Input matrix not square')
else
nm = size(cmat,1)
endif
!------- DM_*/SSL2TP による計算 ------
call dm_cblnc(cmat, nm, nm, dv, info)
if ( info .eq. 30000 ) then
call MessageNotify('E','DM_CEIGEN_SSL2TP','Error in DM_CBLNC')
return
end if
call dm_ches2(cmat, nm, nm, ip, info)
cmatrix = cmat
call dchsqr(cmat, nm, nm, ceig, neig, info)
if ( info .eq. 15000 ) then
call MessageNotify('E','DM_CEIGEN_SSL2TP', 'Eigenvalue not completely obtaind in DCHSQR')
return
else if ( info .ge. 20000 ) then
call MessageNotify('E','DM_CEIGEN_SSL2TP','Error in DM_CHSQR')
return
end if
ind = 1
call dm_chvec(cmatrix, nm, nm, ceig, ind, neig, cvec, zaw, info)
if ( info .ge. 15000 ) then
call MessageNotify('E','DM_CEIGEN_SSL2TP', 'Eigenvector not completely obtaind in DM_CHVEC')
return
else if ( info .ge. 20000 ) then
call MessageNotify('E','DM_CEIGEN_SSL2TP','Error in DM_CHSQR')
return
end if
call dm_chbk2(cvec, nm, nm, ind, neig, cmat, ip, dv, info)
call dm_cnrml(cvec, nm, nm, ind, neig, 1, info)
!------- 固有ベクトル入れ換え -------
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) = cvec(:,j)
enddo
end subroutine dm_ceigen_ssl2tp