Zhbgv Zhbgvd Zhbgvx Zhegv Zhegvd Zhegvx Zhpgv Zhpgvd Zhpgvx

Zhpgvd() Sub Zhpgvd ( IType As  Long, Jobz As  String, Uplo As  String, N As  Long, Ap() As  Complex, Bp() As  Complex, W() As  Double, Z() As  Complex, 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] 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] エルミート行列 A の上または下三角部分. 1次元配列に列ごとに圧縮されており, A のj列が配列 Ap() に次のように格納される Uplo = "U": Ap(i + j*(j + 1)/2) = Aij. ただし, 0 <= i <= j <= N - 1. Uplo = "L": Ap((i + j*(2*N - j - 1)/2) = Aij. ただし, 0 <= j < = i <= N - 1. [out] Ap() の内容は壊される. [in,out] Bp() 配列 Bp(LBp - 1) (LBp >= N(N + 1)/2) [in] 正定値エルミート行列 B の上または下三角部分. 1次元配列に列ごとに圧縮されており, B のj列が配列 Bp() に次のように格納される Uplo = "U": Bp(i + j*(j + 1)/2) = Bij. ただし, 0 <= i <= j <= N - 1. Uplo = "L": Bp((i + j*(2*N - j - 1)/2) = Bij. ただし, 0 <= j < = i <= N - 1. [out] コレスキー分解 B = U^H*U あるいは B = L*L^H の三角行列 U あるいは L が B と同じ格納形式で入る. [out] W() 配列 W(LW - 1) (LW >= N) Info = 0 の場合, 求められた固有値(昇順). [out] Z() 配列 Z(LZ1 - 1, LZ2 - 1) (LZ1 >= N, LZ2 >= N) Jobz = "V": Info = 0 の場合, 固有ベクトルからなる行列 Z を返す. 固有ベクトルは次のように正規化される.   Itype = 1 または 2: Z^H*B*Z = I   Itype = 3: Z^H*inv(B)*Z = I Jobz = "N": Z() は参照されない. [out] Info = 0: 正常終了 = -1: パラメータ Itype の誤り (Itype < 1 or Itype > 3) = -2: パラメータ Jobz の誤り (Jobz <> "V" および "N") = -3: パラメータ Uplo の誤り (Uplo <> "U" および "L") = -4: パラメータ N の誤り (N < 0) = -5: パラメータ Ap() の誤り. = -6: パラメータ Bp() の誤り. = -7: パラメータ W() の誤り. = -8: パラメータ Z() の誤り. = i (0 < i <= N): Zspevd が収束しなかった. 中間の3重対角形の非対角要素のうち i 個が 0 に収束しなかった. = 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_Zhpgvd() Const N = 3 Dim Ap(N * (N + 1) / 2) As Complex, Bp(N * (N + 1) / 2) As Complex Dim W(N - 1) As Double, Z(N - 1, N - 1) As Complex, 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) Call Zhpgvd(1, "V", "L", N, Ap(), Bp(), W(), Z(), 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 "Info =", Info End Sub 実行結果 Eigenvalues = -4.97466704628586 5.03789053905172E-02 0.392451765416646 Eigenvectors = -0.819151009277443 -3.74256105252903E-17 -0.413822310464714 3.74256105252903E-17 0.320545559438499 -0.890210156692715 -0.332130735287233 1.67477550936547E-02 0.935832442766443 -6.14475107289119E-02 -0.126697019588926 0.325735727937686 -0.354935871366059 -8.42076236819031E-17 0.17251768647939 0.196424402514893 -0.125117400789647 -0.228553336347363 Info = 0