|
|
◆ znaupd()
| void znaupd |
( |
int * |
ido, |
|
|
char |
bmat, |
|
|
int |
n, |
|
|
char const * |
which, |
|
|
int |
nev, |
|
|
double |
tol, |
|
|
doublecomplex |
resid[], |
|
|
int |
ncv, |
|
|
int |
ldv, |
|
|
doublecomplex |
v[], |
|
|
int |
iparam[], |
|
|
int |
ipntr[], |
|
|
doublecomplex |
workd[], |
|
|
doublecomplex |
workl[], |
|
|
int |
lworkl, |
|
|
double |
rwork[], |
|
|
int * |
info |
|
) |
| |
アーノルディ分解 (複素疎行列)
- 目的
- リスタート付きアーノルディ法 (Implicitly restarted Arnoldi method (IRAM)) により複素疎行列の固有値・固有ベクトルを求める.
本ルーチンはアルゴリズムの前半部であるアーノルディ分解を行う. 分解結果から求められる近似固有値はリッツ値とよばれ, 対応する近似固有ベクトルはリッツベクトルとよばれる. アルゴリズムの後半部であるリッツ値およびリッツベクトルの計算は, 本ルーチンの出力をもとに zneupd が行う.
本ルーチンは, 半正定値エルミート行列 B によって定義される半内積に関する複素線形演算子「OP」のいくつかの近似固有値対を求める.
本ルーチンは, 以下の問題のどれかを解くために, 収束するまで繰り返し呼び出される.
- 通常モード (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) は 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) |
| [in] | tol | 停止基準: リッツ値の最大相対誤差. tol <= 0 の場合, マシン精度とみなす. |
| [in,out] | resid[] | 配列 resid[lresid] (lresid >= n)
[in] info != 0 の場合, 初期残差ベクトル.
[out] 最終残差ベクトル. |
| [in] | ncv | 行列 V の列数. (nev < ncv <= n)
各反復において生成されるアーノルディ基底ベクトルの数を表す. |
| [in] | ldv | 配列 v の整合寸法. (ldv >= n) |
| [out] | v[] | 配列 v[ldv * lv] (lv >= ncv)
アーノルディ基底ベクトルの最終値. |
| [in,out] | iparam[] | 配列 iparam[liparam] (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, ..., 3)
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] (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中のリッツ値の配列 ritz へのポインター.
ipntr[6]: workl中の(射影された)リッツベクトルの配列 Q へのポインター.
ipntr[7]: workl中の誤差限界の配列へのポインター.
ipntr[8], ..., ipntr[12]: zneupd でのみ使用される.
ipntr[13]: workl中の np 個のシフト値へのポインター. |
| [out] | workd[] | 配列 workd[lworkd] (lworkd >= 3*n)
リバースコミュニケーションに使われる分散作業領域. |
| [out] | workl[] | 配列 workl[lworkl]
ローカル作業領域. |
| [in] | lworkl | 配列 workl[]のサイズ. (lworkl >= 3*ncv^2 + 5*ncv) |
| [out] | rwork[] | 配列 rwork[lrwork] (lrwork >= 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)
|