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

Zhseqr() Sub Zhseqr ( Job As  String, Compz As  String, N As  Long, Ilo As  Long, Ihi As  Long, H() As  Complex, W() As  Complex, Z() As  Complex, Info As  Long  ) ヘッセンベルグ行列の固有値およびシュール分解 (QR法) (複素行列) 目的 本ルーチンは複素ヘッセンベルグ行列 H の固有値, および, 必要によりシュール分解 H = Z T Z^H の行列 T および Z を求める. ここで, T は上準三角行列(シュール形), Z はシュールベクトルからなるユニタリ行列である. 必要により行列 A のシュール分解を与えるために Z を入力されたユニタリ行列 Q に乗算することができる. ただし, A はユニタリ行列 Q により A = Q*H*Q^H = (QZ)*T*(QZ)^H のようにヘッセンベルグ形 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 は通常 Zgebal の呼び出しによって設定され, Zgebal の出力行列がヘッセンベルグ形に変換されるときに Zgehrd に渡される. そうでなければ, これらはそれぞれ 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 が入る. Info = 0 で Job = "E" の場合, H() の内容は不定である. (Info > 0 の場合の H() の出力値については下記 Info の説明を参照せよ.) [out] W() 配列 W(LW - 1) (LW >= N) 求められた固有値. Job = "S" の場合, 固有値は H() に返されるシュール形の対角要素と同じ順に格納され, W(i) = H(i, 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 への変換のための Zgehrd の呼び出し後 Zunghr によって生成されるユニタリ行列である. (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: パラメータ W() の誤り. = -8: パラメータ Z() の誤り. = i > 0: すべての固有値を求めることはできなかった. W() の要素 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.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で固有値を求める. 固有ベクトルも求める例としては Ztrevc3, Zhsein の使用例を参照せよ. Sub Ex_Zgehrd_Zhseqr() Const N = 3 Dim A(N - 1, N - 1) As Complex, Tau(N - 2) As Complex Dim W(N - 1) As Complex, Z() As Complex Dim Ilo As Long, Ihi 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 Call Zhseqr("E", "N", N, Ilo, Ihi, A(), W(), Z(), Info) If Info <> 0 Then Debug.Print "Error in Zhseqr: 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)) End Sub 実行結果 Eigenvalues = -1.15894122423918 -0.50662892448174 1.05593587167591 0.900255855387815 0.21300535256327 1.29637306909393