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

Dhseqr() Sub Dhseqr ( Job As  String, Compz As  String, N As  Long, Ilo As  Long, Ihi As  Long, H() As  Double, Wr() As  Double, Wi() As  Double, Z() As  Double, Info As  Long  ) ヘッセンベルグ行列の固有値およびシュール分解 (QR法) 目的 本ルーチンはヘッセンベルグ行列 H の固有値, および, 必要によりシュール分解 H = Z T Z^T の行列 T および Z を求める. ここで, T は上準三角行列(シュール形), Z はシュールベクトルからなる直交行列である. 必要により行列 A のシュール分解を与えるために Z を入力された直交行列 Q に乗算することができる. ただし, A は直交行列 Q により A = Q*H*Q^T = (QZ)*T*(QZ)^T のようにヘッセンベルグ形 H に変換された行列である. 引数 [in] Job = "E": 固有値のみ求める. = "S": 固有値およびシュール形 T を求める. [in] Compz = "N": シュールベクトルを求めない. = "I": Z() は単位行列に初期化される. H のシュールベクトルからなる行列 Z が返される. = "V": Z() に直交行列 Q を入力しておく. 積 Q*Z が返される. [in] N 行列 A の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Ilo [in] Ihi 行および列 1〜Ilo-1 および Ihi+1〜N において H はすでに上三角形であるとみなす. Ilo および Ihi は通常 Dgebal の呼び出しによって設定され, Dgebal の出力行列がヘッセンベルグ形に変換されるときに Dgehrd に渡される. そうでなければ, これらはそれぞれ 1 および N に設定しなければならない. (1 <= Ilo <= Ihi <= N (N > 0 の場合). Ilo = 1 かつ Ihi = 0 (N = 0 の場合)) [in,out] H() 配列 H(LH1 - 1, LH2 - 1) (LH1 >= N, LH2 >= N) [in] 上ヘッセンベルグ行列 H. [out] Info = 0 で Job = "S" の場合, シュール分解(シュール形)の上準三角行列 T が入る. 2 x 2 対角ブロック(固有値の共役複素対に対応)は標準的な形式 (H(i, i) = H(i+1, i+1), H(i+1, i])*H(i, i+1) < 0) で返される. Info = 0 で Job = "E" の場合, H() の内容は不定である. (Info > 0 の場合の H() の出力値については下記 Info の説明を参照せよ.) [out] Wr() 配列 Wr(LWr - 1) (LWr >= N) [out] Wi() 配列 Wi(LWi - 1) (LWi >= N) 求められた固有値の実数部および虚数部. 2つの固有値が共役複素対の場合, それらは Wr() と Wi() の連続する要素 i 番目および i+1 番目に Wi(i) > 0 および Wi(i+1) < 0 となるように格納される. Job = "S" の場合, 固有値は H() に返されるシュール形の対角要素と同じ順に格納される. すなわち, Wr(i) = H(i, i). また, H(i〜i+1, i〜i+1) が 2 x 2 対角ブロックであれば, Wi(i) = sqrt(-H(i+1, i)*H(i, i+1)), Wi(i+1) = -Wi(i) となる. [in,out] Z() 配列 Z(LZ1 - 1, LZ2 - 1) (LZ1 >= N, LZ2 >= N) Compz = "I": [in] Z() は設定する必要がない. [out] Info = 0 の場合, Z() に H のシュールベクトルからなる直交行列 Z を返す. Compz = "V": [in] Z() に N x N 行列 Q を設定する. ただし, 小行列 Z(Ilo〜Ihi, Ilo〜Ihi) 以外は単位行列とみなす. [out]Info = 0 の場合, Z() に Q*Z を返す. 通常 Q は, ヘッセンベルグ行列 H への変換のための Dgehrd の呼び出し後 Dorghr によって生成される直交行列である. (Info > 0 の場合の Z() の出力値については下記 Info の説明を参照せよ.) Compz = "N": Z() は参照されない. [out] Info = 0: 正常終了. = -1: パラメータ Job の誤り. (Job <> "E" および "S") = -2: パラメータ Compz の誤り. (Compz <> "N", "I" および "V") = -3: パラメータ N の誤り. (N < 0) = -4: パラメータ Ilo の誤り. (Ilo < 1 または Ilo > N) = -5: パラメータ Ihi の誤り. (Ihi < min(Ilo, N) または Ihi > N) = -6: パラメータ H() の誤り. = -7: パラメータ Wr() の誤り. = -8: パラメータ Wi() の誤り. = -9: パラメータ Z() の誤り. = i > 0: すべての固有値を求めることはできなかった. Wr() と Wi() の要素 1〜Ilo-1 および i+1〜N には正常に計算済の固有値が格納される. (失敗することはまれである.) Job = "E" の場合, 収束しなかった残りの固有値は H() の最終出力の上ヘッセンベルグ行列の行および列 Ilo〜i の固有値である. Job = "S" の場合, 次のようになる. (H() の初期値) * U = U * (H() の最終値) --- (*) ただし, U はの直交行列である. H() の最終値は上ヘッセンベルグかつ行および列i+1〜Ihiにおいて準三角行列である. Compz = "V" の場合, 次のようになる. (Z() の最終値) = (Z() の初期値) * U ただし, U は (*) の直交行列である (Job の値に関係ない). Compz = "I" の場合, 次のようになる. (Z() の最終値) = U ただし, U は (*) の直交行列である (Job の値に関係ない). Compz = "N" の場合, Z() は参照されない. 出典 LAPACK 使用例 行列 A の固有値を求める. ただし, ( 0.20 -0.11 -0.93 ) A = ( -0.32 0.81 0.37 ) ( -0.80 -0.92 -0.29 ) とする. 固有ベクトルも求める例としては Dtrevc3, Dhsein の使用例を参照せよ. Sub Ex_Dgehrd_Dhseqr() 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, Z() As Double Dim Ilo As Long, Ihi 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 Call Dhseqr("E", "N", N, Ilo, Ihi, A(), Wr(), Wi(), Z(), Info) If Info <> 0 Then Debug.Print "Error in Dhseqr: 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) End Sub 実行結果 Eigenvalues (r) = -0.904130023851345 0.812065011925673 0.812065011925673 Eigenvalues (i) = 0 0.48915757543818 -0.48915757543818