Dgsev Dgsgv Dnaupd Dneupd Dsaupd Dseupd Dssev Dssgv

Dnaupd() Sub Dnaupd ( Ido As  Long, Bmat As  String, N As  Long, Which As  String, Nev As  Long, Resid() As  Double, Ncv As  Long, V() As  Double, IParam() As  Long, Ipntr() As  Long, Workd() As  Double, Workl() As  Double, Info As  Long, Optional Tol As  Double = 0  ) アーノルディ分解 (一般疎行列) (Arpack) 目的 リスタート付きアーノルディ法 (Implicitly restarted Arnoldi method (IRAM)) により一般疎行列の固有値・固有ベクトルを求める. 本ルーチンはアルゴリズムの前半部であるアーノルディ分解を行う. 分解結果から求められる近似固有値はリッツ値とよばれ, 対応する近似固有ベクトルはリッツベクトルとよばれる. アルゴリズムの後半部であるリッツ値およびリッツベクトルの計算は, 本ルーチンの出力をもとに Dneupd が行う. 本ルーチンは, 半正定値実対称行列 B によって定義される半内積に関する線形演算子「OP」のいくつかの近似固有値対を求める. 注 - 線形演算子「OP」が実数で, 半正定値実対称行列 B に関して対称である場合, つまり B*OP = (OP^T)*B の場合には, 代わりに Dsaupd を使用すべきである. 本ルーチンは, 以下の問題のどれかを解くために, 収束するまで繰り返し呼び出される. 通常モード (mode = 1): A*x = λ*x. (OP = A, B = I) インバースモード (mode = 2): A*x = λ*M*x, M は正定値対称. (OP = M^(-1)*A, B = M) シフト＆インバートモード (mode = 3): A*x = λ*M*x, M は半正定値対称. (OP = (A - σ*M)^(-1)*M の実数部, B = M) シフト＆インバートモード (mode = 4): A*x = λ*M*x, M は半正定値対称. (OP = (A - σ*M)^(-1)*M の虚数部, B = M) 注 - モード (mode) は Iparam(6) で指定される. 引数 [in,out] Ido リバースコミュニケーション用フラグ. [in] 初回呼び出し時 Ido = 0 とすること. それ以降の呼び出し時には値を変更してはならない. [out] 要求する操作を示す値. 呼び出し元に戻ったときにその操作を実行し, 再呼び出しして結果を返すこと. = -1: y = OP * x を計算. (初期化フェーズ用) = 1: y = OP * x を計算. = 2: y = B * x を計算. = 3: np シフトを計算. (注 3) = 99: 終了. 注 1 - x の格納先: Workd(Ipntr(0) - 1), Workd(Ipntr(0)), ..., Workd(Ipntr(0) + n - 2). 注 2 - y の格納先: Workd(Ipntr(1) - 1), Workd(Ipntr(1)), ..., Workd(Ipntr(1) + n - 2). 注 3 - Iparam(0) = 0 かつ Ido = 3 の場合, np (= Iparam(7)) シフトを Workl(Ipntr(10) - 1), Workl(Ipntr(10)), ..., Workl(Ipntr(10) + np - 2) に格納すること. [in] Bmat = "I": 標準固有値問題 A*x = λ*x = "G": 一般化固有値問題 A*x = λ*B*x [in] N 行列の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Which = "LM": 絶対値最大のものから Nev 個の固有値を求める. = "SM": 絶対値最小のものから Nev 個の固有値を求める. = "LR": 実数部が最大のものから Nev 個の固有値を求める. = "SR": 実数部が最小のものから Nev 個の固有値を求める. = "LI": 虚数部が最大のものから Nev 個の固有値を求める. = "SI": 虚数部が最小のものから Nev 個の固有値を求める. [in] Nev 求める固有値の数. (0 < Nev < N - 1) [in] Tol 停止基準: リッツ値の最大相対誤差. tol <= 0 の場合, マシン精度とみなす. [in,out] Resid() 配列 Resid(LResid - 1) (LResid >= N) [in] Info <> 0 の場合, 初期残差ベクトル. [out] 最終残差ベクトル. [in] Ncv 行列 V の列数. (Nev + 2 < Ncv <= N) 各反復において生成されるアーノルディ基底ベクトルの数を表す. [out] V() 配列 V(Ldv - 1, LV - 1) (Ldv >= N, LV >= Ncv) アーノルディ基底ベクトルの最終値. [in,out] Iparam() 配列 Iparam(LIparam - 1) (LIparam >= 11) Iparam(0) (ishift): [in] シフトの選択方法. = 0: シフトをリバースコミュニケーションによりユーザーが与える. = 1: 現在のヘッセンベルグ行列 H に関するシフトを使う. Iparam(1) 使用されない. Iparam(2) (mxiter): [in] アーノルディ反復の最大回数. [out] 実際の反復回数. Iparam(3) (nb): [in] 漸化式で用いられるブロックサイズ. (現バージョンでは nb = 1 でなければならない) Iparam(4) (nconv): [out] 収束条件を満たすリッツ値の数. Iparam(5) 使用されない. Iparam(6) (mode): [in] 固有値問題のタイプ (mode = 1, ..., 4) Iparam(7) (np): [out] Ido = 3 かつ Iparam(0) = 0 の場合, ユーザーが与えるべきシフトの数を返す. (0 < np <= Ncv - Nev) Iparam(8) (numop): [out] OP*x 操作の回数. Iparam(9) (numopb): [out] B*x 操作の回数 (Bmat = "G" の場合). Iparam(10) (numreo): [out] 再直交化ステップ数. [out] Ipntr() 配列 Ipntr(LIpntr - 1) (LIpntr >= 14) アーノルディ反復で用いられる行列/ベクトルの先頭へのポインター(作業配列中でのインデックス(1から始まる))を示す. Ipntr(0): Workd中のベクトル x へのポインター. Ipntr(1): Workd中のベクトル y へのポインター. Ipntr(2): Workd中のベクトル B*x へのポインター (シフト＆インバートモードの場合). Ipntr(3): Workl中の次に使用可能なプログラムにより未使用な場所へのポインター. Ipntr(4): Workl中の Ncv x Ncv 上ヘッセンベルグ行列 H へのポインター. Ipntr(5): Workl中の Ncv リッツ値の実数部の配列 ritzr へのポインター. Ipntr(6): Workl中の Ncv リッツ値の虚数部の配列 ritzi へのポインター. Ipntr(7): Workl中の リッツ値に関連するリッツ推定値へのポインター. Ipntr(8), ..., Ipntr(12): Dseupd でのみ使用される. Ipntr(13): Workl中の np 個のシフト値へのポインター. [out] Workd() 配列 Workd(LWorkd - 1) (LWorkd >= 3*N) リバースコミュニケーションに使われる分散作業領域. [out] Workl() 配列 Workl(LWorkl - 1) (LWorkl >= 3*Ncv^2 + 6*Ncv) ローカル作業領域. [in,out] Info [in] = 0: 初期残差ベクトルとしてランダムベクトルを使う. <> 0: 初期残差ベクトルを Resid() に与える. 前回実行値のままでもよい. [out] リターンコード. = 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = 1: 最大反復回数に達した. Iparam(4) に収束したリッツ値の数を返す. = 3: リスタート付きアーノルディ反復中にシフト値が適用できなかった. Ncv を Nev に比較して大きくすると改善するかもしれない (Ncv >= 2*Nev が推奨される). = 11: 初期残差ベクトルが 0 である. = 12: アーノルディ分解に失敗した. Iparam(4) に現状の分解の大きさを返す. = 13: LAPACK ルーチンで固有値の計算に失敗した. 出典 ARPACK-NG (https://github.com/opencollab/arpack-ng) 使用例 行列 A について最大の固有値とその固有ベクトルを求める. ただし, ( 0.2 -0.11 -0.93 ) A = ( -0.32 0.81 0.37 ) ( -0.8 -0.92 -0.29 ) とする. Sub Ex_Dnaupd() Const N = 3, Nnz = 9, Ncv = N, Nev = Ncv - 2 Const Rvec = 1, Bmat = "I", Which = "LM" Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim Dr(Nev) As Double, Di(Nev) As Double, Z(N - 1, Nev) As Double Dim V(N - 1, Ncv - 1) As Double, Resid(N - 1) As Double Dim Workd(3 * N - 1) As Double, Workl(3 * Ncv * (Ncv + 2) - 1) As Double Dim Workev(3 * Ncv - 1) As Double, XX(N - 1) As Double, YY(N - 1) As Double Dim Selct(Ncv - 1) As Long, IParam(10) As Long, Ipntr(13) As Long Dim Sigmar As Double, Sigmai As Double Dim Ido As Long, Nconv As Long, Info As Long A(0) = 0.2: A(1) = -0.11: A(2) = -0.93: A(3) = -0.32: A(4) = 0.81: A(5) = 0.37: A(6) = -0.8: A(7) = -0.92: A(8) = -0.29 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 IParam(0) = 1: IParam(2) = 300: IParam(3) = 1: IParam(6) = 1 Info = 0 Ido = 0 Do Call Dnaupd(Ido, Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), Info) If Ido = -1 Or Ido = 1 Then Call DCopy(N, Workd(Ipntr(0) - 1), XX(0)) Call CsrDusmv("N", N, N, 1, A(), Ia(), Ja(), XX(), 0, YY()) Call DCopy(N, YY(0), Workd(Ipntr(1) - 1)) End If Loop While (Ido <> 99) Nconv = IParam(4) If Info <> 0 Then Debug.Print "Dnaupd Info =" + Str(Info) Call Dneupd(Rvec, "A", Selct(), Dr(), Di(), Z(), Sigmar, Sigmai, Workev(), Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), Info) If Info <> 0 Then Debug.Print "Dneupd Info =" + Str(Info) Debug.Print "D =", Dr(0), Di(0) Debug.Print "Z =" Debug.Print Z(0, 0), Z(0, 1) Debug.Print Z(1, 0), Z(1, 1) Debug.Print Z(2, 0), Z(2, 1) Debug.Print "Nconv =" + Str(Nconv) End Sub 実行結果 D = 0.812065011925673 0.48915757543818 Z = -0.365053973674333 -0.416255621362224 0.64300369896322 -0.194693188992647 -5.47400351005566E-02 0.488989966998924 Nconv = 2