Zgsev Zgsgv Znaupd Zneupd

Zgsgv() Sub Zgsgv ( Jobz As  String, Uplo2 As  String, N As  Long, Val() As  Complex, Ptr() As  Long, Ind() As  Long, Val2() As  Complex, Ptr2() As  Long, Ind2() 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 Base2 As  Long = -1, Optional Format As  Long = 0, Optional Mode As  Long = 2, Optional Sigmar As  Double = 0, Optional Sigmai As  Double = 0, Optional Maxiter As  Long = 300, Optional Tol As  Double = 0  ) 一般化固有値問題 (複素疎行列) (リスタート付きアーノルディ法 (IRAM)) (Arpack) (ドライバ) 目的 一般化固有値問題 A*z = λ*B*z において, λは固有値, z はその固有ベクトルを表す. A は複素疎行列, B は半正定値エルミート疎行列である. 本ルーチンは, Arpack ルーチン Znaupd および Zneupd を使用して, リスタート付きアーノルディ法 (Implicitly restarted Arnoldi method (IRAM)) により上の問題の固有値 λ および固有ベクトル z を求める. 疎行列 A は CSC 形式または CSR 形式で格納されているものとする. また, エルミート疎行列 B は上三角部分のみ, 下三角部分のみあるいは両方(行列全体)が CSC 形式または CSR 形式で格納されているものとする. 本ルーチンは以下の2つのモードで動作する. (1) 通常インバースモード Mode = 2 であれば, 通常インバースモードで計算を行う. 大きい方から, あるいは, 小さい方から数個の固有値およびその固有ベクトルを求めることができる. Sigmar および Sigmai は参照されない. (2) シフト＆インバートモード Mode = 3 であれば, シフト＆インバートモードで計算を行う. 指定した σ (シフト) の近くの固有値を求めることができる. OP = (A - σ*M)^(-1)*M, M = B として, 別の固有値問題 OP*x = ν*M*x の固有値 ν を求めたとすると, 元の問題の固有値 λ は ν = 1/(λ - σ) より求めることができる. すなわち, 最大の ν を求めれば σ に最も近い λ を求められる. パラメータ Which は求める固有値 ν を指定する. 一般的には, "LM"とすれば σ の最も近くの固有値 λ を求めることができる. 引数 [in] Jobz = 'N': 固有値のみ求める. = 'V': 固有値と固有ベクトルを求める. [in] Uplo2 行列 B の上三角部分あるいは下三角部分のどちらが格納されているかを指定. = 'U': 上三角部分. = 'L': 下三角部分. = 'B': 行列全体(上および下三角部分). 指定された以外の三角部分(対角成分をを除く)は無視される. [in] N 行列の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Val() 配列 Val(LVal - 1) (LVal >= Nnz) 行列 A の非ゼロ要素の値. (Nnz は非ゼロ要素数) [in] Ptr() 配列 Ptr(LPtr - 1) (LPtr >= N + 1 行列 A の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind() 配列 Ind(LInd - 1) (LInd >= Nnz) 行列 A の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). (Nnz は非ゼロ要素数) [in] Val2() 配列 Val2(LVal2 - 1) (LVal2 >= Nnz2) 行列 B の非ゼロ要素の値. (Nnz2 は非ゼロ要素数) [in] Ptr2() 配列 Ptr2(LPtr2 - 1) (LPtr2 >= N + 1 行列 B の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind2() 配列 Ind2(LInd2 - 1) (LInd2 >= Nnz2) 行列 B の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). (Nnz2 は非ゼロ要素数) [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 = λ*B*z の固有値を近似するリッツ値を返す. [out] Z() 配列 Z(Ldz - 1, LZ - 1) (Ldz >= N, LZ >= Nev) 固有値問題 A*z = λ*B*z のリッツ値に対応するリッツベクトル(近似固有ベクトル)を返す. Jobz = "N" の場合, Z() は参照されない. [out] Info リターンコード. = 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = 1: 最大反復回数に達した. = 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] Base2 (省略可) Ptr2() および Ind2() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Ptr2(0) = 1 であれば 1, そうでなければ 0 とみなす) [in] Format (省略可) 行列の格納形式. (省略時 = 0) = 0: CSR 形式. = 1: CSC 形式. [in] Mode (省略可) 動作モード. (省略時 = 2) = 2: 通常インバースモード. = 3: シフト＆インバートモード. [in] Sigmar (省略可) シフト(σ)の実数部. Mode = 2 のときには参照されない. (省略時 = 0) [in] Sigmai (省略可) シフト(σ)の虚数部. Mode = 2 のときには参照されない. (省略時 = 0) [in] Maxiter (省略可) アーノルディ反復の最大回数. (MaxIter > 0) (省略時 = 300) [in] Tol (省略可) 停止基準: リッツ値の最大相対誤差. (省略時 = 0) Tol <= 0 の場合, マシン精度とみなす. 使用例 行列 A および エルミート行列 B について最大の一般化固有値とその固有ベクトルを求める. ただし, ( 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 ) ( 2.31 0.25-0.23i -0.81+0.83i ) B = ( 0.25+0.23i 2.26 -0.56+0.08i ) ( -0.81-0.83i -0.56-0.08i 2.09 ) とする. Sub Ex_Zgsgv() Const N = 3, Nnz = 9, Nnz2 = 6, Ncv = 3, Nev = 2 Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(Nnz2 - 1) As Complex, Ib(N) As Long, Jb(Nnz2 - 1) As Long Dim D(Nev - 1) As Complex, Z(N - 1, Nev - 1) As Complex Dim Nconv As Long, Info 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 B(0) = Cmplx(2.31): B(1) = Cmplx(0.25, 0.23): B(2) = Cmplx(2.26): B(3) = Cmplx(-0.81, -0.83): B(4) = Cmplx(-0.56, -0.08): B(5) = Cmplx(2.09) Ib(0) = 0: Ib(1) = 1: Ib(2) = 3: Ib(3) = 6 Jb(0) = 0: Jb(1) = 0: Jb(2) = 1: Jb(3) = 0: Jb(4) = 1: Jb(5) = 2 Call Zgsgv("V", "L", N, A(), Ia(), Ja(), B(), Ib(), Jb(), "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 = 0.790264815830272 0.731095506039657 0.144386602592976 0.763021818472598 Z = -0.602926109593061 -2.54716468513815E-02 -6.04572102722823E-03 -0.162145611949676 0.174190955405643 -0.156732370014386 0.578522100432737 -5.20895534988314E-02 -0.403225790941503 -0.558909798961031 0.533641052631575 -0.179122219344393 Nconv = 2, Info = 0