NAG Library Routine Document
f11gef
(real_symm_basic_solver)
1
Purpose
f11gef is an iterative solver for a symmetric system of simultaneous linear equations;
f11gef is the second in a suite of three routines, where the first routine,
f11gdf, must be called prior to
f11gef to set up the suite, and 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
Fortran Interface
Integer, Intent (In)  :: 
lwork  Integer, Intent (Inout)  :: 
irevcm,
ifail  Real (Kind=nag_wp), Intent (In)  :: 
wgt(*)  Real (Kind=nag_wp), Intent (Inout)  :: 
u(*),
v(*),
work(lwork) 

C Header Interface
#include nagmk26.h
void 
f11gef_ (
Integer *irevcm,
double u[],
double v[],
const double wgt[],
double work[],
const Integer *lwork,
Integer *ifail) 

3
Description
f11gef solves the symmetric system of linear simultaneous equations
$Ax=b$ using the preconditioned conjugate gradient method (see
Hestenes and Stiefel (1952),
Golub and Van Loan (1996),
Barrett et al. (1994) and
Dias da Cunha and Hopkins (1994)), a preconditioned Lanczos method based upon the algorithm SYMMLQ (see
Paige and Saunders (1975) and
Barrett et al. (1994)), or the MINRES algorithm (see
Paige and Saunders (1975)).
For a general description of the methods employed you are referred to
Section 3 in
f11gdf.
f11gef can solve the system after the first routine in the suite,
f11gdf, has been called to initialize the computation and specify the method of solution. The third routine in the suite,
f11gff, can be used to return additional information generated by the computation during monitoring steps and after
f11gef has completed its tasks.
f11gef uses
reverse communication, i.e.,
f11gef returns repeatedly to the calling program with the argument
irevcm (see
Section 5) set to specified values which require the calling program to carry out a specific task: either to compute the matrixvector product
$v=Au$; to solve the preconditioning equation
$Mv=u$; to notify the completion of the computation; or, to allow the calling program to monitor the solution. Through the argument
irevcm the calling program can cause immediate or tidy termination of the execution. On final exit, the last iterates of the solution and of the residual vectors of the original system of equations are returned.
Reverse communication has the following advantages.
1. 
Maximum flexibility in the representation and storage of sparse matrices: all matrix operations are performed outside the solver routine, thereby avoiding the need for a complicated interface with enough flexibility to cope with all types of storage schemes and sparsity patterns. This applies also to preconditioners. 
2. 
Enhanced user interaction: you can closely monitor the solution and tidy or immediate termination can be requested. This is useful, for example, when alternative termination criteria are to be employed or in case of failure of the external routines used to perform matrix operations. 
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 onenorm 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
Arguments
Note: this routine uses
reverse communication. Its use involves an initial entry, intermediate exits and reentries, and a final exit, as indicated by the argument
irevcm. Between intermediate exits and reentries,
all arguments other than irevcm and v must remain unchanged.
 1: $\mathbf{irevcm}$ – IntegerInput/Output

On initial entry: ${\mathbf{irevcm}}=0$, otherwise an error condition will be raised.
On intermediate reentry: must either be unchanged from its previous exit value, or can have one of the following values.
 ${\mathbf{irevcm}}=5$
 Tidy termination: the computation will terminate at the end of the current iteration. Further reverse communication exits may occur depending on when the termination request is issued. f11gef will then return with the termination code ${\mathbf{irevcm}}=4$. Note that before calling f11gef with ${\mathbf{irevcm}}=5$ the calling program must have performed the tasks required by the value of irevcm returned by the previous call to f11gef, otherwise subsequently returned values may be invalid.
 ${\mathbf{irevcm}}=6$
 Immediate termination: f11gef will return immediately with termination code ${\mathbf{irevcm}}=4$ and with any useful information available. This includes the last iterate of the solution and, for conjugate gradient only, the last iterate of the residual vector. The residual vector is generally not available when the Lanczos method (SYMMLQ) is used. f11gef will then return with the termination code ${\mathbf{irevcm}}=4$.
Immediate termination may be useful, for example, when errors are detected during matrixvector multiplication or during the solution of the preconditioning equation.
Changing
irevcm to any other value between calls will result in an error.
On intermediate exit:
has the following meanings.
 ${\mathbf{irevcm}}=1$
 The calling program must compute the matrixvector product $v=Au$, where $u$ and $v$ are stored in u and v, respectively.
 ${\mathbf{irevcm}}=2$
 The calling program must solve the preconditioning equation $Mv=u$, where $u$ and $v$ are stored in u and v, respectively.
 ${\mathbf{irevcm}}=3$
 Monitoring step: the solution and residual at the current iteration are returned in the arrays u and v, respectively. No action by the calling program is required. To return additional information f11gff can be called at this step.
On final exit: if
${\mathbf{irevcm}}=4$,
f11gef has completed its tasks. The value of
ifail determines whether the iteration has been successfully completed, errors have been detected or the calling program has requested termination.
Constraint:
on initial entry,
${\mathbf{irevcm}}=0$; on reentry, either
irevcm must remain unchanged or be reset to
${\mathbf{irevcm}}=5$ or
$6$.
Note: any values you return to f11gef as part of the reverse communication procedure should not include floatingpoint NaN (Not a Number) or infinity values, since these are not handled by f11gef. If your code inadvertently does return any NaNs or infinities, f11gef is likely to produce unexpected results.
 2: $\mathbf{u}\left(*\right)$ – Real (Kind=nag_wp) arrayInput/Output

Note: the dimension of the array
u
must be at least
$\mathit{n}$.
On initial entry: an initial estimate, ${x}_{0}$, of the solution of the system of equations $Ax=b$.
On intermediate reentry: must remain unchanged.
On intermediate exit:
the returned value of
irevcm determines the contents of
u in the following way.
If
${\mathbf{irevcm}}=1$ or
$2$,
u holds the vector
$u$ on which the operation specified by
irevcm is to be carried out.
If
${\mathbf{irevcm}}=3$,
u holds the current iterate of the solution vector.
On final exit: if
${\mathbf{ifail}}={\mathbf{3}}$ or
${{\mathit{i}}}$, the array
u is unchanged from the initial entry to
f11gef. If
${\mathbf{ifail}}={\mathbf{1}}$, the array
u is unchanged from the last entry to
f11gef. Otherwise,
u holds the last iterate of the solution of the system of equations, for all returned values of
ifail.
 3: $\mathbf{v}\left(*\right)$ – Real (Kind=nag_wp) arrayInput/Output

Note: the dimension of the array
v
must be at least
$\mathit{n}$.
On initial entry: the righthand side $b$ of the system of equations $Ax=b$.
On intermediate reentry: the returned value of
irevcm determines the contents of
v in the following way.
If
${\mathbf{irevcm}}=1$ or
$2$,
v must store the vector
$v$, the result of the operation specified by the value of
irevcm returned by the previous call to
f11gef.
If
${\mathbf{irevcm}}=3$,
v must remain unchanged.
On intermediate exit:
if
${\mathbf{irevcm}}=3$,
v holds the current iterate of the residual vector. Note that this is an approximation to the true residual vector. Otherwise, it does not contain any useful information.
On final exit: if
${\mathbf{ifail}}={\mathbf{3}}$ or
${{\mathit{i}}}$, the array
v is unchanged from the last entry to
f11gef. If
${\mathbf{ifail}}={\mathbf{1}}$, the array
v is unchanged from the initial entry to
f11gef. If
${\mathbf{ifail}}={\mathbf{0}}$ or
${\mathbf{2}}$, the array
v contains the true residual vector of the system of equations (see also
Section 6). Otherwise,
v stores the last iterate of the residual vector unless the Lanczos method (SYMMLQ) was used and
${\mathbf{ifail}}\ge {\mathbf{5}}$, in which case
v is set to
$0.0$.
 4: $\mathbf{wgt}\left(*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the dimension of the array
wgt
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,\mathit{n}\right)$.
On entry: the usersupplied weights, if these are to be used in the computation of the vector norms in the termination criterion (see
Sections 3 and
5 in
f11gdf).
Weights are NOT used in the MINRES algorithm.
Constraint:
if weights are to be used, at least one element of
wgt must be nonzero.
 5: $\mathbf{work}\left({\mathbf{lwork}}\right)$ – Real (Kind=nag_wp) arrayCommunication Array

On initial entry: the array
work as returned by
f11gdf (see also
Section 5 in
f11gdf).
On intermediate reentry: must remain unchanged.
 6: $\mathbf{lwork}$ – IntegerInput

On initial entry: the dimension of the array
work as declared in the (sub)program from which
f11gef is called (see also
Section 3 in
f11gdf).
The required amount of workspace is as follows:
Method 
Requirements 
CG 
${\mathbf{lwork}}=120+5\mathit{n}+p$. 
SYMMLQ 
${\mathbf{lwork}}=120+6\mathit{n}+p$, 
MINRES 
${\mathbf{lwork}}=120+9\mathit{n}$, 
where
 $p=2*\left({\mathbf{maxits}}+1\right)$, when an estimate of ${\sigma}_{1}\left(A\right)$ (sigmax) is computed;
 $p=0$, otherwise.
Constraint:
${\mathbf{lwork}}\ge {\mathbf{lwreq}}$, where
lwreq is returned by
f11gdf.
 7: $\mathbf{ifail}$ – IntegerInput/Output

On initial entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, because for this routine the values of the output arguments may be useful even if
${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit, the recommended value is
$1$.
When the value $\mathbf{1}\text{ or}1$ is used it is essential to test the value of ifail on exit.
On intermediate exit:
the value of
ifail is meaningless and should be ignored.
On final exit: (i.e., when
${\mathbf{irevcm}}=4$)
${\mathbf{ifail}}={\mathbf{0}}$, unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

f11gef has already completed its tasks. You need to set a new problem.
 ${\mathbf{ifail}}=2$

The required accuracy could not be obtained. However, a reasonable accuracy may have been achieved.
Userrequested termination: the required accuracy could not be obtained. However, a reasonable accuracy may have been achieved.
 ${\mathbf{ifail}}=3$

Either
f11gdf was not called before calling this routine or it has returned an error.
 ${\mathbf{ifail}}=4$

Userrequested tidy termination. The solution has not converged after $\u2329\mathit{\text{value}}\u232a$ iterations.
 ${\mathbf{ifail}}=5$

The solution has not converged after $\u2329\mathit{\text{value}}\u232a$ iterations.
 ${\mathbf{ifail}}=6$

The preconditioner appears not to be positive definite. The computation cannot continue.
 ${\mathbf{ifail}}=7$

The matrix of the coefficients $A$ appears not to be positive definite. The computation cannot continue.
 ${\mathbf{ifail}}=8$

Userrequested immediate termination.
 ${\mathbf{ifail}}=9$

The matrix of the coefficients $A$ appears to be singular. The computation cannot continue.
 ${\mathbf{ifail}}=10$

The weights in array
wgt are all zero.
 ${\mathbf{ifail}}=1$

On initial entry, ${\mathbf{irevcm}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{irevcm}}=0$.
On intermediate reentry,
${\mathbf{irevcm}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: either
irevcm must be unchanged from its previous exit value or
${\mathbf{irevcm}}=5$ or
$6$.
 ${\mathbf{ifail}}=6$

On entry,
${\mathbf{lwork}}=\u2329\mathit{\text{value}}\u232a$.
Constraint:
${\mathbf{lwork}}\ge {\mathbf{lwreq}}$, where
lwreq is returned by
f11gdf.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
On completion, i.e.,
${\mathbf{irevcm}}=4$ on exit, the arrays
u and
v will return the solution and residual vectors,
${x}_{k}$ and
${r}_{k}=bA{x}_{k}$, respectively, at the
$k$th iteration, the last iteration performed, unless an immediate termination was requested and the Lanczos method (SYMMLQ) was used.
On successful completion, the termination criterion is satisfied to within the userspecified tolerance, as described in
Section 3 in
f11gdf. The computed values of the left and righthand sides of the termination criterion selected can be obtained by a call to
f11gff.
8
Parallelism and Performance
f11gef is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f11gef 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 implementationspecific information.
The number of operations carried out by f11gef for each iteration is likely to be principally determined by the computation of the matrixvector products $v=Au$ and by the solution of the preconditioning equation $Mv=u$ in the calling program. Each of these operations is carried out once every iteration.
The number of the remaining operations in f11gef for each iteration is approximately proportional to $\mathit{n}$. Note that the Lanczos method (SYMMLQ) requires a slightly larger number of operations than the conjugate gradient method.
The number of iterations required to achieve a prescribed accuracy cannot be easily determined at the onset, as it can depend dramatically on the conditioning and spectrum of the preconditioned matrix of the coefficients $\stackrel{}{A}={E}^{1}A{E}^{\mathrm{T}}$.
Additional matrixvector products are required for the computation of
${\Vert A\Vert}_{1}={\Vert A\Vert}_{\infty}$, when this has not been supplied to
f11gdf and is required by the termination criterion employed.
The number of operations required to compute
${\sigma}_{1}\left(\stackrel{}{A}\right)$ is negligible for reasonable values of
sigtol and
maxits (see
Sections 5 and
9 in
f11gdf).
If the termination criterion
${\Vert {r}_{k}\Vert}_{p}\le \tau \left({\Vert b\Vert}_{p}+{\Vert A\Vert}_{p}\times {\Vert {x}_{k}\Vert}_{p}\right)$ is used (see
Section 3 in
f11gdf) and
$\Vert {x}_{0}\Vert \gg \Vert {x}_{k}\Vert $, so that because of loss of significant digits the required accuracy could not be obtained, the iteration is restarted automatically at some suitable point:
f11gef sets
${x}_{0}={x}_{k}$ and the computation begins again. For particularly badly scaled problems, more than one restart may be necessary. Naturally, restarting adds to computational costs: it is recommended that the iteration should start from a value
${x}_{0}$ which is as close to the true solution
$\stackrel{~}{x}$ as can be estimated. Otherwise, the iteration should start from
${x}_{0}=0$.
10
Example
See
Section 10 in
f11gdf.