◆ Dgeesx_r()

Sub Dgeesx_r	(	Jobvs As	String,
		Sort As	String,
		Sense As	String,
		N As	Long,
		A() As	Double,
		Sdim As	Long,
		Wr() As	Double,
		Wi() As	Double,
		Vs() As	Double,
		RConde As	Double,
		RCondv As	Double,
		Info As	Long,
		IRev As	Long,
		Selct() As	Long
	)

(エキスパートドライバ) シュール分解 (一般行列) (リバースコミュニケーション版)

目的: 本ルーチンはn×n一般実行列Aの固有値, 実シュール形T, および, 必要によりシュールベクトルからなる行列Zを求める. これらにより, シュール分解 A = Z*T*Z^T が与えられる.

また, 必要により, 選択された固有値が左上にくるように固有値を実シュール形の対角要素上に並べ, 選択された固有値の平均の条件数の逆数(RConde), および, 選択された固有値に対応する右不変部分空間の条件数の逆数(RCondv)を求める. Zの先行する側の列はこの不変部分空間の直交基底を形成する.

条件数の逆数RCondeおよびRCondvについてのさらなる説明は LAPACK Users' Guide Third Edition Section 4.8.1 を参照のこと(これらはそれぞれsおよびsepと呼ばれている).

行列が1×1および2×2ブロックよりなる準上三角行列であれば実シュール形である. 2×2ブロックは次の形式に標準化できる.
[ a b ]

[ c a ]

ここで, b*c < 0 である. このようなブロックの固有値は a±sqrt(bc) である.

Dgeesx_rはDgeesxのリバースコミュニケーション版である.

引数

[in]	Jobvs	= "N": シュールベクトルを求めない. = "V": シュールベクトルを求める.
[in]	Sort	シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する. = "N": 固有値の並べ替えを行わない. = "S": 固有値の並べ替えを行う(Selctを参照のこと).
[in]	Sense	どの条件数の逆数を計算するかを決定する. = 'N': 条件数を計算しない. = 'E': 選択された固有値の平均について条件数を計算する. = 'V': 選択された右不変部分空間について条件数を計算する. = 'B': 両方について条件数を計算する. Sense = "E", "V"あるいは"B"の場合、Sort = "S"でなければならない.
[in]	N	行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in,out]	A()	配列 A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= N) [in] N×N行列 A. [out] A()はその実シュール形Tで上書きされる.
[out]	Sdim	Sort = "N": Sdim = 0. Sort = "S": Sdim = Selct(i)がtrueであった固有値(並べ替え後)の数. (どちらかの固有値についてSelct(i)がtrueであった複素共役対は2とカウントする)
[out]	Wr()	配列 Wr(LWr - 1) (LWr >= N)
[out]	Wi()	配列 Wi(LWi - 1) (LWi >= N) 求められた固有値がシュール形Tの出力の対角要素に現れるのと同じ順で, Wr()およびWi()に求められた固有値のそれぞれ実数部および虚数部が入る. 複素共役対の固有値は, 虚数部が正の固有値を先頭に, 続けて現れる.
[out]	Vs()	配列 Vs(LVs1 - 1, LVs2 - 1) (LVs1 >= N, LVs2 >= N) Jobvs = "V": Vs()にシュールベクトルからなる直交行列Zが入る. Jobvs = "N": Vs()は参照されない.
[out]	RConde	Sense = 'E'または'B': RCondeには選択された固有値の平均の条件数の逆数が入る. Sense = 'N'または'V': 参照されない.
[out]	RCondv	Sense = 'V'または'B': RCondvには選択された右不変部分空間の条件数の逆数が入る. Sense = 'N'または'E': 参照されない.
[out]	Info	= 0: 正常終了. = -1: パラメータ Jobvs の誤り. (Jobvs <> "V"および"N") = -2: パラメータ Sort の誤り. (Sort <> "S"および"N") = -3: パラメータ Sense の誤り. (Sense <> "E", "V", "B"および"N") = -4: パラメータ N の誤り. (N < 0) = -5: パラメータ A() の誤り. = -6: パラメータ Wr() の誤り. = -7: パラメータ Wi() の誤り. = -9: パラメータ Vs() の誤り. = -14: パラメータ Selct() の誤り. = i (0 < i <= N): QRアルゴリズムが失敗しすべての固有値を求めることはできなかった. Wr()およびWi()の要素 0〜Ilo-2 および i〜N-1 には収束した固有値が入る. Jobvs = "V"であれば, 部分的に収束したシュール形にAを変換する行列がVs()に入る. = N+1: 一部の固有値が近すぎて分離できないため固有値の並べ替えができなかった(問題が非常に悪条件である). = N+2: 並べ替え後, 丸め誤差により値が変わった複素固有値があり, シュール形の先頭のいくつかの固有値が Selct(i) = true を満たさなくなった. これはスケーリングによるアンダーフローでも起こりうる.
[in,out]	IRev	リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 結果はInfoをチェックすること. = 1: Sort = "S"の場合, シュール形の左上にくるように並べ替えを行う固有値の選択のための判定値を Selct(i) (i = 0 To N-1) に設定する. Wr(i)とWi(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の固有値, 実シュール形T, および, シュールベクトルを求める. ただし,
( 0.20 -0.11 -0.93 )

A = ( -0.32 0.81 0.37 )

( -0.80 -0.92 -0.29 )

とする.
Sub Ex_Dgeesx_r()

Const N = 3

Dim A(N - 1, N - 1) As Double, Wr(N - 1) As Double, Wi(N - 1) As Double

Dim Sdim As Long, Vs(N - 1, N - 1) As Double

Dim IRev As Long, Selct(N - 1) As Long

Dim RConde As Double, RCondv 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

IRev = 0

Do

Call Dgeesx_r("V", "S", "B", N, A(), Sdim, Wr(), Wi(), Vs(), RConde, RCondv, Info, IRev, Selct())

If IRev = 1 Then Call Selct_r(Wr(), Wi(), Selct())

Loop While IRev <> 0

Debug.Print "Eigenvalues (r) =", Wr(0), Wr(1), Wr(2)

Debug.Print "Eigenvalues (i) =", Wi(0), Wi(1), Wi(2)

Debug.Print "Schur form T ="

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 vectors ="

Debug.Print Vs(0, 0), Vs(0, 1), Vs(0, 2)

Debug.Print Vs(1, 0), Vs(1, 1), Vs(1, 2)

Debug.Print Vs(2, 0), Vs(2, 1), Vs(2, 2)

Debug.Print "Rconde =", RConde, "Rcondv =", RCondv

Debug.Print "Sdim =", Sdim, "Info =", Info

End Sub

Sub Selct_r(Wr() As Double, Wi() As Double, Selct() As Long)

Const N = 3

Dim I As Long

For I = 0 To N - 1

Selct(I) = 0

If Wi(I) <> 0 Then Selct(I) = 1

Next

End Sub

実行結果: Eigenvalues (r) = 0.812065011925672 0.812065011925672 -0.904130023851345

Eigenvalues (i) = 0.48915757543818 -0.48915757543818 0

Schur form T =

0.812065011925672 0.540472392276116 0.684902131153596

-0.442714812131086 0.812065011925672 -0.537914483752962

0 0 -0.904130023851345

Schur vectors =

-0.492366426634308 -0.620320586279211 0.610555216308549

0.867251026977165 -0.290139488675864 0.404592057902722

-7.38306042939846E-02 0.728712184164039 0.680828608770563

Rconde = 0.894479266547508 Rcondv = 1.47509658058209

Sdim = 2 Info = 0