<< Spilu Spilu spilu_ilu0 >>

Spilu >> Spilu > Overview

Overview

An overview of the Sparse Incomplete LU toolbox.

Purpose

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.

Quick start

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.

[Lfull,Ufull]=lu(full(A));
nnz(sparse(Lfull)) + nnz(sparse(Ufull)) // 6778

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.

spilu_ilutpplot(A);

Choosing a method

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.

Depending on the sparsity pattern of the L and U matrix, the algorithms can be sorted in two sets.

Algorithms based on the levels of fill can cause some difficulties for realistic problems that arise in many applications. On the other hand, with threshold strategies, the zero pattern is determined dynamically.

Depending on the pivoting strategy, the algorithms can be sorted in two sets.

Algorithms without pivoting may fail for many of the matrices that arise from real applications, for one of the following reasons.

In this case, the solution is to use pivoting.

History

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.

Authors

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

Acknowledgements

Michael Baudin thanks Thierry Clopeau for his advices on this topic.

Licence

This toolbox is distributed under the CeCILL.

Bibliography

"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).


Report an issue
<< Spilu Spilu spilu_ilu0 >>