◆ ZBicg_s()

Sub ZBicg_s	(	N As	Long,
		Matvec As	LongPtr,
		MatvecTrans As	LongPtr,
		Psolve As	LongPtr,
		PsolveTrans 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`
	)

Solution of linear system Ax = b using bi-conjugate gradient (BICG) method (Complex matrices) (Subroutine version)

Purpose: This routine solves the linear system Ax = b using the bi-conjugate gradient (BICG) iterative method with preconditioning.

Parameters

[in]	N	Dimension of the matrix A. (N >= 0) (if N = 0, returns without computation)
[in]	Matvec	User supplied subroutine which calculates the product of matrix A and vector x as follows. Sub Matvec(N As Long, X() As Complex, Y() As Complex) Compute A*x, and return in Y(). End Sub
[in]	Psolve	User supplied subroutine which perform the preconditioner solve routine for the linear system Mx = b as follows, where M is a preconditioner metrix. Sub Psolve(N As Long, B() As Complex, X() As Complex) Solve Mx = b, and return solution in X(). End Sub
[in]	ChkConv	User supplied subroutine which is called on every iteration for the convergence test as follows, where X() is the current approximate solution, Res is the current residual norm norm(b - A*x), and Iter is the current number of iterations. This routine can also be used to output the intermediate results. Sub ChkConv(N As Long, X() As Complex, Res As Double, Iter As Long, IChk As Long) Set IChk = 1 if converged. Otherwise, set IChk = 0. End Sub
[in]	B()	Array B(LB - 1) (LB >= N) Right hand side vector b.
[in,out]	X()	Array X(LX - 1) (LX >= N) [in] Initial guess of solution. [out] Obtained approximate solution.
[out]	Info	(Optional) = 0: Successful exit = i < 0: The (-i)-th argument is invalid. = 11: Maximum number of iterations exceeded. = 12: Breakdown occurred.
[out]	Iter	(Optional) Actual number of iterations performed for convergence.
[out]	Res	(Optional) Final residual norm value norm(b - A*x).
[in]	MaxIter	(Optional) Maximum number of iterations. (MaxIter > 0) (default = 500)

Example Program: Solve the system of linear equations Ax = B, where
( 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 MatvecTrans(N As Long, X() As Complex, Y() As Complex)

Call CsrZusmv("C", 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 PsolveTrans(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_ZBicg_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 ZBicg_s(N, AddressOf Matvec, AddressOf MatvecTrans, AddressOf Psolve, AddressOf PsolveTrans, 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

Example Results: X =

(0.589999999999999,-0.28)

(-0.200000000000002,-4.00000000000004E-02)

(0.24,-0.489999999999999)

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