◆ Zgees_r()

Sub Zgees_r	(	Jobvs As	String,
		Sort As	String,
		N As	Long,
		A() As	Complex,
		Sdim As	Long,
		W() As	Complex,
		Vs() As	Complex,
		Info As	Long,
		IRev As	Long,
		Selct() As	Long
	)

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

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

また, 必要により, 選択された固有値が左上にくるように固有値をシュール形の対角要素上に並べる. その結果Zの先行する側の列は選択された固有値に対応する不変部分空間の直交基底を形成する.

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

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

引数

[in]	Jobvs	= "N": シュールベクトルを求めない. = "V": シュールベクトルを求める.
[in]	Sort	シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する. = "N": 固有値の並べ替えを行わない. = "S": 固有値の並べ替えを行う(IRevを参照のこと).
[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]	Info	= 0: 正常終了. = -1: パラメータ Jobvs の誤り. (Jobvs <> "V"および"N") = -2: パラメータ Sort の誤り. (Sort <> "S"および"N") = -3: パラメータ N の誤り. (N < 0) = -4: パラメータ A() の誤り. = -6: パラメータ W() の誤り. = -7: パラメータ Vs() の誤り. = -10: パラメータ 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_Zgees_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, 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 Zgees_r("V", "S", N, A(), Sdim, W(), Vs(), 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 "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

Sdim = 3 Info = 0