Zgges Zgges_r Zggesx Zggesx_r Zggev Zggevx

Zggesx() Sub Zggesx ( Jobvsl As  String, Jobvsr As  String, Sort As  String, Selctg As  LongPtr, Sense As  String, N As  Long, A() As  Complex, B() As  Complex, Sdim As  Long, Alpha() As  Complex, Beta() As  Complex, Vsl() As  Complex, Vsr() As  Complex, RConde() As  Double, RCondv() As  Double, Info As  Long  ) (エキスパートドライバ) 一般化シュール分解 (複素行列) 目的 本ルーチンはn×n複素行列のペア(A, B)の一般化固有値, 一般化複素シュール形(S, T), および, 必要により左および／または右シュールベクトルからなる行列(VSLおよびVSR)を求める. これらにより, 一般化シュール分解が次のように与えられる. (A, B) = ((VSL)*S*(VSR)^H, (VSL)*T*(VSR)^H) (VSR)^HはVSRの共役転置である. また, 必要により, 選択された固有値が上三角行列Sおよび上三角行列Tの左上の対角ブロックにくるように固有値を並べ替え, 選択された固有値の平均の条件数の逆数(RConde), および, 選択された固有値に対応する右および左部分空間の条件数の逆数(RCondv)を求める. その結果VSLおよびVSRの先行する側の列は, 対応する左および右固有空間(部分空間)のユニタリ基底を形成する. 行列のペア(A, B)の一般化固有値は, A - wB が特異となるようなスカラーwあるいは比α/β = wである. β = 0 あるいは両方共0の場合に適当な解釈ができるため, これは通常ペア(α, β)で表される. Tが非負の対角要素を持つ上三角行列でSが上三角行列であれば, 行列のペア(S, T)は一般化複素シュール形である. 引数 [in] Jobvsl = "N": 左シュールベクトルを求めない. = "V": 左シュールベクトルを求める. [in] Jobvsr = "N": 右シュールベクトルを求めない. = "V": 右シュールベクトルを求める. [in] Sort 一般化シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する. = "N": 固有値の並べ替えを行わない. = "S": 固有値の並べ替えを行う(Selctgを参照のこと). [in] Selctg Sort = "S": シュール形の左上にくるように並べ替える固有値を選択するためにSelctgが使われる.   Selctg(Alpha(j), Beta(j))がtrue (= 1)であれば, 固有値 Alpha(j)/Beta(j) が選択される.   並べ替え後, 選択された複素固有値は Selctg(Alpha(j), Beta(j)) = true を満たさなくなるかもしれないことに注意せよ. これは, 並べ替えによって複素固有値の値が変わる可能性があるからである(特に固有値が悪条件の場合). この場合, Info = N+2 に設定される(下のInfoを参照のこと). Sort = "N": Selctgは参照されない. [in] Sense どの条件数を計算するかを決定する. = "N": 条件数を計算しない. = "E": 選択された固有値の平均について条件数を計算する. = "V": 選択された部分空間について条件数を計算する. = "B": 両方について条件数を計算する. Sense = "E", "V"あるいは"B"の場合, Sort = "S"でなければならない. [in] N 行列A, B, VSLおよびVSRの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= N) [in] 行列ペアの1つ目の行列. [out] A()は一般化シュール形Sで上書きされる. [in,out] B() 配列 B(LB1 - 1, LB2 - 1) (LB1 >= N, LB2 >= N) [in] 行列ペアの2つ目の行列. [out] B()は一般化シュール形Tで上書きされる. [in] Sdim Sort = "N": Sdim = 0. Sort = "S": Sdim = Selctgがtrueであった固有値(並べ替え後)の数. [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), j = 0, ..., N-1 は複素シュール形(S, T)の対角要素である. Beta(j) は非負の実数である. 注: 商Alpha(j)/Beta(j)はオーバーフローあるいはアンダーフローし易く, Beta(j)が0になることさえある. 従って, 比α/βを単純に計算することは避けなければならない. しかし, Alphaは常にnorm(A)より小さく, 通常は同程度の大きさである. また, Betaは常にnorm(B)より小さく, 通常は同程度の大きさである. [out] Vsl() 配列 Vsl(LVsl1 - 1, LVsl2 - 1) (LVsl1 >= N, LVsl2 >= N) Jobvsl = "V": Vsl()に左シュールベクトルが入る. Jobvsl = "N": 参照されない. [out] Vsr() 配列 Vsr(LVsr1 - 1, LVsr2 - 1) (LVsr1 >= N, LVsr2 >= N) Jobvsr = "V": Vsr()に右シュールベクトルが入る. Jobvsr = "N": 参照されない. [out] RConde() 配列 RConde(LRConde - 1) (LRConde >= 2) Sense = "E"または"B": RConde(0)およびRConde(1)には選択された固有値の平均についての条件数の逆数が入る. Sense = "N"または"V": 参照されない. [out] RCondv() 配列 RCondv(LRCondv - 1) (LRCondv >= 2) Sense = "V"または"B": RCondv(0)およびRCondv(1)には選択された部分空間についての条件数の逆数が入る. Sense = "N"または"E": 参照されない. [out] Info = 0: 正常終了. = -1: パラメータ Jobvl の誤り. (Jobvl <> "V"および"N") = -2: パラメータ Jobvr の誤り. (Jobvr <> "V"および"N") = -3: パラメータ Sort の誤り. (Sort <> "S"および"N") = -5: パラメータ Sense の誤り. (Sense <> "E", "V", "B"および"N") = -6: パラメータ N の誤り. (N < 0) = -7: パラメータ A() の誤り. = -8: パラメータ B() の誤り. = -10: パラメータ Alpha() の誤り. = -11: パラメータ Beta() の誤り. = -12: パラメータ Vsl() の誤り. = -13: パラメータ Vsr() の誤り. = -14: パラメータ RConde() の誤り. = -15: パラメータ RCondv() の誤り. = i (0 < i <= N): QZ反復が失敗した. (A, B)はシュール形ではないが, Alpha(j)およびBeta(j) (j = i〜N-1) は正しい. = N+1: ZhgeqzにおいてQZ反復の失敗以外のエラーが起きた. = N+2: 並べ替え後, 丸め誤差により値が変化した複素固有値があり, シュール形の先頭側のいくつかの固有値が Selctg = true を満たさなくなった. これは, スケーリングによっても起こりうる. = N+3: Ztgsenにおいて並べ替えに失敗した. 出典 LAPACK 使用例 行列のペア(A, B)の一般化固有値, 一般化シュール形(S, T), および, 左および右シュールベクトルを求める. ただし, ( 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_Zggesx() 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, Sdim As Long Dim Vsl(N - 1, N - 1) As Complex, Vsr(N - 1, N - 1) As Complex Dim RConde(1) As Double, RCondv(1) As Double, Info As Long 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 Zggesx("V", "V", "S", AddressOf Selctg, "B", N, A(), B(), Sdim, Alpha(), Beta(), Vsl(), Vsr(), 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 "Schur form S =" Debug.Print Creal(A(0, 0)), Cimag(A(0, 0)), Creal(A(0, 1)), Cimag(A(0, 1)) Debug.Print Creal(A(1, 0)), Cimag(A(1, 0)), Creal(A(1, 1)), Cimag(A(1, 1)) Debug.Print Creal(A(2, 0)), Cimag(A(2, 0)), Creal(A(2, 1)), Cimag(A(2, 1)) Debug.Print Creal(A(0, 2)), Cimag(A(0, 2)) Debug.Print Creal(A(1, 2)), Cimag(A(1, 2)) Debug.Print Creal(A(2, 2)), Cimag(A(2, 2)) Debug.Print "Schur form T =" Debug.Print Creal(B(0, 0)), Cimag(B(0, 0)), Creal(B(0, 1)), Cimag(B(0, 1)) Debug.Print Creal(B(1, 0)), Cimag(B(1, 0)), Creal(B(1, 1)), Cimag(B(1, 1)) Debug.Print Creal(B(2, 0)), Cimag(B(2, 0)), Creal(B(2, 1)), Cimag(B(2, 1)) Debug.Print Creal(B(0, 2)), Cimag(B(0, 2)) Debug.Print Creal(B(1, 2)), Cimag(B(1, 2)) Debug.Print Creal(B(2, 2)), Cimag(B(2, 2)) Debug.Print "Left Schur vectors =" Debug.Print Creal(Vsl(0, 0)), Cimag(Vsl(0, 0)), Creal(Vsl(0, 1)), Cimag(Vsl(0, 1)) Debug.Print Creal(Vsl(1, 0)), Cimag(Vsl(1, 0)), Creal(Vsl(1, 1)), Cimag(Vsl(1, 1)) Debug.Print Creal(Vsl(2, 0)), Cimag(Vsl(2, 0)), Creal(Vsl(2, 1)), Cimag(Vsl(2, 1)) Debug.Print Creal(Vsl(0, 2)), Cimag(Vsl(0, 2)) Debug.Print Creal(Vsl(1, 2)), Cimag(Vsl(1, 2)) Debug.Print Creal(Vsl(2, 2)), Cimag(Vsl(2, 2)) Debug.Print "Right Schur vectors =" Debug.Print Creal(Vsr(0, 0)), Cimag(Vsr(0, 0)), Creal(Vsr(0, 1)), Cimag(Vsr(0, 1)) Debug.Print Creal(Vsr(1, 0)), Cimag(Vsr(1, 0)), Creal(Vsr(1, 1)), Cimag(Vsr(1, 1)) Debug.Print Creal(Vsr(2, 0)), Cimag(Vsr(2, 0)), Creal(Vsr(2, 1)), Cimag(Vsr(2, 1)) Debug.Print Creal(Vsr(0, 2)), Cimag(Vsr(0, 2)) Debug.Print Creal(Vsr(1, 2)), Cimag(Vsr(1, 2)) Debug.Print Creal(Vsr(2, 2)), Cimag(Vsr(2, 2)) Debug.Print "Rconde =", RConde(0), RConde(1) Debug.Print "Rcondv =", RCondv(0), RCondv(1) Debug.Print "Sdim =", Sdim, "Info =", Info End Sub Function Selctg(Alpha As Complex, Beta As Complex) As Long Selctg = 0 If Cimag(Alpha) <> 0 Then Selctg = 1 End Function 実行結果 Eigenvalues = -0.4784814787767 -0.640760182519056 1.62433774044009 1.10642263894432 -0.149178317709927 5.63446258373312 Schur form S = -0.72797081435173 -0.974864717993083 -0.273608952637649 0.236841131861002 0 0 1.50621819846239 1.02596515027546 0 0 0 0 1.07441516363149E-02 4.41155734877746E-02 -0.689834642053671 -0.367063708571792 -2.75225804680911E-02 1.03952740743988 Schur form T = 1.5214190029108 0 -0.330233470430214 -1.05504057337832 0 0 0.927281415042602 0 0 0 0 0 0.359567173199973 1.92492599356072 -0.911905206658519 -8.20406969051636E-02 0.184494508924602 0 Left Schur vectors = -0.129055100819468 0.380797390644917 0.459273707935736 -0.209026938555656 -0.693189782060827 0.176375779993619 -0.368487931998788 -0.564073489521078 0.149492866673543 0.551696946993847 0.512426331896747 -0.16980105101744 0.759605142163285 8.19362946843842E-02 -0.150340091122171 0.108697330557172 -0.614629176696597 6.40327579983544E-02 Right Schur vectors = -0.609843105250231 -0.453993250889564 0.531526779080892 0.119433675168335 -0.264988296049601 -9.90695388688886E-02 -0.740613399363287 -0.14310787368587 0.289272174294141 -0.508202276391652 -5.56309308836406E-02 -0.362121539068532 -0.265832927214418 -0.233514989432741 -0.344926967779475 -0.481667671217571 -0.51570410282413 0.507813473877466 Rconde = 1 1 Rcondv = 3.96614422329799 3.96614422329799 Sdim = 3 Info = 0