|
|
◆ Cg_s()
| Sub Cg_s |
( |
N As |
Long, |
|
|
Matvec As |
LongPtr, |
|
|
Psolve 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 |
|
) |
| |
共役勾配(CG)法による連立一次方程式 Ax = b の解 (正定値対称行列) (サブルーチン形式)
- 目的
- 前処理付き反復法(共役勾配(CG)法)により正定値対称な係数行列の連立一次方程式 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] | Psolve | 前処理行列の適用, すなわち方程式 M*x = b の解をを求めるユーザーサブルーチンで, 次のように定義すること. ここで, M は前処理行列である.
Sub Psolve(N As Long, B() As Double, X() As Double)
M*x = b の解を求め 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: 正常終了.
< 0: (-Info)番目の入力パラメータの誤り.
= 1: (警告) 行列 A が正定値でない(計算は続行する).
= 2: (警告) 前処理行列 M が正定値でない(計算は続行する).
= 11: 最大反復回数を超えた.
= 12: 行列 A が特異である(対角要素が 0). |
| [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 )
Const N = 3, Nnz = 6, 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 SsrDusmv("L", N, 1, A(), Ia(), Ja(), X(), 0, Y()) End Sub
Sub Psolve(N As Long, B() As Double, X() As Double)
Dim I As Long
For I = 0 To N - 1
X(I) = B(I)
Next
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_Cg_s()
Dim B(N - 1) As Double, X(N - 1) As Double
Dim Iter As Long, Res As Double, 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
Call Cg_s(N, AddressOf Matvec, AddressOf Psolve, 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
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) (対称行列) Sub Cg_s(N As Long, Matvec As LongPtr, Psolve 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) 共役勾配(CG)法による連立一次方程式 Ax = b の解 (正定値対称行列) (サブルーチン形式)
- 実行結果
X = -0.8 -0.92 -0.29
Iter = 3, Res = 1.03670172979888E-16, Info = 0
|
1.9.7