Documentation 
Sparse incomplete LU factorization
ilu(A,setup)
[L,U] = ilu(A,setup)
[L,U,P] = ilu(A,setup)
ilu produces a unit lower triangular matrix, an upper triangular matrix, and a permutation matrix.
ilu(A,setup) computes the incomplete LU factorization of A. setup is an input structure with up to five setup 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 

type  Type of factorization. Values for type include:
If type is not specified, the ILU factorization with 0 level of fill in is performed. Pivoting is only performed with type set to 'ilutp'. 
droptol  Drop tolerance of the incomplete LU factorization. droptol is a nonnegative scalar. The default value is 0, which produces the complete LU factorization. The nonzero entries of U satisfy abs(U(i,j)) >= droptol*norm(A(:,j)), with the exception of the diagonal entries, which are retained regardless of satisfying the criterion. The entries of L are 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). 
milu  Modified incomplete LU factorization. Values for milu include:

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. 
ilu(A,setup) returns L+Uspeye(size(A)), where L is a unit lower triangular matrix and U is an upper triangular matrix.
[L,U] = ilu(A,setup) returns a unit lower triangular matrix in L and an upper triangular matrix in U.
[L,U,P] = ilu(A,setup) returns a unit lower triangular matrix in L, an upper triangular matrix in U, and a permutation matrix in P.
Start with a sparse matrix and compute the LU factorization.
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*eL*U*e) ans = 1.4251e014
This shows that A and L*U, where L and U are given by the modified Crout ILU, have the same rowsum.
Start with a sparse matrix and compute the LU factorization.
A = gallery('neumann', 1600) + speye(1600); setup.type = 'nofill'; nnz(A) ans = 7840 nnz(lu(A)) ans = 126478 nnz(ilu(A,setup)) ans = 7840
This shows that A has 7840 nonzeros, the complete LU factorization has 126478 nonzeros, and the incomplete LU factorization, with 0 level of fillin, has 7840 nonzeros, the same amount as A.