|
|
◆ csr_ilu()
| void csr_ilu |
( |
int |
n, |
|
|
const double |
val[], |
|
|
const int |
rowptr[], |
|
|
const int |
colind[], |
|
|
int |
base, |
|
|
int |
p, |
|
|
int |
md, |
|
|
int |
nnz2, |
|
|
double |
val2[], |
|
|
int |
rowptr2[], |
|
|
int |
colind2[], |
|
|
int |
base2, |
|
|
double |
d[], |
|
|
double |
work[], |
|
|
int |
iwork[], |
|
|
int * |
info |
|
) |
| |
不完全LU分解(レベル指定) (ILU(p)) (CSR)
- 目的
- 連立一次方程式の疎な係数行列 A の不完全LU分解を求める. ここで, Rは完全なLU分解との差分であるが, Rが小さいものとしてこの分解を使って連立一次方程式を解くことにすると次の前処理行列が得られる. フィルインを抑えるために, 分解により生じる新たな非ゼロ要素のうち指定された値 p よりも大きなレベル値のものを廃棄する. これを ILU(p) と呼ぶ.
まず, レベル値の初期値を次のように設定する. lev(a(i,j)) = 0 (a(i,j) != 0 または i = j の場合),
= ∞ (a(i,j) = 0 の場合)
LU分解においては消去の過程で a(i,j) = a(i,j) - a(i,k)*a(k,j) という演算を行うが, このとき, a(i,j)のレベル値を次のように更新する. lev(a(i,j)) = min{ lev(a(i,j)), lev(a(i,k)) + lev(a(k,j)) + 1 }
このようにすると, 分解前の A の非ゼロ要素の場所は最後までレベル0のままになる. 分解後にレベル値が p よりも大きい要素は廃棄される.
val2[], colind2[] および rowptr2[] に下三角行列 L と上三角行列 U を出力する. また, d[] に U の対角要素をコピーする. val2[], colind2[], rowptr2[] および d[] を csr_ilu_solve() が使用する.
- 引数
-
| [in] | n | 行列 A の次数. (n >= 0) (n = 0 の場合, 処理を行わずに戻る) |
| [in] | val[] | 配列 val[lval] (lval >= nnz) (nnz は行列 A の非ゼロ要素数)
行列 A の非ゼロ要素の値. |
| [in] | rowptr[] | 配列 rowptr[lrowptr] (lrowptr >= n + 1)
行列 A の行ポインタ. |
| [in] | colind[] | 配列 colind[lcolind] (lcolind >= nnz)
行列 A の列インデクス. |
| [in] | base | rowptr[] および colind[] のインデクス形式.
= 0: 0-ベース(C形式): 開始インデクス値が 0.
= 1: 1-ベース(Fortran形式): 開始インデクス値が 1. |
| [in] | p | ILU(p) アルゴリズムのレベル p の値. (p >= 0)
注: p = 0 とすることができるが, csr_ilu0()を使用した方が高速である. |
| [in] | md | 廃棄した要素の分を対角要素に修正を加えるかどうか指定する. このような修正を行う方法を修正ILU分解(Modified ILU (MILU)分解)という.
= 0: ILU (修正なし).
= 1: MILU (対角要素の修正あり).
(上記以外の値であれば md = 0 とみなす)) |
| [in] | nnz2 | 行列 L*U の非ゼロ要素数の最大値. 分解中に非ゼロ要素数がこれを超えるとエラー(info = -8)で終了する. |
| [out] | val2[] | 配列 val2[lval2] (lval2 >= nnz2)
下三角行列 L と上三角行列 U (行列 L*U) の非ゼロ要素の値. |
| [out] | rowptr2[] | 配列 rowptr2[lrowptr2] (lrowptr2 >= n + 1)
行列 L*U の列ポインタ. |
| [out] | colind2[] | 配列 colind2[lcolind2] (lcolind2 >= nnz2)
行列 L*U の行インデクス. |
| [in] | base2 | rowptr2[] および colind2[] のインデクス形式.
= 0: 0-ベース(C形式): 開始インデクス値が 0.
= 1: 1-ベース(Fortran形式): 開始インデクス値が 1. |
| [out] | d[] | 配列 d[ld] (ld >= n)
上三角行列 U の対角要素. |
| [out] | work[] | 配列 work[lwork] (lwork >= n)
作業領域. |
| [out] | iwork[] | 配列 iwork[liwork] (liwork >= nnz2 + 2*n)
整数作業領域. |
| [out] | info | = 0: 正常終了.
= i < 0: (-i)番目の入力パラメータの誤り.
= j > 0: 行列が特異である(j番目の対角要素が0). |
|
1.9.7