# HG changeset patch # User dbateman # Date 1113518448 0 # Node ID bb224e6c26a73dfbd820720eb05aa1b97434293a # Parent 5bdc3f24cd5f9371698aa129e178d6b95a69a67e [project @ 2005-04-14 22:40:48 by dbateman] diff -r 5bdc3f24cd5f -r bb224e6c26a7 src/DLD-FUNCTIONS/luinc.cc --- a/src/DLD-FUNCTIONS/luinc.cc Thu Apr 14 22:17:27 2005 +0000 +++ b/src/DLD-FUNCTIONS/luinc.cc Thu Apr 14 22:40:48 2005 +0000 @@ -265,100 +265,3 @@ ;;; mode: C++ *** ;;; End: *** */ - -/* - LUINC Sparse Incomplete LU factorization. - LUINC produces two different kinds of incomplete LU factorizations -- the - drop tolerance and the 0 level of fill-in factorizations. These factors - may be useful as preconditioners for a system of linear equations being - solved by an iterative method such as BICG (BiConjugate Gradients). - - LUINC(X,DROPTOL) performs the incomplete LU factorization of - X with drop tolerance DROPTOL. - - LUINC(X,OPTS) allows additional options to the incomplete LU - factorization. OPTS is a structure with up to four fields: - droptol --- the drop tolerance of incomplete LU - milu --- modified incomplete LU - udiag --- replace zeros on the diagonal of U - thresh --- the pivot threshold (see also LU) - - Only the fields of interest need to be set. - - droptol is a non-negative scalar used as the drop - tolerance for the incomplete LU factorization. This factorization - is computed in the same (column-oriented) manner as the - LU factorization except after each column of L and U has - been calculated, all entries in that column which are smaller - in magnitude than the local drop tolerance, which is - droptol * NORM of the column of X, are "dropped" from L or U. - The only exception to this dropping rule is the diagonal of the - upper triangular factor U which remains even if it is too small. - 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 is either 0 (unmodified, the default) or 1 (modified). - Instead of discarding those entries from the newly-formed - column of the factors, they are subtracted from the diagonal - of the upper triangular factor, U. - - 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. - - Example: - - load west0479 - A = west0479; - nnz(A) - nnz(lu(A)) - nnz(luinc(A,1e-6)) - - This shows that A has 1887 nonzeros, its complete LU factorization - has 16777 nonzeros, and its incomplete LU factorization with a - drop tolerance of 1e-6 has 10311 nonzeros. - - - [L,U,P] = LUINC(X,'0') produces the incomplete LU factors of a sparse - matrix with 0 level of fill-in (i.e. no fill-in). L is unit lower - trianglar, U is upper triangular and P is a permutation matrix. U has the - same sparsity pattern as triu(P*X). L has the same sparsity pattern as - tril(P*X), except for 1's on the diagonal of L where P*X may be zero. Both - L and U may have a zero because of cancellation where P*X is nonzero. L*U - differs from P*X only outside of the sparsity pattern of P*X. - - [L,U] = LUINC(X,'0') produces upper triangular U and L is a permutation of - unit lower triangular matrix. Thus no comparison can be made between the - sparsity patterns of L,U and X, although nnz(L) + nnz(U) = nnz(X) + n. L*U - differs from X only outside of its sparsity pattern. - - LU = LUINC(X,'0') returns "L+U-I", where L is unit lower triangular, U is - upper triangular and the permutation information is lost. - - Example: - - load west0479 - A = west0479; - [L,U,P] = luinc(A,'0'); - isequal(spones(U),spones(triu(P*A))) - spones(L) ~= spones(tril(P*A)) - D = (L*U) .* spones(P*A) - P*A - - spones(L) differs from spones(tril(P*A)) at some positions on the - diagonal and at one position in L where cancellation zeroed out a - nonzero element of P*A. The entries of D are of the order of eps. - - LUINC works only for sparse matrices. - - See also LU, CHOLINC, BICG. - -*/