` MGStep`, ` MGIter`, ` NestedMGIter`,
` MGPCGIter`,
` BPXPrecond`, ` BPXPCGIter`
-- multilevel solvers

#include <laspack/mlsolv.h> Vector *MGStep(int NoLevels, QMatrix *A, Vector *x, Vector *b, Matrix *R, Matrix *P, int Level, int Gamma, IterProcType SmoothProc, int Nu1, int Nu2, PrecondProcType PrecondProc, double Omega, IterProcType SolvProc, int NuC, PrecondProcType PrecondProcC, double OmegaC); Vector *MGIter(int NoLevels, QMatrix *A, Vector *x, Vector *b, Matrix *R, Matrix *P, int MaxIter, int Gamma, IterProcType SmoothProc, int Nu1, int Nu2, PrecondProcType PrecondProc, double Omega, IterProcType SolvProc, int NuC, PrecondProcType PrecondProcC, double OmegaC); Vector *NestedMGIter(int NoLevels, QMatrix *A, Vector *x, Vector *b, Matrix *R, Matrix *P, int Gamma, IterProcType SmoothProc, int Nu1, int Nu2, PrecondProcType PrecondProc, double Omega, IterProcType SolvProc, int NuC, PrecondProcType PrecondProcC, double OmegaC); Vector *MGPCGIter(int NoLevels, QMatrix *A, Vector *x, Vector *b, Matrix *R, Matrix *P, int MaxIter, int NoMGIter, int Gamma, IterProcType SmoothProc, int Nu1, int Nu2, PrecondProcType PrecondProc, double Omega, IterProcType SolvProc, int NuC, PrecondProcType PrecondProcC, double OmegaC); Vector *BPXPrecond(int NoLevels, QMatrix *A, Vector *y, Vector *c, Matrix *R, Matrix *P, int Level, IterProcType SmoothProc, int Nu, PrecondProcType PrecondProc, double Omega, IterProcType SmoothProcC, int NuC, PrecondProcType PrecondProcC, double OmegaC); Vector *BPXPCGIter(int NoLevels, QMatrix *A, Vector *x, Vector *b, Matrix *R, Matrix *P, int MaxIter, IterProcType SmoothProc, int Nu, PrecondProcType PrecondProc, double Omega, IterProcType SmoothProcC, int NuC, PrecondProcType PrecondProcC, double OmegaC);

The number of grid levels is set by the parameter ` NoLevels`.
Matrices ` A[i]` and vectors ` x[i]` and ` b[i]`
should be defined at each level,
from the coarsest level ` (i = 0)`
to
the finest level
` (i = NoLevels - 1).`
Furthermore, the initial approach ` x[Level]` of the solution vector
and the right hand side vector ` b[Level]`
are required at the actual level ` Level`.

Matrices ` R[i]` ` (0 <= i <= NoLevels - 2)`
should define restriction operators
which map a vector from level ` i` + 1 to level ` i`.
Prolongation operators are set
by matrices ` P[i]` (` 1 <= i <= NoLevels - 1)`
which have to map a vector from level level ` i` - 1 to level ` i`.
If no prolongation is specified (i.e. if ` P == NULL`),
the transpose of ` R[i - 1]` is used instead of ` P[i]`.

The type of the multigrid iteration is defined by the parameter ` Gamma`:
1 corresponds to a V-cycle, 2 to a W-cycle.

The parameters ` SmoothProc` and ` SolvProc` should
address routines with the prototype:

Vector *(*IterProcType)(QMatrix *, Vector *, Vector *, int, PrecondProcType, double)

The procedure ` SmoothProc` is applied as smoothing method.
For pre-smoothing, ` Nu1` iterations are carried out,
for post-smoothing, ` Nu2` iterations.
Preconditioning can be specified by the procedure ` PrecondProc`
with the prototype:

Vector *(*PrecondProcType)(QMatrix *, Vector *, Vector *, double)The parameter

The system of equations on the coarsest grid is solved by the
procedure ` SolvProc` with parameters ` NuC`, ` PrecondProcC`,
and ` OmegaC`.

If any preconditioner was specified by the procedures
` PrecondProc` and ` PrecondProcC`
(i.e. if ` PrecondProc != NULL` and ` PrecondProcC != NULL`,
respectively),
the preconditioned variant of the corresponding iterative methods
is used (if it is available).
The parameters ` Omega` and ` OmegaC`
are passed as relaxation parameters to the preconditioner
in this case.

In all above solvers,
the solution process is terminated either after ` MaxIter` iteration steps
or if the residual on the finest grid satisfies the condition

|| r ||_2 = || b - A x ||_2 <= eps || b ||_2,

where ` eps`
is the accuracy defined
by ` LASPack`
termination control (module ** RTC**).

For systems with singular matrices, orthogonalization of the solution vector to the null space of the matrix (which should be specified) is applied to ensure unique solution as well as convergence in the resulting subspace.

The multigrid algorithms and their theoretical foundation are described in

W. Hackbusch: Multi-Grid Methods and Applications, Springer-Verlag, Berlin, 1985.The implementation of the BPX preconditioner is based on the theory given inS. F. McCormick: Multigrid Methods, SIAM, Philadelphia, 1987.

J. H. Bramble, J. E. Pasciak, J. Xu: Parallel multilevel preconditioners, Mathematics of Computations, 55 (1990), pp. 1-22.

` mlsolv.h ... ` header file

` mlsolv.c ... ` source file

In the following code fragment, a simple Poisson boundary value problem

- u''(x) = 1 for 0 <= x <= 1, u(0) = 0, u(1) = 0

should be solved by means of two-grid V-cycle multigrid method. For simplicity, equidistant grids with 7 and 3 inner points are used.

As smoother, SSOR method with one pre-smoothing and one post-smoothing step is applied. The relaxation parameter is set to 1.2. Solution on the coarse grid should be carried out by CG method with diagonal preconditioning. The maximum number of iterations is set to 10.

The multigrid iterations should be terminated if the accuracy of 1e-5 was reached but at least after 100 iterations.

For generation of the matrices ` L[0]` and ` L[1]`,
and of the restriction operator ` R[0]`,
look at examples of modules ** QMATRIX** and
** MATRIX**.

QMatrix A[2]; Vector x[2], b[2]; Matrix R[2]; size_t Row, Ind; double h; Q_Constr(&A[0], "A[0]", 3, True, Rowws, Normal, True); Q_Constr(&A[1], "A[1]", 7, True, Rowws, Normal, True); V_Constr(&x[0], "x[0]", 3, Normal, True); V_Constr(&x[1], "x[1]", 7, Normal, True); V_Constr(&b[0], "b[0]", 3, Normal, True); V_Constr(&b[1], "b[1]", 7, Normal, True); M_Constr(&R[0], "R[0]", 3, 7, Rowws, Normal, True); /* generation of the matrix A[0] for the coarse grid */ h = 1.0 / 4.0; ... /* generation of the matrix A[1] for the fine grid */ h = 1.0 / 8.0; ... /* generation of the right hand side vector for the fine grid */ h = 1.0 / 8.0; V_SetAllCmp(&b[1], h); /* generation of the restriction operator R[0] */ ... /* solution of the system of equations by the multigrid solver */ SetRTCAccuracy(1e-5); V_SetAllCmp(&x, 0.0); MGIter(2, &A, &x, &b, &R, NULL, 100, 1, SSORIter, 1, 1, NULL, 1.2, CGIter, 10, JacobiPrecond, 1.0); /* output of the solution vector */ for (Ind = 1; Ind <= Dim; Ind++) printf("%d %12.5e\n", Ind, V_GetCmp(&x[1], Ind)); Q_Destr(&A[0]); Q_Destr(&A[1]); V_Destr(&x[0]); V_Destr(&x[1]); V_Destr(&b[0]); V_Destr(&b[1]); V_Destr(&R[0]);

qmatrix(3LAS), vector(3LAS), matrix(3LAS), itersolv(3LAS), rtc(3LAS)

Tomas Skalicky (skalicky@msmfs1.mw.tu-dresden.de)