level-k of the sparse ILU factorization with dual truncation
[L,U]=spilu_iluk(A) [L,U]=spilu_iluk(A,lfil)
a n-by-n sparse real matrix of doubles
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]
.
a n-by-n sparse real matrix of doubles, lower triangular
a n-by-n sparse real matrix of doubles, upper triangular
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.
// 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); | ![]() | ![]() |
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