NAG Library Routine Document
F11GDF
1 Purpose
F11GDF is a setup routine, the first in a suite of three routines for the iterative solution of a symmetric system of simultaneous linear equations. F11GDF must be called before the iterative solver,
F11GEF. The third routine in the suite,
F11GFF, can be used to return additional information about the computation.
These three routines are suitable for the solution of large sparse symmetric systems of equations.
2 Specification
SUBROUTINE F11GDF ( |
METHOD, PRECON, SIGCMP, NORM, WEIGHT, ITERM, N, TOL, MAXITN, ANORM, SIGMAX, SIGTOL, MAXITS, MONIT, LWREQ, WORK, LWORK, IFAIL) |
INTEGER |
ITERM, N, MAXITN, MAXITS, MONIT, LWREQ, LWORK, IFAIL |
REAL (KIND=nag_wp) |
TOL, ANORM, SIGMAX, SIGTOL, WORK(LWORK) |
CHARACTER(*) |
METHOD |
CHARACTER(1) |
PRECON, SIGCMP, NORM, WEIGHT |
|
3 Description
The suite consisting of the routines F11GDF,
F11GEF and
F11GFF is designed to solve the symmetric system of simultaneous linear equations
of order
, where
is large and the matrix of the coefficients
is sparse.
F11GDF is a setup routine which must be called before
F11GEF, the iterative solver. The third routine in the suite,
F11GFF can be used to return additional information about the computation. One of the following methods can be used:
- Conjugate Gradient Method (CG)
For this method (see
Hestenes and Stiefel (1952),
Golub and Van Loan (1996),
Barrett et al. (1994) and
Dias da Cunha and Hopkins (1994)), the matrix
should ideally be positive definite. The application of the Conjugate Gradient method to indefinite matrices may lead to failure or to lack of convergence.
- Lanczos Method (SYMMLQ)
This method, based upon the algorithm SYMMLQ (see
Paige and Saunders (1975) and
Barrett et al. (1994)), is suitable for both positive definite and indefinite matrices. It is more robust than the Conjugate Gradient method but less efficient when
is positive definite.
- Minimum Residual Method (MINRES)
This method may be used when the matrix is indefinite. It seeks to reduce the norm of the residual at each iteration and often takes fewer iterations than the other methods. It does however require slightly more memory.
The CG and SYMMLQ methods start from the residual
, where
is an initial estimate for the solution (often
), and generate an orthogonal basis for the Krylov subspace
, for
, by means of three-term recurrence relations (see
Golub and Van Loan (1996)). A sequence of symmetric tridiagonal matrices
is also generated. Here and in the following, the index
denotes the iteration count. The resulting symmetric tridiagonal systems of equations are usually more easily solved than the original problem. A sequence of solution iterates
is thus generated such that the sequence of the norms of the residuals
converges to a required tolerance. Note that, in general, the convergence is not monotonic.
In exact arithmetic, after iterations, this process is equivalent to an orthogonal reduction of to symmetric tridiagonal form, ; the solution would thus achieve exact convergence. In finite-precision arithmetic, cancellation and round-off errors accumulate causing loss of orthogonality. These methods must therefore be viewed as genuinely iterative methods, able to converge to a solution within a prescribed tolerance.
The orthogonal basis is not formed explicitly in either method. The basic difference between the Conjugate Gradient and Lanczos methods lies in the method of solution of the resulting symmetric tridiagonal systems of equations: the conjugate gradient method is equivalent to carrying out an (Cholesky) factorization whereas the Lanczos method (SYMMLQ) uses an factorization.
Faster convergence for all the methods can be achieved using a
preconditioner (see
Golub and Van Loan (1996) and
Barrett et al. (1994)). A preconditioner maps the original system of equations onto a different system, say
with, hopefully, better characteristics with respect to its speed of convergence: for example, the condition number of the matrix of the coefficients can be improved or eigenvalues in its spectrum can be made to coalesce. An orthogonal basis for the Krylov subspace
, for
, is generated and the solution proceeds as outlined above. The algorithms used are such that the solution and residual iterates of the original system are produced, not their preconditioned counterparts. Note that an unsuitable preconditioner or no preconditioning at all may result in a very slow rate, or lack, of convergence. However, preconditioning involves a trade-off between the reduction in the number of iterations required for convergence and the additional computational costs per iteration. Also, setting up a preconditioner may involve non-negligible overheads.
A preconditioner must be
symmetric and positive definite, i.e., representable by
, where
is nonsingular, and such that
in
(1), where
is the identity matrix of order
. Also, we can define
and
. These are formal definitions, used only in the design of the algorithms; in practice, only the means to compute the matrix-vector products
and to solve the preconditioning equations
are required, that is, explicit information about
,
or their inverses is not required at any stage.
The first termination criterion
is available for both conjugate gradient and Lanczos (SYMMLQ) methods. In
(2),
and
denotes a user-specified tolerance subject to
, where
is the
machine precision. Facilities are provided for the estimation of the norm of the matrix of the coefficients
, when this is not known in advance, used in
(2), by applying Higham's method (see
Higham (1988)). Note that
cannot be estimated internally. This criterion uses an error bound derived from
backward error analysis to ensure that the computed solution is the exact solution of a problem as close to the original as the termination tolerance requires. Termination criteria employing bounds derived from
forward error analysis could be used, but any such criteria would require information about the condition number
which is not easily obtainable.
The second termination criterion
is available only for the Lanczos method (SYMMLQ). In
(3),
is the largest singular value of the (preconditioned) iteration matrix
. This termination criterion monitors the progress of the solution of the preconditioned system of equations and is less expensive to apply than criterion
(2). When
is not supplied, facilities are provided for its estimation by
. The interlacing property
and Gerschgorin's theorem provide lower and upper bounds from which
can be easily computed by bisection. Alternatively, the less expensive estimate
can be used, where
by Gerschgorin's theorem. Note that only order of magnitude estimates are required by the termination criterion.
Termination criterion
(2) is the recommended choice, despite its (small) additional costs per iteration when using the Lanczos method (SYMMLQ). Also, if the norm of the initial estimate is much larger than the norm of the solution, that is, if
, a dramatic loss of significant digits could result in complete lack of convergence. The use of criterion
(2) will enable the detection of such a situation, and the iteration will be restarted at a suitable point. No such restart facilities are provided for criterion
(3).
Optionally, a vector
of user-specified weights can be used in the computation of the vector norms in termination criterion
(2), i.e.,
, where
, for
. Note that the use of weights increases the computational costs.
The MINRES algorithm terminates when the norm of the residual of the preconditioned system , , where is the preconditioned matrix.
The sequence of calls to the routines comprising the suite is enforced: first, the setup routine F11GDF must be called, followed by the solver
F11GEF. The diagnostic routine
F11GFF can be called either when
F11GEF is carrying out a monitoring step or after
F11GEF has completed its tasks. Incorrect sequencing will raise an error condition.
4 References
Barrett R, Berry M, Chan T F, Demmel J, Donato J, Dongarra J, Eijkhout V, Pozo R, Romine C and Van der Vorst H (1994) Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods SIAM, Philadelphia
Dias da Cunha R and Hopkins T (1994) PIM 1.1 — the parallel iterative method package for systems of linear equations user's guide — Fortran 77 version Technical Report Computing Laboratory, University of Kent at Canterbury, Kent, UK
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Hestenes M and Stiefel E (1952) Methods of conjugate gradients for solving linear systems J. Res. Nat. Bur. Stand. 49 409–436
Higham N J (1988) FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation ACM Trans. Math. Software 14 381–396
Paige C C and Saunders M A (1975) Solution of sparse indefinite systems of linear equations SIAM J. Numer. Anal. 12 617–629
5 Parameters
- 1: METHOD – CHARACTER(*)Input
On entry: the iterative method to be used.
- Conjugate gradient method (CG).
- Lanczos method (SYMMLQ).
- Minimum residual method (MINRES).
Constraint:
, or .
- 2: PRECON – CHARACTER(1)Input
On entry: determines whether preconditioning is used.
- No preconditioning.
- Preconditioning.
Constraint:
or .
- 3: SIGCMP – CHARACTER(1)Input
On entry: determines whether an estimate of
, the largest singular value of the preconditioned matrix of the coefficients, is to be computed using the bisection method on the sequence of tridiagonal matrices
generated during the iteration. Note that
when a preconditioner is not used.
If
(see below), i.e., when
is supplied, the value of
SIGCMP is ignored.
- is to be computed using the bisection method.
- The bisection method is not used.
If the termination criterion
(3) is used, requiring
, an inexpensive estimate is computed and used (see
Section 3).
It is not used if .
Suggested value:
.
Constraint:
or .
- 4: NORM – CHARACTER(1)Input
On entry: if
or
,
NORM defines the matrix and vector norm to be used in the termination criteria.
- Use the norm.
- Use the norm.
- Use the norm.
It has no effect if .
Suggested values:
- if , ;
- if , .
Constraints:
- if , , or ;
- if , .
- 5: WEIGHT – CHARACTER(1)Input
On entry: specifies whether a vector
of user-supplied weights is to be used in the vector norms used in the computation of termination criterion
(2) (
):
, where
, for
. The suffix
denotes the vector norm used, as specified by the parameter
NORM. Note that weights cannot be used when
, i.e., when criterion
(3) is used.
- User-supplied weights are to be used and must be supplied on initial entry to F11GEF.
- All weights are implicitly set equal to one. Weights do not need to be supplied on initial entry to F11GEF.
It has no effect if .
Suggested value:
.
Constraints:
- if , or ;
- if , .
- 6: ITERM – INTEGERInput
On entry: defines the termination criterion to be used.
- Use the termination criterion defined in (2) (both conjugate gradient and Lanczos (SYMMLQ) methods).
- Use the termination criterion defined in (3) (Lanczos method (SYMMLQ) only).
It has no effect if .
Suggested value:
.
Constraints:
- if , ;
- if , or .
- 7: N – INTEGERInput
On entry: , the order of the matrix .
Constraint:
.
- 8: TOL – REAL (KIND=nag_wp)Input
On entry: the tolerance
for the termination criterion.
If , is used, where is the machine precision.
Otherwise is used.
Constraint:
.
- 9: MAXITN – INTEGERInput
On entry: the maximum number of iterations.
Constraint:
.
- 10: ANORM – REAL (KIND=nag_wp)Input
On entry: if
, the value of
to be used in the termination criterion
(2) (
).
If
,
and
or
, then
is estimated internally by
F11GEF.
If
, then
ANORM is not referenced.
It has no effect if .
Constraint:
if and , .
- 11: SIGMAX – REAL (KIND=nag_wp)Input
On entry: if
, the value of
.
If
,
is estimated by
F11GEF when either
or termination criterion
(3) (
) is employed, though it will be used only in the latter case.
Otherwise, or if
,
SIGMAX is not referenced.
- 12: SIGTOL – REAL (KIND=nag_wp)Input
On entry: the tolerance used in assessing the convergence of the estimate of
when the bisection method is used.
If , the default value is used. The actual value used is .
If
or
, then
SIGTOL is not referenced.
It has no effect if .
Suggested value:
should be sufficient in most cases.
Constraint:
if and , .
- 13: MAXITS – INTEGERInput
On entry: the maximum iteration number
for which
is computed by bisection (see also
Section 3). If
or
, or if
, then
MAXITS is not referenced.
Suggested value:
when
SIGTOL is of the order of its default value
.
Constraint:
if and , .
- 14: MONIT – INTEGERInput
On entry: if
, the frequency at which a monitoring step is executed by
F11GEF: the current solution and residual iterates will be returned by
F11GEF and a call to
F11GFF made possible every
MONIT iterations, starting from the (
MONIT)th. Otherwise, no monitoring takes place.
There are some additional computational costs involved in monitoring the solution and residual vectors when the Lanczos method (SYMMLQ) is used.
Constraint:
.
- 15: LWREQ – INTEGEROutput
On exit: the minimum amount of workspace required by
F11GEF. (See also
Section 5 in F11GEF.)
- 16: WORK(LWORK) – REAL (KIND=nag_wp) arrayCommunication Array
On exit: the array
WORK is initialized by F11GDF. It must
not be modified before calling the next routine in the suite, namely
F11GEF.
- 17: LWORK – INTEGERInput
On entry: the dimension of the array
WORK as declared in the (sub)program from which F11GDF is called.
Constraint:
.
Note: although the minimum value of
LWORK ensures the correct functioning of F11GDF, a larger value is required by the other routines in the suite, namely
F11GEF and
F11GFF. The required value is as follows:
Method |
Requirements |
CG |
. |
SYMMLQ |
, |
MINRES |
, |
where
- , when an estimate of (SIGMAX) is computed;
- , otherwise.
- 18: IFAIL – INTEGERInput/Output
-
On entry:
IFAIL must be set to
,
. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction 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 parameter, 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, the th argument had an illegal value.
F11GDF has been called out of sequence.
7 Accuracy
Not applicable.
When
is not supplied (
) but it is required, it is estimated by
F11GEF using either of the two methods described in
Section 3, as specified by the parameter
SIGCMP. In particular, if
, then the computation of
is deemed to have converged when the differences between three successive values of
differ, in a relative sense, by less than the tolerance
SIGTOL, i.e., when
The computation of
is also terminated when the iteration count exceeds the maximum value allowed, i.e.,
.
Bisection is increasingly expensive with increasing iteration count. A reasonably large value of
SIGTOL, of the order of the suggested value, is recommended and an excessive value of
MAXITS should be avoided. Under these conditions,
usually converges within very few iterations.
9 Example
This example solves a symmetric system of simultaneous linear equations using the conjugate gradient method, where the matrix of the coefficients
, has a random sparsity pattern. An incomplete Cholesky preconditioner is used (
F11JAF and
F11JBF).
9.1 Program Text
Program Text (f11gdfe.f90)
9.2 Program Data
Program Data (f11gdfe.d)
9.3 Program Results
Program Results (f11gdfe.r)