Zgsev Zgsgv Znaupd Zneupd

Zgsev() Sub Zgsev ( Jobz As  String, N As  Long, Val() As  Complex, Ptr() As  Long, Ind() As  Long, Which As  String, Nev As  Long, Ncv As  Long, D() As  Complex, Z() As  Complex, Optional Info As  Long, Optional Nconv As  Long, Optional Niter As  Long, Optional Base As  Long = -1, Optional Format As  Long = 0, Optional Mode As  Long = 1, Optional Sigmar As  Double = 0, Optional Sigmai As  Double = 0, Optional Maxiter As  Long = 300, Optional Tol As  Double = 0  ) 固有値・固有ベクトル (複素疎行列) (リスタート付きアーノルディ法 (IRAM)) (Arpack) (ドライバ) 目的 固有値問題 A*z = λ*z において, λは固有値, z はその固有ベクトルを表す. A は複素疎行列である. 本ルーチンは, Arpack ルーチン Znaupd および Zneupd を使用して, リスタート付きアーノルディ法 (Implicitly restarted Arnoldi method (IRAM)) により上の問題の固有値 λ および固有ベクトル z を求める. 疎行列 A は CSC 形式または CSR 形式で格納されているものとする. 本ルーチンは以下の2つのモードで動作する. (1) 通常モード Mode = 1 であれば通常モードで計算を行う. 大きい方から, あるいは, 小さい方から数個の固有値およびその固有ベクトルを求めることができる. Sigmar および Sigmai は参照されない. (2) シフト＆インバートモード Mode = 3 であればシフト＆インバートモードで計算を行う. 指定した σ (シフト) の近くの固有値を求めることができる. OP = (A - σ*I)^(-1) として, 別の固有値問題 OP*x = ν*x の固有値 ν を求めたとすると, 元の問題の固有値 λ は ν = 1/(λ - σ) より求めることができる. すなわち, 最大の ν を求めれば σ に最も近い λ を求められる. パラメータ Which は求める固有値 ν を指定する. 一般的には, "LM"とすれば σ の最も近くの固有値 λ を求めることができる. 引数 [in] Jobz = 'N': 固有値のみ求める. = 'V': 固有値と固有ベクトルを求める. [in] N 行列の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Val() 配列 Val(LVal - 1) (LVal >= Nnz) 入力行列の非ゼロ要素の値. (Nnz は非ゼロ要素数) [in] Ptr() 配列 Ptr(LPtr - 1) (LPtr >= N + 1 入力行列の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind() 配列 Ind(LInd - 1) (LInd >= Nnz) 入力行列の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). (Nnz は非ゼロ要素数) [in] Which = "LM": 絶対値最大のものから Nev 個の固有値を求める. = "SM": 絶対値最小のものから Nev 個の固有値を求める. = "LR": 実数部が最大のものから Nev 個の固有値を求める. = "SR": 実数部が最小のものから Nev 個の固有値を求める. = "LI": 虚数部が最大のものから Nev 個の固有値を求める. = "SI": 虚数部が最小のものから Nev 個の固有値を求める. [in] Nev 求める固有値の数. (0 < Nev < N) [in] Ncv 行列 V の列数. (Nev < Ncv <= N) 各反復において生成されるアーノルディ基底ベクトルの数 (Ncv >= 2*Nev が推奨される). [out] D() 配列 D(LD - 1) (LD >= Nev) A*z = λ*z の固有値を近似するリッツ値を返す. [out] Z() 配列 Z(Ldz - 1, LZ - 1) (Ldz >= N, LZ >= Nev) 固有値問題 A*z = λ*z のリッツ値に対応するリッツベクトル(近似固有ベクトル)を返す. Jobz = "N" の場合, Z() は参照されない. [out] Info リターンコード. = 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = 1: 最大反復回数に達した. iparam[4] に収束したリッツ値の数を返す. = 3: リスタート付きアーノルディ反復中にシフト値が適用できなかった. Ncv を Nev に比較して大きくすると改善するかもしれない (Ncv >= 2*Nev が推奨される). = 11: 初期残差ベクトルが 0 である. = 12: アーノルディ分解に失敗した. = 13: LAPACK ルーチンで固有値の計算に失敗した. = 14: Znaupd は十分な精度の固有値を得られなかった. = 15 〜 19: 内部処理エラー. [out] Nconv (省略可) 収束条件を満たしたリッツ値の数. [out] Niter (省略可) アーノルディ反復回数. [in] Base (省略可) Ptr() および Ind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Ptr(0) = 1 であれば 1, そうでなければ 0 とみなす) [in] Format (省略可) 行列の格納形式. (省略時 = 0) = 0: CSR 形式. = 1: CSC 形式. [in] Mode (省略可) 動作モード. (省略時 = 1) = 1: 通常モード. = 3: シフト＆インバートモード. [in] Sigmar (省略可) シフト(σ)の実数部. Mode = 1 のときには参照されない. (省略時 = 0) [in] Sigmai (省略可) シフト(σ)の虚数部. Mode = 1 のときには参照されない. (省略時 = 0) [in] Maxiter (省略可) ランチョス反復の最大回数. (MaxIter > 0) (省略時 = 300) [in] Tol (省略可) 停止基準: リッツ値の最大相対誤差. (省略時 = 0) Tol <= 0 の場合, マシン精度とみなす. 使用例 行列 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 ) とする. Sub Ex_Zgsev() Const N = 3, Nnz = 9, Ncv = 3, Nev = 2 Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long Dim D(Nev - 1) As Complex, Z(N - 1, Nev - 1) As Complex Dim Nconv As Long, Info As Long, I As Long A(0) = Cmplx(0.2, -0.11): A(1) = Cmplx(-0.93, -0.32): A(2) = Cmplx(0.81, 0.37): A(3) = Cmplx(-0.8, -0.92): A(4) = Cmplx(-0.29, 0.86): A(5) = Cmplx(0.64, 0.51): A(6) = Cmplx(0.71, 0.59): A(7) = Cmplx(-0.15, 0.19): A(8) = Cmplx(0.2, 0.94) Ia(0) = 0: Ia(1) = 3: Ia(2) = 6: Ia(3) = 9 Ja(0) = 0: Ja(1) = 1: Ja(2) = 2: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2: Ja(6) = 0: Ja(7) = 1: Ja(8) = 2 Call Zgsev("V", N, A(), Ia(), Ja(), "LM", Nev, Ncv, D(), Z(), Info, Nconv) Debug.Print "D =" Debug.Print Creal(D(0)), Cimag(D(0)), Creal(D(1)), Cimag(D(1)) Debug.Print "Z =" 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 "Nconv =" + Str(Nconv) + ", Info =" + Str(Info) End Sub 実行結果 D = 1.05593587167591 0.900255855387814 0.213005352563271 1.29637306909393 Z = -0.626960983365491 0.172271418863805 7.9679684322005E-03 0.150485895565997 0.268918491094688 -5.24016295711404E-02 0.3557286693982 0.561279089020254 -0.662645215037481 -0.251158748471167 0.241195085578284 0.691041486022676 Nconv = 2, Info = 0