Accelerating the pace of engineering and science

• Trials

# luinc

Sparse incomplete LU factorization

 Note:   luinc has been removed. Use ilu instead.

## Syntax

luinc(A,'0')
luinc(A,droptol)
luinc(A,options)
[L,U] = luinc(A,'0')
[L,U] = luinc(A,options)
[L,U,P] = luinc(...)

## Description

luinc produces a unit lower triangular matrix, an upper triangular matrix, and a permutation matrix.

luinc(A,'0') computes the incomplete LU factorization of level 0 of a square sparse matrix. The triangular factors have the same sparsity pattern as the permutation of the original sparse matrix A, and their product agrees with the permuted A over its sparsity pattern. luinc(A,'0') returns the strict lower triangular part of the factor and the upper triangular factor embedded within the same matrix. The permutation information is lost, but nnz(luinc(A,'0')) = nnz(A), with the possible exception of some zeros due to cancellation.

luinc(A,droptol) computes the incomplete LU factorization of any sparse matrix using the drop tolerance specified by the non-negative scalar droptol. The result is an approximation of the complete LU factors returned by lu(A). For increasingly smaller values of the drop tolerance, this approximation improves until the drop tolerance is 0, at which time the complete LU factorization is produced, as in lu(A).

As each column j of the triangular incomplete factors is being computed, the entries smaller in magnitude than the local drop tolerance (the product of the drop tolerance and the norm of the corresponding column of A)

`  droptol*norm(A(:,j))`

are dropped from the appropriate factor.

The only exceptions to this dropping rule are the diagonal entries of the upper triangular factor, which are preserved to avoid a singular factor.

luinc(A,options) computes the factorization with up to four options. These options are specified by fields of the input structure options. The fields must be named exactly as shown in the table below. You can include any number of these fields in the structure and define them in any order. Any additional fields are ignored.

Field Name

Description

droptol

Drop tolerance of the incomplete factorization.

milu

If milu is 1, luinc produces the modified incomplete LU factorization that subtracts the dropped elements in any column from the diagonal element of the upper triangular factor. The default value is 0.

udiag

If udiag is 1, any zeros on the diagonal of the upper triangular factor are replaced by the local drop tolerance. The default is 0.

thresh

Pivot threshold between 0 (forces diagonal pivoting) and 1, the default, which always chooses the maximum magnitude entry in the column to be the pivot. thresh is described in greater detail in the lu reference page.

luinc(A,options) is the same as luinc(A,droptol) if options has droptol as its only field.

[L,U] = luinc(A,'0') returns the product of permutation matrices and a unit lower triangular matrix in L and an upper triangular matrix in U. The exact sparsity patterns of L, U, and A are not comparable but the number of nonzeros is maintained with the possible exception of some zeros in L and U due to cancellation:

`  nnz(L)+nnz(U) = nnz(A)+n, where A is n-by-n.`

The product L*U agrees with A over its sparsity pattern. (L*U).*spones(A)-A has entries of the order of eps.

[L,U] = luinc(A,options) returns a permutation of a unit lower triangular matrix in L and an upper triangular matrix in U. The product L*U is an approximation to A. luinc(A,options) returns the strict lower triangular part of the factor and the upper triangular factor embedded within the same matrix. The permutation information is lost.

[L,U,P] = luinc(...) returns a unit lower triangular matrix in L, an upper triangular matrix in U, and a permutation matrix in P.

[L,U,P] = luinc(A,'0') returns a unit lower triangular matrix in L, an upper triangular matrix in U and a permutation matrix in P. L has the same sparsity pattern as the lower triangle of permuted A

`  spones(L) = spones(tril(P*A))`

with the possible exceptions of 1s on the diagonal of L where P*A may be zero, and zeros in L due to cancellation where P*A may be nonzero. U has the same sparsity pattern as the upper triangle of P*A

`  spones(U) = spones(triu(P*A))`

with the possible exceptions of zeros in U due to cancellation where P*A may be nonzero. The product L*U agrees within rounding error with the permuted matrix P*A over its sparsity pattern. (L*U).*spones(P*A)-P*A has entries of the order of eps.

[L,U,P] = luinc(A,options) returns a unit lower triangular matrix in L, an upper triangular matrix in U, and a permutation matrix in P. The nonzero entries of U satisfy

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

with the possible exception of the diagonal entries, which were retained despite not satisfying the criterion. The entries of L were tested against the local drop tolerance before being scaled by the pivot, so for nonzeros in L

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

The product L*U is an approximation to the permuted P*A.

## Limitations

luinc(X,'0') works on square matrices only.

## Examples

```load west0479;
S = west0479;
[L,U] = lu(S);```

Compute the incomplete LU factorization of level 0.

```[L,U,P] = luinc(S,'0');
D = (L*U).*spones(P*S)-P*S;```

spones(U) and spones(triu(P*S)) are identical.

spones(L) and spones(tril(P*S)) disagree at 73 places on the diagonal, where L is 1 and P*S is 0, and also at position (206,113), where L is 0 due to cancellation, and P*S is -1. D has entries of the order of eps.

```[IL0,IU0,IP0] = luinc(S,0);
[IL1,IU1,IP1] = luinc(S,1e-10);
.
.
.```

A drop tolerance of 0 produces the complete LU factorization. Increasing the drop tolerance increases the sparsity of the factors (decreases the number of nonzeros) but also increases the error in the factors, as seen in the plot of drop tolerance versus norm(L*U-P*S,1)/norm(S,1) in the second figure below.

expand all

### Tips

These incomplete factorizations may be useful as preconditioners for solving large sparse systems of linear equations. The lower triangular factors all have 1s along the main diagonal but a single 0 on the diagonal of the upper triangular factor makes it singular. The incomplete factorization with a drop tolerance prints a warning message if the upper triangular factor has zeros on the diagonal. Similarly, using the udiag option to replace a zero diagonal only gets rid of the symptoms of the problem but does not solve it. The preconditioner may not be singular, but it probably is not useful and a warning message is printed.

### Algorithms

luinc(A,'0') is based on the "KJI" variant of the LU factorization with partial pivoting. Updates are made only to positions which are nonzero in A.

luinc(A,droptol) and luinc(A,options) are based on the column-oriented lu for sparse matrices.

## References

[1] Saad, Yousef, Iterative Methods for Sparse Linear Systems, PWS Publishing Company, 1996, Chapter 10 - Preconditioning Techniques.