◆ CsrIlu0()

Sub CsrIlu0	(	N As	Long,
		Val() As	Double,
		Rowptr() As	Long,
		Colind() As	Long,
		Val2() As	Double,
		D() As	Double,
		Optional Info As	Long,
		Optional Base As	Long = `-1`
	)

不完全LU分解(フィルインなし)(ILU0)前処理のための初期化 (CSR)

目的: 連立一次方程式の疎な係数行列 A の不完全LU分解(フィルインなし)を求める.
A = L * U + R

ここで, Rは完全なLU分解との差分であるが, Rが小さいものとしてこの分解を使って連立一次方程式を解くことにすると次の前処理行列が得られる.
M = L * U

Val2() に下三角行列 L と上三角行列 U を出力する. また, D() に U の対角要素をコピーする. Val2() および 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 の列インデクス.
[out]	Val2()	配列 Val2(LVal2 - 1) (LVal2 >= Nnz) 下三角行列 L と上三角行列 U の非ゼロ要素の値. (A の下三角および上三角要素と同じ場所に書き込まれる)
[out]	D()	配列 D(LD) (LD >= N) 上三角行列 U の対角要素.
[out]	Info	(省略可) = 0: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = j > 0: 行列が特異である(j番目の対角要素が0).
[in]	Base	(省略可) Rowptr() および Colind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Rowptr(0) = 1 であれば 1, そうでなければ 0 とみなす)

使用例: 連立一次方程式 Ax = B を DS 前処理付き FGMRES 法で解く. ただし,
( 0.2 -0.11 -0.93 ) ( -0.3727 )

A = ( -0.32 0.81 0.37 ), B = ( 0.4319 )

( -0.8 -0.92 -0.29 ) ( -1.4247 )

とする.
Sub Ex_Fgmres_Ilu0_Csr()

Const N = 3, Nnz = N * N, Tol = 0.0000000001 '1.0e-10

Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long

Dim B(N - 1) As Double, X(N - 1) As Double

Dim XX(N - 1) As Double, YY(N - 1) As Double

Dim Iter As Long, Res As Double, IRev 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

B(0) = -0.3727: B(1) = 0.4319: B(2) = -1.4247

Dim M(Nnz - 1) As Double, D(N - 1) As Double

Call CsrIlu0(N, A(), Ia(), Ja(), M(), D(), Info)

If Info <> 0 Then Debug.Print "Ilu0 Info =" + Str(Info)

IRev = 0

Do

Call Fgmres_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res)

If IRev = 1 Then '- Matvec

Call CsrDusmv("N", N, N, 1, A(), Ia(), Ja(), XX(), 0, YY())

ElseIf IRev = 3 Then '- Psolve

Call CsrIluSolve("N", N, M(), Ia(), Ja(), 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 =", X(0), X(1), X(2)

Debug.Print "Iter = " + CStr(Iter) + ", Res = " + CStr(Res) + ", Info = " + CStr(Info)

End Sub

CsrDusmv
Sub CsrDusmv(Trans As String, M As Long, N As Long, Alpha As Double, Val() As Double, Rowptr() As Long, Colind() As Long, X() As Double, Beta As Double, Y() As Double, 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 (CSR)

Fgmres_r
Sub Fgmres_r(N As Long, B() As Double, X() As Double, Info As Long, XX() As Double, YY() As Double, IRev As Long, Optional Iter As Long, Optional Res As Double, Optional M As Long=0, Optional MaxIter As Long=500)
最小残差(FGMRES)法による連立一次方程式 Ax = b の解 (リバースコミュニケーション版)

CsrIlu0
Sub CsrIlu0(N As Long, Val() As Double, Rowptr() As Long, Colind() As Long, Val2() As Double, D() As Double, Optional Info As Long, Optional Base As Long=-1)
不完全LU分解(フィルインなし)(ILU0)前処理のための初期化 (CSR)

CsrIluSolve
Sub CsrIluSolve(ByVal Trans As String, N As Long, Val() As Double, Rowptr() As Long, Colind() As Long, D() As Double, B() As Double, X() As Double, Optional Info As Long, Optional Base As Long=-1)
不完全LU分解による前処理 (ILU) (CSR)

実行結果: X = 0.86 0.64 0.51

Iter = 1, Res = 2.00597268669339E-16, Info = 0