◆ 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`
	)

Solution of linear system Ax = b using bi-conjugate gradient (BICG) method (Complex matrices) (Reverse communication 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]	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	= 0: Successful exit = i < 0: The (-i)-th argument is invalid. = 11: Maximum number of iterations exceeded. = 12: Breakdown occurred.
[in,out]	XX()	Array XX(LXX - 1) (LXX >= N) Vector XX for Matvec, MatvecTrans, Psolve and PsolveTrans operations.
[in,out]	YY()	Array YY(LYY - 1) (LYY >= N) Vector YY for Matvec, MatvecTrans, Psolve and PsolveTrans operations.
[in,out]	IRev	Control variable for reverse communication. [in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered (except if converged). [out] If IRev is not zero, complete the following process and call this routine again. = 0: Computation finished. See return code in info. = 1: Matvec operation. User should set AXX in YY. Do not alter any other variables. = 3: Psolve operation. User should set solution of MXX = YY in XX. Do not alter any other variables. = 10: To be returned for the convergence test on every iteration . Set IRev = 11 if converged. Do not alter IRev otherwise. The latest values in X(), Iter and Res can be used to decide the convergence. Further, these values may be used to output the intermediate results.
[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 )

Sub Ex_ZBicg_r()

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

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

Dim XX(N - 1) As Complex, YY(N - 1) As Complex

Dim Iter As Long, Res As Double, IRev As Long, Info As Long, I 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)

IRev = 0

Do

Call ZBicg_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res)

If IRev = 1 Then '- Matvec

Call CsrZusmv("N", N, N, Cmplx(1), A(), Ia(), Ja(), XX(), Cmplx(0), YY())

ElseIf IRev = 2 Then '- MatvecTrans

Call CsrZusmv("C", N, N, Cmplx(1), A(), Ia(), Ja(), XX(), Cmplx(0), YY())

ElseIf IRev = 3 Or IRev = 4 Then '- Psolve, PsolveTrans

For I = 0 To N - 1

XX(I) = YY(I)

Next

ElseIf IRev = 10 Then '- Check convergence

If Res < Tol Then IRev = 11

End If

Loop While IRev <> 0

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
Building complex number

Cimag
Function Cimag(A As Complex) As Double
Imaginary part of complex number

Creal
Function Creal(A As Complex) As Double
Real part of complex number

CsrZusmv
Sub CsrZusmv(Trans As String, M As Long, N As Long, Alpha As Complex, Val() As Complex, Rowptr() As Long, Colind() As Long, X() As Complex, Beta As Complex, Y() As Complex, 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 or y <- αAHx + βy (Complex matrices) (CSR)

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)
Solution of linear system Ax = b using bi-conjugate gradient (BICG) method (Complex matrices) (Revers...

Example Results: X =

(0.589999999999999,-0.28)

(-0.200000000000002,-4.00000000000004E-02)

(0.24,-0.489999999999999)

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