Zgges Zgges_r Zggesx Zggesx_r Zggev Zggevx

Zggevx() Sub Zggevx ( Balanc As  String, Jobvl As  String, Jobvr As  String, Sense As  String, N As  Long, A() As  Complex, B() As  Complex, Alpha() As  Complex, Beta() As  Complex, Vl() As  Complex, Vr() As  Complex, Ilo As  Long, Ihi As  Long, LScale() As  Double, RScale() As  Double, AbNrm As  Double, BbNrm As  Double, RConde() As  Double, RCondv() As  Double, Info As  Long  ) (エキスパートドライバ) 一般化固有値問題 (複素行列) 目的 本ルーチンはn×n複素行列のペア(A, B)の一般化固有値, および, 必要により左および／または右一般化固有ベクトルを求める. また, 必要により固有値および固有ベクトルの条件を改善するための均衡化の変換(Ilo, Ihi, Lscale, Rscale, AbNrmおよびBbNrm), 固有値の条件数の逆数(RConde), および, 右固有ベクトルの条件数の逆数(RCondv)を求める. 行列のペア(A, B)の一般化固有値は, A - λB が特異となるようなスカラーλあるいは比α/β=λである. β = 0 あるいは両方共0の場合に適当な解釈ができるため, これは通常ペア(α, β)で表される. (A, B)の一般化固有値λ(j)に対応する右一般化固有ベクトルv(j)は次式を満たす. A * v(j) = λ(j) * B * v(j) (A, B)の一般化固有値λ(j)に対応する左一般化固有ベクトルu(j)は次式を満たす. u(j)^H * A = λ(j) * u(j)^H * B ここで, u(j)^Hはu(j)の共役転置である. 引数 [in] Balanc 行うべき均衡化の方法を指定する. = "N": 対角スケーリングおよび入替を行わない. = "P": 入替えのみを行う. = "S": 対角スケーリングのみを行う. = "B": 対角スケーリングおよび入替を行う. 条件数の逆数は入替および／またはスケーリング後の行列について求められる. 入替により条件数は変わらない(丸め誤差がなければ)が, スケーリングにより変わる. [in] Jobvl = "N": 左一般化固有ベクトルを求めない. = "V": 左一般化固有ベクトルを求める. [in] Jobvr = "N": 右一般化固有ベクトルを求めない. = "V": 右一般化固有ベクトルを求める. [in] Sense どの条件数の逆数を計算するかを決定する. = "N": 計算しない. = "E": 固有値についてのみ計算する. = "V": 固有ベクトルについてのみ計算する. = "B": 固有値および固有ベクトルについて計算する. [in] N 行列A, B, VLおよびVRの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= N) [in] ペア(A, B)の行列 A. [out] A()は上書きされる. Jobvl = "V"あるいはJobvr = "V"あるいは両方の場合, 均衡化された入力AおよびBの複素シュール形の最初の部分が入る. [in,out] B() 配列 B(LB1 - 1, LB2 - 1) (LB1 >= N, LB2 >= N) [in] ペア(A, B)の行列 B. [out] B()は上書きされる. Jobvl = "V"あるいはJobvr = "V"あるいは両方の場合, 均衡化された入力AおよびBの複素シュール形の続きの部分が入る. [out] Alpha() 配列 Alpha(LAlpha - 1) (LAlpha >= N) [out] Beta() 配列 Beta(LBeta - 1) (LBeta >= N) Alpha(j)/Beta(j) (j = 0〜N-1) は一般化固有値である. 注: 商Alpha(j)/Beta(j)はオーバーフローあるいはアンダーフローし易く, Beta(j)が0になることさえある. 従って, 比α/βを単純に計算することは避けなければならない. しかし, Alphaは常にnorm(A)より小さく, 通常は同程度の大きさである. また, Betaは常にnorm(B)より小さく, 通常は同程度の大きさである. [out] Vl() 配列 Vl(LVl1 - 1, LVl2 - 1) (LVl1 >= N, LVl2 >= N) Jobvl = "V": 左一般化固有ベクトルu(j)が固有値と同じ順にVl()の列に入る. 固有ベクトルは, 最大要素が |実数部| + |虚数部| = 1 となるように正規化される. Jobvl = "N": 参照されない. [out] Vr() 配列 Vr(LVr1 - 1, LVr2 - 1) (LVr1 >= N, LVr2 >= N) Jobvr = "V": 右一般化固有ベクトルv(j)が固有値と同じ順にVr()の列に入る. 固有ベクトルは, 最大要素が |実数部| + |虚数部| = 1 となるように正規化される. Jobvr = "N": 参照されない. [out] Ilo [out] Ihi IloとIhiは, i > j かつ j = 0〜Ilo-2 または i = Ihi〜N-1 に対して, A(i, j) = 0, B(i, j) = 0 となるような整数値である. balanc = "N"あるいは"S"の場合, Ilo = 1, Ihi = N となる. [out] LScale() 配列 LScale(LLScale - 1) (LLScale >= N) 左からAとBに適用される入替とスケーリング因子の詳細. Pl(j)を行jと入替を行った行のインデックス, Dl(j)を行jに適用したスケーリング因子としたとき, 次のようになる.   LScale(j) = Pl(j) (j = 0 〜 Ilo-2)     = Dl(j) (j = Ilo-1 〜 Ihi-1)     = Pl(j) (j = Ihi 〜 N-1) 入替の順番はN-1からIhi, そして0からIlo-2である. [out] RScale() 配列 RScale(LRScale - 1) (LRScale >= N)右からAとBに適用される入替とスケーリング因子の詳細. Pr(j)を列jと入替を行った列のインデックス, Dr(j)を列jに適用したスケーリング因子としたとき, 次のようになる.   RScale(j) = Pr(j) (j = 0 〜 Ilo-2)     = Dr(j) (j = Ilo-1 〜 Ihi-1)     = Pr(j) (j = Ihi 〜 N-1) 入替の順番はN-1からIhi, そして0からIlo-2である. [out] AbNrm 均衡化された行列Aの1ノルム. [out] BbNrm 均衡化された行列Bの1ノルム. [out] RConde() 配列 RConde(LRConde - 1) (LRConde >= N) Sense = "E"または"B": 固有値の条件数の逆数が配列の要素に順に格納される. Sense = "N"または"V": RConde()は参照されない. [out] RCondv() 配列 RCondv(LRCondv - 1) (LRCondv >= N) Sense = "V"または"B": 固有ベクトルの条件数の逆数の推定値が配列の要素に順に格納される. RCondv(j)を求めるために固有値を並べ替えることができない場合, RCondv(j)は0に設定される. これは, 真の値が非常に小さい場合にのみ起こりうる. Sense = "N"または"E": RCondv()は参照されない. [out] Info = 0: 正常終了. = -1: パラメータ Balanc の誤り. (Balanc <> "N", "P", "S"および"B") = -2: パラメータ Jobvl の誤り. (Jobvl <> "V"および"N") = -3: パラメータ Jobvr の誤り. (Jobvr <> "V"および"N") = -4: パラメータ Sense の誤り. (Sense <> "N", "E", "V"および"B") = -5: パラメータ N の誤り. (N < 0) = -6: パラメータ A() の誤り. = -7: パラメータ B() の誤り. = -8: パラメータ Alpha() の誤り. = -9: パラメータ Beta() の誤り. = -10: パラメータ Vl() の誤り. = -11: パラメータ Vr() の誤り. = -14: パラメータ LScale() の誤り. = -15: パラメータ RScale() の誤り. = -18: パラメータ RConde() の誤り. = -19: パラメータ RCondv() の誤り. = i (0 < i <= N): QZ反復が収束しなかった. 固有ベクトルは計算されないが, Alpha(j), Beta(j) (j = i〜N-1) は正しい値である. = N+1: ZhgeqzにおいてQZ反復の収束以外のエラーが起きた. = N+2: Ztgevcにおいてエラーがあった. 出典 LAPACK 使用例 行列のペア(A, B)の一般化固有値および左および右一般化固有ベクトルを求める. ただし, ( 0.2-0.11i -0.93-0.32i 0.81+0.37i ) A = ( -0.8-0.92i -0.29+0.86i 0.64+0.51i ) ( 0.71+0.59i -0.15+0.19i 0.2+0.94i ) ( 0.57-0.91i -0.28-0.45i 0.25+0.91i ) B = ( 0.83-0.46i 0.63-0.19i -0.69+0.09i ) ( 0.24-1.33i -0.56-0.67i 0.9+1.25i ) とする. Sub Ex_Zggevx() Const N = 3 Dim A(N - 1, N - 1) As Complex, B(N - 1, N - 1) As Complex Dim Alpha(N - 1) As Complex, Beta(N - 1) As Complex Dim Vl(N - 1, N - 1) As Complex, Vr(N - 1, N - 1) As Complex, Info As Long Dim Ilo As Long, Ihi As Long, LScale(N - 1) As Double, RScale(N - 1) As Double Dim AbNrm As Double, BbNrm As Double Dim RConde(N - 1) As Double, RCondv(N - 1) As Double A(0, 0) = Cmplx(0.2, -0.11): A(0, 1) = Cmplx(-0.93, -0.32): A(0, 2) = Cmplx(0.81, 0.37) A(1, 0) = Cmplx(-0.8, -0.92): A(1, 1) = Cmplx(-0.29, 0.86): A(1, 2) = Cmplx(0.64, 0.51) A(2, 0) = Cmplx(0.71, 0.59): A(2, 1) = Cmplx(-0.15, 0.19): A(2, 2) = Cmplx(0.2, 0.94) B(0, 0) = Cmplx(0.57, -0.91): B(0, 1) = Cmplx(-0.28, -0.45): B(0, 2) = Cmplx(0.25, 0.91) B(1, 0) = Cmplx(0.83, -0.46): B(1, 1) = Cmplx(0.63, -0.19): B(1, 2) = Cmplx(-0.69, 0.09) B(2, 0) = Cmplx(0.24, -1.33): B(2, 1) = Cmplx(-0.56, -0.67): B(2, 2) = Cmplx(0.9, 1.25) Call Zggevx("N", "V", "V", "N", N, A(), B(), Alpha(), Beta(), Vl(), Vr(), Ilo, Ihi, LScale(), RScale(), AbNrm, BbNrm, RConde(), RCondv(), Info) Debug.Print "Eigenvalues =" Debug.Print Creal(Cdiv(Alpha(0), Beta(0))), Cimag(Cdiv(Alpha(0), Beta(0))), Debug.Print Creal(Cdiv(Alpha(1), Beta(1))), Cimag(Cdiv(Alpha(1), Beta(1))) Debug.Print Creal(Cdiv(Alpha(2), Beta(2))), Cimag(Cdiv(Alpha(2), Beta(2))) Debug.Print "Eigenvectors (L) =" Debug.Print Creal(Vl(0, 0)), Cimag(Vl(0, 0)), Creal(Vl(0, 1)), Cimag(Vl(0, 1)) Debug.Print Creal(Vl(1, 0)), Cimag(Vl(1, 0)), Creal(Vl(1, 1)), Cimag(Vl(1, 1)) Debug.Print Creal(Vl(2, 0)), Cimag(Vl(2, 0)), Creal(Vl(2, 1)), Cimag(Vl(2, 1)) Debug.Print Creal(Vl(0, 2)), Cimag(Vl(0, 2)) Debug.Print Creal(Vl(1, 2)), Cimag(Vl(1, 2)) Debug.Print Creal(Vl(2, 2)), Cimag(Vl(2, 2)) Debug.Print "Eigenvectors (R) =" Debug.Print Creal(Vr(0, 0)), Cimag(Vr(0, 0)), Creal(Vr(0, 1)), Cimag(Vr(0, 1)) Debug.Print Creal(Vr(1, 0)), Cimag(Vr(1, 0)), Creal(Vr(1, 1)), Cimag(Vr(1, 1)) Debug.Print Creal(Vr(2, 0)), Cimag(Vr(2, 0)), Creal(Vr(2, 1)), Cimag(Vr(2, 1)) Debug.Print Creal(Vr(0, 2)), Cimag(Vr(0, 2)) Debug.Print Creal(Vr(1, 2)), Cimag(Vr(1, 2)) Debug.Print Creal(Vr(2, 2)), Cimag(Vr(2, 2)) Debug.Print "Info =", Info End Sub 実行結果 Eigenvalues = -0.4784814787767 -0.640760182519056 1.62433774044009 1.10642263894432 -0.149178317709927 5.63446258373312 Eigenvectors (L) = -0.488182728920006 0.511817271079994 0.112214636564846 -0.78086108170552 -0.493389319144031 1.33311637129139E-02 -0.132533171028514 -0.347190877300084 0.339735049171677 -6.82051937200691E-02 0.672543628409034 0.327456371590966 0.902635460243872 9.73645397561283E-02 -0.178648471173719 0.129164561360569 -0.730361156069673 7.60898456031052E-02 Eigenvectors (R) = -0.573248979253811 -0.426751020746189 0.644105768156504 -0.298230968832631 -0.249087460228498 -9.31247915124554E-02 -0.701765691258405 -0.298234308741595 0.271914164828682 -0.477707190075455 0.274425774522266 -0.280815941661082 2.08352227991442E-02 2.31726206010056E-02 -0.737220272842649 -0.262779727157351 -0.474367740415068 -4.82054262401725E-02 Info = 0