Zhbev Zhbevd Zhbevx Zheev Zheevd Zheevr Zheevx Zhpev Zhpevd Zhpevx

Zheevr() Sub Zheevr ( Jobz As  String, Range As  String, Uplo As  String, N As  Long, A() 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, Isuppz() As  Long, Info As  Long  ) (MRRR法ドライバ) 固有値・固有ベクトル (エルミート行列) 目的 本ルーチンはエルミート行列 A のすべての, あるいは, 選択された固有値, および, 必要により固有ベクトルを求める. 必要な固有値の範囲あるいは番号の範囲を指定することにより, 求める固有値・固有ベクトルを選択することができる. 本ルーチンはまず行列 A を3重対角形 T に変換する. 次に, すべての固有値および固有ベクトルを求める場合にはMRRR法(*)を使用する. また, すべての固有値のみを求める場合にはQL法またはQR法を使用する. 選択された固有値・固有ベクトルを求める場合には二分法および逆反復法を使用する. (*) MRRR法: 固有値は dqds アルゴリズムにより求められる. 固有ベクトルは適切なシフト付きLDL^T分解(RRR(Relatively Robust Representations)という)より求められる. 引数 [in] Jobz = "N": 固有値のみ求める. = "V": 固有値と固有ベクトルを求める. [in] Range = "A": すべての固有値を求める. = "V": 半開区間 (Vl, Vu] のすべての固有値を求める. = "I": Il番目からIu番目までの固有値を求める. Range = "V" または "I" かつ Iu - Il < N - 1 であれば Dstebz および Zstein を呼び出す. [in] Uplo = "U": A の上三角部分を格納する. = "L": A の下三角部分を格納する. [in] N 行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= N) [in] エルミート行列 A. Uplo = "U" の場合, A() の N x N 上三角部分に行列 A の上三角部分を格納する. Uplo = "L" の場合, A() の N x N 下三角部分に行列 A の下三角部分を格納する. [out] A() の下三角部分(Uplo = "L"の場合)あるいは上三角部分(Uplo = "U"の場合)は対角部分を含め壊される. [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ノルムである. Demmel および Kahan, "Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy", LAPACK Working Note #3. を参照せよ. 高い相対精度が重要な場合は, Abstol を Dlamch("S") に設定せよ. それにより, 将来のリリースで可能な場合, 固有値が高い相対精度で計算されることが保証される. 現在のコードは高い相対精度を保証するものではないが, 将来のリリースでは保証される. どの行列が固有値を高い相対精度で定義するかの議論については J. Barlow および J. Demmel, "Computing Accurate Eigensystems of Scaled Diagonally Dominant Matrices", LAPACK Working Note #7 を参照せよ. [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() 配列 A(LZ1 - 1, LZ2 - 1) (LZ1 >= N, LZ2 >= M) Jobz = "V": Info = 0 の場合, 求められた固有値に対応して Z() の最初の M 列に行列 A の正規直交固有ベクトルが入る. W(i) に関連する固有ベクトルが Z() の i 列に入る. Jobz = "N": Z() は参照されない. 注: 配列 Z() は少なくても max(1, M) 列を割り当てること. Range = "V" の場合, M の値をあらかじめ知ることはできないが上限値を使用すること. [out] Isuppz() 配列 Isuppz(LIsuppz - 1) (LIsuppz >= 2*max(1, M)) Z() の固有ベクトルの非ゼロ要素のインデックスを表す. i(= 1 〜 M)番目の固有ベクトルにおいては Isuppz(2*(i-1)) 〜 Isuppz(2*(i-1) + 1) 番目の要素のみが非ゼロである. これは Zstemr の出力である(三重対角行列). A の固有ベクトルに関しては Zunmtr によって適用されるユニタリ変換のため通常 1 〜 N である. Range = "A" の場合, または, Range = "I" で Iu - Il = N - 1 の場合に有効である. [out] Info = 0: 正常終了. = -1: パラメータ Jobz の誤り. (Jobz <> "N" および "V") = -2: パラメータ Range の誤り. (Range <> "A", "V" および "I") = -3: パラメータ Uplo の誤り. (Uplo <> "U" および "L") = -4: パラメータ N の誤り. (N < 0) = -5: パラメータ A() の誤り. = -7: パラメータ Vu の誤り. (Vu < Vl) = -8: パラメータ Il の誤り. (Il < 1 または Il > N) = -9: パラメータ Iu の誤り. (Iu < min(N, Il) または Iu > N) = -12: パラメータ W() の誤り. = -13: パラメータ Z() の誤り. = -14: パラメータ Isuppz() の誤り. = > 0: 内部エラー. 出典 LAPACK 使用例 エルミート行列 A の固有値・固有ベクトルを求める. ただし, ( 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 ) とする. Sub Ex_Zheevr() Const N = 3 Dim A(N - 1, N - 1) As Complex, 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, Isuppz(2 * N - 1) As Long, Info As Long A(0, 0) = Cmplx(0.2, 0) A(1, 0) = Cmplx(-0.11, -0.93): A(1, 1) = Cmplx(-0.32, 0) A(2, 0) = Cmplx(0.81, 0.37): A(2, 1) = Cmplx(-0.8, -0.92): A(2, 2) = Cmplx(-0.29, 0) AbsTol = 0 Call Zheevr("V", "A", "L", N, A(), Vl, Vu, Il, Iu, AbsTol, M, W(), Z(), Isuppz(), 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.05348849668514 0.124622388617309 1.51886610806783 Eigenvectors = -0.449276526719114 -0 0.654793596518192 0 0.227247885813611 -0.597641779578735 0.519997178670922 -3.19846835072559E-02 0.621236109316912 -5.83009495222977E-02 0.204907317474215 -0.507777757881847 0.607779522934083 0 -0.392237107311198 -0.407323787101333 0.23846608290599 0.503959683819116 M = 3 Info = 0