◆ ZTfqmr_s()

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

転置不要疑似最小残差(TFQMR)法による連立一次方程式 Ax = b の解 (複素行列) (サブルーチン形式)

目的: 前処理付き反復法(転置不要疑似最小残差(TFQMR)法)により連立一次方程式 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	前処理行列の適用, すなわち方程式 Mx = b の解をを求めるユーザーサブルーチンで, 次のように定義すること. ここで, M は前処理行列である. Sub Psolve(N As Long, B() As Double, X() As Double) Mx = 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]	MaxIter	(省略可) 最大反復回数. (MaxIter > 0) (省略時 = 500)

使用例: 連立一次方程式 Ax = B を解く. ただし,
( 0.78+0.16i -0.9-1.46i 0.48-1.08i )

A = ( 0.73+0.63i 1.58-1.24 -0.41-0.91i )

( 0.23-1.37i 0.79+0.64i -0.73-1.5i )

( 0.2126-0.2904i )

B = ( -0.3028+0.3346i )

( -1.2905-1.0346i )

とする.
Const N = 3, Nnz = N * N, Tol = 0.0000000001 '1.0e-10

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

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

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

End Sub

Sub Psolve(N As Long, B() As Complex, X() As Complex)

Dim I As Long

For I = 0 To N - 1

X(I) = B(I)

Next

End Sub

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

If Res < Tol Then

IChk = 1

Else

IChk = 0

End If

End Sub

Sub Ex_ZTfqmr_s()

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

Dim Iter As Long, Res As Double, Info As Long

A(0) = Cmplx(0.78, 0.16): A(1) = Cmplx(-0.9, -1.46): A(2) = Cmplx(0.48, -1.08): A(3) = Cmplx(0.73, 0.63): A(4) = Cmplx(1.58, -1.24): A(5) = Cmplx(-0.41, -0.91): A(6) = Cmplx(0.23, -1.37): A(7) = Cmplx(0.79, 0.64): A(8) = Cmplx(-0.73, -1.5)

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) = Cmplx(0.2126, -0.2904): B(1) = Cmplx(-0.3028, 0.3346): B(2) = Cmplx(-1.2905, -1.0346)

Call ZTfqmr_s(N, AddressOf Matvec, AddressOf Psolve, AddressOf ChkConv, B(), X(), Info, Iter, Res)

Debug.Print "X ="

Debug.Print "(" + CStr(Creal(X(0))) + "," + CStr(Cimag(X(0))) + ")"

Debug.Print "(" + CStr(Creal(X(1))) + "," + CStr(Cimag(X(1))) + ")"

Debug.Print "(" + CStr(Creal(X(2))) + "," + CStr(Cimag(X(2))) + ")"

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

End Sub

実行結果: X =

(0.59,-0.28)

(-0.2,-4.00000000000002E-02)

(0.24,-0.49)

Iter = 3, Res = 2.36560039920782E-15, Info = 0