Zgges Zgges_r Zggesx Zggesx_r Zggev Zggevx

Zggesx_r() Sub Zggesx_r ( Jobvsl As  String, Jobvsr As  String, Sort As  String, 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, 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の左上の対角ブロックにくるように固有値を並べ替え, 選択された固有値の平均の条件数の逆数(RConde), および, 選択された固有値に対応する右および左部分空間の条件数の逆数(RCondv)を求める. その結果VSLおよびVSRの先行する側の列は, 対応する左および右固有空間(部分空間)のユニタリ基底を形成する. 行列のペア(A, B)の一般化固有値は, A - wB が特異となるようなスカラーwあるいは比α/β = wである. β = 0 あるいは両方共0の場合に適当な解釈ができるため, これは通常ペア(α, β)で表される. Tが非負の対角要素を持つ上三角行列でSが上三角行列であれば, 行列のペア(S, T)は一般化複素シュール形である. Zggesx_rはZggesxのリバースコミュニケーション版である. 引数 [in] Jobvsl = "N": 左シュールベクトルを求めない. = "V": 左シュールベクトルを求める. [in] Jobvsr = "N": 右シュールベクトルを求めない. = "V": 右シュールベクトルを求める. [in] Sort 一般化シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する. = "N": 固有値の並べ替えを行わない. = "S": 固有値の並べ替えを行う(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 = Selct(i)が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") = -4: パラメータ Sense の誤り. (Sense <> "E", "V", "B"および"N") = -5: パラメータ N の誤り. (N < 0) = -6: パラメータ A() の誤り. = -7: パラメータ B() の誤り. = -9: パラメータ Alpha() の誤り. = -10: パラメータ Beta() の誤り. = -11: パラメータ Vsl() の誤り. = -12: パラメータ Vsr() の誤り. = -13: パラメータ RConde() の誤り. = -14: パラメータ RCondv() の誤り. = -17: パラメータ 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_Zggesx_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 RConde(1) As Double, RCondv(1) As Double, Info As Long Dim IRev As Long, Selct(N - 1) 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 Zggesx_r("V", "V", "S", "B", N, A(), B(), Sdim, Alpha(), Beta(), Vsl(), Vsr(), RConde(), RCondv(), 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 "Rconde =", RConde(0), RConde(1) Debug.Print "Rcondv =", RCondv(0), RCondv(1) 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 Rconde = 1 1 Rcondv = 3.96614422329799 3.96614422329799 Sdim = 3 Info = 0