◆ Zgees()

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

(シンプルドライバ) シュール分解 (複素行列)

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

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

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

引数

[in]	Jobvs	= "N": シュールベクトルを求めない. = "V": シュールベクトルを求める.
[in]	Sort	シュール形の対角要素上の固有値の並べ替えを行うかどうか指定する. = "N": 固有値の並べ替えを行わない. = "S": 固有値の並べ替えを行う(Selctを参照のこと).
[in]	Selct	Sort = "S": シュール形の左上にくるように並べ替える固有値を選択するためにSelctが使われる. Selct(W(j))が true (= 1) であれば, 固有値 W(j) が選択される. Sort = "N": Selctは参照されない.
[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が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") = -4: パラメータ N の誤り. (N < 0) = -5: パラメータ A() の誤り. = -7: パラメータ W() の誤り. = -8: パラメータ Vs() の誤り. = i (0 < i <= N): QRアルゴリズムが失敗しすべての固有値を求めることはできなかった. W()の要素 0〜Ilo-2 および i〜N-1 には収束した固有値が入る. Jobvs = "V"であれば, 部分的に収束したシュール形にAを変換する行列がVs()に入る. = N+1: 一部の固有値が近すぎて分離できないため固有値の並べ替えができなかった(問題が非常に悪条件である). = N+2: 並べ替え後, 丸め誤差により値が変わった複素固有値があり, シュール形の先頭のいくつかの固有値が Selct = true を満たさなくなった. これはスケーリングによるアンダーフローでも起こりうる.

出典: 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()

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, 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)

Call Zgees("V", "S", AddressOf Selct, N, A(), Sdim, W(), Vs(), Info)

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

Function Selct(W As Complex) As Long

Selct = 0

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

End Function

実行結果: 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