Zgebak Zgebal Zgehrd Zhsein Zhseqr Ztrevc3 Ztrexc Ztrsen Ztrsna Ztrsyl Zunghr Zunmhr

Zhsein() Sub Zhsein ( Side As  String, Eigsrc As  String, Initv As  String, Selct() As  Boolean, N As  Long, H() As  Complex, W() As  Complex, Vl() As  Complex, Vr() As  Complex, 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 W() に入力された固有値の求め方を指定する. = "Q": Zhseqr を使用して固有値を求めた. H の副対角要素に 0 のものがある, すなわち, ブロック三角行列であれば, j 番目の固有値は j 番目の行／列を含むブロックの固有値であるとみなすことができる. これにより一つの対角ブロックにだけ逆反復法を適用すればよくなる. = "N": 固有値と対角ブロックの間に対応関係を仮定しない. この場合, 常に行列 H 全体を用いて逆反復法を行わなければならない. [in] Initv = "N": 初期ベクトルは入力されない. = "U": 配列 Vl() および／または Vr() にユーザーにより初期ベクトルが入力される. [in] Selct() 配列 Selct(LSelct - 1) (LSelct >= N) 求める固有ベクトルを指定する. 固有値 W(j) に対応する固有ベクトルを選択するためには Selct(j) を True に設定しなければならない. [in] N 行列 H の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] H() 配列 H(LH1 - 1, LH2 - 1) (LH1 >= N, LH2 >= N) 上ヘッセンベルグ行列 H. H() の中で NaN が検出された場合, Info = -6 を返す. [in] W() 配列 W(LW - 1) (LW >= N) [in] H の固有値. [out] 独立な固有ベクトルを探すときに近接固有値には小さな摂動を加えるため, W() の実数部が変更されることがある. [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() の連続する列に対応する固有値と同じ順に格納される. 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() の連続する列に対応する固有値と同じ順に格納される. Side = "L" の場合, Vr() は参照されない. [in] MM 配列 Vl() および／または Vr() の中の列数. (MM >= M) [out] M 固有ベクトルを格納するのに必要な配列 Vl() および／または Vr() の中の列数 (= Selct() 中の True の数). [out] Ifaill() 配列 Ifaill(LIfaill - 1) (LIfaill >= MM) Side = "L" または "B" の場合, Vl() の第 i 列の左固有ベクトル(固有値 W(j) に対応)が収束しなかったならば Ifaill(i) = j > 0 に設定される. 固有ベクトルが正常に収束したならば Ifaill(i) = 0 に設定される. Side = "R" の場合, Ifaill() は参照されない. [out] Ifailr() 配列 Ifailr(LIfailr - 1) (LIfailr >= MM) Side = "R" または "B" の場合, Vr() の第 i 列の右固有ベクトル(固有値 W(j) に対応)が収束しなかったならば Ifailr(i) = j > 0 に設定される. 固有ベクトルが正常に収束したならば Ifailr(i) = 0 に設定される. 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: パラメータ W() の誤り. = -8: パラメータ Vl() の誤り. = -9: パラメータ Vr() の誤り. = -10: パラメータ MM の誤り. (MM < M) = -12: パラメータ Ifaill() の誤り. = -13: パラメータ Ifailr() の誤り. = i > 0: i 個の固有ベクトルが収束しなかった. 詳細は Ifaill() および Ifailr() を参照せよ. 出典 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で固有値を, ZhseinおよびZunmhrでその固有ベクトルを求める. Sub Ex_Zgehrd_Zhseqr_Zhsein() Const N = 3 Dim A(N - 1, N - 1) As Complex, Tau(N - 2) As Complex, H(N - 1, N - 1) As Complex Dim W(N - 1) As Complex, Z() As Complex, Selct(N - 1) 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, Ifaill(N - 1) As Long, Ifailr(N - 1) As Long Dim 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 H(I, J) = A(I, J) Next Next Call Zhseqr("E", "N", N, Ilo, Ihi, H(), W(), Z(), Info) If Info <> 0 Then Debug.Print "Error in Zhseqr: Info =", Info Exit Sub End If For I = 0 To N - 1 Selct(I) = True Next Call Zhsein("B", "Q", "N", Selct(), N, A(), W(), Vl(), Vr(), N, M, Ifaill(), Ifailr(), Info) If Info <> 0 Then Debug.Print "Error in Zhsein: Info =", Info Exit Sub End If Call Zunmhr("L", "N", N, N, Ilo, Ihi, A(), Tau(), Vr(), Info) If Info <> 0 Then Debug.Print "Error in Zunmhr: Info =", Info Exit Sub End If Call Zunmhr("L", "N", N, N, Ilo, Ihi, A(), Tau(), Vl(), Info) If Info <> 0 Then Debug.Print "Error in Zunmhr: 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 -0.871083495245762 0.2 -0.8 0.267351466048342 -0.480647278057965 -0.03692440974696 0.616773162002434 -0.285413896331544 0.542078734362481 0.166244018883372 -0.368066071419648 -2.22091099794888E-02 0.344241204682076 7.27399270176405E-02 -0.5924933721717 0.347217078108291 -0.4230226508624 Eigenvectors (R) = 0.690958738203183 -0.309041261796817 -0.507806593830091 0.492193406169909 0.634809967097297 -0.421369311149223 0.229049609768182 -0.190619676572438 -0.363524659234702 0.113190794408639 -0.763294421913407 0.107049283139747 -0.111157316336775 -2.52863433046769E-02 -0.481933704163733 0.142931109108698 -0.552705878022003 3.27787741584283E-02 M = 3