Cg Cg1 Cg_r Cg_s Cr Cr_r Cr_s

Cg_r() Sub Cg_r ( N As  Long, B() As  Double, X() As  Double, Info As  Long, XX() As  Double, YY() As  Double, 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 (symmetric positive definite matrix) (Reverse communication version) Purpose This routine solves the linear system Ax = b with symmetric 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 ( 2.2 -0.11 -0.32 ) ( -1.566 ) A = ( -0.11 2.93 0.81 ), B = ( -2.8425 ) ( -0.32 0.81 2.37 ) ( -1.1765 ) Sub Ex_Cg_r() Const N = 3, Nnz = 6, Tol = 0.0000000001 '1.0e-10 Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(N - 1) As Double, X(N - 1) As Double Dim XX(N - 1) As Double, YY(N - 1) As Double Dim Iter As Long, Res As Double, IRev As Long, Info As Long A(0) = 2.2: A(1) = -0.11: A(2) = 2.93: A(3) = -0.32: A(4) = 0.81: A(5) = 2.37 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) = -1.566: B(1) = -2.8425: B(2) = -1.1765 IRev = 0 Do Call Cg_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res) If IRev = 1 Then '- Matvec Call SsrDusmv("L", N, 1, A(), Ia(), Ja(), XX(), 0, YY()) ElseIf IRev = 3 Then '- Psolve Call Dcopy(N, YY(0), XX(0)) ElseIf IRev = 10 Then '- Check convergence If Res < Tol Then IRev = 11 End If Loop While IRev <> 0 Debug.Print "X =", X(0), X(1), X(2) Debug.Print "Iter = " + CStr(Iter) + ", Res = " + CStr(Res) + ", Info = " + CStr(Info) End Sub Example Results X = -0.8 -0.92 -0.29 Iter = 3, Res = 1.03670172979888E-16, Info = 0