Zgges Zgges_r Zggesx Zggesx_r Zggev Zggevx

Zgges_r() Sub Zgges_r ( Jobvsl As  String, Jobvsr As  String, Sort 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, Info As  Long, IRev As  Long, Selct() 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の左上の対角ブロックにくるように固有値を並べ替える. その結果, VSLおよびVSRの先行する側の列は, 対応する左および右固有空間(部分空間)のユニタリ基底を形成する. (一般化固有値のみ必要ならば代わりにより高速なZggevを使用せよ) 行列のペア(A, B)の一般化固有値は, A - wB が特異となるようなスカラーwあるいは比α/β = wである. β = 0 あるいは両方共0の場合に適当な解釈ができるため, これは通常ペア(α, β)で表される. SおよびTが上三角行列で, さらに, Tの対角要素が非負の実数であれば, 行列のペア(S, T)は一般化複素シュール形である. Zgges_rはZggesのリバースコミュニケーション版である. 引数 [in] Jobvsl = "N": 左シュールベクトルを求めない. = "V": 左シュールベクトルを求める. [in] Jobvsr = "N": 右シュールベクトルを求めない. = "V": 右シュールベクトルを求める. [in] Sort 一般化シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する. = "N": 固有値の並べ替えを行わない. = "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 = Selct(i)がtrueであった固有値(並べ替え後)の数. [out] Alphar() 配列 Alphar(LAlphar - 1) (LAlphar >= N) [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] Info = 0: 正常終了. = -1: パラメータ Jobvl の誤り. (Jobvl <> "V"および"N") = -2: パラメータ Jobvr の誤り. (Jobvr <> "V"および"N") = -3: パラメータ Sort の誤り. (Sort <> "S"および"N") = -4: パラメータ N の誤り. (N < 0) = -5: パラメータ A() の誤り. = -6: パラメータ B() の誤り. = -8: パラメータ Alpha() の誤り. = -9: パラメータ Beta() の誤り. = -10: パラメータ Vsl() の誤り. = -11: パラメータ Vsr() の誤り. = -14: パラメータ Selct() の誤り. = i (0 < i <= N): QZ反復が失敗した. (A, B)はシュール形ではないが, Alpha(j), Beta(j) (j = i〜N-1) は正しい. = N+1: ZhgeqzにおいてQZ反復の失敗以外のエラーが起きた. = N+2: 並べ替え後, 丸め誤差により値が変わった複素固有値があり, シュール形の先頭のいくつかの固有値が Selct(i) = true を満たさなくなった. これはスケーリングによっても起こりうる. = N+3: Ztgsenにおいて並べ替えに失敗した. [in,out] IRev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 結果はInfoをチェックすること. = 1: Sort = "S"の場合, シュール形の左上にくるように並べ替えを行う固有値の選択のための判定値を Selct(i) (i = 0 To N-1) に設定する. Alpha(i)およびBeta(i)(Alpha(i)/Beta(i)は固有値)を参照して判定を行い, 選択する場合 Selct(i) = true (1), 選択しない場合 Selct(i) = false (0) とする. Selct()以外の変数を変更してはならない.   Sort = "N"の場合, 常に IRev = 0 を返す. [in] Selct() 配列 Selct(LSelct - 1) (LSelct >= N) IRev = 1 の場合, 並べ替えを行う固有値の選択のための判定値を入力する. 出典 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_Zgges_r() 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 IRev As Long, Selct(N - 1) As Long, 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) IRev = 0 Do Call Zgges_r("V", "V", "S", N, A(), B(), Sdim, Alpha(), Beta(), Vsl(), Vsr(), Info, IRev, Selct()) If IRev = 1 Then Call Selctg_r(Alpha(), Beta(), Selct()) Loop While IRev <> 0 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 "Sdim =", Sdim, "Info =", Info End Sub Sub Selctg_r(Alpha() As Complex, Beta() As Complex, Selct() As Long) Const N = 3 Dim I As Long For I = 0 To N - 1 Selct(I) = 0 If Cimag(Alpha(I)) <> 0 Then Selct(I) = 1 Next End Sub 実行結果 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 Sdim = 3 Info = 0