Dgsev Dgsgv Dnaupd Dneupd Dsaupd Dseupd Dssev Dssgv

Dsaupd() Sub Dsaupd ( 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 Lanczos method: IRLM) により実対称疎行列の固有値・固有ベクトルを求める. 本ルーチンはアルゴリズムの前半部であるランチョス分解を行う. ランチョス分解から求められる近似固有値はリッツ値とよばれ, 対応する近似固有ベクトルはリッツベクトルとよばれる. アルゴリズムの後半部であるリッツ値およびリッツベクトルの計算は, 本ルーチンの出力をもとに Dseupd が行う. 本ルーチンは, 半正定値の実対称行列 B について対称な実線形演算子「OP」のいくつかの近似固有値対を求める. B*OP = (OP^T)*B または <x, OPy> = <OPx, y> ただし <z, w> = z^T*B*w. 本ルーチンは, 以下の問題のどれかを解くために, 収束するまで繰り返し呼び出される. 通常モード (mode = 1): A*x = λ*x, A は対称行列. (OP = A, B = I) インバースモード (mode = 2): A*x = λ*M*x, A は対称行列, M は正定値対称行列. (OP = M^(-1)*A, B = M) シフト＆インバートモード (mode = 3): K*x = λ*M*x, K は対称行列, M は半正定値対称行列. (OP = (K - σ*M)^(-1)*M, B = M) バックリングモード (mode = 4): K*x = λ*KG*x, K は半正定値対称行列, KG は不定値対称行列. (OP = (K - σ*KG)^(-1)*K, B = K) ケーリー変換モード (mode = 5): A*x = λ*M*x, A は対称行列, M は半正定値対称行列. (OP = (A - σ*M)^(-1)*(A + σ*M), B = M) 注1 - モード (mode) は Iparam(6) で指定される. 注2 - M がコレスキー分解できる場合 (M = L*L^T), mode = 2 を選ばずに, 代わりに OP = L^(-1)*A*L^(-T) として mode = 1 を使用するとよい. 収束後, 元の問題の固有ベクトル z は L^T*z = x を解くことにより求めることができる (x は OP のリッツベクトルとする). 引数 [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 = "LA": 最大のものから Nev 個の固有値を求める. = "SA": 最小のものから Nev 個の固有値を求める. = "LM": 絶対値最大のものから Nev 個の固有値を求める. = "SM": 絶対値最小のものから Nev 個の固有値を求める. = "BE": スペクトルの両端から半分づつ, 合計 Nev 個の固有値を求める. Nev が奇数の場合, 大きい側を1個多く求める. [in] Nev 求める固有値の数. (0 < Nev < N) [in] Tol 停止基準: リッツ値の最大相対誤差. Tol <= 0 の場合, マシン精度とみなす. [in,out] Resid() 配列 Resid(LResid - 1) (LResid >= N) [in] Info <> 0 の場合, 初期残差ベクトル. [out] 最終残差ベクトル. [in] Ncv 行列 V の列数. (Nev < Ncv <= N) 各反復において生成されるランチョス基底ベクトルの数を表す. [out] V() 配列 V(Ldv - 1, LV - 1) (Ldv >= N, LV >= Ncv) Ncv 本のランチョス基底ベクトルからなる行列 V. [in,out] Iparam() 配列 Iparam(LIparam - 1) (LIparam >= 11) Iparam(0) (ishift): [in] シフトの選択方法. = 0: シフトをリバースコミュニケーションによりユーザーが与える. = 1: 縮小された三重対角行列 T に関するシフトを使う. 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, ..., 5) 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() Array Ipntr(LIpntr - 1) (LIpntr >= 11) ランチョス反復で用いられる行列/ベクトルの先頭へのポインター(作業配列中でのインデックス(1から始まる))を示す. Ipntr(0): Workd中のベクトル x へのポインター. Ipntr(1): Workd中のベクトル y へのポインター. Ipntr(2): Workd中のベクトル B*x へのポインター (シフト＆インバートモードの場合). Ipntr(3): Workl中の次に使用可能なプログラムにより未使用な場所へのポインター. Ipntr(4): Workl中の Ncv x 2 三重対角行列 T へのポインター. Ipntr(5): Workl中の Ncv リッツ値の配列へのポインター. Ipntr(6): Workl中の リッツ値に関連するリッツ推定値へのポインター. Ipntr(7), Ipntr(8), Ipntr(9): Dseupd でのみ使用される. Ipntr(10): Workl中の np 個のシフト値へのポインター. [out] Workd() 配列 Workd(LWorkd - 1) (LWorkd >= 3*N) リバースコミュニケーションに使われる分散作業領域. [out] Workl() 配列 Workl(LWorkl - 1) (LWorkl >= Ncv^2 + 8*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 ルーチン DSTEQR). 出典 ARPACK-NG (https://github.com/opencollab/arpack-ng) 使用例 対称行列 A について大きい方から2個の固有値とその固有ベクトルを求める. ただし, ( 2.2 -0.11 -0.32 ) A = ( -0.11 2.93 0.81 ) ( -0.32 0.81 2.37 ) とする. Sub Ex_Dsaupd() Const N = 3, Nnz = 6, Ncv = N, Nev = Ncv - 1 Const Rvec = 1, Bmat = "I", Which = "LM" Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim D(Nev - 1) As Double, Z(N - 1, Nev - 1) As Double Dim V(N - 1, Ncv - 1) As Double, Resid(N - 1) As Double Dim Workd(3 * N - 1) As Double, Workl(Ncv * (Ncv + 8) - 1) As Double Dim XX(N - 1) As Double, YY(N - 1) As Double Dim Selct(Ncv - 1) As Long, IParam(10) As Long, Ipntr(10) As Long, Sigma As Double Dim Ido As Long, Nconv As Long, Info As Long A(0) = 2.2: A(1) = -0.11: A(2) = 2.93: A(3) = -0.32: A(4) = 0.81: A(5) = 2.37 Ia(0) = 0: Ia(1) = 1: Ia(2) = 3: Ia(3) = 6 Ja(0) = 0: Ja(1) = 0: Ja(2) = 1: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2 IParam(0) = 1: IParam(2) = 300: IParam(3) = 1: IParam(6) = 1 Info = 0 Ido = 0 Do Call Dsaupd(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 SsrDusmv("L", 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 "Dsaupd Info =" + Str(Info) Call Dseupd(Rvec, "A", Selct(), D(), Z(), Sigma, Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), Info) If Info <> 0 Then Debug.Print "Dseupd Info =" + Str(Info) Debug.Print "D =", D(0), D(1) 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 = 2.22943643244226 3.56350401844728 Z = -0.894521385341403 -0.200931256445799 -0.390994588677271 0.78468897753835 0.216690384632057 0.586421212707156 Nconv = 2