Navigation

Operators and Keywords

Function List:

C++ API

: ilu (A)
: ilu (A, opts)
: [L, U] = ilu (…)
: [L, U, P] = ilu (…)

Compute the incomplete LU factorization of the sparse square matrix A.

ilu returns a unit lower triangular matrix L, an upper triangular matrix U, and optionally a permutation matrix P, such that L*U approximates P*A.

The factors given by this routine may be useful as preconditioners for a system of linear equations being solved by iterative methods such as BICG (BiConjugate Gradients) or GMRES (Generalized Minimum Residual Method).

The factorization may be modified by passing options in a structure opts. The option name is a field of the structure and the setting is the value of field. Names and specifiers are case sensitive.

type

Type of factorization.

"nofill"

ILU factorization with no fill-in (ILU(0)).

Additional supported options: milu.

"crout"

Crout version of ILU factorization (ILUC).

Additional supported options: milu, droptol.

"ilutp" (default)

ILU factorization with threshold and pivoting.

Additional supported options: milu, droptol, udiag, thresh.

droptol

A non-negative scalar specifying the drop tolerance for factorization. The default value is 0 which produces the complete LU factorization.

Non-diagonal entries of U are set to 0 unless

abs (U(i,j)) >= droptol * norm (A(:,j)).

Non-diagonal entries of L are set to 0 unless

abs (L(i,j)) >= droptol * norm (A(:,j))/U(j,j).

milu

Modified incomplete LU factorization:

"row"

Row-sum modified incomplete LU factorization. The factorization preserves row sums: A * e = L * U * e, where e is a vector of ones.

"col"

Column-sum modified incomplete LU factorization. The factorization preserves column sums: e' * A = e' * L * U.

"off" (default)

Row and column sums are not necessarily preserved.

udiag

If true, any zeros on the diagonal of the upper triangular factor are replaced by the local drop tolerance droptol * norm (A(:,j))/U(j,j). The default is false.

thresh

Pivot threshold for factorization. It can range between 0 (diagonal pivoting) and 1 (default), where the maximum magnitude entry in the column is chosen to be the pivot.

If ilu is called with just one output, the returned matrix is L + U - speye (size (A)), where L is unit lower triangular and U is upper triangular.

With two outputs, ilu returns a unit lower triangular matrix L and an upper triangular matrix U. For opts.type == "ilutp", one of the factors is permuted based on the value of opts.milu. When opts.milu == "row", U is a column permuted upper triangular factor. Otherwise, L is a row-permuted unit lower triangular factor.

If there are three named outputs and opts.milu != "row", P is returned such that L and U are incomplete factors of P*A. When opts.milu == "row", P is returned such that L and U are incomplete factors of A*P.

EXAMPLES 

A = gallery ("neumann", 1600) + speye (1600);
opts.type = "nofill";
nnz (A)
ans = 7840

nnz (lu (A))
ans = 126478

nnz (ilu (A, opts))
ans = 7840

This shows that A has 7,840 nonzeros, the complete LU factorization has 126,478 nonzeros, and the incomplete LU factorization, with 0 level of fill-in, has 7,840 nonzeros, the same amount as A. Taken from: http://www.mathworks.com/help/matlab/ref/ilu.html 

A = gallery ("wathen", 10, 10);
b = sum (A, 2);
tol = 1e-8;
maxit = 50;
opts.type = "crout";
opts.droptol = 1e-4;
[L, U] = ilu (A, opts);
x = bicg (A, b, tol, maxit, L, U);
norm (A * x - b, inf)

This example uses ILU as preconditioner for a random FEM-Matrix, which has a large condition number. Without L and U BICG would not converge.

See also: lu, ichol, bicg, gmres.

Package: octave