Cg Cg1 Cg_r Cg_s Cr Cr_r Cr_s

Cg1() Sub Cg1 ( N As  Long, Val() As  Double, Rowptr() As  Long, Colind() As  Long, B() As  Double, X() As  Double, Optional Info As  Long, Optional Iter As  Long, Optional Res As  Double, Optional MaxIter As  Long = 500, Optional Tol As  Double = 1.0E-10, Optional Uplo As  String = "F"  ) Solution of linear system Ax = b using conjugate gradient (CG) method (Symmetric positive definite matrix) (Simple driver) Purpose This routine solves the linear system Ax = b with symmetric positive definite coefficient matrix using the conjugate gradient (CG) method. Matrix A is represented in CSR format. Parameters [in] N Dimension of the matrix A. (N >= 0) (if N = 0, returns without computation) [in] Val() Array Val(LVal - 1) (LVal >= Nnz) (Nnz is the number of nonzero elements of the matrix A) Values of nonzero elements of the matrix A. [in] Rowptr() Array Rowptr(LRowptr - 1) (LRowptr >= N + 1) Row pointers of the matrix A. [in] Colind() Array Colind(LColind - 1) (LColind >= Nnz) Column indices of the matrix A. [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. = 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). [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) [in] Tol (Optional) Tolerance for convergence test. (default = 1.0e-10) Assumed to be converged if norm(b - A*x) <= Tol*norm(b). If Tol < eps (machine epsilon), Tol = eps is assumed. [in] Uplo (Optional) Specifies whether the upper or lower triangular part of matrix A is stored. (default = "F") = "U": Upper triangular part is stored. = "L": Lower triangular part is stored. = "F": Full matrix (both upper and lower triangular parts) are stored. Note Indexing of Rowptr() and Ind() is assumed to be zero-based (C style) if Rowptr(0) = 1, or one-based (Fortran style) otherwise. 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_Cg1() Const N = 3, Nnz = 6 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 Iter As Long, Res As Double, 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 Call Cg1(N, A(), Ia(), Ja(), B(), X(), Info, Iter, Res) 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