<< spilu_ilut Spilu spilu_milu0 >>

Spilu >> Spilu > spilu_ilutp

spilu_ilutp

sparse ILU factorization with truncation and pivoting

Calling Sequence

[L,U,perm]=spilu_ilutp(A)
[L,U,perm]=spilu_ilutp(A,lfil)
[L,U,perm]=spilu_ilutp(A,lfil,drop)
[L,U,perm]=spilu_ilutp(A,lfil,drop,permtol)
[L,U,perm]=spilu_ilutp(A,lfil,drop,permtol,bloc)

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

permtol

a 1-by-1 matrix of doubles, positive, tolerance ratio used to determine whether or not to permute two columns (default: permtol=0.5)

bloc

a 1-by-1 matrix of doubles, positive, permuting can be done within the diagonal blocks of size bloc (default: bloc=n). bloc 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

perm

a 1-by-n dense real matrix of doubles, integer value, the permutation arrays

Description

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

        P*A ≈ L*U
      

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 .ge. 0.

The permtol parameter is the tolerance ratio used to determine whether or not to permute two columns. At step i, columns i and j are permuted when

        abs(a(i,j))*permtol > abs(a(i,i)).
      

If permtol==0, then never permute. Recommended values of permtol are from 0.1 to 0.01.

If desired, permuting can be done only within the diagonal blocks of size bloc. The bloc parameter is useful for PDE problems with several degrees of freedom. To discard the feature, set bloc=n.

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,perm]=spilu_ilutp(A)
// Check the quality of the decomposition
P = spilu_permVecToMat(perm);
[norm(P*A-L*U,"inf")/norm(A,"inf") nnz(L)+nnz(U)]
// Configure lfil
[L,U,perm]=spilu_ilutp(A,1);
// Configure drop
[L,U,perm]=spilu_ilutp(A,[],1);
// Configure permtol
[L,U,perm]=spilu_ilutp(A,[],[],0);
// Configure bloc
[L,U,perm]=spilu_ilutp(A,[],[],[],2);

// Decompose 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;
drop=0;
permtol=0;
bloc=3;
[L,U,perm]=spilu_ilutp(A,lfil,drop,permtol,bloc);
P = spilu_permVecToMat(perm);
x=U\(L\(P*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


Report an issue
<< spilu_ilut Spilu spilu_milu0 >>