Ddisna Dopgtr Dopmtr Dorgtr Dormtr Dpteqr Dsbtrd Dsptrd Dstebz Dstedc Dstein Dstemr Dsteqr Dsterf Dsytrd

Dstebz() Sub Dstebz ( Range As  String, Order As  String, N As  Long, Vl As  Double, Vu As  Double, Il As  Long, Iu As  Long, AbsTol As  Double, D() As  Double, E() As  Double, M As  Long, Nsplit As  Long, W() As  Double, Iblock() As  Long, Isplit() As  Long, Info As  Long  ) 対称三重対角行列の固有値 (二分法) 目的 本ルーチンは実対称三重対角行列 T の固有値を求める. 全固有値, 半開区間 (Vl, Vu] 内の全固有値, または, Il 〜 Iu番目の固有値を求めることができる. オーバーフローを避けるために行列は, 最大要素が絶対値で overflow^(1/2) * underflow^(1/4) よりも大きくないように, スケーリングされていなければならない. また, 最大精度を得るためにはこれより小さすぎないようにしなければならない. 引数 [in] Range = "A": すべての固有値を求める. = "V": 半開区間 (Vl, Vu] のすべての固有値を求める. = "I": Il 番目から Iu 番目までの固有値を求める. [in] Order = "B": 分割されたブロック (Iblock(), Isplit() を参照) ごとに固有値がグループ分けされ, ブロック内で最小から最大の順に並べられる. = "E": 行列全体の固有値が最小から最大の順に並べられる. [in] N 三重対角行列の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Vl Range = "V": 固有値を求める区間の下端. Vl 以下または Vu より大きい固有値は返されない. (Vl < Vu) Range = "A"または"I": 参照されない. [in] Vu Range = "V": 固有値を求める区間の上端. Vl 以下または Vu より大きい固有値は返されない. (Vl < Vu) Range = "A"または"I": 参照されない. [in] Il Range = "I": 求める最小固有値の番号. (1 <= Il <= Iu <= N (N > 0 の場合). Il = 1, Iu = 0 (N = 0 の場合)) Range = "A"または"V": 参照されない. [in] Iu Range = "I": 求める最大固有値の番号. (1 <= Il <= Iu <= N (N > 0 の場合). Il = 1, Iu = 0 (N = 0 の場合)) Range = "A"または"V": 参照されない. [in] Abstol 固有値の絶対誤差限界. 固有値(またはクラスター)は, Abstol 以下の幅の区間に入っていると判断された場合にそこにあると見なされる. Abstol がゼロ以下の場合, ulp*|T| が使用される. ただし, |T| は T の1ノルムを表す. Abstol がゼロではなくアンダーフローしきい値の 2倍(2*Dlamch("S")) に設定された場合に, 固有値が最も正確に計算される. [in] D() 配列 D(LD - 1) (LD >= N) 三重対角行列 T の N 個の対角要素. [in] E() 配列 E(LE - 1) (LE >= N - 1) 三重対角行列 T の N - 1 個の副対角要素. [out] M 実際に求められた固有値の数 (0 <= M <= N) (Info = 2, 3 の説明も参照せよ.) [out] Nsplit 行列 T の対角ブロック数. (1 <= Nsplit <= N) [out] W() 配列 W(LW - 1) (LW >= N) M 個の求められた固有値が W() の先頭から入る. (本ルーチンは残りの N - M 要素を作業領域として使うことがある.) [out] Iblock() 配列 Iblock(LIblock - 1) (LIblock >= N) E(j) がゼロまたは小さいような行／列 j それぞれにおいて, 行列 T はブロック対角行列に分割されると考えられる. 終了時 Info = 0 であれば Iblock(i) は固有値 W(i) がどのブロック(1 〜 ブロック数)に属するかを示す. (本ルーチンは残りの N - M 要素を作業領域として使うことがある.) [out] Isplit() 配列 Isplit(LIsplit - 1) (LIsplit >= N) T が小行列に分割する分割点. 最初の小行列は行／列 1 〜 Isplit(0), ２番目は行／列 Isplit(0) + 1 〜 Isplit(1), ... よりなる. Nsplit番目は行／列 Isplit(Nsplit - 2) + 1 〜 Isplit(Nsplit - 1) = N よりなる. (最初の Nsplit 要素のみが実際に使われるが, Nsplit の値を前もって知ることはできないので Isplit() には N 要素分を用意しておかなければならない.) [out] Info = 0: 正常終了. = -1: パラメータ Range の誤り. (Range <> "A", "V" および "I") = -2: パラメータ Order の誤り. (Order <> "B" および "E") = -3: パラメータ N の誤り. (N < 0) = -5: パラメータ Vu の誤り. (Vu < Vl) = -6: パラメータ Il の誤り. (Il < 1 または Il > N) = -7: パラメータ Iu の誤り. (Iu < min(N, Il) または Iu > N) = -9: パラメータ D() の誤り. = -10: パラメータ E() の誤り. = -13: パラメータ W() の誤り. = -14: パラメータ Iblock() の誤り. = -15: パラメータ Isplit() の誤り. = > 0: 一部または全部の固有値が収束しなかった, あるいは, 計算されなかった.   = 1 or 3: 二分法がいくつかの固有値について収束しなかった. それらの固有値は負のブロック番号で示される.   = 2 or 3: 固有値 Il 〜 Iu の全部を求めることができなかった. M < Iu + 1 - Il となる. (Range = "I")   = 4: ゲルシュゴリン区間の初期値が小さすぎる. 固有値は計算されなかった. (Range = "I") 出典 LAPACK 使用例 対称行列 A の固有値・固有ベクトルを求める. ただし, ( 2.20 -0.11 -0.32 ) A = ( -0.11 2.93 0.81 ) ( -0.32 0.81 2.37 ) とする. Dsytrdで対称三重対角形に変換したのち, Dstebzで固有値を, DsteinおよびDormtrでその固有ベクトルを求める. Sub Ex_Dsytrd_Dstebz_Dstein() Const N = 3 Dim A(N - 1, N - 1) As Double, W(N - 1) As Double, Z(N - 1, N - 1) As Double Dim D(N - 1) As Double, E(N - 2) As Double, Tau(N - 2) As Double 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) = 2.2 A(1, 0) = -0.11: A(1, 1) = 2.93 A(2, 0) = -0.32: A(2, 1) = 0.81: A(2, 2) = 2.37 Call Dsytrd("L", N, A(), D(), E(), Tau(), Info) If Info <> 0 Then Debug.Print "Error in Dsytrd: 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 Dstein(N, D(), E(), M, W(), Iblock(), Isplit(), Z(), Ifail(), Info) If Info <> 0 Then Debug.Print "Error in Dstein: Info =", Info Exit Sub End If Call Dormtr("L", "L", "N", N, N, A(), Tau(), Z(), Info) If Info <> 0 Then Debug.Print "Error in Dormtr: Info =", Info Exit Sub End If Debug.Print "Eigenvalues =", W(0), W(1), W(2) Debug.Print "Eigenvectors =" Debug.Print Z(0, 0), Z(0, 1), Z(0, 2) Debug.Print Z(1, 0), Z(1, 1), Z(1, 2) Debug.Print Z(2, 0), Z(2, 1), Z(2, 2) Debug.Print "M =", M, "Nsplit =", Nsplit End Sub 実行結果 Eigenvalues = 1.70705954911046 2.22943643244226 3.56350401844728 Eigenvectors = 0.399322077213382 0.894521385341403 0.200931256445799 -0.481026444340548 0.390994588677271 -0.78468897753835 0.780484105216166 -0.216690384632057 -0.586421212707156 M = 3 Nsplit = 1