◆ ZBicg_r()

Sub ZBicg_r	(	N As	Long,
		B() As	Complex,
		X() As	Complex,
		Info As	Long,
		XX() As	Complex,
		YY() As	Complex,
		IRev As	Long,
		Optional Iter As	Long,
		Optional Res As	Double,
		Optional MaxIter As	Long = `500`
	)

双共役勾配(BICG)法による連立一次方程式 Ax = b の解 (複素行列) (リバースコミュニケーション版)

目的: 前処理付き反復法(双共役勾配(BICG)法)により連立一次方程式 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: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = 11: 最大反復回数を超えた. = 12: ブレークダウンが発生した.
[in,out]	XX()	配列 XX(LXX - 1) (LXX >= N) Matvec, MatvecTrans, Psolve および PsolveTrans演算のためのベクトル XX.
[in,out]	YY()	配列 YY(LYY - 1) (LYY >= N) Matvec, MatvecTrans, Psolve および PsolveTrans演算のためのベクトル 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 を解く. ただし,
( 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 )

とする.
Sub Ex_ZCg()

Const N = 3, Nnz = 6

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

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

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

A(0) = Cmplx(1.4): A(1) = Cmplx(-1.5, -0.46): A(2) = Cmplx(1.44): A(3) = Cmplx(0.16, -0.23): A(4) = Cmplx(-0.12, -0.04): A(5) = Cmplx(0.05)

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) = Cmplx(-2.3215, -1.1316): B(1) = Cmplx(1.7972, 2.0692): B(2) = Cmplx(-0.4042, -0.0049)

Call ZCg(N, A(), Ia(), Ja(), B(), X(), Info, Iter, Res, Uplo:="L")

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

Cmplx
Function Cmplx(R As Double, Optional I As Double=0) As Complex
複素数の作成

Cimag
Function Cimag(A As Complex) As Double
複素数の虚数部

Creal
Function Creal(A As Complex) As Double
複素数の実数部

ZCg
Sub ZCg(N As Long, Val() As Complex, Ptr() As Long, Ind() As Long, B() As Complex, X() As Complex, Optional Info As Long, Optional Iter As Long, Optional Res As Double, Optional Format As Long=0, Optional MaxIter As Long=500, Optional Tol As Double=1.0E-10, Optional Uplo As String="F", Optional Base As Long=-1, Optional Precon As Long=0, Optional Omega As Double=1.5)
共役勾配(CG)法による連立一次方程式 Ax = b の解 (正定値エルミート行列) (ドライバ)

実行結果: X =

(0.589999999999999,-0.28)

(-0.200000000000002,-4.00000000000004E-02)

(0.24,-0.489999999999999)

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