SscIc0 SscIcSolve SscSsorSolve SsrIc0 SsrIcSolve SsrSsorSolve

SsrIc0() Sub SsrIc0 ( N As  Long, Val() As  Double, Rowptr() As  Long, Colind() As  Long, Val2() As  Double, Idiag() As  Long, Optional Info As  Long, Optional ByVal Uplo As  String = "L", Optional Base As  Long = -1  ) Initialize incomplete Cholesky decomposition without fill-in (IC0) preconditioner (symmetric positive definite matrix) (CSR) Purpose This routine computes the incomplete Cholesky decomposition, without fill-in, of the symmetric positive definite coefficient matrix of the sparse linear equations. A = L * D * L^T + R or A = U^T * D * U + R where L is the lower triangular matrix, U is the upper triangular matrix, and D is the diagonal matrix. It is called that fill-in occurs when a position that is zero in coefficient matrix A but it becomes non-zero after Cholesky decomposition is performed. To avoid this, such positions in L or U in the above equation are forcibly set to zeros. R shows the differences from the complete Cholesky decomposition arising from this modification. Assuming that R is not large, it is expected that this incomplete decomposition can be used to obtain a good approximate solution of linear equations. This decomposition is used as the preconditioner matrix M. M = L * D * L^T or M = U^T * D * U This routine outputs L or U and D in Val2(), and indices of diagonal elements in Idiag(). Val2() and Idiag() will be used by SsrIcSolve(). 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 number of nonzero elements of matrix A) Values of non-zero elements of matrix A. (Only the diagonal and lower triangular elements are referenced) [in] Rowptr() Array Rowptr(LRowptr - 1) (LRowptr >= N + 1) Row pointers of matrix A. [in] Colind() Array Colind(LColind - 1) (LColind >= Nnz) Column indices of matrix A. [out] Val2() Array Val2(LVal2 - 1) (LVal2 >= Nnz) Values of non-zero elements of preconditioner matrix M. (Values to be stored in the same location of diagonal and lower triangular elements of A) [out] Idiag() Array Idiag(LIdiag) (LIdiag >= N) Indices of diagonal elements. [out] Info (Optional) = 0: Successful exit. = i < 0: The (-i)-th argument is invalid. = j > 0: Matrix is singular (j-th diagonal element is zero). [in] Uplo (Optional) Specifies whether the upper or lower triangular part of input matrix A is stored. (default = "L") = "U": Upper triangular part (to be factorized to U^T*D*U). = "L": Lower triangular part (to be factorized to L*D*L^T). [in] Base (Optional) Indexing of Rowptr() and Colind(). = 0: Zero-based (C style) indexing: Starting index is 0. = 1: One-based (Fortran style) indexing: Starting index is 1. (default: Assumes 1 if Rowptr(0) = 1, 0 otherwise) Example Program Solve the system of linear equations Ax = B by CG method with IC(0) preconditioner, 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_Ic0_Ssr() 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 Dim M(Nnz - 1) As Double, Id(N - 1) As Long Call SsrIc0(N, A(), Ia(), Ja(), M(), Id(), Info) If Info <> 0 Then Debug.Print "Ic0 Info =" + Str(Info) 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 SsrIcSolve(N, M(), Ia(), Ja(), Id(), YY(), XX(), Info) If Info <> 0 Then Debug.Print "Ic0Solve Info =" + Str(Info) 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 = 1, Res = 2.22044604925031E-16, Info = 0