|
|
◆ Zhbgvx()
| Sub Zhbgvx |
( |
Jobz As |
String, |
|
|
Range As |
String, |
|
|
Uplo As |
String, |
|
|
N As |
Long, |
|
|
Ka As |
Long, |
|
|
Kb As |
Long, |
|
|
Ab() As |
Complex, |
|
|
Bb() As |
Complex, |
|
|
Q() As |
Complex, |
|
|
Vl As |
Double, |
|
|
Vu As |
Double, |
|
|
Il As |
Long, |
|
|
Iu As |
Long, |
|
|
AbsTol As |
Double, |
|
|
M As |
Long, |
|
|
W() As |
Double, |
|
|
Z() As |
Complex, |
|
|
IFail() As |
Long, |
|
|
Info As |
Long |
|
) |
| |
(エキスパートドライバ) 一般化固有値問題 (エルミート帯行列)
- 目的
- 本ルーチンはエルミート帯行列の一般化固有値問題 の選択された固有値, および, 必要により固有ベクトルを求める. ここで, AとBはエルミート行列, さらにBは正定値である.
全固有値, または, 必要な固有値の範囲あるいは番号の範囲を指定することにより, 求める固有値・固有ベクトルを選択することができる.
- 引数
-
| [in] | Jobz | = "N": 固有値のみ求める.
= "V": 固有値と固有ベクトルを求める. |
| [in] | Range | = "A": すべての固有値を求める.
= "V": 半開区間(Vl, Vu]のすべての固有値を求める.
= "I": Il番目からIu番目までの固有値を求める. |
| [in] | Uplo | = "U": AおよびBの上三角部分を格納する.
= "L": AおよびBの下三角部分を格納する. |
| [in] | N | 行列AおよびBの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) |
| [in] | Ka | 行列Aの上帯幅(Uplo = "U"の場合)または下帯幅(Uplo = "L"の場合). (ka >= 0) |
| [in] | Kb | 行列Bの上帯幅(Uplo = "U"の場合)または下帯幅(Uplo = "L"の場合). (kb >= 0) |
| [in,out] | Ab() | 配列 Ab(LAb1 - 1, LAb2 - 1) (LAb1 >= Ka + 1, LAb2 >= N)
[in] Ka+1×N対称帯行列形式のN×Nエルミート帯行列 A. Uploに従って上または下三角部分を格納する.
[out] Ab()の内容は壊される. |
| [in,out] | Bb() | 配列 Bb(LBb1 - 1, LBb2 - 1) (LBb1 >= Kb + 1, LBb2 >= N)
[in] Kb+1×N対称帯行列形式のN×N正定値エルミート帯行列 B. Uploに従って上または下三角部分を格納する.
[out] Zpbstfにより求められたスプリットコレスキー分解 B = S^H*S のSが入る. |
| [out] | Q() | 配列 Q(LQ1 - 1, LQ2 - 1) (LQ1 >= N, LQ2 >= N)
Jobz = "V": Ax = λBx を標準形 Cx = λx に変換(Cは3重対角行列になる)する際に使われたN×N行列. Jobz = "N": 配列Q()は参照されない. |
| [in] | Vl | Range = "V": 固有値を求める区間の下端. (Vl < Vu)
Range = "A"または"I": 参照されない. |
| [in] | Vu | Range = "V": 固有値を求める区間の上端. (Vl < Vu)
Range = "A"または"I": 参照されない. |
| [in] | Il | Range = "I": 求める最小固有値の番号. (1 <= Il <= Iu <= N (N > 0 の場合), Il = 1, Iu = 0 (N = 0 の場合))
Range = "A"または"V": 参照されない. |
| [in] | Iu | Range = "I": 求める最大固有値の番号. (1 <= Il <= Iu <= N (N > 0 の場合), Il = 1, Iu = 0 (N = 0 の場合))
Range = "A"または"V": 参照されない. |
| [in] | AbsTol | 固有値の絶対誤差限界.
固有値の近似値は区間[a, b]に入っているときに収束したものとみなされる. この区間の幅は AbsTol + eps*max(|a|, |b|) に等しいかこれより小さい. ここでepsはマシンイプシロンである. AbsTol <= 0 の場合, eps*|T|が代わりに使用される. ここで, |T|は行列Aを変換して得られた3重対角行列の1ノルムである.
AbsTolを0ではなくアンダフロー限界の2倍(2*Dlamch("S"))に設定したときに固有値が最も正確に求められる. Info > 0 (固有ベクトルのいくつかが収束しなかったことを示す)で戻ったときには, AbsTolを2*Dlamch("S")に設定してみるとよい. |
| [out] | M | 求められた固有値の数. (0 <= M <= N)
Range = "A"の場合 M = N, Range = "I"の場合 M = Iu - Il + 1 となる. |
| [out] | W() | 配列 W(LW - 1) (LW >= N)
正常終了時, M個の求められた固有値が先頭から昇順に入る. |
| [out] | Z() | 配列 Z(LZ1 - 1, LZ2 - 1) (LZ1 >= N, LZ2 >= M)
Jobz = "V": Info = の場合, 求められた固有値に対応してZ()の最初のM列に行列Aの正規直交固有ベクトルが入る. W(i)に関連する固有ベクトルがZ()のi列に入る.
固有ベクトルは Z^H*B*Z = I となるように正規化される.
固有ベクトルの収束に失敗した場合, Z()のその列には固有ベクトルの最終近似が入り, IFail()に固有ベクトルの番号が入る.
Jobz = "N": Z()は参照されない.
注: 配列Z()は少なくてもmax(1, M)列を割り当てること. Range = "V"の場合, Mの値をあらかじめ知ることはできないが上限値を使用すること. |
| [out] | ifail[] | 配列 ifail[lifail] (lifail >= n)
jobz = "V": info = 0の場合, ifail[]の最初のm要素が0に設定される. info > 0 の場合, 収束しなかった固有ベクトルの番号がifail[]に入る.
jobz = "N": ifail[]は参照されない. |
| [out] | info | = 0: 正常終了.
= -1: パラメータ Jobz の誤り. (Jobz <> "V"および"N")
= -2: パラメータ Range の誤り. (Range <> "A", "V"および"I")
= -3: パラメータ Uplo の誤り. (Uplo <> "U"および"L")
= -4: パラメータ N の誤り. (N < 0)
= -5: パラメータ Ka の誤り. (Ka < 0)
= -6: パラメータ Kb の誤り. (Kb < 0)
= -7: パラメータ Ab() の誤り.
= -8: パラメータ Bb() の誤り.
= -9: パラメータ Q() の誤り.
= -11: パラメータ Vu の誤り. (Vu <= Vl)
= -12: パラメータ Il の誤り. (Il < 1 または Il > N)
= -13: パラメータ Iu の誤り. (Iu < min(N, Il) または Iu > N)
= -16: パラメータ W() の誤り.
= -17: パラメータ Z() の誤り.
= -18: パラメータ IFail() の誤り.
= i (0 < i <= N): i個の固有ベクトルが収束しなかった. IFail()にその番号が入る.
= i (i > N): ZpbstfがInfo = i-N を返した. Bの i-N次小行列が正定値でない. Bの分解が完了できず, 固有値・固有ベクトルは計算されなかった. |
- 出典
- LAPACK
- 使用例
- 一般化固有値問題 Ax = λBx の固有値および固有ベクトルを求める. ここで, Aはエルミート帯行列, Bは正定値エルミート帯行列である. ただし,
( 0.20 -0.11+0.93i 0 )
A = ( -0.11-0.93i -0.32 -0.80+0.92i )
( 0 -0.80-0.92i -0.29 )
( 2.20 -0.11+0.93i 0 )
B = ( -0.11-0.93i 2.32 -0.80+0.92i )
( 0 -0.80-0.92i 2.29 )
Sub Ex_Zhbgvx()
Const N = 3, Ka = 1, Kb = 1
Dim Ab(Ka, N - 1) As Complex, Bb(Kb, N - 1) As Complex
Dim W(N - 1) As Double, Q(N - 1, N - 1) As Complex
Dim Vl As Double, Vu As Double, Il As Long, Iu As Long, AbsTol As Double
Dim M As Long, Z(N - 1, N - 1) As Complex, IFail(N - 1) As Long, Info As Long
Ab(0, 0) = Cmplx(0.2): Ab(0, 1) = Cmplx(-0.32): Ab(0, 2) = Cmplx(-0.29)
Ab(1, 0) = Cmplx(-0.11, -0.93): Ab(1, 1) = Cmplx(-0.8, -0.92)
Bb(0, 0) = Cmplx(2.2): Bb(0, 1) = Cmplx(2.32): Bb(0, 2) = Cmplx(2.29)
Bb(1, 0) = Cmplx(-0.11, -0.93): Bb(1, 1) = Cmplx(-0.8, -0.92)
Call Zhbgvx("V", "A", "L", N, Ka, Kb, Ab(), Bb(), Q(), Vl, Vu, Il, Iu, AbsTol, M, W(), Z(), IFail(), Info)
Debug.Print "Eigenvalues =", W(0), W(1), W(2)
Debug.Print "Eigenvectors ="
Debug.Print Creal(Z(0, 0)), Cimag(Z(0, 0)), Creal(Z(0, 1)), Cimag(Z(0, 1))
Debug.Print Creal(Z(1, 0)), Cimag(Z(1, 0)), Creal(Z(1, 1)), Cimag(Z(1, 1))
Debug.Print Creal(Z(2, 0)), Cimag(Z(2, 0)), Creal(Z(2, 1)), Cimag(Z(2, 1))
Debug.Print Creal(Z(0, 2)), Cimag(Z(0, 2))
Debug.Print Creal(Z(1, 2)), Cimag(Z(1, 2))
Debug.Print Creal(Z(2, 2)), Cimag(Z(2, 2))
Debug.Print "M =", M, "Info =", Info
End Sub
- 実行結果
Eigenvalues = -2.33466883563969 4.53015245998351E-03 0.35977634455495
Eigenvectors =
0.47826114359381 -6.39634991743266E-18 0.507837765924615 -8.10231908112231E-19
9.59938458849389E-02 0.811584333390847 1.21596344669029E-02 0.102804182311088
-0.441771484124976 0.486432071889564 0.281208537371776 -0.309637123223848
0.332163577636529 2.56485488839366E-18
-3.84923101537441E-02 -0.325434985845291
-0.154385502431493 0.169992999811227
M = 3 Info = 0
|
1.9.6