<< spilu_iludp Spilu spilu_ilut >>

Spilu >> Spilu > spilu_iluk

spilu_iluk

level-k of the sparse ILU factorization with dual truncation

Calling Sequence

[L,U]=spilu_iluk(A)
[L,U]=spilu_iluk(A,lfil)

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). lfil must be in the range [0,n].

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
      

All the diagonal elements of the input matrix must be nonzero.

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.

It may happen that the memory internally allocated by the toolbox for the LU factors is too small. This is because the number of entries in the LU factors is not predictable. Hence, the toolbox makes some assumption about the size of the produced LU. If this assumption is wrong, the "not enough memory for matrix L" error is produced (see below for an example). This is not a bug: this is a consequence of the algorithm. In this case, we may use another type of decomposition.

Examples

// 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_iluk(A)
// Check the quality of the decomposition
[norm(A-L*U,"inf")/norm(A,"inf") nnz(L)+nnz(U)]
//
// Here, the default value of elts is :
// 25 = 5+6+2+2+5+5
// Indeed, for the row #1, the bandwidth is 5-1+1=5 and 
// for the row #3, the bandwidth is 4-3+1=2.
//
// Here, the default value of lfil is :
n = size(A,"r") // n=6
lfil=floor(1+nnz(A)/n) // lfil = 3
// Configure lfil
[L,U]=spilu_iluk(A,1)

// A 225-by-225 sparse matrix
// nnz(A)=1065
path = spilu_getpath (  );
filename = fullfile(path,"tests","matrices","pde225.mtx");
A=mmread(filename);
n=size(A,1);
b=ones(n,1);
lfil=10;
[L,U]=spilu_iluk(A,lfil);
x=U\(L\b);
norm(A*x-b)

// In the following case, there is 
// not enough memory for matrix L.
// An error is generated.
A = [
    56.  -12.    3.     5.   -1.     7.   
  -12.    24.    0.     0.     0.     0.   
    3.     0.     24.    0.     9.     0.   
    5.     0.     0.     10.    0.     0.   
  -1.     0.     9.     0.     20.    0.   
    7.     0.     0.     0.     0.     14.  
];
A = sparse(A);
[L,U]=spilu_iluk(A);

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


Report an issue
<< spilu_iludp Spilu spilu_ilut >>