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

Dhsein() Sub Dhsein ( Side As  String, Eigsrc As  String, Initv As  String, Selct() As  Boolean, N As  Long, H() As  Double, Wr() As  Double, Wi() As  Double, Vl() As  Double, Vr() As  Double, Mm As  Long, M As  Long, Ifaill() As  Long, Ifailr() As  Long, Info As  Long  ) ヘッセンベルグ行列の固有ベクトル (逆反復法) 目的 本ルーチンは逆反復法を使用して上ヘッセンベルグ行列 H の右および／または左固有ベクトルを求める. 行列 H の固有値 w に対応する右固有ベクトル x および左固有ベクトル y は次のように定義される. H * x = w * x, y^H * H = w * y^H ただし, y^H はベクトル y の共役転置を表す. 引数 [in] Side = "R": 右固有ベクトルのみ求める. = "L": 左固有ベクトルのみ求める. = "B": 右および左固有ベクトルの両方を求める. [in] Eigsrc (Wr(), Wi()) に入力された固有値の求め方を指定する. = "Q": 固有値が Dhseqr を使用して求められたとすると, H の副対角要素に 0 のものがある, すなわち, ブロック三角行列であれば, j 番目の固有値は j 番目の行／列を含むブロックの固有値であるとみなすことができる. これにより一つの対角ブロックにだけ逆反復法を適用すればよい. = "N": 固有値と対角ブロックの間に対応関係を仮定しない. この場合, 常に行列 H 全体を用いて逆反復法を行わなければならない. [in] Initv = "N": 初期ベクトルは入力されない. = "U": 配列 Vl() および／または Vr() にユーザーにより初期ベクトルが入力される. [in,out] Selct() 配列 Selct(LSelct - 1) (LSelct >= N) [in] 求める固有ベクトルを指定する. 実数固有値 W(j) に対応する固有ベクトルを選択するためには Selct(j) を True に設定しなければならない. 共役複素対の固有値 (Wr(j), Wi(j)) および (Wr(j+1), Wi(j+1)) に対応する固有ベクトルを選択するためには Selct(j) または Selct(j+1) のどちらかまたは両方を True に設定しなければならない. [out] 共役複素対の固有値に対応する固有ベクトルを選択した場合, Selct(j) は True に, Selct(j+1) は False に設定される. [in] N 行列 H の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] H() 配列 H(LH1 - 1, LH2 - 1) (LH1 >= N, LH2 >= N) 上ヘッセンベルグ行列 H. H() の中で NaN が検出された場合, Info = -7 を返す. [in,out] Wr() 配列 Wr(LWr - 1) (LWr >= N) [in] Wi() 配列 Wi(LWi - 1) (LWi >= N) [in] H の固有値の実数部および虚数部. 複素共役対をなす固有値は, Wr() と Wi() の連続した要素に格納されなければならない. [out] 独立な固有ベクトルを探すときに近接固有値には小さな摂動を加えるため, Wr() が変更されることがある. [in,out] Vl() 配列 Vl(LVl1 - 1, LVl2 - 1) (LVl1 >= N, LVl2 >= MM) [in] Initv = "U" かつ Side = "L" または "B" の場合, Vl() に左固有ベクトルの逆反復のための初期ベクトルを入れておかなければならない. 各固有ベクトルのための初期ベクトルは, それぞれの固有ベクトルが格納されるのと同じ列に入れておかなければならない. [out] Side = "L" または "B" の場合, Selct() により指定された左固有ベクトルが Vl() の連続する列に対応する固有値と同じ順に格納される. 複素固有値に対応する複素固有ベクトルは連続する2列に格納される(1列目には実数部, 2列目には虚数部が入る). Side = "R" の場合, Vl() は参照されない. [in,out] Vr() 配列 Vr(LVr1 - 1, LVr2 - 1) (LVr1 >= N, LVr2 >= MM) [in] Initv = "U" かつ Side = "R" または "B" の場合, Vr() に右固有ベクトルの逆反復のための初期ベクトルを入れておかなければならない. 各固有ベクトルのための初期ベクトルは, それぞれの固有ベクトルが格納されるのと同じ列に入れておかなければならない. [out] Side = "R" または "B" の場合, Selct() により指定された右固有ベクトルが Vr() の連続する列に対応する固有値と同じ順に格納される. 複素固有値に対応する複素固有ベクトルは連続する2列に格納される(1列目には実数部, 2列目には虚数部が入る). Side = "L" の場合, Vr() は参照されない. [in] MM 配列 Vl() および／または Vr() の中の列数. (MM >= M) [out] M 固有ベクトルを格納するのに必要な配列 Vl() および／または Vr() の中の列数. 選択された実固有ベクトルそれぞれには 1 列, また, 選択された複素固有ベクトルそれぞれには 2 列が必要である. [out] Ifaill() 配列 Ifaill(LIfaill - 1) (LIfaill >= MM) Side = "L" または "B" の場合, Vl() の第 i 列の左固有ベクトル(固有値 W(j) に対応)が収束しなかったならば Ifaill(i) = j > 0 に設定される. 固有ベクトルが正常に収束したならば Ifaill(i) = 0 に設定される. Vl() の第 i および i+1 列に複素固有ベクトルが入っていれば Ifaill(i) と Ifaill(i+1) は同じ値に設定される. Side = "R" の場合, Ifaill() は参照されない. [out] Ifailr() 配列 Ifailr(LIfailr - 1) (LIfailr >= MM) Side = "R" または "B" の場合, Vr() の第 i 列の右固有ベクトル(固有値 W(j) に対応)が収束しなかったならば Ifailr(i) = j > 0 に設定される. 固有ベクトルが正常に収束したならば Ifailr(i) = 0 に設定される. Vr() の第 i および i+1 列に複素固有ベクトルが入っていれば Ifailr(i) と Ifailr(i+1) は同じ値に設定される. Side = "L" の場合, Ifailr() は参照されない. [out] Info = 0: 正常終了. = -1: パラメータ Side の誤り. (Side <> "R", "L" および "B") = -2: パラメータ Eigsrc の誤り. (Eigsrc <> "Q" および "N") = -3: パラメータ Initv の誤り. (Initv <> "N" および "U") = -4: パラメータ Selct() の誤り. = -5: パラメータ N の誤り. (N < 0) = -6: パラメータ H() の誤り. = -7: パラメータ Wr() の誤り. = -8: パラメータ Wi() の誤り. = -9: パラメータ Vl() の誤り. = -10: パラメータ Vr() の誤り. = -11: パラメータ MM の誤り. (MM < M) = -13: パラメータ Ifaill() の誤り. = -14: パラメータ Ifailr() の誤り. = i > 0: i 個の固有ベクトルが収束しなかった. 詳細は Ifaill() および Ifailr() を参照せよ. 出典 LAPACK 使用例 行列Aの固有値・固有ベクトルを求める. ただし, ( 0.20 -0.11 -0.93 ) A = ( -0.32 0.81 0.37 ) ( -0.80 -0.92 -0.29 ) とする. Dgehrdでヘッセンベルグ形に変換したのち, Dhseqrで固有値を, DhseinおよびDormhrでその固有ベクトルを求める. Sub Ex_Dgehrd_Dhseqr_Dhsein() Const N = 3 Dim A(N - 1, N - 1) As Double, Tau(N - 2) As Double Dim H(N - 1, N - 1) As Double, Z() As Double Dim Wr(N - 1) As Double, Wi(N - 1) As Double, Selct(N - 1) 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, Ifaill(N - 1) As Long, Ifailr(N - 1) As Long Dim 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 H(I, J) = A(I, J) Next Next Call Dhseqr("E", "N", N, Ilo, Ihi, H(), Wr(), Wi(), Z(), Info) If Info <> 0 Then Debug.Print "Error in Dhseqr: Info =", Info Exit Sub End If For I = 0 To N - 1 Selct(I) = True Next Call Dhsein("B", "Q", "N", Selct(), N, A(), Wr(), Wi(), Vl(), Vr(), N, M, Ifaill(), Ifailr(), Info) If Info <> 0 Then Debug.Print "Error in Dhsein: Info =", Info Exit Sub End If Call Dormhr("L", "N", N, N, Ilo, Ihi, A(), Tau(), Vr(), Info) If Info <> 0 Then Debug.Print "Error in Dormhr: Info =", Info Exit Sub End If Call Dormhr("L", "N", N, N, Ilo, Ihi, A(), Tau(), Vl(), Info) If Info <> 0 Then Debug.Print "Error in Dormhr: 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.780366759827849 0.442486777319186 -0.353989421855349 -0.517119803163016 0.539408759469125 0.533652513982974 -0.870185040161694 -0.337856087926769 0.327926867989764 Eigenvectors (R) = -0.922812932547959 -0.385121132500322 -0.406438554278919 6.58962255629428E-02 0.64221680239353 -0.221790002576845 -1.10339145165208 -3.62765327766942E-02 0.49628885836462 M = 3