◆ Sor_s()

Sub Sor_s	(	N As	Long,
		Matvec As	LongPtr,
		Matsol As	LongPtr,
		ChkConv As	LongPtr,
		B() As	Double,
		X() As	Double,
		Optional Info As	Long,
		Optional Iter As	Long,
		Optional Res As	Double,
		Optional MaxIter As	Long = `500`
	)

逐次的過剰緩和(SOR)法による連立一次方程式 Ax = b の解 (サブルーチン形式)

目的: 反復法(逐次的過剰緩和(SOR)法)により連立一次方程式 Ax = b の解を求める.

引数

[in]	N	行列 A の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in]	Matvec	行列 A とベクトル x の積をを求めるユーザーサブルーチンで, 次のように定義すること. Sub Matvec(N As Long, X() As Double, Y() As Double) A*x を計算し Y() に入れる. End Sub
[in]	Matsol	行列 A の対角および下三角部分を使って (D/ω + L)^(-1)b を求めるユーザーサブルーチンで, 次のように定義すること. Sub Matsol(N As Long, B() As Double, X() As Double) x = (D/ω + L)^(-1)b を求め (あるいは (D/ω + L)*x = b を解き) 解 x を X() に返す, End Sub
[in]	ChkConv	反復ごとに呼び出され収束判定を行うユーザーサブルーチンで, 次のように定義すること. ここで, X() は現在の近似解, Res は現在の残差ノルム norm(b - A*x), Iterは現在の反復回数である. 本ルーチンは中間結果を出力するために使用することもできる. Sub ChkConv(N As Long, X() As Double, Res As Double, Iter As Long, IChk As Long) 収束であれば IChk = 1, そうでなければ IChk = 0 に設定する. End Sub
[in]	B()	配列 B(LB - 1) (LB >= N) 右辺ベクトル b.
[in,out]	X()	配列 X(LX - 1) (LX >= N) [in] 解の初期推定値. [out] 求められた近似解.
[out]	Info	(省略可) = 0: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = 11: 最大反復回数を超えた. = 12: 行列が特異である(対角要素が0).
[out]	Iter	(省略可) 収束時の反復回数.
[out]	Res	(省略可) 最終的な残差ノルム norm(b - A*x) の値.
[in]	MaxIter	(省略可) 最大反復回数. (MaxIter > 0) (省略時 = 500)

参考: (D/ω + L)^(-1)*b の計算は CscDussv または CsrDussv を用いて行うことができる.

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

とする.
Const N = 3, Nnz = N * N, Omega = 0.4, Tol = 0.0000000001

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

Sub Matvec(N As Long, X() As Double, Y() As Double)

Call CsrDusmv("N", N, N, 1, A(), Ia(), Ja(), X(), 0, Y())

End Sub

Sub Matsol(N As Long, B() As Double, X() As Double)

Dim I As Long

For I = 0 To N - 1

X(I) = B(I)

Next

Call CsrDussv("L", "N", "N", N, A(), Ia(), Ja(), X(), Omega:=Omega)

End Sub

Sub ChkConv(N As Long, X() As Double, Res As Double, Iter As Long, IChk As Long)

If Res < Tol Then

IChk = 1

Else

IChk = 0

End If

End Sub

Sub Ex_Sor_s()

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

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

Call Sor_s(N, AddressOf Matvec, AddressOf Matsol, AddressOf ChkConv, B(), X(), Info, Iter, Res)

Debug.Print "X =", X(0), X(1), X(2)

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

End Sub

CsrDussv
Sub CsrDussv(Uplo As String, Trans As String, Diag As String, N As Long, Val() As Double, Rowptr() As Long, Colind() As Long, X() As Double, Optional Info As Long, Optional Base As Long=-1, Optional IncX As Long=1, Optional Omega As Double=1)
Ax = b または ATx = b の解 (三角行列) (CSR)

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)

Sor_s
Sub Sor_s(N As Long, Matvec As LongPtr, Matsol As LongPtr, ChkConv As LongPtr, B() As Double, X() As Double, Optional Info As Long, Optional Iter As Long, Optional Res As Double, Optional MaxIter As Long=500)
逐次的過剰緩和(SOR)法による連立一次方程式 Ax = b の解 (サブルーチン形式)

実行結果: X = 0.859999999965658 0.640000000003731 0.509999999994652

Iter = 57, Res = 2.83730635879748E-11, Info = 0