An overview of the Sparse Incomplete LU toolbox.
Spilu is a Scilab toolbox which provides preconditioners based on Incomplete LU (ILU) factorizations. This module is based on a set of Fortran routines from the Sparskit module by Yousef Saaf. More specifically, this module provides some of the preconditioners from the ITSOL sub-module of Sparskit.
The preconditioners which are provided in this toolbox may be used in preconditioned iterative algorithms for solving sparse linear systems of equations. According to Y. Saad, "roughly speaking, a preconditioner is any form of implicit or explicit modification of an original linear system which makes it easier to solve by a given iterative method." Examples of preconditioned iterative algorithms are the Generalized Minimum Residual Method (GMRES) or the Preconditioned Conjugate Gradient (PCG). Hence, the Spilu toolbox is the companion of the Imsl toolbox, which provides these iterative methods.
The need for incomplete factorizations is a consequence of a property of the general LU decomposition of a sparse matrix A. In general, the L and U matrices from the exact (i.e. complete) LU decomposition have much more nonzeros than A.
For example (see below for details), consider the 225-by-225 sparse matrix in the pde225 file provided with this toolbox. The matrix A has 1065 nonzeros, while the exact LU decomposition produces a L matrix with 3389 nonzeros. The matrix U also has 3389 nonzeros.
The fact that the exact LU decomposition increases the number of nonzeros is known as the "fill-in". This fill-in should be avoided for at least two reasons. First, the fill-in increases the memory required to store the decomposition. Second, and more importantly, the exact decomposition of A is not required to provide a good preconditionner to the iterative solvers such as GMRES or PCG. Indeed, even an incomplete decomposition can be used to provide good preconditioners to iteratives solvers for sparse linear systems of equations.
The spilu_ilutp
function computes a lower triangular matrix L
,
an upper triangular matrix U
and a permutation matrix P
such that
P*A ≈ L*U
In the following example, we decompose a 225-by-225 sparse matrix
with 1065 nonzeros.
First, we read the matrix data file with the Matrix Market function mmread
.
Then, we plot the sparsity pattern of the matrix A
with the
PlotSparse
function.
path = spilu_getpath ( ); filename = fullfile(path,"tests","matrices","pde225.mtx"); A=mmread(filename); nnz(A) // 1065 scf(); PlotSparse(A); xtitle("Sparsity pattern of A"); | ![]() | ![]() |
The previous script produces the following output.
-->nnz(A) ans = 1065.
Then we decompose the sparse matrix A
and plot the
sparsity pattern of the lower triangular matrix L
and the
upper triangular matrix U
.
[L,U,perm]=spilu_ilutp(A); P = spilu_permVecToMat(perm); nnz(L) + nnz(U) // 2316 norm(P*A-L*U,"inf")/norm(A,"inf") // 0.0395390 scf(); PlotSparse(L); xtitle("Sparsity pattern of L"); scf(); PlotSparse(U); xtitle("Sparsity pattern of U"); | ![]() | ![]() |
We notice that the number of nonzeros in L
and U
(2316) is larger than the number
of nonzeros in A
(1065).
But if we consider the full LU decomposition of A
, the number of
nonzeros is much larger, as shown in the following session.
In the following script, we plot the sparsity pattern of the matrix Lfull
.
scf(); PlotSparse(sparse(Lfull)); xtitle("Sparsity pattern of Lfull"); | ![]() | ![]() |
In the following script, we create a plot of the sensitivity of the ILUTP decomposition with respect to its parameters.
In this section, we give some general descriptions and advices on the choice of a decomposition algorithm.
This section is mainly based on the comments made by Y. Saad in "Iterative Methods for Sparse Linear Systems, Second Edition".
If no other information is provided, the spilu_ilutp
function
may be used.
This is because it is based on a threshold strategy, which take into account
for the magnitude of the elements.
Moreover, the pivoting in the ILUTP algorithm may prevent from encountering a
zero pivot or may limit the growth of the entries in the factors.
In the following list, we review all the algorithms provided in this toolbox and analyze their main properties.
spilu_ilu0
- ILU(0) preconditioning.
The Incomplete LU factorization technique with no fill-in, denoted by ILU(0), takes the zero pattern P to be precisely the zero pattern of A. The accuracy of the ILU(0) incomplete factorization may be insufficient to yield an adequate rate of convergence of a preconditioned iterative method.
spilu_ilud
- sparse ILU factorization with single dropping and diagonal compensation.
This algorithm computes the ILU factorization with a standard threshold dropping.
This function is controled by two parameters: drop
and alph
.
The drop
parameter is the threshold parameter for dropping small terms in the factorization.
Hence, a smaller value of drop
(e.g. drop=0.
) leads to a better
decomposition of the matrix, but also leads to a larger number of nonzeros
in the matrices L
and U
.
The parameter alph
is the diagonal compensation parameter, which modifies the
diagonal entries of U.
spilu_iludp
- ILUD with column pivoting.
This algorithm has four parameters: drop
and alph
(the same parameters as in ILUD) and permtol
and bloc
,
which control the pivoting strategy.
The permtol
parameter is a tolerance ratio used for determining whether to permute two columns.
The bloc
parameter controls the size of the diagonal blocks within which the
pivoting is done.
spilu_iluk
- level-k of the sparse ILU factorization with dual truncation.
This algorithm has one main parameter: lfil
is the
level of fill.
Each row of L and each row of U will have a maximum of lfil
elements
(excluding the diagonal element).
Hence, ILUK differ from ILU(0) by allowing some fill-in.
spilu_ilut
- sparse ILU factorization with dual truncation.
This function includes a set of rules for dropping small elements.
The ILUT algorithm has two parameters: lfil
and
drop
.
spilu_ilutp
- sparse ILU factorization with truncation and pivoting.
This algorithm has four parameters: lfil
and drop
(the same parameters as in ILUT) and permtol
and bloc
,
which control the pivoting strategy.
spilu_milu0
- MILU(0) preconditioning.
The algorithm is the same than in the ILU(0) algorithm, but with diagonal compensation in U at the end.
Depending on the sparsity pattern of the L and U matrix, the algorithms can be sorted in two sets.
based on the same sparsity pattern as A:
spilu_ilu0
, spilu_milu0
,
based on the levels of fill: spilu_iluk
,
based on a threshold strategy: spilu_ilud
, spilu_iludp
,
spilu_ilut
, spilu_ilutp
.
Depending on the pivoting strategy, the algorithms can be sorted in two sets.
without pivoting:
spilu_ilu0
, spilu_iluk
,
spilu_ilud
, spilu_ilut
,
spilu_milu0
,
with pivoting: spilu_iludp
,
spilu_ilutp
.
Algorithms without pivoting may fail for many of the matrices that arise from real applications, for one of the following reasons.
The procedure encounters a zero pivot.
The procedure encounters an overflow or underflow condition, because of an exponential growth of the entries of the factors.
The Sparskit toolkit was developped in 1993-1996 by Y.Saad, University of Minnesota.
In 2005, this module was developed by Aladin Group (IRISA-INRIA) in order to provide a Scilab connection to the Sparskit ILU decomposition functions.
In 2010 - 2011, Michael Baudin (Digitéo) and Benoit Goepfert (NII internship) updated this module for Scilab 5. Additionnally, bugs were fixed, the gateways were updated, help pages were updated and unit test were created. Finally, new functions were created to help the use of the toolbox.
Copyright (C) 1993 - Univ. of Tennessee and Oak Ridge National Laboratory
Copyright (C) 2000 - 2001 - INRIA - Aladin Group
Copyright (C) 2010 - 2011 - DIGITEO - Michael Baudin
Copyright (C) 2011 - NII - Benoit Goepfert
Michael Baudin thanks Thierry Clopeau for his advices on this topic.
This toolbox is distributed under the CeCILL.
"SPARSKIT: A basic tool-kit for sparse matrix computations", Youcef Saad, 1994
"Iterative Methods for Sparse Linear Systems, Second Edition", Saad, SIAM Publications, 2003 (ftp ftp.cs.umn.edu; cd dept/users/saad/PS; get all_ps.zip).