A=rand(4);
A=A*A';
A
A =
2.3967 1.2645 1.6027 1.5817
1.2645 1.0480 1.0544 1.1012
1.6027 1.0544 1.3409 1.1634
1.5817 1.1012 1.1634 1.8520
R=chol(A);
R'*R
ans =
2.3967 1.2645 1.6027 1.5817
1.2645 1.0480 1.0544 1.1012
1.6027 1.0544 1.3409 1.1634
1.5817 1.1012 1.1634 1.8520
R*R'
ans =
5.1795 1.2959 0.3017 0.7985
1.2959 0.6821 0.0885 0.3378
0.3017 0.0885 0.1653 -0.0807
0.7985 0.3378 -0.0807 0.6108
H=ichol(A);
{Error using ichol
First argument must be sparse.
}
H=ichol(sparse(A));
H'*H
ans =
(1,1) 5.1795
(2,1) 1.2959
(3,1) 0.3017
(4,1) 0.7985
(1,2) 1.2959
(2,2) 0.6821
(3,2) 0.0885
(4,2) 0.3378
(1,3) 0.3017
(2,3) 0.0885
(3,3) 0.1653
(4,3) -0.0807
(1,4) 0.7985
(2,4) 0.3378
(3,4) -0.0807
(4,4) 0.6108
full(H'*H)
ans =
5.1795 1.2959 0.3017 0.7985
1.2959 0.6821 0.0885 0.3378
0.3017 0.0885 0.1653 -0.0807
0.7985 0.3378 -0.0807 0.6108
A
A =
2.3967 1.2645 1.6027 1.5817
1.2645 1.0480 1.0544 1.1012
1.6027 1.0544 1.3409 1.1634
1.5817 1.1012 1.1634 1.8520
full(H*H')
ans =
2.3967 1.2645 1.6027 1.5817
1.2645 1.0480 1.0544 1.1012
1.6027 1.0544 1.3409 1.1634
1.5817 1.1012 1.1634 1.8520
fem2d
Preconditioner? (0 = No, 1 = Yes) 1
Reordering of the points? (0 = No, 1 = Yes) 0
norm of HH'-A
ans =
2.0399
pcg converged at iteration 32 to a solution with relative residual 5.8e-05.
ans =
8.5971e-05
normest(A)
ans =
4.5036e+15
fem2d
Preconditioner? (0 = No, 1 = Yes) 1
Reordering of the points? (0 = No, 1 = Yes) 1
norm of HH'-A
ans =
0.7111
pcg converged at iteration 19 to a solution with relative residual 3.9e-05.
ans =
8.5943e-05
help ilu
ilu Sparse Incomplete LU factorization
The factors given by this factorization may be useful as
preconditioners for a system of linear equations being solved by
iterative methods such as BICG (BiConjugate Gradients) and GMRES
(Generalized Minimum Residual Method).
ilu(A) computes the incomplete LU factorization of A with zero level of
fill in.
ilu(A,SETUP) performs the incomplete LU factorization of A. SETUP is
a structure with up to five fields:
type --- type of factorization
droptol --- the drop tolerance of incomplete LU
milu --- modified incomplete LU
udiag --- replace zeros on the diagonal of U
thresh --- the pivot threshold
type may be 'nofill' which is the ilu factorization with zero level of
fill in, known as ilu(0), 'crout' which is the Crout Version of ilu,
known as ILUC, or 'ilutp' which is the ilu factorization with
threshold and pivoting. If type is not specified the ilu factorization
with 0 level of fill in will be performed. Pivoting is never
performed with type 'nofill' and with type 'crout'.
droptol is a non-negative scalar used as the drop tolerance which
means that all entries which are smaller in magnitude than the local
drop tolerance, which is droptol * NORM of the column of A for the
column and droptol * NORM of the row of A for the row, are "dropped"
from L or U. The only exception to this dropping rule is the diagonal
of the upper triangular factor U which is never dropped. Note that
entries of the lower triangular factor L are tested before being
scaled by the pivot. Setting droptol = 0 produces the complete LU
factorization, which is the default.
milu stands for modified incomplete LU factorization. Its value can
be 'row' (row-sum), 'col' (column-sum), or 'off'. When milu is equal
to 'row', the diagonal element of the upper triangular factor U is
compensated in such a way as to preserve row sums. That is, the
product A*e is equal to L*U*e, where e is the vector of ones. When
milu is equal to 'col', the diagonal of the upper triangular factor U
is compensated so that column sums are preserved. That is, the
product e'*A is equal to the product e'*L*U. The default is 'off'.
udiag is either 0 or 1. If it is 1, any zero diagonal entries of the
upper triangular factor U are replaced by the local drop tolerance in
an attempt to avoid a singular factor. The default is 0.
thresh is a pivot threshold in [0,1]. Pivoting occurs when the
diagonal entry in a column has magnitude less than thresh times the
magnitude of any sub-diagonal entry in that column. thresh = 0 forces
diagonal pivoting. thresh = 1 is the default.
For SETUP.type == 'nofill', only SETUP.milu is used; all other fields
are ignored. For SETUP.type == 'crout', only SETUP.droptol and
SETUP.milu are used; all other fields are ignored.
W = ilu(A,SETUP) returns "L+U-speye(size(A))" where L is unit lower
triangular and U is upper triangular. If SETUP.type == 'ilutp', the
permutation information is lost.
[L,U] = ilu(A,SETUP) returns unit lower triangular L and upper
triangular U when SETUP.type == 'nofill' or when SETUP.type ==
'crout'. For SETUP.type == 'ilutp', one of the factors is permuted
based on the value of SETUP.milu. When SETUP.milu == 'row', U is a
column permuted upper triangular factor. Otherwise, L is a
row-permuted unit lower triangular factor.
[L,U,P] = ilu(A,SETUP) returns unit lower triangular L, upper
triangular U and a permutation matrix P. When SETUP.type == 'nofill'
or when SETUP.type == 'crout', P is always an identity matrix as
neither of these methods performs pivoting. When SETUP.type =
'ilutp', the role of P is determined by the value of SETUP.milu. When
SETUP.milu ~= 'row', P is returned such that L and U are incomplete
factors of P*A. When SETUP.milu == 'row', P is returned such that L
and U are incomplete factors of A*P.
Example:
A = gallery('neumann', 1600) + speye(1600);
setup.type = 'nofill';
nnz(A)
nnz(lu(A))
nnz(ilu(A,setup))
This shows that A has 7840 nonzeros, its complete LU factorization has
126478 nonzeros, and its incomplete LU factorization with 0 level of
fill-in has 7840 nonzeros, the same amount as A.
Example:
A = gallery('neumann', 1600) + speye(1600);
setup.type = 'crout';
setup.milu = 'row';
setup.droptol = 0.1;
[L,U] = ilu(A,setup);
e = ones(size(A,2),1);
norm(A*e-L*U*e)
This shows that A and L*U, where L and U are given by the modified
Crout ilu, have the same row-sum.
ilu works only for sparse matrices.
See also ichol, lu, gmres, bicg.
Reference page in Help browser
doc ilu
diary off