◆ Ztrevc3()

Sub Ztrevc3	(	Side As	String,
		Howmny As	String,
		Selct() As	Boolean,
		N As	Long,
		T() As	Complex,
		Vl() As	Complex,
		Vr() As	Complex,
		Mm As	Long,
		M As	Long,
		Info As	Long
	)

シュール分解の三角行列の固有ベクトル (複素行列)

目的: 本ルーチンは上三角複素行列 T の一部または全部の右または左固有ベクトルを求める. このタイプの行列は一般複素行列のシュール分解 A = Q*T*Q^H (Zhseqr により計算) により生成される.

T の固有値 w に対応する右固有ベクトル x および左固有ベクトル y は次のように定義される.
T*x = w*x, (y^H)*T = w*(y^H)

ここで, y^H はベクトル y の共役転置を表す. 固有値は本ルーチンに入力されないが, T の対角ブロックより直接読み込まれる.

本ルーチンは T の右および左固有ベクトル x および／または y を返す. あるいは, 入力された行列 Q との積 Q*X および／または Q*Y を返す. Q が行列 A をシュール形 T に変換するユニタリ行列であれば Q*X および Q*Y は A の右および左固有ベクトルである.

本ルーチンは逆変換の Level 3 BLAS 版を使用する.

引数

[in]	Side	= "R": 右固有ベクトルのみ求める. = "L": 左固有ベクトルのみ求める. = "B": 右および左固有ベクトルの両方を求める.
[in]	Howmny	= "A": すべての右および／または左固有ベクトルを求める. = "B": Vl() および／または Vr() に入っている行列により逆変換されたすべての右および／または左固有ベクトルを求める. = "S": 論理型配列 Selct() により指定された一部の右および／または左固有ベクトルを求める.
[in,out]	Selct()	配列 Selct(LSelct - 1) (LSelct >= N) Howmny = "S" の場合, Selct() は求める固有ベクトルを指定する. Selct(j) = True であれば j 番目の固有値に対応する固有ベクトルが求められる. Howmny = "A" または "B" の場合, 参照されない.
[in]	N	行列 A の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in]	T()	配列 T(LT1 - 1, LT2 - 1) (LT1 >= N, LT2 >= N) 上三角行列 T. T() は変更されるが終了時に復元される.
[in,out]	Vl()	配列 Vl(LVl1 - 1, LVl2 - 1) (LVl1 >= N, LVl2 >= MM) [in] Side = "L" または "B" かつ Howmny = "B" の場合, Vl() に N x N 行列 Q を入れる (通常は Zhseqr により返されたシュールベクトルからなるユニタリ行列 Q). [out] Side = "L" または "B" の場合, Vl() に以下を返す. Howmny = "A" の場合, T の左固有ベクトルからなる行列 Y. Howmny = "B" の場合, 行列 Q*Y. Howmny = "S" の場合, Selct() により指定された T の左固有ベクトルが, 固有値と同じ順に Vl() の列に続けて格納される. Side = "R" の場合, 参照されない.
[in,out]	Vr()	配列 Vr(LVr1 - 1, LVr2 - 1) (LVr1 >= N, LVr2 >= MM) [in] Side = "R" または "B" かつ Howmny = "B" の場合, Vr() に N x N 行列 Q を入れる (通常は Zhseqr により返されたシュールベクトルからなるユニタリ行列 Q). [out] Side = "R" または "B" の場合, Vr() に以下を返す. Howmny = "A" の場合, T の右固有ベクトルからなる行列 X. Howmny = "B" の場合, 行列 Q*X. Howmny = "S" の場合, Selct() により指定された T の右固有ベクトルが, 固有値と同じ順に Vr() の列に続けて格納される. Side = "L" の場合, 参照されない.
[in]	MM	配列 Vl() および／または Vr() の列数. (MM >= M)
[out]	M	配列 Vl() および／または Vr() において実際に固有ベクトルを格納するのに使われた列数. Howmny = "A" または "B" の場合, M = N と設定される. 選択された固有ベクトルはそれぞれ 1 列を占める.
[out]	Info	= 0: 正常終了. = -1: パラメータ Side の誤り. (Side <> "R", "L" および "B") = -2: パラメータ Howmny の誤り. (Howmny <> "A", "B" および "S") = -3: パラメータ Selct() の誤り. = -4: パラメータ N の誤り. (N < 0) = -5: パラメータ T() の誤り. = -6: パラメータ Vl() の誤り. = -7: パラメータ Vr() の誤り. = -8: パラメータ MM の誤り. (MM < M)

詳細: 本プログラムで使用されているアルゴリズムは基本的には後退(前進)代入である. 起こりうるオーバーフローに対して強くするためにスケーリングを行う.

固有ベクトルはそれぞれ最大要素の大きさが 1 になるように正規化される. ここで, 複素数 (x, y) の大きさは |x| + |y| とする.

出典: LAPACK

使用例: 行列Aの固有値・固有ベクトルを求める. ただし,
( 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 )

とする.
Zgehrdでヘッセンベルグ形に変換したのち, ZhseqrおよびZunghrで固有値およびシュール形を求め, Ztrevc3でシュール形より固有ベクトルを求める.
Sub Ex_Zgehrd_Zhseqr_Ztrevc3()

Const N = 3

Dim A(N - 1, N - 1) As Complex, Tau(N - 2) As Complex

Dim W(N - 1) As Complex, Selct() As Boolean

Dim Vl(N - 1, N - 1) As Complex, Vr(N - 1, N - 1) As Complex, M As Long

Dim Ilo As Long, Ihi As Long, I As Long, J 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)

Ilo = 1: Ihi = N

Call Zgehrd(N, Ilo, Ihi, A(), Tau(), Info)

If Info <> 0 Then

Debug.Print "Error in Zgehrd: Info =", Info

Exit Sub

End If

For I = 0 To N - 1

For J = 0 To N - 1

Vr(I, J) = A(I, J)

Next

Next

Call Zunghr(N, Ilo, Ihi, Vr(), Tau(), Info)

If Info <> 0 Then

Debug.Print "Error in Zunghr: Info =", Info

Exit Sub

End If

Call Zhseqr("S", "V", N, Ilo, Ihi, A(), W(), Vr(), Info)

If Info <> 0 Then

Debug.Print "Error in Zhseqr: Info =", Info

Exit Sub

End If

For I = 0 To N - 1

For J = 0 To N - 1

Vl(I, J) = Vr(I, J)

Next

Next

Call Ztrevc3("B", "B", Selct(), N, A(), Vl(), Vr(), N, M, Info)

If Info <> 0 Then

Debug.Print "Error in Dtrevc3: Info =", Info

Exit Sub

End If

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 "Eigenvectors (L) ="

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

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

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

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

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

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

Debug.Print "Eigenvectors (R) ="

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

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

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

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

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

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

Debug.Print "M =", M

End Sub

実行結果: Eigenvalues =

-1.15894122423918 -0.50662892448174 1.05593587167591 0.900255855387815

0.21300535256327 1.29637306909393

Eigenvectors (L) =

-0.36286784799871 0.63713215200129 -0.144564860055844 0.855435139944156

-0.395771083048871 0.240186600504172 -1.15502382771539E-02 -0.649951770460259

0.434572937991399 -0.277594817887723 -0.144369992549788 0.399622008123156

1.19529397618877E-02 0.453615415642407

2.46422882501396E-02 -0.784859171023358

0.404424115184401 -0.595575884815599

Eigenvectors (R) =

-0.505480633843522 0.494519366156478 4.45520906484036E-03 0.657348688695353

-0.419909506510805 0.574781768215114 1.90776326854836E-02 -0.276338131107688

0.28193442840661 -0.216881788717049 -0.418424064044568 0.581575935955432

0.105185244310007 -0.110887166536361

0.164796997542067 -0.65350282278796

0.340224561176925 -0.659775438823075

M = 3