General description of incomplete factorizations

The IFPACK incomplete factorizations (IC, ICT, ILU, ILUT) compute a "Relaxed" LDL^T or ILU factorization of a given Epetra_RowMatrix. The factorization that is produced is a function of several parameters.

1. The pattern of the matrix (for Ifpack_IC and Ifpack_ILU) - All fill is derived from the original matrix nonzero structure. Level zero fill is defined as the original matrix pattern (nonzero structure), even if the matrix value at an entry is stored as a zero. (Thus it is possible to add entries to the ILU factors by adding zero entries the original matrix.)

2. Level of fill (for Ifpack_IC and Ifpack_ILU) - (integer value) Starting with the original matrix pattern as level fill of zero, the next level of fill is determined by analyzing the graph of the previous level and determining nonzero fill that is a result of combining entries that were from previous level only (not the current level). This rule limits fill to entries that are direct decendents from the previous level graph. Fill for level k is determined by applying this rule recursively. For sufficiently large values of k, the fill would eventually be complete and an exact LU factorization would be computed. object.

3. Level of fill (for Ifpack_ICT and Ifpack_ILUT) - (double value) During the factorization, a complete row of the factors is computed, then only the fill-in-value * nonzeros-in-row elements with largest magnitude are kept (fill-in-value is the specified level of fill, and nonzeros-in-row is the number of nonzeros in the original matrix nonzero structure).

4. Drop tolerance (for Ifpack_ICT and Ifpack_ILUT) - (double value) Entries with magnitude less than the drop tolerance will be dropped. Ifpack may perform this dropping on the original matrix A instead of the factors L and U. Thus, Ifpack_ILUT will give a slightly different preconditioner than the standard ILUT.

5. Fraction of relaxation - Ifpack_ILU computes the ILU factorization row-by-row. As entries at a given row are computed, some number of them will be dropped because they do match the prescribed sparsity pattern. The relaxation factor determines how these dropped values will be handled. If the RelaxValue is zero, then these extra entries will by dropped. This is a classical ILU approach. If the RelaxValue is 1, then the sum of the extra entries will be added to the diagonal. This is a classical Modified ILU (MILU) approach. If RelaxValue is between 0 and 1, then RelaxValue times the sum of extra entries will be added to the diagonal.

   For most situations, RelaxValue should be set to zero. For certain kinds of problems, e.g., reservoir modeling, there is a conservation principle involved such that any operator should obey a zero row-sum property. MILU was designed for these cases and you should set the RelaxValue to 1. For other situations, setting RelaxValue to some nonzero value may improve the stability of factorization, and can be used if the computed ILU factors are poorly conditioned.

6. Diagonal perturbation - Prior to computing the factorization, it is possible to modify the diagonal entries of the matrix for which the factorization will be computing. If the absolute and relative perturbation values are zero and one, respectively, the factorization will be compute for the original user matrix A. Otherwise, the factorization will computed for a matrix that differs from the original user matrix in the diagonal values only. More details can be found below.

   Given the above method to estimate the conditioning of the incomplete factors, if we detect that our factorization is too ill-conditioned we can improve the conditioning by perturbing the matrix diagonal and restarting the factorization using this more diagonally dominant matrix. In order to apply perturbation, prior to starting the factorization, we compute a diagonal perturbation of our matrix \(A\) and perform the factorization on this perturbed matrix. The overhead cost of perturbing the diagonal is minimal since the first step in computing the incomplete factors is to copy the matrix \(A\) into the memory space for the incomplete factors. We simply compute the perturbed diagonal at this point. The actual perturbation values we use are the diagonal values \((d_1, d_2, \ldots, d_n)\) with \(d_i = \rho sgn(d_i)\alpha + d_i\rho\), \(i=1, 2, \ldots, n\), where \(n\) is the matrix dimension and \(sgn(d_i)\) returns the sign of the diagonal entry. This has the effect of forcing the diagonal values to have minimal magnitude of \(\alpha\) and to increase each by an amount proportional to \(\rho\), and still keep the sign of the original diagonal entry.