|
|
◆ Dggesx()
| Sub Dggesx |
( |
Jobvsl As |
String, |
|
|
Jobvsr As |
String, |
|
|
Sort As |
String, |
|
|
Selctg As |
LongPtr, |
|
|
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 |
|
) |
| |
(エキスパートドライバ) 一般化シュール分解 (一般行列)
- 目的
- 本ルーチンは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の対応する要素が次の形式になるように"標準化"される. また, SおよびTの対応する2×2ブロックのペアは複素共役対の一般化固有値である.
- 引数
-
| [in] | Jobvsl | = "N": 左シュールベクトルを求めない.
= "V": 左シュールベクトルを求める. |
| [in] | Jobvsr | = "N": 右シュールベクトルを求めない.
= "V": 右シュールベクトルを求める. |
| [in] | Sort | 一般化シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する.
= "N": 固有値の並べ替えを行わない.
= "S": 固有値の並べ替えを行う(Selctgを参照のこと). |
| [in] | Selctg | Sort = "S": シュール形の左上にくるように並べ替える固有値を選択するためにSelctgが使われる.
Selctg(Alphar(j), Alphai(j), Beta(j))がtrue (= 1)であれば, 固有値 (Alphar(j)± Alphai(j)*i)/Beta(j) が選択される. すなわち, 固有値の複素共役対のうち一つが選択されれば両方の複素固有値が選択される.
並べ替え後, 選択された複素固有値は Selctg(Alphar(j), Alphai(j), Beta(j)) = true を満たさなくなるかもしれないことに注意せよ. これは, 並べ替えによって複素固有値の値が変わる可能性があるからである(特に, 固有値が悪条件の場合). この場合, Info = N+2 に設定される.
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であった固有値(並べ替え後)の数. (どちらかの固有値についてSelctgが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")
= -5: パラメータ Sense の誤り. (Sense <> "E", "V", "B"および"N")
= -6: パラメータ N の誤り. (N < 0)
= -7: パラメータ A() の誤り.
= -8: パラメータ B() の誤り.
= -10: パラメータ Alphar() の誤り.
= -11: パラメータ Alphai() の誤り.
= -12: パラメータ Beta() の誤り.
= -13: パラメータ Vsl() の誤り.
= -14: パラメータ Vsr() の誤り.
= -15: パラメータ RConde() の誤り.
= -16: パラメータ RCondv() の誤り.
= i (0 < i <= N): QZ反復が失敗した. (A, B)はシュール形ではないが, Alphar(j), Alphai(j), Beta(j) (j = i〜N-1) は正しい.
= N+1: DhgeqzにおいてQZ反復の失敗以外のエラーが起きた.
= N+2: 並べ替え後, 丸め誤差により値が変化した複素固有値があり, シュール形の先頭側のいくつかの固有値が Selctg = true を満たさなくなった. これは, スケーリングによっても起こりうる.
= N+3: Dtgsenにおいて並べ替えに失敗した. |
- 出典
- 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()
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
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
Call Dggesx("V", "V", "S", AddressOf Selctg, "B", N, A(), B(), Sdim, Alphar(), Alphai(), Beta(), Vsl(), Vsr(), RConde(), RCondv(), Info)
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
Function Selctg(Alphar As Double, Alphai As Double, Beta As Double) As Long
Selctg = 0
If Alphai <> 0 Then Selctg = 1
End Function
Function Beta(A As Double, B As Double, Optional Info As Long) As Double ベータ関数 B(a, b)
Sub Dggesx(Jobvsl As String, Jobvsr As String, Sort As String, Selctg As LongPtr, 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) (エキスパートドライバ) 一般化シュール分解 (一般行列)
- 実行結果
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
|