ILUFactor -- incomplete factorization of quadratic matrices
#include <laspack/factor.h> QMatrix *ILUFactor(QMatrix *Q);
Q = (D + L) D^{-1} (D + U) + R
is performed, where D, U, and L are certain diagonal, upper, and lower triangular matrices, respectively. The remainder matrix R contains fill elements, which have been ignored during the factorization process. For symmetric matrices, the incomplete Cholesky factorization with L = U^T is applied. In order to be able to perform the incomplete factorization for singular matrices too, regularization be means of increasing of diagonal entries is applied.
In the current implementation, matrices L and U are laid down by the same position of non-zero elements as in matrix Q.
Matrices are stored D, U, and L in connection with Q as matrix D + L + U of type QMatrix. This one is also the one returned by the procedure ILUFactor. How to extract the particular matrices from it, is shown in the example below.
The incomplete factorization is comprehensively describe and analyzed e.g. in:
W. Hackbusch: Iterative Solution of Large Sparse Systems of Equations, Springer-Verlag, Berlin, 1994.
factor.h ... header file
factor.c ... source file
The following example shown the usage of the procedure ILUFactor in implementation of an ILU preconditioner. This have to solve the system of equations
W y = c
which arises during the solution of preconditioned systems
W^{-1} A x = W^{-1} b
at which
W = (D + L) D^{-1} (D + U).
The corresponding LASPack routine could be build as follows:
Vector *ILUPrecond(QMatrix *A, Vector *y, Vector *c, double Omega) { Q_Lock(A); V_Lock(y); V_Lock(c); Asgn_VV(y, MulInv_QV(Add_QQ(Diag_Q(ILUFactor(A)), Upper_Q(ILUFactor(A))), Mul_QV(Diag_Q(ILUFactor(A)), MulInv_QV(Add_QQ(Diag_Q(ILUFactor(A)), Lower_Q(ILUFactor(A))), c)))); Q_Unlock(A); V_Unlock(y); V_Unlock(c); return(y); }
qmatrix(3LAS), operats(3LAS), errhandl(3LAS)
In the current implementation, it is assumed that the matrix Q has symmetric structure with regard to non-zero elements. Because during discretization of differential equations even such matrices arise, this restriction is for many applications not grave. For matrices which do not have suffice to the above condition, the error LASILUStructErr is raised.