◆ 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 の解 (正定値対称行列) (リバースコミュニケーション版)

目的: 前処理付き反復法(共役勾配(CG)法)により正定値対称な係数行列の連立一次方程式 Ax = b の解を求める.

引数

[in]	N	行列 A の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in]	B()	配列 B(LB - 1) (LB >= N) 右辺ベクトル b.
[in,out]	X()	配列 X(LX - 1) (LX >= N) [in] 解の初期推定値. [out] 求められた近似解.
[out]	Info	= 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = 1: (警告) 行列 A が正定値でない(計算は続行する). = 2: (警告) 前処理行列 M が正定値でない(計算は続行する). = 11: 最大反復回数を超えた. = 12: 行列 A が特異である(対角要素が 0).
[in,out]	XX()	配列 XX(LXX - 1) (LXX >= N) Matvec および Psolve演算のためのベクトル XX.
[in,out]	YY()	配列 YY(LYY - 1) (LYY >= N) Matvec および Psolve演算のためのベクトル YY.
[in,out]	IRev	リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に IRev = 0 に設定しておくこと. それ以降の呼び出し時にはIRevの値を変更してはならない(収束時を除く). [out] 0 以外のときには下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかは Info をチェックすること. = 1: Matvec演算. AXX を YY に設定すること. 他の変数を変更してはならない. = 3: Psolve演算: MXX = YY の解を XX に設定すること(Mは前処理行列). 他の変数を変更してはならない. = 10: 収束判定を行うために反復ごとに戻る. 収束時には IRev = 11 として, それ以外のときは IRev を変更せずに再度呼び出すこと. X(), Iter および Res にはその反復時の最新の値が入っているので収束判定のために使用してよい. また, 中間結果出力のために使うこともできる.
[out]	Iter	(省略可) 収束時の反復回数.
[out]	Res	(省略可) 最終的な残差ノルム norm(b - A*x) の値.
[in]	MaxIter	(省略可) 最大反復回数. (MaxIter > 0) (省略時 = 500)

使用例: 連立一次方程式 Ax = B を解く. ただし,
( 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_r()

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

IRev = 0

Do

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

If IRev = 1 Then '- Matvec

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

ElseIf IRev = 3 Then '- Psolve

Call Dcopy(N, YY(0), XX(0))

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

SsrDusmv
Sub SsrDusmv(Uplo As String, 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 (CSR) (対称行列)

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 の解 (正定値対称行列) (リバースコミュニケーション版)

Dcopy
Sub Dcopy(N As Long, X_I As Double, Y_I As Double, Optional IncX As Long=1, Optional IncY As Long=1)
y <- x (BLAS 1)

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

Iter = 3, Res = 1.03670172979888E-16, Info = 0