MLSOLV(3LAS)

Next: OPERATS(3LAS) Up: Manual Pages Previous: MATRIX(3LAS)

MLSOLV(3LAS)

NAME

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

SYNOPSIS

#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);

DESCRIPTION

The procedure MGStep performs one multigrid step based on a recursively use of a two-grid iteration. It forms the basic algorithm from which all other multigrid solvers are derived.

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 Omega is used as relaxation 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.

The procedure MGIter carries out multigrid iterations of a system of NoLevels grids. Meaning of the parameters agrees with that of MGStep.

The procedure NestedMGIter performs nested multigrid iterations (also referred to as Full Multigrid Method). The parameters correspond to those of MGStep.

The procedure MGPCGIter is a multilevel solver based on the preconditioned conjugate gradient method. For preconditioning, NoMGIter multigrid iterations are applied to NoLevels grid levels. Meaning of the parameters agrees with that of MGStep.

The procedure BPXPCGIter carries out conjugate gradient method in connection with the BPX multilevel preconditioner which is implemented as procedure BPXPrecond. At the preconditioning step, residuals are smoothed at each level. The procedure SmoothProc with parameters Nu, PrecondProc, and Omega is used as smoothing method. In order to keep efficiency even if the system of equations on coarsest level has fairly high dimension, the smoother is for this level given separately by the parameters SmoothProcC, NuC, PrecondProcC, and OmegaC. All remaining parameters correspond to those of MGStep.

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.

REFERENCES

The multigrid algorithms and their theoretical foundation are described in

W. Hackbusch: Multi-Grid Methods and Applications, Springer-Verlag, Berlin, 1985.
S. F. McCormick: Multigrid Methods, SIAM, Philadelphia, 1987.

The implementation of the BPX preconditioner is based on the theory given in

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

FILES

mlsolv.h ... header file
mlsolv.c ... source file

EXAMPLE

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]);