◆ SscIc0()

Sub SscIc0	(	N As	Long,
		Val() As	Double,
		Colptr() As	Long,
		Rowind() As	Long,
		Val2() As	Double,
		Idiag() As	Long,
		Optional Info As	Long,
		Optional ByVal Uplo As	String = `"U"`,
		Optional Base As	Long = `-1`
	)

不完全コレスキー分解(フィルインなし)(IC0)前処理のための初期化 (正定値対称行列) (CSC)

目的: 連立一次方程式の正定値対称で疎な係数行列の不完全コレスキー分解(フィルインなし)を求める.
A = L * D * L^T + R または A = U^T * D * U + R

ここで, L は下三角行列, U は上三角行列, D は対角行列である.
係数行列 A の要素が 0 であった場所がコレスキー分解後に 0 でなくなることをフィルインが生じるというが, それを避けるために上式の L または U ではそのような場所を強制的に 0 にすることにする. R はこの操作による完全なコレスキー分解との差分を表すが, R が大きくはないものとすると, この不完全分解を使って連立一次方程式を解いてもよい近似解になっていることが期待される. これを前処理行列 M として使用する.
M = L * D * L^T または M = U^T * D * U

本ルーチンは Val2() に L または U および D を出力する. また, Idiag() に対角要素のインデックスを出力する. Val2() および Idiag() を SscIcSolve() が使用する.

引数

[in]	N	行列 A の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in]	Val()	配列 Val(LVal - 1) (LVal >= Nnz) (Nnz は行列 A の非ゼロ要素数) 行列 A の非ゼロ要素の値. (対角要素と下三角要素のみ参照する)
[in]	Colptr()	配列 Colptr(LColptr - 1) (LColptr >= N + 1) 行列 A の列ポインタ.
[in]	Rowind()	配列 Rowind(LRowind - 1) (LRowind >= Nnz) 行列 A の行インデクス.
[out]	Val2()	配列 Val2(LVal2 - 1) (LVal2 >= Nnz) 前処理行列 M の非ゼロ要素の値. (A の対角要素と下三角要素と同じ場所に書き込まれる)
[out]	Idiag()	配列 Idiag(LIdiag) (LIdiag >= N) 対角要素のインデックス.
[out]	Info	(省略可) = 0: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = j > 0: 行列が特異である(j番目の対角要素が0).
[in]	Uplo	(省略可) 入力行列 A が上三角部分あるいは下三角部分のどちらを格納しているか指定する. (省略時 = "U") = "U": 上三角部分 (U^TDU と分解する). = "L": 下三角部分 (LDL^T と分解する).
[in]	Base	(省略可) Colptr() および Rowind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Colptr(0) = 1 であれば 1, そうでなければ 0 とみなす)

使用例: 連立一次方程式 Ax = B を IC(0) 前処理付き CG 法で解く. ただし,
( 2.2 -0.11 -0.32 ) ( -1.566 )

A = ( -0.11 2.93 0.81 ), B = ( -2.8425 )

( -0.32 0.81 2.37 ) ( -1.1765 )

とする.
Sub Ex_Cg_Ic0_Ssc()

Const N = 3, Nnz = 6, 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) = 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

B(0) = -1.566: B(1) = -2.8425: B(2) = -1.1765

Dim M(Nnz - 1) As Double, Id(N - 1) As Long

Call SscIc0(N, A(), Ia(), Ja(), M(), Id(), Info)

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

IRev = 0

Do

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

If IRev = 1 Then '- Matvec

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

ElseIf IRev = 3 Then '- Psolve

Call SscIcSolve(N, M(), Ia(), Ja(), Id(), YY(), XX(), Info)

If Info <> 0 Then Debug.Print "Ic0Solve 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

SscDusmv
Sub SscDusmv(Uplo As String, N As Long, Alpha As Double, Val() As Double, Colptr() As Long, Rowind() 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 (CSC) (対称行列)

Cg_r
Sub Cg_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 MaxIter As Long=500)
共役勾配(CG)法による連立一次方程式 Ax = b の解 (正定値対称行列) (リバースコミュニケーション版)

SscIc0
Sub SscIc0(N As Long, Val() As Double, Colptr() As Long, Rowind() As Long, Val2() As Double, Idiag() As Long, Optional Info As Long, Optional ByVal Uplo As String="U", Optional Base As Long=-1)
不完全コレスキー分解(フィルインなし)(IC0)前処理のための初期化 (正定値対称行列) (CSC)

SscIcSolve
Sub SscIcSolve(N As Long, Val() As Double, Colptr() As Long, Rowind() As Long, Idiag() As Long, B() As Double, X() As Double, Optional Info As Long, Optional ByVal Uplo As String="U", Optional Base As Long=-1)
不完全コレスキー分解による前処理 (IC) (正定値対称行列) (CSC)

実行結果: X = -0.8 -0.92 -0.29

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