ZCg ZCg_r ZCg_s ZCr ZCr_r ZCr_s

ZCg_r() Sub ZCg_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 conjugate gradient (CG) method (Hermitian positive definite matrices) (Reverse communication version) Purpose This routine solves the linear system Ax = b with Hermitian positive definite coefficient matrix using the conjugate gradient (CG) 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 < 0: The (-Info)-th argument is invalid. = 1: (Warning) Matrix A is not positive definite (computation continued). = 2: (Warning) Preconditioner matrix M is not positive definite (computation continued). = 11: Maximum number of iterations exceeded. = 12: Matrix A is singular (zero diagonal element). [in,out] XX() Array XX(LXX - 1) (LXX >= N) Vector XX for Matvec and Psolve operations. [in,out] YY() Array YY(LYY - 1) (LYY >= N) Vector YY for Matvec and Psolve 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 A*XX in YY. Do not alter any other variables. = 3: Psolve operation. User should set solution of M*XX = 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 ( 1.4 -1.5+0.46i 0.16+0.23i ) A = ( -1.5-0.46i 1.44 -0.12+0.04i ) ( 0.16-0.23i -0.12-0.04i 0.05 ) ( -2.3215-1.1316i ) B = ( 1.7972+2.0692i ) ( -0.4042-0.0049i ) Sub Ex_ZCg_r() Const N = 3, Nnz = 6, 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(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) IRev = 0 Do Call ZCg_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res) If IRev = 1 Then '- Matvec Call HsrZusmv("L", N, Cmplx(1), A(), Ia(), Ja(), XX(), Cmplx(0), YY()) ElseIf IRev = 3 Then '- Psolve 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 Example Results X = (-0.819999999999983,-0.939999999999992) (0.739999999999985,0.199999999999988) (0.480000000000001,0.209999999999999) Iter = 3, Res = 8.08057086495744E-14, Info = 1