|
|
◆ Gpbicg_s()
| Sub Gpbicg_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 Mode As |
Long = 0, |
|
|
Optional MaxIter As |
Long = 500 |
|
) |
| |
積型双共役勾配(GPBICG)法, 安定化双共役勾配(BICGSTAB)法 または BICGSTAB2法による連立一次方程式 Ax = b の解 (サブルーチン形式)
- 目的
- 前処理付き反復法(積型双共役勾配(GPBICG)法, 安定化双共役勾配(BICGSTAB)法 または BICGSTAB2法)により連立一次方程式 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: 正常終了.
= i < 0: (-i)番目の入力パラメータの誤り.
= 11: 最大反復回数を超えた.
= 12: ブレークダウンが発生した. |
| [out] | Iter | (省略可)
収束時の反復回数. |
| [out] | Res | (省略可)
最終的な残差ノルム norm(b - A*x) の値. |
| [in] | Mode | (省略可)
解法の選択. (省略時 = 0)
= 0: GPBICG.
= 1: BICGSTAB.
= 2: BICGSTAB2. |
| [in] | MaxIter | (省略可)
最大反復回数. (MaxIter > 0) (省略時 = 500) |
- 使用例
- 連立一次方程式 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, 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 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_Gpbicg_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 Gpbicg_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 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) Sub Gpbicg_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 Mode As Long=0, Optional MaxIter As Long=500) 積型双共役勾配(GPBICG)法, 安定化双共役勾配(BICGSTAB)法 または BICGSTAB2法による連立一次方程式 Ax = b の解 (サブルーチン形式)
- 実行結果
X = 0.859999999999996 0.640000000000003 0.509999999999995
Iter = 3, Res = 4.09907879275695E-15, Info = 0
|
1.9.7