NAG Library Routine Document
F11BEF
1 Purpose
F11BEF is an iterative solver for a real general (nonsymmetric) system of simultaneous linear equations; F11BEF is the second in a suite of three routines, where the first routine,
F11BDF, must be called prior to F11BEF to set up the suite, and the third routine in the suite,
F11BFF, can be used to return additional information about the computation.
These three routines are suitable for the solution of large sparse general (nonsymmetric) systems of equations.
2 Specification
INTEGER 
IREVCM, LWORK, IFAIL 
REAL (KIND=nag_wp) 
U(*), V(*), WGT(*), WORK(LWORK) 

3 Description
F11BEF solves the general (nonsymmetric) system of linear simultaneous equations
$Ax=b$ of order
$\mathit{n}$, where
$\mathit{n}$ is large and the coefficient matrix
$A$ is sparse, using one of four available methods: RGMRES, the preconditioned restarted generalized minimum residual method (see
Saad and Schultz (1986)); CGS, the preconditioned conjugate gradient squared method (see
Sonneveld (1989)); BiCGSTAB(
$\ell $), the biconjugate gradient stabilized method of order
$\ell $ (see
Van der Vorst (1989) and
Sleijpen and Fokkema (1993)); or TFQMR, the transposefree quasiminimal residual method (see
Freund and Nachtigal (1991) and
Freund (1993)).
For a general description of the methods employed you are referred to
Section 3 in F11BDF.
F11BEF can solve the system after the first routine in the suite,
F11BDF, has been called to initialize the computation and specify the method of solution. The third routine in the suite,
F11BFF, can be used to return additional information generated by the computation, during monitoring steps and after F11BEF has completed its tasks.
F11BEF uses
reverse communication, i.e., it returns repeatedly to the calling program with the parameter
IREVCM (see
Section 5) set to specified values which require the calling program to carry out one of the following tasks:
– 
compute the matrixvector product $v=Au$ or $v={A}^{\mathrm{T}}u$ (the four methods require the matrix transposevector product only if ${\Vert A\Vert}_{1}$ or ${\Vert A\Vert}_{\infty}$ is estimated internally by Higham's method (see Higham (1988))); 
– 
solve the preconditioning equation $Mv=u$; 
– 
notify the completion of the computation; 
– 
allow the calling program to monitor the solution. 
Through the parameter
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.
 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.
 Enhanced user interaction: you can closely monitor the progress of 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
Freund R W (1993) A transposefree quasiminimal residual algorithm for nonHermitian linear systems SIAM J. Sci. Comput. 14 470–482
Freund R W and Nachtigal N (1991) QMR: a QuasiMinimal Residual Method for NonHermitian Linear Systems Numer. Math. 60 315–339
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
Saad Y and Schultz M (1986) GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems SIAM J. Sci. Statist. Comput. 7 856–869
Sleijpen G L G and Fokkema D R (1993) BiCGSTAB$\left(\ell \right)$ for linear equations involving matrices with complex spectrum ETNA 1 11–32
Sonneveld P (1989) CGS, a fast Lanczostype solver for nonsymmetric linear systems SIAM J. Sci. Statist. Comput. 10 36–52
Van der Vorst H (1989) BiCGSTAB, a fast and smoothly converging variant of BiCG for the solution of nonsymmetric linear systems SIAM J. Sci. Statist. Comput. 13 631–644
5 Parameters
Note: this routine uses
reverse communication. Its use involves an initial entry, intermediate exits and reentries, and a final exit, as indicated by the parameter
IREVCM. Between intermediate exits and reentries,
all parameters other than IREVCM and V must remain unchanged.
 1: 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. F11BEF will then return with the termination code ${\mathbf{IREVCM}}=4$. Note that before calling F11BEF with ${\mathbf{IREVCM}}=5$ the calling program must have performed the tasks required by the value of IREVCM returned by the previous call to F11BEF, otherwise subsequently returned values may be invalid.
 ${\mathbf{IREVCM}}=6$
 Immediate termination: F11BEF will return immediately with termination code ${\mathbf{IREVCM}}=4$ and with any useful information available. This includes the last iterate of the solution. The residual vector is generally not available.
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={A}^{\mathrm{T}}u$, where $u$ and $v$ are stored in U and V, respectively; RGMRES, CGS and BiCGSTAB($\ell $) methods return ${\mathbf{IREVCM}}=1$ only if the matrix norm ${\Vert A\Vert}_{1}$ or ${\Vert A\Vert}_{\infty}$ is estimated internally using Higham's method. This can only happen if ${\mathbf{ITERM}}=1$ in F11BDF.
 ${\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. F11BFF can be called at this step to return additional information.
On final exit:
${\mathbf{IREVCM}}=4$: F11BEF 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
$5$ or
$6$.
 2: U($*$) – 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$, $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
${\mathbf{IFAIL}}<{\mathbf{0}}$, the array
U is unchanged from the initial entry to F11BEF.
If
${\mathbf{IFAIL}}={\mathbf{1}}$, the array
U is unchanged from the last entry to F11BEF.
Otherwise,
U holds the last available iterate of the solution of the system of equations, for all returned values of
IFAIL.
 3: V($*$) – 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$, $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 F11BEF;
 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
${\mathbf{IFAIL}}<{\mathbf{0}}$, the array
V is unchanged from the initial entry to F11BEF.
If
${\mathbf{IFAIL}}={\mathbf{1}}$, the array
V is unchanged from the last entry to F11BEF.
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 available iterate of the residual vector unless
${\mathbf{IFAIL}}={\mathbf{8}}$ is returned on last entry, in which case
V is set to
$0.0$.
 4: WGT($*$) – 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 F11BDF).
Constraint:
if weights are to be used, at least one element of
WGT must be nonzero.
 5: WORK(LWORK) – REAL (KIND=nag_wp) arrayCommunication Array
On initial entry: the array
WORK as returned by
F11BDF (see also
Section 5 in F11BDF).
On intermediate reentry: must remain unchanged.
 6: LWORK – INTEGERInput
On initial entry: the dimension of the array
WORK as declared in the (sub)program from which F11BEF is called (see also
Sections 3 and
5 in F11BDF).
The required amount of workspace is as follows:
Method 
Requirements 
RGMRES 
${\mathbf{LWORK}}=100+\mathit{n}\left(m+3\right)+m\left(m+5\right)+1$, where $m$ is the dimension of the basis. 
CGS 
${\mathbf{LWORK}}=100+7\mathit{n}$. 
BiCGSTAB($\ell $) 
${\mathbf{LWORK}}=100+\left(2\mathit{n}+\ell \right)\left(\ell +2\right)+p$, where $\ell $ is the order of the method. 
TFQMR 
${\mathbf{LWORK}}=100+10\mathit{n}$, 
where
$p=2\mathit{n}$ 
if $\ell >1$ and ${\mathbf{ITERM}}=2$ was supplied. 
$p=\mathit{n}$ 
if $\ell >1$ and a preconditioner is used or ${\mathbf{ITERM}}=2$ was supplied. 
$p=0$ 
otherwise. 
Constraint:
${\mathbf{LWORK}}\ge {\mathbf{LWREQ}}$, where
LWREQ is returned by
F11BDF.
 7: IFAIL – INTEGERInput/Output
On initial entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. 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
$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 parameters 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 final exit:
${\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}}={\mathbf{0}}$ or
${{\mathbf{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}}=i$
On entry, the $i$th argument had an illegal value.
 ${\mathbf{IFAIL}}=1$
F11BEF has been called again after returning the termination code
${\mathbf{IREVCM}}=4$. No further computation has been carried out and all input data and data stored for access by
F11BFF have remained unchanged.
 ${\mathbf{IFAIL}}=2$
The required accuracy could not be obtained. However, F11BEF has terminated with reasonable accuracy: the last iterate of the residual satisfied the termination criterion but the exact residual
$r=bAx$, did not. After the first occurrence of this situation, the iteration was restarted once, but F11BEF could not improve on the accuracy. This error code usually implies that your problem has been fully and satisfactorily solved to within or close to the accuracy available on your system. Further iterations are unlikely to improve on this situation. You should call
F11BFF to check the values of the left and righthand sides of the termination condition.
 ${\mathbf{IFAIL}}=3$
F11BDF was either not called before calling F11BEF or it returned an error. The arguments
U and
V remain unchanged.
 ${\mathbf{IFAIL}}=4$
The calling program requested a tidy termination before the solution had converged. The arrays
U and
V return the last iterates available of the solution and of the residual vector, respectively.
 ${\mathbf{IFAIL}}=5$
The solution did not converge within the maximum number of iterations allowed. The arrays
U and
V return the last iterates available of the solution and of the residual vector, respectively.
 ${\mathbf{IFAIL}}=6$
Algorithmic breakdown. The last available iterates of the solution and residuals are returned, although it is possible that they are completely inaccurate.
 ${\mathbf{IFAIL}}=8$
The calling program requested an immediate termination. However, the array
U returns the last iterate of the solution, the array
V returns the last iterate of the residual vector, for the CGS and TFQMR methods only.
 ${\mathbf{IFAIL}}=10$
Usersupplied weights are to be used, but all elements of the array
WGT are zero.
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.
On successful completion, the termination criterion is satisfied to within the userspecified tolerance, as described in
Section 3 in F11BDF. The computed values of the left and righthand sides of the termination criterion selected can be obtained by a call to
F11BFF.
The number of operations carried out by F11BEF 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 F11BEF for each iteration is approximately proportional to $\mathit{n}$.
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}={M}^{1}A$ (RGMRES, CGS and TFQMR methods) or $\stackrel{}{A}=A{M}^{1}$ (BiCGSTAB($\ell $) method).
Additional matrixvector products are required for the computation of
${\Vert A\Vert}_{1}$ or
${\Vert A\Vert}_{\infty}$, when this has not been supplied to
F11BDF and is required by the termination criterion employed.
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 F11BDF) and
$\Vert {x}_{0}\Vert \gg \Vert {x}_{k}\Vert $, then the required accuracy cannot be obtained due to loss of significant digits. The iteration is restarted automatically at some suitable point: F11BEF sets
${x}_{0}={x}_{k}$ and the computation begins again. For particularly badly scaled problems, more than one restart may be necessary. This does not apply to the RGMRES method which, by its own nature, selfrestarts every superiteration. 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$.
9 Example