|
|
◆ ZCsrIlu()
| Sub ZCsrIlu |
( |
N As |
Long, |
|
|
Val() As |
Complex, |
|
|
Rowptr() As |
Long, |
|
|
Colind() As |
Long, |
|
|
P As |
Long, |
|
|
Val2() As |
Complex, |
|
|
Rowptr2() As |
Long, |
|
|
Colind2() As |
Long, |
|
|
D() As |
Complex, |
|
|
Optional Info As |
Long, |
|
|
Optional Md As |
Long = 0, |
|
|
Optional Base As |
Long = -1, |
|
|
Optional Base2 As |
Long = 0 |
|
) |
| |
不完全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 の場合)
lev(a(i,j)) = min{ lev(a(i,j)), lev(a(i,k)) + lev(a(k,j)) + 1 }
Val2(), Rowptr2() および Colind2() に下三角行列 L と上三角行列 U を出力する. また, D() に U の対角要素をコピーする. Val2(), Rowptr2(), Colind2() および D() を CsrIluSolve() が使用する.
- 引数
-
| [in] | N | 行列 A の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) |
| [in] | Val() | 配列 Val(LVal - 1) (LVal >= Nnz) (Nnz は非ゼロ要素数)
行列 A の非ゼロ要素の値. |
| [in] | Rowptr() | 配列 Rowptr(LRowptr - 1) (LRowptr >= N + 1)
行列 A の行ポインタ. |
| [in] | Colind() | 配列 Colind(LColind - 1) (LColind >= Nnz)
行列 A の列インデクス. |
| [in] | P | ILU(p) アルゴリズムのレベル p の値. (P >= 0)
注: P = 0 とすることができるが, CsrIlu0()を使用した方が高速である. |
| [out] | Val2() | 配列 Val2(LVal2 - 1) (LVal2 > 0) 下三角行列 L と上三角行列 U (行列 L*U) の非ゼロ要素の値.
Nnz2 = min(LVal2, LColind2) とすると, Nnz2 は行列 L*U の非ゼロ要素数の最大値である. 分解中に非ゼロ要素数がこれを超えるとエラー番号 -6 で終了する. |
| [out] | Rowptr2() | 配列 Rowptr2(LRowptr2 - 1) (LRowptr2 >= N + 1)
行列 L*U の行ポインタ. |
| [out] | Colind2() | 配列 Colind2(LColind2 - 1) (LColind2 > 0)
行列 L*U の列インデクス. |
| [out] | D() | 配列 D(LD) (LD >= N)
上三角行列 U の対角要素. |
| [out] | Info | (省略可)
= 0: 正常終了.
= i < 0: (-i)番目の入力パラメータの誤り.
= j > 0: 行列が特異である(j番目の対角要素が0). |
| [in] | Md | (省略可)
廃棄した要素の分を対角要素に修正を加えるかどうか指定する. このような修正を行う方法を修正ILU分解(Modified ILU (MILU)分解)という. (省略時 = 0)
= 0: ILU (修正なし).
= 1: MILU (対角要素の修正あり). |
| [in] | Base | (省略可)
Rowptr() および Colind() のインデクス形式.
= 0: 0-ベース(C形式): 開始インデクス値が 0.
= 1: 1-ベース(Fortran形式): 開始インデクス値が 1.
(省略時: Rowptr(0) = 1 であれば 1, そうでなければ 0 とみなす) |
| [in] | Base2 | (省略可)
Rowptr2() および Colind2() のインデクス形式. (省略時 = 0)
= 0: 0-ベース(C形式): 開始インデクス値が 0.
= 1: 1-ベース(Fortran形式): 開始インデクス値が 1. |
- 使用例
- 連立一次方程式 Ax = B を ILU(p) 前処理付き FGMRES 法で解く. ただし,
( 0.78+0.16i -0.9-1.46i 0.48-1.08i )
A = ( 0.73+0.63i 1.58-1.24 -0.41-0.91i )
( 0.23-1.37i 0.79+0.64i -0.73-1.5i )
( 0.2126-0.2904i )
B = ( -0.3028+0.3346i )
( -1.2905-1.0346i )
Sub Ex_ZFgmres_Ilu_Csr()
Const N = 3, Nnz = N * N, NnzM = Nnz, Tol = 0.0000000001 '1.0e-10
Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long
Dim B(N - 1) As Complex, X(N - 1) As Complex
Dim XX(N - 1) As Complex, YY(N - 1) As Complex
Dim Iter As Long, Res As Double, IRev As Long, Info As Long, I As Long
A(0) = Cmplx(0.78, 0.16): A(1) = Cmplx(-0.9, -1.46): A(2) = Cmplx(0.48, -1.08): A(3) = Cmplx(0.73, 0.63): A(4) = Cmplx(1.58, -1.24): A(5) = Cmplx(-0.41, -0.91): A(6) = Cmplx(0.23, -1.37): A(7) = Cmplx(0.79, 0.64): A(8) = Cmplx(-0.73, -1.5) 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(0.2126, -0.2904): B(1) = Cmplx(-0.3028, 0.3346): B(2) = Cmplx(-1.2905, -1.0346) Dim M(NnzM - 1) As Complex, Im(N) As Long, Jm(NnzM - 1) As Long
Dim P As Long, D(N - 1) As Complex
P = 3
Call ZCsrIlu(N, A(), Ia(), Ja(), P, M(), Im(), Jm(), D(), Info) If Info <> 0 Then Debug.Print "Ilu Info =" + Str(Info)
IRev = 0
Do
Call ZFgmres_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res) If IRev = 1 Then '- Matvec
ElseIf IRev = 3 Then '- Psolve
Call ZCsrIluSolve("N", N, M(), Im(), Jm(), D(), YY(), XX(), Info) If Info <> 0 Then Debug.Print "IluSolve Info =" + Str(Info)
ElseIf IRev = 10 Then '- Check convergence
If Res < Tol Then IRev = 11
End If
Loop While IRev <> 0
Debug.Print "X ="
Debug.Print "(" + CStr( Creal(X(0))) + "," + CStr( Cimag(X(0))) + ")" Debug.Print "(" + CStr( Creal(X(1))) + "," + CStr( Cimag(X(1))) + ")" Debug.Print "(" + CStr( Creal(X(2))) + "," + CStr( Cimag(X(2))) + ")" Debug.Print "Iter =" + Str(Iter) + ", Res =" + Str(Res) + ", Info =" + Str(Info)
End Sub
Function Cmplx(R As Double, Optional I As Double=0) As Complex 複素数の作成 Function Cimag(A As Complex) As Double 複素数の虚数部 Function Creal(A As Complex) As Double 複素数の実数部 Sub CsrZusmv(Trans As String, M As Long, N As Long, Alpha As Complex, Val() As Complex, Rowptr() As Long, Colind() As Long, X() As Complex, Beta As Complex, Y() As Complex, Optional Info As Long, Optional Base As Long=-1, Optional IncX As Long=1, Optional IncY As Long=1) y <- αAx + βy, y <- αATx + βy または y <- αAHx + βy (複素行列) (CSR) Sub ZFgmres_r(N As Long, B() As Complex, X() As Complex, Info As Long, XX() As Complex, YY() As Complex, IRev As Long, Optional Iter As Long, Optional Res As Double, Optional M As Long=0, Optional MaxIter As Long=500) 最小残差(FGMRES)法による連立一次方程式 Ax = b の解 (複素行列) (リバースコミュニケーション版) Sub ZCsrIlu(N As Long, Val() As Complex, Rowptr() As Long, Colind() As Long, P As Long, Val2() As Complex, Rowptr2() As Long, Colind2() As Long, D() As Complex, Optional Info As Long, Optional Md As Long=0, Optional Base As Long=-1, Optional Base2 As Long=0) 不完全LU分解(レベル指定)(ILU(p))前処理のための初期化 (複素行列) (CSR) Sub ZCsrIluSolve(ByVal Trans As String, N As Long, Val() As Complex, Rowptr() As Long, Colind() As Long, D() As Complex, B() As Complex, X() As Complex, Optional Info As Long, Optional Base As Long=-1) 不完全LU分解(ILU)前処理 (複素行列) (CSR)
- 実行結果
X =
(0.59,-0.28)
(-0.2,-3.99999999999999E-02)
(0.24,-0.49)
Iter = 1, Res = 1.83807799103727E-16, Info = 0
|
1.9.7