Zhbtrd Zhetrd Zhptrd Zpteqr Zstedc Zstein Zstemr Zsteqr Zungtr Zunmtr Zupgtr Zupmtr

Zstein() Sub Zstein ( N As  Long, D() As  Double, E() As  Double, M As  Long, W() As  Double, Iblock() As  Long, Isplit() As  Long, Z() As  Complex, Ifail() As  Long, Info As  Long  ) エルミート行列を変換した対称三重対角行列の固有ベクトル (逆反復法) 目的 本ルーチンはエルミート行列を変換した実対称三重対角行列 T の指定された固有値に対応する固有ベクトルを逆反復法により求める. 各固有ベクトルの計算に許される最大反復回数は内部パラメータ Maxits (現在 5 に設定されている) で指定される. 引数 [in] N 三重対角行列の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] D() 配列 D(LD - 1) (LD >= N) 三重対角行列 T の N 個の対角要素. [in] E() 配列 E(LE - 1) (LE >= N - 1) 三重対角行列 T の N - 1 個の副対角要素. [in] M 求める固有ベクトルの数 (0 <= M <= N) [in] W() 配列 W(LW - 1) (LW >= N) W() の先頭 M 個にその固有ベクトルを求める固有値を入れる. 固有値は分割されたブロックごとにグループ分けされ, ブロック内で最小から最大の順に並べられていること. (Order = "B" とした Dstebz の出力である配列 W() が想定されている.) [in] Iblock() 配列 Iblock(LIblock - 1) (LIblock >= N) W() 内の固有値に対応する小行列のインデックス. 固有値 W(i) が先頭の小行列に属するとき Iblock(i) = 1, W(i) が2つ目の小行列に属するときは 2, ... である. (Dstebz の出力である配列 Iblock() が想定されている.) [in] Isplit() 配列 Isplit(LIsplit - 1) (LIsplit >= N) T が小行列に分割する分割点. 最初の小行列は行／列 1 〜 Isplit(0), 2番目は行／列 Isplit(0) + 1 〜 Isplit(1), ... よりなる. (Dstebz の出力である配列 Isplit() が想定されている.) [out] Z() 配列 Z(LZ1 - 1, LZ2 - 1) (LZ1 >= N, LZ2 >= M) 求められた固有ベクトル. 固有値 W(i) に関連する固有ベクトルが Z() の i 列に入る. 収束に失敗したときは Maxits 回の反復後時点のベクトルが設定される. [out] Ifail() 配列 Ifail(LIfail - 1) (LIfail >= M) 正常終了時, Ifail() のすべての要素はゼロである. 1つあるいは複数の固有ベクトルが Maxits 回の反復で収束しなかったならばそれらのインデックスが配列 Ifail() に格納される. [out] Info = 0: 正常終了. = -1: パラメータ N の誤り. (N < 0) = -2: パラメータ D() の誤り. = -3: パラメータ E() の誤り. = -4: パラメータ M の誤り. (M < 0 または M > N) = -5: パラメータ W() の誤り. = -6: パラメータ Iblock() の誤り. = -7: パラメータ Isplit() の誤り. = -8: パラメータ Z() の誤り. = -9: パラメータ Ifail() の誤り. = i > 0: i 個の固有ベクトルが Maxits 回の反復で収束しなかった. そのインデックスが Ifail() に格納される. 出典 LAPACK 使用例 エルミート行列Aの固有値・固有ベクトルを求める. ただし, ( 0.20 -0.11+0.93i 0.81-0.37i ) A = ( -0.11-0.93i -0.32 -0.80+0.92i ) ( 0.81+0.37i -0.80-0.92i -0.29 ) とする. Zhetrdで実対称三重対角形に変換したのち, Dstebzで固有値を, ZsteinおよびZunmtrでその固有ベクトルを求める. Sub Ex_Zhetrd_Dstebz_Zstein() Const N = 3 Dim A(N - 1, N - 1) As Complex, W(N - 1) As Double, Z(N - 1, N - 1) As Complex Dim D(N - 1) As Double, E(N - 2) As Double, Tau(N - 2) As Complex Dim Vl As Double, Vu As Double, Il As Long, Iu As Long, Abstol As Double Dim Iblock(N - 1) As Long, Isplit(N - 1) As Long, Ifail(N - 1) As Long Dim M As Long, Nsplit As Long, Info As Long A(0, 0) = Cmplx(0.2, 0) A(1, 0) = Cmplx(-0.11, -0.93): A(1, 1) = Cmplx(-0.32, 0) A(2, 0) = Cmplx(0.81, 0.37): A(2, 1) = Cmplx(-0.8, -0.92): A(2, 2) = Cmplx(-0.29, 0) Call Zhetrd("L", N, A(), D(), E(), Tau(), Info) If Info <> 0 Then Debug.Print "Error in Zhetrd: Info =", Info Exit Sub End If Abstol = 0 Call Dstebz("A", "E", N, Vl, Vu, Il, Iu, Abstol, D(), E(), M, Nsplit, W(), Iblock(), Isplit(), Info) If Info <> 0 Then Debug.Print "Error in Dstebz: Info =", Info Exit Sub End If Call Zstein(N, D(), E(), M, W(), Iblock(), Isplit(), Z(), Ifail(), Info) If Info <> 0 Then Debug.Print "Error in Zstein: Info =", Info Exit Sub End If Call Zunmtr("L", "L", "N", N, N, A(), Tau(), Z(), Info) If Info <> 0 Then Debug.Print "Error in Zunmtr: Info =", Info Exit Sub End If Debug.Print "Eigenvalues =", W(0), W(1), W(2) Debug.Print "Eigenvectors =" Debug.Print Creal(Z(0, 0)), Cimag(Z(0, 0)), Creal(Z(0, 1)), Cimag(Z(0, 1)) Debug.Print Creal(Z(1, 0)), Cimag(Z(1, 0)), Creal(Z(1, 1)), Cimag(Z(1, 1)) Debug.Print Creal(Z(2, 0)), Cimag(Z(2, 0)), Creal(Z(2, 1)), Cimag(Z(2, 1)) Debug.Print Creal(Z(0, 2)), Cimag(Z(0, 2)) Debug.Print Creal(Z(1, 2)), Cimag(Z(1, 2)) Debug.Print Creal(Z(2, 2)), Cimag(Z(2, 2)) Debug.Print "M =", M, "Nsplit =", Nsplit End Sub 実行結果 Eigenvalues = -2.05348849668514 0.124622388617308 1.51886610806783 Eigenvectors = -0.449276526719113 0 -0.654793596518192 0 0.227247885813611 -0.597641779578735 -0.519997178670921 3.19846835072554E-02 0.621236109316912 -5.83009495222983E-02 -0.204907317474214 0.507777757881846 0.607779522934083 0 -0.392237107311198 -0.407323787101333 0.23846608290599 0.503959683819116 M = 3 Nsplit = 1