Dgebak Dgebal Dgehrd Dhsein Dhseqr Dorghr Dormhr Dtrevc3 Dtrexc Dtrsen Dtrsna Dtrsyl

Dtrevc3() Sub Dtrevc3 ( Side As  String, Howmny As  String, Selct() As  Boolean, N As  Long, T() As  Double, Vl() As  Double, Vr() As  Double, Mm As  Long, M As  Long, Info As  Long  ) シュール分解の準三角行列の固有ベクトル 目的 本ルーチンは実上準三角行列 T の一部または全部の右または左固有ベクトルを求める. このタイプの行列は一般実行列のシュール分解 A = Q*T*Q^T (Dhseqr により計算) により生成される. T の固有値 w に対応する右固有ベクトル x および左固有ベクトル y は次のように定義される. T*x = w*x, (y^T)*T = w*(y^T) ここで, y^T はベクトル 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() は求める固有ベクトルを指定する. W(j) が実数固有値の場合, Selct(j) が True であれば対応する実数固有ベクトルが求められる. W(j) と W(j+1) が複素固有値の実数部と虚数部の場合, Selct(j) あるいは Selct(j+1) のどちらかが True であれば対応する複素固有ベクトルが求められ, 終了時 Selct(j) は True に, また, Selct(j+1) は False に設定される. Howmny = "A" または "B" の場合, 参照されない. [in] N 行列 A の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] T() 配列 T(LT1 - 1, LT2 - 1) (LT1 >= N, LT2 >= N) シュール正準形の上準三角行列 T. [in,out] Vl() 配列 Vl(LVl1 - 1, LVl2 - 1) (LVl1 >= N, LVl2 >= MM) [in] Side = "L" または "B" かつ Howmny = "B" の場合, Vl() に N x N 行列 Q を入れる (通常は Dhseqr により返されたシュールベクトルからなる直交行列 Q). [out] Side = "L" または "B" の場合, Vl() に以下を返す. Howmny = "A" の場合, T の左固有ベクトルからなる行列 Y. Howmny = "B" の場合, 行列 Q*Y. Howmny = "S" の場合, Selct() により指定された T の左固有ベクトルが, 固有値と同じ順に Vl() の列に続けて格納される. 複素固有値に対応する複素固有ベクトルは, 最初は実数部, 2番目は虚数部というように2つの列に続けて格納される. 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 を入れる (通常は Dhseqr により返されたシュールベクトルからなる直交行列 Q). [out] Side = "R" または "B" の場合, Vr() に以下を返す. Howmny = "A" の場合, T の右固有ベクトルからなる行列 X. Howmny = "B" の場合, 行列 Q*X. Howmny = "S" の場合, Selct() により指定された T の右固有ベクトルが, 固有値と同じ順に Vr() の列に続けて格納される. 複素固有値に対応する複素固有ベクトルは, 最初は実数部, 2番目は虚数部というように2つの列に続けて格納される. Side = "L" の場合, 参照されない. [in] MM 配列 Vl() および／または Vr() の列数. (MM >= M) [out] M 配列 Vl() および／または Vr() において実際に固有ベクトルを格納するのに使われた列数. Howmny = "A" または "B" の場合, M = N と設定される. 選択された実数固有ベクトルはそれぞれ1列を占め, 選択された複素固有ベクトルはそれぞれ2列を占める. [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.11 -0.93 ) A = ( -0.32 0.81 0.37 ) ( -0.80 -0.92 -0.29 ) とする. Dgehrdでヘッセンベルグ形に変換したのち, DhseqrおよびDorghrで固有値およびシュール形を求め, Dtrevc3でシュール形より固有ベクトルを求める. Sub Ex_Dgehrd_Dhseqr_Dtrevc3() Const N = 3 Dim A(N - 1, N - 1) As Double, Tau(N - 2) As Double Dim Wr(N - 1) As Double, Wi(N - 1) As Double, Selct() As Boolean Dim Vl(N - 1, N - 1) As Double, Vr(N - 1, N - 1) As Double, M As Long Dim Ilo As Long, Ihi As Long, I As Long, J As Long, 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 Ilo = 1: Ihi = N Call Dgehrd(N, Ilo, Ihi, A(), Tau(), Info) If Info <> 0 Then Debug.Print "Error in Dgehrd: 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 Dorghr(N, Ilo, Ihi, Vr(), Tau(), Info) If Info <> 0 Then Debug.Print "Error in Dorghr: Info =", Info Exit Sub End If Call Dhseqr("S", "V", N, Ilo, Ihi, A(), Wr(), Wi(), Vr(), Info) If Info <> 0 Then Debug.Print "Error in Dhseqr: 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 Dtrevc3("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 (r) =", Wr(0), Wr(1), Wr(2) Debug.Print "Eigenvalues (i) =", Wi(0), Wi(1), Wi(2) Debug.Print "Eigenvectors (L) =" Debug.Print Vl(0, 0), Vl(0, 1), Vl(0, 2) Debug.Print Vl(1, 0), Vl(1, 1), Vl(1, 2) Debug.Print Vl(2, 0), Vl(2, 1), Vl(2, 2) Debug.Print "Eigenvectors (R) =" Debug.Print Vr(0, 0), Vr(0, 1), Vr(0, 2) Debug.Print Vr(1, 0), Vr(1, 1), Vr(1, 2) Debug.Print Vr(2, 0), Vr(2, 1), Vr(2, 2) Debug.Print "M =", M End Sub 実行結果 Eigenvalues (r) = -0.904130023851345 0.812065011925673 0.812065011925673 Eigenvalues (i) = 0 0.48915757543818 -0.48915757543818 Eigenvectors (L) = -0.89678255062032 -0.426130142157969 -0.33620818638403 -0.594264184393386 0.381115511928175 -0.618884488071825 -1 0.379151460023587 0.244224395790644 Eigenvectors (R) = -0.83634228919053 -0.578201236867243 -0.117907376795369 5.97215298924774E-02 0.420188120674838 -0.579811879325162 -1 0.268818029145007 0.450298711668266 M = 3