Dgges Dgges_r Dggesx Dggesx_r Dggev Dggevx

Dgges_r() Sub Dgges_r ( Jobvsl As  String, Jobvsr As  String, Sort As  String, N As  Long, A() As  Double, B() As  Double, Sdim As  Long, Alphar() As  Double, Alphai() As  Double, Beta() As  Double, Vsl() As  Double, Vsr() As  Double, Info As  Long, IRev As  Long, Selct() As  Long  ) (シンプルドライバ) 一般化シュール分解 (一般行列) (リバースコミュニケーション版) 目的 本ルーチンはn×n一般実行列のペア(A, B)の一般化固有値, 一般化実シュール形(S, T), および, 必要により左および／または右シュールベクトルからなる行列(VSLおよびVSR)を求める. これらにより, 一般化シュール分解が次のように与えられる. (A, B) = ((VSL)*S*(VSR)^T, (VSL)*T*(VSR)^T) また, 必要により, 選択された固有値が上準三角行列Sおよび上三角行列Tの左上の対角ブロックにくるように固有値を並べ替える. その結果, VSLおよびVSRの先行する側の列は, 対応する左および右固有空間(部分空間)の正規直交基底を形成する. (一般化固有値のみ必要ならば代わりにより高速なDggevを使用せよ.) 行列のペア(A, B)の一般化固有値は, A - wB が特異となるようなスカラーwあるいは比α/β = wである. β = 0 あるいは両方共0の場合に適当な解釈ができるため, これは通常ペア(α, β)で表される. Tが非負の対角成分を持つ上三角行列で, Sが1×1および2×2ブロックよりなるブロック上三角行列であれば, 行列のペア(S, T)は一般化実シュール形である. 1×1ブロックは実一般化固有値に対応し, 一方, Sの2×2ブロックはTの対応する要素が次の形式になるように"標準化"される. [ a 0 ] [ 0 b ] また, SおよびTの対応する2×2ブロックの対は複素共役対の一般化固有値である. Dgges_rはDggesのリバースコミュニケーション版である. 引数 [in] Jobvsl = "N": 左シュールベクトルを求めない. = "V": 左シュールベクトルを求める. [in] Jobvsr = "N": 右シュールベクトルを求めない. = "V": 右シュールベクトルを求める. [in] Sort 一般化シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する. = "N": 固有値の並べ替えを行わない. = "S": 固有値の並べ替えを行う(Selctgを参照のこと). [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であった固有値(並べ替え後)の数. (どちらかの固有値についてSelct(i)がtrueであった複素共役対は2とカウントする.) [out] Alphar() 配列 Alphar(LAlphar - 1) (LAlphar >= N) [out] Alphai() 配列 Alphai(LAlphai - 1) (LAlphai >= N) [out] Beta() 配列 Beta(LBeta - 1) (LBeta >= N) (Alphar(j) + Alphai(j)*i)/Beta(j) (j = 0〜N-1) は一般化固有値である. Alphar(j) + Alphai(j)*i および Beta(j), j=0, ..., N-1 は複素シュール形(S, T)の対角成分である. (S, T)は, (A, B)の実シュール形の2×2対角ブロックを2×2複素ユニタリ変換によりさらに三角型に変換した場合に得られる. Alphai(j)が0の場合, j番目の固有値は実数である. 正数の場合, j番目とj+1番目の固有値は複素共役対であり, Alphai(j+1)は負である. 注: 商Alphar(j)/Beta(j)およびAlphai(j)/Beta(j)はオーバーフローあるいはアンダーフローし易く, Beta(j)が0になることさえある. 従って, 比α/βを単純に計算することは避けなければならない. しかし, AlpharとAlphaiは常に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: パラメータ Alphar() の誤り. = -9: パラメータ Alphai() の誤り. = -10: パラメータ Beta() の誤り. = -11: パラメータ Vsl() の誤り. = -12: パラメータ Vsr() の誤り. = -15: パラメータ Selct() の誤り. = i (0 < i <= N): QZ反復が失敗した. (A, B)はシュール形ではないが, Alphar(j), Alphai(j), Beta(j) (j = i〜N-1) は正しい. = N+1: DhgeqzにおいてQZ反復の失敗以外のエラーが起きた. = N+2: 並べ替え後, 丸め誤差により値が変わった複素固有値があり, シュール形の先頭のいくつかの固有値が Selct(i) = true を満たさなくなった. これはスケーリングによっても起こりうる. = N+3: Dtgsenにおいて並べ替えに失敗した. [in,out] IRev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 結果はInfoをチェックすること. = 1: Sort = "S"の場合, シュール形の左上にくるように並べ替えを行う固有値の選択のための判定値を Selct(i) (i = 0 To N-1) に設定する. Alphar(i), Alphai(i)およびBeta(i) (Alphar(i)/Beta(i)とAlphai(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.20 -0.11 -0.93 ) ( -0.58 -0.79 0.82 ) A = ( -0.32 0.81 0.37 ), B = ( 0.77 0.71 -0.55 ) ( -0.80 -0.92 -0.29 ) ( -1.36 -1.22 1.66 ) とする. Sub Ex_Dgges_r() Const N = 3 Dim A(N - 1, N - 1) As Double, B(N - 1, N - 1) As Double Dim Alphar(N - 1) As Double, Alphai(N - 1) As Double, Beta(N - 1) As Double Dim Sdim As Long, Vsl(N - 1, N - 1) As Double, Vsr(N - 1, N - 1) As Double Dim IRev As Long, Selct(N - 1) As Long, Info As Long A(0, 0) = 0.2: A(0, 1) = -0.11: A(0, 2) = -0.93 A(1, 0) = -0.32: A(1, 1) = 0.81: A(1, 2) = 0.37 A(2, 0) = -0.8: A(2, 1) = -0.92: A(2, 2) = -0.29 B(0, 0) = -0.58: B(0, 1) = -0.79: B(0, 2) = 0.82 B(1, 0) = 0.77: B(1, 1) = 0.71: B(1, 2) = -0.55 B(2, 0) = -1.36: B(2, 1) = -1.22: B(2, 2) = 1.66 IRev = 0 Do Call Dgges_r("V", "V", "S", N, A(), B(), Sdim, Alphar(), Alphai(), Beta(), Vsl(), Vsr(), Info, IRev, Selct()) If IRev = 1 Then Call Selctg_r(Alphar(), Alphai(), Beta(), Selct()) Loop While IRev <> 0 Debug.Print "Eigenvalues =" Debug.Print " (r)", Alphar(0) / Beta(0), Alphar(1) / Beta(1), Alphar(2) / Beta(2) Debug.Print " (i)", Alphai(0) / Beta(0), Alphai(1) / Beta(1), Alphai(2) / Beta(2) Debug.Print "Schur form S =" Debug.Print A(0, 0), A(0, 1), A(0, 2) Debug.Print A(1, 0), A(1, 1), A(1, 2) Debug.Print A(2, 0), A(2, 1), A(2, 2) Debug.Print "Schur form T =" Debug.Print B(0, 0), B(0, 1), B(0, 2) Debug.Print B(1, 0), B(1, 1), B(1, 2) Debug.Print B(2, 0), B(2, 1), B(2, 2) Debug.Print "Left Schur vectors =" Debug.Print Vsl(0, 0), Vsl(0, 1), Vsl(0, 2) Debug.Print Vsl(1, 0), Vsl(1, 1), Vsl(1, 2) Debug.Print Vsl(2, 0), Vsl(2, 1), Vsl(2, 2) Debug.Print "Right Schur vectors =" Debug.Print Vsr(0, 0), Vsr(0, 1), Vsr(0, 2) Debug.Print Vsr(1, 0), Vsr(1, 1), Vsr(1, 2) Debug.Print Vsr(2, 0), Vsr(2, 1), Vsr(2, 2) Debug.Print "Sdim =", Sdim, "Info =", Info End Sub Sub Selctg_r(Alphar() As Double, Alphai() As Double, Beta() As Double, Selct() As Long) Const N = 3 Dim I As Long For I = 0 To N - 1 Selct(I) = 0 If Alphai(I) <> 0 Then Selct(I) = 1 Next End Sub 実行結果 Eigenvalues = (r) 1.10765548065266 1.10765548065266 -1.74329621491312 (i) 1.40252005479606 -1.40252005479606 0 Schur form S = 0.61463908695189 -1.30352251020676 6.87282197451471E-02 0.715639689206334 0.422919612615492 0.382300104551802 0 0 -0.681226655778666 Schur form T = 1.62147306790067 -0 2.50284121079862 0 0.230317162252267 0.163903374557831 0 0 0.390769308136549 Left Schur vectors = -0.352844167536935 0.427271157509602 -0.832430388318189 0.480836332469667 -0.680393780608842 -0.553046765369735 -0.802681390597968 -0.595402100491098 3.46254807692483E-02 Right Schur vectors = 0.988795256019328 0.1471244329068 2.52654490465408E-02 0.106663051199989 -0.814736341679618 0.569936564062371 0.104436273307598 -0.560855680891587 -0.821300170479256 Sdim = 2 Info = 0