ILUD with column pivoting
[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)
a n-by-n sparse real matrix of doubles
a 1-by-1 matrix of doubles, positive, the diagonal compensation
parameter (default: alph=0.5
)
a 1-by-1 matrix of doubles, positive, threshold for dropping
elements in L
and U
(default: drop=0.001*max(abs(A))
)
a 1-by-1 matrix of doubles, positive, tolerance ratio used for
determining whether to permute two columns (default: permtol=0.5
)
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]
.
a n-by-n sparse real matrix of doubles, lower triangular
a n-by-n sparse real matrix of doubles, upper triangular
a 1-by-n dense real matrix of doubles, integer value, the permutation arrays
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:
alph
= 0 : ILU with threshold,
alph
= 1 : MILU with threshold.
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.
// 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) | ![]() | ![]() |
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