◆ Zgeesx_r()

Sub Zgeesx_r	(	Jobvs As	String,
		Sort As	String,
		Sense As	String,
		N As	Long,
		A() As	Complex,
		Sdim As	Long,
		W() As	Complex,
		Vs() As	Complex,
		RConde As	Double,
		RCondv As	Double,
		Info As	Long,
		IRev As	Long,
		Selct() As	Long
	)

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

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

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

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

複素行列が上三角行列であればシュール形である.

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

引数

[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であった固有値の数.
[out]	W()	配列 W(LW - 1) (LW >= N) 求められた固有値がシュール形Tの出力の対角要素に現れるのと同じ順でW()に入る.
[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() の誤り. = -7: パラメータ W() の誤り. = -8: パラメータ Vs() の誤り. = -13: パラメータ Selct() の誤り. = i (0 < i <= N): QRアルゴリズムが失敗しすべての固有値を求めることはできなかった. W()の要素 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.11i -0.93-0.32i 0.81+0.37i )

A = ( -0.80-0.92i -0.29+0.86i 0.64+0.51i )

( 0.71+0.59i -0.15+0.19i 0.20+0.94i )

とする.
Sub Ex_Zgeesx_r()

Const N = 3

Dim A(N - 1, N - 1) As Complex, W(N - 1) As Complex

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

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

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

IRev = 0

Do

Call Zgeesx_r("V", "S", "B", N, A(), Sdim, W(), Vs(), RConde, RCondv, Info, IRev, Selct())

If IRev = 1 Then Call Selct_r(W(), Selct())

Loop While IRev <> 0

Debug.Print "Eigenvalues ="

Debug.Print Creal(W(0)), Cimag(W(0)), Creal(W(1)), Cimag(W(1))

Debug.Print Creal(W(2)), Cimag(W(2))

Debug.Print "Schur form T ="

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

Debug.Print Creal(Vs(0, 0)), Cimag(Vs(0, 0)), Creal(Vs(0, 1)), Cimag(Vs(0, 1))

Debug.Print Creal(Vs(1, 0)), Cimag(Vs(1, 0)), Creal(Vs(1, 1)), Cimag(Vs(1, 1))

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

Debug.Print Creal(Vs(0, 2)), Cimag(Vs(0, 2))

Debug.Print Creal(Vs(1, 2)), Cimag(Vs(1, 2))

Debug.Print Creal(Vs(2, 2)), Cimag(Vs(2, 2))

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

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

End Sub

Sub Selct_r(W() As Complex, Selct() As Long)

Const N = 3

Dim I As Long

For I = 0 To N - 1

Selct(I) = 0

If Cimag(W(I)) <> 0 Then Selct(I) = 1

Next

End Sub

実行結果: Eigenvalues =

-1.15894122423918 -0.50662892448174 1.05593587167591 0.900255855387815

0.21300535256327 1.29637306909393

Schur form T =

-1.15894122423918 -0.50662892448174 2.31358655627147E-02 -0.442808752291521

0 0 1.05593587167591 0.900255855387815

0 0 0 0

-0.644484909898823 -0.733094546116673

0.379301413619788 -0.286520422490734

0.21300535256327 1.29637306909393

Schur vectors =

-0.474826946245658 0.464530398931381 -0.103329315153953 0.628092734605153

-0.39444507925059 0.539925475846714 -9.30826818180923E-02 -0.290418461635993

0.264837176182039 -0.203729501266495 -0.367753983153771 0.605452152301986

1.03227229675376E-02 0.391748503946406

2.12814186266949E-02 -0.67781516115713

0.349266221164391 -0.514347515136026

Rconde = 1 Rcondv = 2.76522085257846

Sdim = 3 Info = 0