ZCscIlu ZCscIlu0 ZCscIluSolve ZCscSsorSolve ZCsrIlu ZCsrIlu0 ZCsrIluSolve ZCsrSsorSolve ZCsxDs ZCsxDsSolve ZCsxSsor ZSscSsorSolve ZSsrSsorSolve

ZCsrIlu() Sub ZCsrIlu ( N As  Long, Val() As  Complex, Rowptr() As  Long, Colind() As  Long, P As  Long, Val2() As  Complex, Rowptr2() As  Long, Colind2() As  Long, D() As  Complex, Optional Info As  Long, Optional Md As  Long = 0, Optional Base As  Long = -1, Optional Base2 As  Long = 0  ) Initialize incomplete LU decomposition with level (ILU(p)) preconditioner (Complex matrices) (CSR) Purpose This routine computes the incomplete LU decomposition of the coefficient matrix A of the sparse linear equations. A = L * U + R where R is the difference from the complete LU decomposition. Assuming that R is small, the following preconditioner matrix is obtained by solving the equations using this decomposition. M = L * U In order to suppress fill-ins, new non-zero elements generated during decomposition with level values larger than the specified value p are dropped. This is called as ILU(p). Level values are initialized as follows. lev(a(i,j)) = 0 (if a(i,j) != 0 or i = j), = ∞ (if a(i,j) = 0) In the elimination process of LU decomposition, the operation a(i,j) = a(i,j) - a(i,k)*a(k,j) is required. At the time, the level value of a(i,j) is updated as follows. lev(a(i,j)) = min{ lev(a(i,j)), lev(a(i,k)) + lev(a(k,j)) + 1 } Due to this algorithm, the level of the position of the non-zero element of A before decomposition remains 0 until the last. When decomposition is completed, elements with level values greater than p are dropped. This routine outputs the lower triangular matrix L and the upper triangular matrix U in Val2(), Rowptr2() and Colind2(). The diagonal elements of U are copied to D(). Val2(), Rowptr2(), Colind2() and D() will be used by CsrIluSolve(). 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 non-zero elements) Values of non-zero elements of matrix A. [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. [in] P The value of level p of ILU(p) algorithm. (P >= 0) Note: it is possible to set P = 0. However, it is recommended to use CsrIlu0() since it is faster. [out] Val2() Array Val2(LVal2 - 1) (LVal2 > 0) Values of non-zero elements of the lower triangular matrix L and the upper triangular matrix U (matrix L*U). Nnz2 = min(LVal2, LColind2) is the maximum number of non-zero elements of matrix L*U. If the number of non-zero elements exceeds Nnz2 during the factorization, the routine will terminate with error number -6. [out] Rowptr2() Array Rowptr2(LRowptr2 - 1) (LRowptr2 >= N + 1) Row pointers of matrix L*U. [out] Colind2() Array Colind2(LColind2 - 1) (LColind2 > 0) Column indices of matrix L*U. [out] D() Array D(LD) (LD >= N) The diagonal elements of upper triangular matrix U. [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] Md (Optional) Specify whether the diagonal compensation of the dropped elements are required. ILU with this compensation is called as MILU (modified ILU). (default = 0) = 0: ILU (no compensation). = 1: MILU (with diagonal compensation). [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) [in] Base2 (Optional) Indexing of Rowptr2() and Colind2(). (default = 0) = 0: Zero-based (C style) indexing: Starting index is 0. = 1: One-based (Fortran style) indexing: Starting index is 1. Example Program Solve the system of linear equations Ax = B by FGMRES method with ILU(p) preconditioner, 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_ZFgmres_Ilu_Csr() Const N = 3, Nnz = N * N, NnzM = Nnz, 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) Dim M(NnzM - 1) As Complex, Im(N) As Long, Jm(NnzM - 1) As Long Dim P As Long, D(N - 1) As Complex P = 3 Call ZCsrIlu(N, A(), Ia(), Ja(), P, M(), Im(), Jm(), D(), Info) If Info <> 0 Then Debug.Print "Ilu Info =" + Str(Info) IRev = 0 Do Call ZFgmres_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 = 3 Then '- Psolve Call ZCsrIluSolve("N", N, M(), Im(), Jm(), D(), YY(), XX(), Info) If Info <> 0 Then Debug.Print "IluSolve Info =" + Str(Info) 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.59,-0.28) (-0.2,-3.99999999999999E-02) (0.24,-0.49) Iter = 1, Res = 1.83807799103727E-16, Info = 0