<< spilu_iluk Spilu spilu_ilutp >>

Spilu >> Spilu > spilu_ilut

spilu_ilut

sparse ILU factorization with dual truncation

Calling Sequence

[L,U]=spilu_ilut(A)
[L,U]=spilu_ilut(A,lfil)
[L,U]=spilu_ilut(A,lfil,drop)

Parameters

A

a n-by-n sparse real matrix of doubles

lfil

a 1-by-1 matrix of doubles, positive, integer value, the fill-in parameter (default: lfil=floor(1+nnz(A)/n), i.e. the average number of nonzero elements of the matrix A by line). Each row of L and each row of U will have a maximum of lfil elements (excluding the diagonal element). lfil must be in the range [0,n].

drop

a 1-by-1 matrix of doubles, positive, the threshold for dropping small terms in the factorization (default: drop=0.001*max(abs(A)))

L

a n-by-n sparse real matrix of doubles, lower triangular

U

a n-by-n sparse real matrix of doubles, upper triangular

Description

Builds an incomplete LU factorization of the sparse matrix A, that is, computes a lower triangular matrix L and an upper triangular matrix U such that

        A ≈ L*U
      

The diagonal elements of the input matrix must be nonzero (at least 'structurally').

Any optional input argument equal to the empty matrix [] is replaced by its default value.

Each row of L and each row of U will have a maximum of lfil elements (excluding the diagonal element). lfil must be greater or equal than 0.

The dual drop strategy works as follows.

Flexibility: one can use drop=0 to get a strategy based on keeping the largest elements in each row of L and U. Taking drop different from 0 but lfil=n, will give threshold strategy (however, fill-in is then inpredictible).

Examples

// Decompose a 6-by-6 sparse matrix
// nnz(A) = 16
A=[
-1.    3.    0.    0.    4.    0.  
2.   -5.    0.    0.    0.    1.  
0.    0.   -2.    3.    0.    0.  
0.    0.    7.   -1.    0.    0.  
-3.    0.    0.    4.    6.    0.
0.    5.    0.    0.   -7.    8. 
];
A=sparse(A);
[L,U]=spilu_ilut(A)
// Check the quality of the decomposition
[norm(A-L*U,"inf")/norm(A,"inf") nnz(L)+nnz(U)]
// Configure lfil
[L,U]=spilu_ilut(A,1)
// Configure drop
[L,U]=spilu_ilut(A,[],0)
[L,U]=spilu_ilut(A,[],1)

// Try various values of lfil
mprintf("lfil \t norm(A-L*U) \t nnz(L)+nnz(U)\n");
for lfil = 1:6
  [L,U]=spilu_ilut(A,lfil);
  mprintf("%4d \t %11d \t %13d\n",lfil,..
    norm(A-L*U,"inf")/norm(A,"inf"),nnz(L)+nnz(U));
end

// Decompose a 237-by-237 sparse matrix
// nnz(A)=1017
path = spilu_getpath (  );
filename = fullfile(path,"tests","matrices","nos1.mtx");
A=mmread(filename);
n=size(A,1);
b=ones(n,1);
lfil=10;
drop=0;
[L,U]=spilu_ilut(A,lfil,drop);
x=U\(L\b);
norm(A*x-b)

Authors

Copyright (C) 2011 - DIGITEO - Michael Baudin

Copyright (C) 2011 - NII - Benoit Goepfert

Copyright (C) 2005 - INRIA - Sage Group (IRISA)

Copyright (C) 1993 - Univ. of Tennessee and Oak Ridge National Laboratory - Youcef Saad


<< spilu_iluk Spilu spilu_ilutp >>