Dsbgv Dsbgvd Dsbgvx Dspgv Dspgvd Dspgvx Dsygv Dsygvd Dsygvx

Dspgvd() Sub Dspgvd ( IType As  Long, Jobz As  String, Uplo As  String, N As  Long, Ap() As  Double, Bp() As  Double, W() As  Double, Z() As  Double, 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^T*U あるいは B = L*L^T の三角行列 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^T*B*Z = I   Itype = 3: Z^T*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): Dspevd が収束しなかった. 中間の3重対角形の非対角要素のうち i 個が 0 に収束しなかった. = i (i > N): B の i-N 次小行列が正定値でない. B の分解が完了できず, 固有値・固有ベクトルは計算されなかった. 出典 LAPACK 使用例 一般化固有値問題 Ax = λBx の固有値および固有ベクトルを求める. ここで, Aは対称行列, Bは正定値対称行列である. ただし, ( 0.54 -0.90 -0.94 ) ( 1.18 0.54 -1.22 ) A = ( -0.90 0.70 1.04 ) B = ( 0.54 0.60 -0.71 ) ( -0.94 1.04 1.65 ) ( -1.22 -0.71 1.66 ) とする. Sub Ex_Dspgvd() Const N = 3 Dim Ap(N * (N + 1) / 2) As Double, Bp(N * (N + 1) / 2) As Double Dim W(N - 1) As Double, Z(N - 1, N - 1) As Double Dim Info As Long Ap(0) = 0.54 Ap(1) = -0.9: Ap(3) = 0.7 Ap(2) = -0.94: Ap(4) = 1.04: Ap(5) = 1.65 Bp(0) = 1.18 Bp(1) = 0.54: Bp(3) = 0.6 Bp(2) = -1.22: Bp(4) = -0.71: Bp(5) = 1.66 Call Dspgvd(1, "V", "L", N, Ap(), Bp(), W(), Z(), Info) Debug.Print "Eigenvalues =", W(0), W(1), W(2) Debug.Print "Eigenvectors =" Debug.Print Z(0, 0), Z(0, 1), Z(0, 2) Debug.Print Z(1, 0), Z(1, 1), Z(1, 2) Debug.Print Z(2, 0), Z(2, 1), Z(2, 2) Debug.Print "Info =", Info End Sub 実行結果 Eigenvalues = -0.297963342573455 0.510423243055614 7.37370278804149 Eigenvectors = 1.17970064313729 1.46384155786189 9.13803098094182E-02 0.497345303675336 -0.442269445774235 -1.71612038711754 0.524910948248221 1.35130854246605 -0.947378600715302 Info = 0