NAG FL Interface
d03edf (dim2_ellip_mgrid)
1
Purpose
d03edf solves seven-diagonal systems of linear equations which arise from the discretization of an elliptic partial differential equation on a rectangular region. This routine uses a multigrid technique.
2
Specification
Fortran Interface
Subroutine d03edf ( |
ngx, ngy, lda, a, rhs, ub, maxit, acc, us, u, iout, numit, ifail) |
Integer, Intent (In) |
:: |
ngx, ngy, lda, maxit, iout |
Integer, Intent (Inout) |
:: |
ifail |
Integer, Intent (Out) |
:: |
numit |
Real (Kind=nag_wp), Intent (In) |
:: |
acc |
Real (Kind=nag_wp), Intent (Inout) |
:: |
a(lda,7), rhs(lda), ub(ngx*ngy) |
Real (Kind=nag_wp), Intent (Out) |
:: |
us(lda), u(lda) |
|
C Header Interface
#include <nag.h>
void |
d03edf_ (const Integer *ngx, const Integer *ngy, const Integer *lda, double a[], double rhs[], double ub[], const Integer *maxit, const double *acc, double us[], double u[], const Integer *iout, Integer *numit, Integer *ifail) |
|
C++ Header Interface
#include <nag.h> extern "C" {
void |
d03edf_ (const Integer &ngx, const Integer &ngy, const Integer &lda, double a[], double rhs[], double ub[], const Integer &maxit, const double &acc, double us[], double u[], const Integer &iout, Integer &numit, Integer &ifail) |
}
|
The routine may be called by the names d03edf or nagf_pde_dim2_ellip_mgrid.
3
Description
d03edf solves, by multigrid iteration, the seven-point scheme
which arises from the discretization of an elliptic partial differential equation of the form
and its boundary conditions, defined on a rectangular region. This we write in matrix form as
The algorithm is described in separate reports by
Wesseling (1982a),
Wesseling (1982b) and
McCarthy (1983).
Systems of linear equations, matching the seven-point stencil defined above, are solved by a multigrid iteration. An initial estimate of the solution must be provided by you. A zero guess may be supplied if no better approximation is available.
A ‘smoother’ based on incomplete Crout decomposition is used to eliminate the high frequency components of the error. A restriction operator is then used to map the system on to a sequence of coarser grids. The errors are then smoothed and prolongated (mapped onto successively finer grids). When the finest cycle is reached, the approximation to the solution is corrected. The cycle is repeated for
maxit iterations or until the required accuracy,
acc, is reached.
d03edf will automatically determine the number
of possible coarse grids, ‘levels’ of the multigrid scheme, for a particular problem. In other words,
d03edf determines the maximum integer
so that
and
can be expressed in the form
It should be noted that the rate of convergence improves significantly with the number of levels used (see
McCarthy (1983)), so that
and
should be carefully chosen so that
and
have factors of the form
, with
as large as possible. For good convergence the integer
should be at least
.
d03edf has been found to be robust in application, but being an iterative method the problem of divergence can arise. For a strictly diagonally dominant matrix
no such problem is foreseen. The diagonal dominance of
is not a necessary condition, but should this condition be strongly violated then divergence may occur. The quickest test is to try the routine.
4
References
McCarthy G J (1983) Investigation into the multigrid code MGD1 Report AERE-R 10889 Harwell
Wesseling P (1982a) MGD1 – a robust and efficient multigrid method Multigrid Methods. Lecture Notes in Mathematics 960 614–630 Springer–Verlag
Wesseling P (1982b) Theoretical aspects of a multigrid method SIAM J. Sci. Statist. Comput. 3 387–407
5
Arguments
-
1:
– Integer
Input
-
On entry: the number of interior grid points in the -direction, . should preferably be divisible by as high a power of as possible.
Constraint:
.
-
2:
– Integer
Input
-
On entry: the number of interior grid points in the -direction, . should preferably be divisible by as high a power of as possible.
Constraint:
.
-
3:
– Integer
Input
-
On entry: the first dimension of the array
a, which must also be a lower bound for the dimension of the arrays
rhs,
us and
u as declared in the (sub)program from which
d03edf is called. It is always sufficient to set
, but slightly smaller values may be permitted, depending on the values of
ngx and
ngy. If on entry,
lda is too small, an error message gives the minimum permitted value. (
lda must be large enough to allow space for the coarse-grid approximations.)
-
4:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: must be set to , for , and .
On exit: is overwritten.
-
5:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: must be set to , for and .
On exit: the first elements are unchanged and the rest of the array is used as workspace.
-
6:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: must be set to the initial estimate for the solution .
On exit: the corresponding component of the residual .
-
7:
– Integer
Input
-
On entry: the maximum permitted number of multigrid iterations. If
, no multigrid iterations are performed, but the coarse-grid approximations and incomplete Crout decompositions are computed, and may be output if
iout is set accordingly.
Constraint:
.
-
8:
– Real (Kind=nag_wp)
Input
-
On entry: the required tolerance for convergence of the residual
-norm:
where
and
is the computed solution. Note that the norm is not scaled by the number of equations. The routine will stop after fewer than
maxit iterations if the residual
-norm is less than the specified tolerance. (If
, at least one iteration is always performed.)
If on entry
, the
machine precision is used as a default value for the tolerance; if
, but
acc is less than the
machine precision, the routine will stop when the residual
-norm is less than the
machine precision and
ifail will be set to
.
Constraint:
.
-
9:
– Real (Kind=nag_wp) array
Output
-
On exit: the residual -norm, stored in element .
-
10:
– Real (Kind=nag_wp) array
Output
-
On exit: the computed solution is returned in , for and .
-
11:
– Integer
Input
-
On entry: controls the output of printed information to the current advisory message unit (see
x04abf):
- No output.
- The solution
, for and .
- The residual -norm after each iteration, with the reduction factor over the previous iteration.
- As for and .
- As for , plus the final residual (as returned in ub).
- As for , plus the initial elements of a and rhs.
- As for , plus the Galerkin coarse grid approximations.
- As for , plus the incomplete Crout decompositions.
- As for , plus the residual after each iteration.
The elements
, the Galerkin coarse grid approximations and the incomplete Crout decompositions are output in the format:
- Y-index
- X-index
- where
, for and .
The vectors
,
,
are output in matrix form with
ngy rows and
ngx columns. Where
, the
ngx values for a given
value are produced in rows of
. Values of
may yield considerable amounts of output.
Constraint:
.
-
12:
– Integer
Output
-
On exit: the number of iterations performed.
-
13:
– Integer
Input/Output
-
On entry:
ifail must be set to
,
. If you are unfamiliar with this argument you should refer to
Section 4 in the Introduction to the NAG Library FL Interface for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
.
When the value is used it is essential to test the value of ifail on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
-
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry,
.
Constraint:
lda must be at least
.
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
-
After
maxit iterations the residual norm is not less than the tolerance
, residual norm
, tolerance
. The residual norm has decreased at each iteration after the first.
-
After
maxit iterations the residual norm is not less than the tolerance
, residual norm
, tolerance
. The residual norm increased at one or more iterations after the first.
-
On entry,
acc is less than
machine precision. The routine terminated because the residual norm is less than
machine precision. residual norm
,
machine precision and
.
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 7 in the Introduction to the NAG Library FL Interface for further information.
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library FL Interface for further information.
Dynamic memory allocation failed.
See
Section 9 in the Introduction to the NAG Library FL Interface for further information.
7
Accuracy
8
Parallelism and Performance
d03edf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
The rate of convergence of this routine is strongly dependent upon the number of levels,
, in the multigrid scheme, and thus the choice of
ngx and
ngy is very important. You are advised to experiment with different values of
ngx and
ngy to see the effect they have on the rate of convergence; for example, using a value such as
(
) followed by
(for which
).
10
Example
The program solves the elliptic partial differential equation
on the unit square
, with boundary conditions
For the equation to be elliptic,
must be less than
.
The equation is discretized on a square grid with mesh spacing
in both directions using the following approximations:
Thus the following equations are solved:
10.1
Program Text
10.2
Program Data
10.3
Program Results