<< spilu_ilud Spilu spilu_iluk >>

Spilu >> Spilu > spilu_iludp

spilu_iludp

ILUD with column pivoting

Calling Sequence

[L,U,perm]=spilu_iludp(A)
[L,U,perm]=spilu_iludp(A,alph)
[L,U,perm]=spilu_iludp(A,alph,drop)
[L,U,perm]=spilu_iludp(A,alph,drop,permtol)
[L,U,perm]=spilu_iludp(A,alph,drop,permtol,bloc)

Parameters

A

a n-by-n sparse real matrix of doubles

alph

a 1-by-1 matrix of doubles, positive, the diagonal compensation parameter (default: alph=0.5)

drop

a 1-by-1 matrix of doubles, positive, threshold for dropping elements in L and U (default: drop=0.001*max(abs(A)))

permtol

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

bloc

a 1-by-1 matrix of doubles, positive, integer value, 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 the 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
      

The matrix A must be square.

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

The term

      alph*(sum of all dropped out elements in a given row)
      

is added to the diagonal element of U of the factorization. Thus:

The parameter permtol is a tolerance ratio used for determning whether to permute two columns. Two columns are permuted only 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_iludp(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 alph
[L,U,perm]=spilu_iludp(A,0);
// Configure drop
[L,U,perm]=spilu_iludp(A,[],1);
// Configure permtol
[L,U,perm]=spilu_iludp(A,[],[],0);
// Configure bloc
[L,U,perm]=spilu_iludp(A,[],[],[],2);

// 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);
alph=0;
drop=0;
permtol=0;
bloc=3;
[L,U,iperm]=spilu_iludp(A,alph,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_ilud Spilu spilu_iluk >>