Dgges Dgges_r Dggesx Dggesx_r Dggev Dggevx

Dggesx_r() Sub Dggesx_r ( Jobvsl As  String, Jobvsr As  String, Sort As  String, Sense 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, 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)^T, (VSL)*T*(VSR)^T) また, 必要により, 選択された固有値が上準三角行列Sおよび上三角行列Tの左上の対角ブロックにくるように固有値を並べ替え, 選択された固有値の平均の条件数の逆数(rconde), および, 選択された固有値に対応する右および左部分空間の条件数の逆数(rcondv)を求める. その結果VSLおよびVSRの先行する側の列は, 対応する左および右固有空間(部分空間)の正規直交基底を形成する. 行列のペア(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ブロックのペアは複素共役対の一般化固有値である. Dggesx_rはDggesxのリバースコミュニケーション版である. 引数 [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であった固有値(並べ替え後)の数. (どちらかの固有値について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] 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: パラメータ Alphar() の誤り. = -10: パラメータ Alphai() の誤り. = -11: パラメータ Beta() の誤り. = -12: パラメータ Vsl() の誤り. = -13: パラメータ Vsr() の誤り. = -14: パラメータ RConde() の誤り. = -15: パラメータ RCondv() の誤り. = -18: パラメータ 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_Dggesx_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 RConde(1) As Double, RCondv(1) As Double, Info As Long Dim IRev As Long, Selct(N - 1) 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 Dggesx_r("V", "V", "S", "B", N, A(), B(), Sdim, Alphar(), Alphai(), Beta(), Vsl(), Vsr(), RConde(), RCondv(), 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 "Rconde =", RConde(0), RConde(1) Debug.Print "Rcondv =", RCondv(0), RCondv(1) 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 Rconde = 0.633219382771826 0.547225257604206 Rcondv = 0.306813246166952 0.400117332298039 Sdim = 2 Info = 0