Zhbgv Zhbgvd Zhbgvx Zhegv Zhegvd Zhegvx Zhpgv Zhpgvd Zhpgvx

Zhpgvx() Sub Zhpgvx ( IType As  Long, Jobz As  String, Range As  String, Uplo As  String, N As  Long, Ap() As  Complex, Bp() 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  ) (エキスパートドライバ) 一般化固有値問題 (エルミート行列) (圧縮形式) 目的 本ルーチンはエルミート行列の一般化固有値問題 Ax = λBx, ABx = λx または BAx = λx の選択された固有値, および, 必要により固有ベクトルを求める. ここで, AとBはエルミート行列(圧縮形式)で, さらにBは正定値である. 必要な固有値の範囲あるいは番号の範囲を指定することにより, 求める固有値・固有ベクトルを選択することができる. 引数 [in] IType 解くべき問題のタイプを指定. = 1: Ax = λBx. = 2: ABx = λx. = 3: BAx = λx. [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,out] Ap() 配列 Ap(LAp - 1) (LAp >= N(N + 1)/2) [in] 圧縮形式のN×Nエルミート行列 A. Uploに従い上三角部分あるいは下三角部分を格納する. [out] Ap()の内容は壊される. [in,out] Bp() 配列 Bp(LBp - 1) (LBp >= N(N + 1)/2) [in] 圧縮形式のN×N正定値エルミート行列 B. Uploに従い上三角部分あるいは下三角部分を格納する. [out] コレスキー分解 B = U^H*U あるいは B = L*L^H の三角行列UあるいはLがBと同じ格納形式で入る. [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列に入る.   固有ベクトルは次のように正規化される.     IType = 1 または 2: Z^H*B*Z = I     IType = 3: Z^H*inv(B)*Z = I   固有ベクトルの収束に失敗した場合, Z()のその列には固有ベクトルの最終近似が入り, IFail()に固有ベクトルの番号が入る. Jobz = "N": Z()は参照されない. 注: 配列Z()は少なくてもmax(1, M)列を割り当てること. Range = "V"の場合, Mの値をあらかじめ知ることはできないが上限値を使用すること. [out] IFail() 配列 IFail(LIFail - 1) (LIFail >= N) Jobz = "V": Info = 0の場合, IFail()の最初のM要素が0に設定される. Info > 0 の場合, 収束しなかった固有ベクトルの番号がIFail()に入る. Jobz = "N": IFail()は参照されない. [out] Info = 0: 正常終了. = -1: パラメータ Itype の誤り. (Itype < 1 または Itype > 3) = -2: パラメータ Jobz の誤り. (Jobz <> "V"および"N") = -3: パラメータ Range の誤り. (Range <> "A", "V"および"I") = -4: パラメータ Uplo の誤り. (Uplo <> "U"および"L") = -5: パラメータ N の誤り. (N < 0) = -6: パラメータ Ap() の誤り. = -7: パラメータ Bp() の誤り. = -9: パラメータ Vu の誤り. (Vu <= Vl) = -10: パラメータ Il の誤り. (Il < 1 または Il > N) = -11: パラメータ Iu の誤り. (Iu < min(N, Il) または Iu > N) = -14: パラメータ W() の誤り. = -15: パラメータ Z() の誤り. = -16: パラメータ IFail() の誤り. = i (0 < i <= N): Zhpevxが収束しなかった. i個の固有ベクトルが収束しなかった. 配列IFail()にその番号が入る. = i (i > N): Bの i-N次小行列が正定値でない. Bの分解が完了できず, 固有値・固有ベクトルは計算されなかった. 出典 LAPACK 使用例 一般化固有値問題 Ax = λBx の固有値および固有ベクトルを求める. ここで, Aはエルミート行列, Bは正定値エルミート行列である. ただし, ( 0.20 -0.11+0.93i 0.81-0.37i ) A = ( -0.11-0.93i -0.32 -0.80+0.92i ) ( 0.81+0.37i -0.80-0.92i -0.29 ) ( 2.20 -0.11+0.93i 0.81-0.37i ) B = ( -0.11-0.93i 2.32 -0.80+0.92i ) ( 0.81+0.37i -0.80-0.92i 2.29 ) とする. Sub Ex_Zhpgvx() Const N = 3 Dim Ap(N * (N + 1) / 2) As Complex, Bp(N * (N + 1) / 2) As Complex Dim W(N - 1) As Double 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 Ap(0) = Cmplx(0.2, 0) Ap(1) = Cmplx(-0.11, -0.93): Ap(3) = Cmplx(-0.32, 0) Ap(2) = Cmplx(0.81, 0.37): Ap(4) = Cmplx(-0.8, -0.92): Ap(5) = Cmplx(-0.29, 0) Bp(0) = Cmplx(2.2, 0) Bp(1) = Cmplx(-0.11, -0.93): Bp(3) = Cmplx(2.32, 0) Bp(2) = Cmplx(0.81, 0.37): Bp(4) = Cmplx(-0.8, -0.92): Bp(5) = Cmplx(2.29, 0) AbsTol = 0 Call Zhpgvx(1, "V", "A", "L", N, Ap(), Bp(), 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 = -4.97466704628586 5.03789053905172E-02 0.392451765416646 Eigenvectors = -0.819151009277443 0 -0.413822310464714 9.35640263132256E-17 0.320545559438499 -0.890210156692715 -0.332130735287233 1.67477550936547E-02 0.935832442766443 -6.14475107289119E-02 -0.126697019588926 0.325735727937686 -0.354935871366058 -3.74256105252903E-17 0.17251768647939 0.196424402514893 -0.125117400789647 -0.228553336347363 M = 3 Info = 0