JacobiPrecond, SSORPrecond, ILUPrecond -- pre-defined preconditioners
#include <laspack/precond.h> typedef Vector *(*PrecondProcType)(QMatrix *, Vector *, Vector *, double) Vector *JacobiPrecond(QMatrix *A, Vector *y, Vector *c, double Omega); Vector *SSORPrecond(QMatrix *A, Vector *y, Vector *c, double Omega); Vector *ILUPrecond(QMatrix *A, Vector *y, Vector *c, double Omega);
In LASPack , three preconditioners are currently available. They will be applied to the system
W y = c
which arises during the solution of the preconditioned systems of equations
W^{-1} A x = W^{-1} b.
The matrix W is in all three cases derived from the matrix A . The parameter Omega is used as relaxation parameter.
W = 1 / Omega . diag(A)
so that
y = Omega diag(A)^{-1} c.
A = D + L + U.
The procedure SSORPrecond uses then
W = 1 / (2 - Omega) . (D / Omega + L) (D / Omega)^{-1} (D / Omega + U).
Hence follows
y = (2 - Omega) / Omega . (D / Omega + U)^{-1} D (D / Omega + L)^{-1} c.
A = W + R = (D + L) D^{-1} (D + U) + R
Here 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.
The matrix W can be inverted exactly so that
y = (D + U)^{-1} D (D + L)^{-1} c.
All three preconditioners are applicable for both, symmetric and non-symmetric matrices. For symmetric matrices, the symmetry condition L = U^T is taken into account.
For description and theoretical foundation of above preconditioners see e.g.:
W. Hackbusch: Iterative Solution of Large Sparse Systems of Equations, Springer-Verlag, Berlin, 1994.
precond.h ... header file
precond.c ... source file
Preconditioners as well as iterative solvers are defined in LASPack by means of matrix-vector operations, especially the procedures Mul_QV and MulInv_QV from module OPERATS. This keep they independently on the matrix storage format. The Jacobi preconditioner is e.g. implemented by the following one-line statement:
Vector *JacobiPrecond(QMatrix *A, Vector *y, Vector *c, double Omega) { Q_Lock(A); V_Lock(y); V_Lock(c); Asgn_VV(y, Mul_SV(Omega, MulInv_QV(Diag_Q(A), c))); Q_Unlock(A); V_Unlock(y); V_Unlock(c); return(y); }
qmatrix(3LAS), vector(3LAS), factor(3LAS), itersolv(3LAS), mgsolv(3LAS), operats(3LAS)