Integer type:** int32**** int64**** nag_int** show int32 show int32 show int64 show int64 show nag_int show nag_int

nag_sparse_complex_gen_basic_setup (f11br) is a setup function, the first in a suite of three functions for the iterative solution of a complex general (non-Hermitian) system of simultaneous linear equations. nag_sparse_complex_gen_basic_setup (f11br) must be called before nag_sparse_complex_gen_basic_solver (f11bs), the iterative solver. The third function in the suite, nag_sparse_complex_gen_basic_diag (f11bt), can be used to return additional information about the computation.

These three functions are suitable for the solution of large sparse general (non-Hermitian) systems of equations.

The suite consisting of the functions nag_sparse_complex_gen_basic_setup (f11br), nag_sparse_complex_gen_basic_solver (f11bs) and nag_sparse_complex_gen_basic_diag (f11bt) is designed to solve the general (non-Hermitian) system of simultaneous linear equations Ax = b$Ax=b$ of order n$n$, where n$n$ is large and the coefficient matrix A$A$ is sparse.

nag_sparse_complex_gen_basic_setup (f11br) is a setup function which must be called before nag_sparse_complex_gen_basic_solver (f11bs), the iterative solver. The third function in the suite, nag_sparse_complex_gen_basic_diag (f11bt), can be used to return additional information about the computation. A choice of methods is available:

- restarted generalized minimum residual method (RGMRES);
- conjugate gradient squared method (CGS);
- bi-conjugate gradient stabilized (ℓ$\ell $) method (Bi-CGSTAB(ℓ$\ell $));
- transpose-free quasi-minimal residual method (TFQMR).

The restarted generalized minimum residual method (RGMRES) (see Saad and Schultz (1986), Barrett et al. (1994) and Dias da Cunha and Hopkins (1994)) starts from the residual r_{0} = b − Ax_{0}${r}_{0}=b-A{x}_{0}$, where x_{0}${x}_{0}$ is an initial estimate for the solution (often x_{0} = 0${x}_{0}=0$). An orthogonal basis for the Krylov subspace span{A^{k}r_{0}}$\mathrm{span}\left\{{A}^{\mathit{k}}{r}_{0}\right\}$, for k = 0,1, … $\mathit{k}=0,1,\dots $, is generated explicitly: this is referred to as Arnoldi's method (see Arnoldi (1951)). The solution is then expanded onto the orthogonal basis so as to minimize the residual norm ‖b − Ax‖_{2}${\Vert b-Ax\Vert}_{2}$. The lack of symmetry of A$A$ implies that the orthogonal basis is generated by applying a ‘long’ recurrence relation, whose length increases linearly with the iteration count. For all but the most trivial problems, computational and storage costs can quickly become prohibitive as the iteration count increases. RGMRES limits these costs by employing a restart strategy: every m$m$ iterations at most, the Arnoldi process is restarted from r_{l} = b − Ax_{l}${r}_{l}=b-A{x}_{l}$, where the subscript l$l$ denotes the last available iterate. Each group of m$m$ iterations is referred to as a ‘super-iteration’. The value of m$m$ is chosen in advance and is fixed throughout the computation. Unfortunately, an optimum value of m$m$ cannot easily be predicted.

The conjugate gradient squared method (CGS) (see Sonneveld (1989), Barrett et al. (1994) and Dias da Cunha and Hopkins (1994)) is a development of the bi-conjugate gradient method where the nonsymmetric Lanczos method is applied to reduce the coefficients matrix to tridiagonal form: two bi-orthogonal sequences of vectors are generated starting from the residual r_{0} = b − Ax_{0}${r}_{0}=b-A{x}_{0}$, where x_{0}${x}_{0}$ is an initial estimate for the solution (often x_{0} = 0${x}_{0}=0$) and from the
shadow residual r̂_{0}${\hat{r}}_{0}$ corresponding to the arbitrary problem A^{H}x̂ = b̂${A}^{\mathrm{H}}\hat{x}=\hat{b}$, where b̂$\hat{b}$ can be any vector, but in practice is chosen so that r_{0} = r̂_{0}${r}_{0}={\hat{r}}_{0}$. In the course of the iteration, the residual and shadow residual r_{i} = P_{i}(A)
r_{0}${r}_{i}={P}_{i}\left(A\right){r}_{0}$ and r̂_{i} = P_{i}(A^{H})r̂_{0}${\hat{r}}_{i}={P}_{i}\left({A}^{\mathrm{H}}\right){\hat{r}}_{0}$ are generated, where P_{i}${P}_{i}$ is a polynomial of order i$i$, and bi-orthogonality is exploited by computing the vector product ρ_{i} = (r̂_{i},r_{i}) = (P_{i}(A^{H})r̂_{0},
P_{i}(A)
r_{0}) = (r̂_{0},P_{i}^{2}(A)r_{0})${\rho}_{i}=({\hat{r}}_{i},{r}_{i})=({P}_{i}\left({A}^{\mathrm{H}}\right){\hat{r}}_{0},{P}_{i}\left(A\right){r}_{0})=({\hat{r}}_{0},{P}_{i}^{2}\left(A\right){r}_{0})$. Applying the ‘contraction’ operator P_{i}(A)${P}_{i}\left(A\right)$ twice, the iteration coefficients can still be recovered without advancing the solution of the shadow problem, which is of no interest. The CGS method often provides fast convergence; however, there is no reason why the contraction operator should also reduce the once reduced vector P_{i}(A)r_{0}${P}_{i}\left(A\right){r}_{0}$: this may well lead to a highly irregular convergence which may result in large cancellation errors.

The bi-conjugate gradient stabilized (ℓ$\ell $) method
(Bi-CGSTAB(ℓ$\ell $)) (see Van der Vorst (1989), Sleijpen and Fokkema (1993) and Dias da Cunha and Hopkins (1994)) is similar to the CGS method above. However, instead of generating the sequence {P_{i}^{2}(A)r_{0}}$\left\{{P}_{i}^{2}\left(A\right){r}_{0}\right\}$, it generates the sequence {Q_{i}(A)P_{i}(A)r_{0}}$\left\{{Q}_{i}\left(A\right){P}_{i}\left(A\right){r}_{0}\right\}$, where the Q_{i}(A)${Q}_{i}\left(A\right)$ are polynomials chosen to minimize the residual after the application of the contraction operator P_{i}(A)${P}_{i}\left(A\right)$. Two main steps can be identified for each iteration: an OR (Orthogonal Residuals) step where a basis of order ℓ$\ell $ is generated by a Bi-CG iteration and an MR (Minimum Residuals) step where the residual is minimized over the basis generated, by a method akin to GMRES. For ℓ = 1$\ell =1$, the method corresponds to the Bi-CGSTAB method of Van der Vorst (1989). For ℓ > 1$\ell >1$, more information about complex eigenvalues of the iteration matrix can be taken into account, and this may lead to improved convergence and robustness. However, as ℓ$\ell $ increases, numerical instabilities may arise. For this reason, a maximum value of ℓ = 10$\ell =10$ is imposed, but probably ℓ = 4$\ell =4$ is sufficient in most cases.

The transpose-free quasi-minimal residual method (TFQMR) (see Freund and Nachtigal (1991) and Freund (1993)) is conceptually derived from the CGS method. The residual is minimized over the space of the residual vectors generated by the CGS iterations under the simplifying assumption that residuals are almost orthogonal. In practice, this is not the case but theoretical analysis has proved the validity of the method. This has the effect of remedying the rather irregular convergence behaviour with wild oscillations in the residual norm that can degrade the numerical performance and robustness of the CGS method. In general, the TFQMR method can be expected to converge at least as fast as the CGS method, in terms of number of iterations, although each iteration involves a higher operation count. When the CGS method exhibits irregular convergence, the TFQMR method can produce much smoother, almost monotonic convergence curves. However, the close relationship between the CGS and TFQMR method implies that the overall speed of convergence is similar for both methods. In some cases, the TFQMR method may converge faster than the CGS method.

For each method, a sequence of solution iterates {x_{i}}$\left\{{x}_{i}\right\}$ is generated such that, hopefully, the sequence of the residual norms {‖r_{i}‖}$\left\{\Vert {r}_{i}\Vert \right\}$ converges to a required tolerance. Note that, in general, convergence, when it occurs, is not monotonic.

In the RGMRES and Bi-CGSTAB(ℓ$\ell $) methods above, your program must provide the **maximum** number of basis vectors used, m$m$ or ℓ$\ell $, respectively; however, a **smaller** number of basis vectors may be generated and used when the stability of the solution process requires this (see Section [Further Comments]).

Faster convergence 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 coefficients matrix can be improved or eigenvalues in its spectrum can be made to coalesce. An orthogonal basis for the Krylov subspace span{A^{k}r_{0}}$\mathrm{span}\left\{{\stackrel{-}{A}}^{\mathit{k}}{\stackrel{-}{r}}_{0}\right\}$, for k = 0,1, … $\mathit{k}=0,1,\dots $, 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.

Ax = b,
$$\stackrel{-}{A}\stackrel{-}{x}=\stackrel{-}{b}\text{,}$$ | (1) |

A left preconditioner M^{ − 1}${M}^{-1}$ can be used by the RGMRES, CGS and TFQMR methods, such that A = M^{ − 1}A ∼ I_{n}$\stackrel{-}{A}={M}^{-1}A\sim {I}_{n}$ in (1), where I_{n}${I}_{n}$ is the identity matrix of order n$n$; a right preconditioner M^{ − 1}${M}^{-1}$ can be used by the Bi-CGSTAB(ℓ$\ell $) method, such that A = AM^{ − 1} ∼ I_{n}$\stackrel{-}{A}=A{M}^{-1}\sim {I}_{n}$. These are formal definitions, used only in the design of the algorithms; in practice, only the means to compute the matrix–vector products v = Au$v=Au$ and v = A^{H}u$v={A}^{\mathrm{H}}u$ (the latter only being required when an estimate of ‖A‖_{1}${\Vert A\Vert}_{1}$ or ‖A‖_{∞}${\Vert A\Vert}_{\infty}$ is computed internally), and to solve the preconditioning equations Mv = u$Mv=u$ are required, i.e., explicit information about M$M$, or its inverse is not required at any stage.

The first termination criterion

is available for all four methods. In (2), p = 1$p=1$, ∞ or 2$\infty \text{ or}2$ and τ$\tau $ denotes a user-specified tolerance subject to max (10,sqrt(n))$\mathrm{max}\phantom{\rule{0.125em}{0ex}}(10,\sqrt{n})$, ε ≤ τ < 1$\epsilon \le \tau <1$, where ε$\epsilon $ is the machine precision. Facilities are provided for the estimation of the norm of the coefficients matrix ‖A‖_{1}${\Vert A\Vert}_{1}$ or ‖A‖_{∞}${\Vert A\Vert}_{\infty}$, when this is not known in advance, by applying Higham's method (see Higham (1988)). Note that ‖A‖_{2}${\Vert A\Vert}_{2}$ 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 are not used because any such criteria would require information about the condition number κ(A)$\kappa \left(A\right)$ which is not easily obtainable.

‖r _{k}‖_{p}
≤
τ
(‖b‖_{p}
+
‖A‖_{p} × ‖x_{k}‖_{p})
$${\Vert {r}_{k}\Vert}_{p}\le \tau ({\Vert b\Vert}_{p}+{\Vert A\Vert}_{p}\times {\Vert {x}_{k}\Vert}_{p})$$ | (2) |

The second termination criterion

is available for all methods except TFQMR. In (3), σ_{1}(A) = ‖A‖_{2}${\sigma}_{1}\left(\stackrel{-}{A}\right)={\Vert \stackrel{-}{A}\Vert}_{2}$ is the largest singular value of the (preconditioned) iteration matrix A$\stackrel{-}{A}$. This termination criterion monitors the progress of the solution of the preconditioned system of equations and is less expensive to apply than criterion (2) for the Bi-CGSTAB(ℓ$\ell $) method with ℓ > 1$\ell >1$. Only the RGMRES method provides facilities to estimate σ_{1}(A)${\sigma}_{1}\left(\stackrel{-}{A}\right)$ internally, when this is not supplied (see Section [Further Comments]).

‖r _{k}‖_{2}
≤
τ
(
‖r_{0}‖_{2}
+
σ_{1}
(A)
× ‖Δx_{k}‖_{2})
$${\Vert {\stackrel{-}{r}}_{k}\Vert}_{2}\le \tau ({\Vert {\stackrel{-}{r}}_{0}\Vert}_{2}+{\sigma}_{1}\left(\stackrel{-}{A}\right)\times {\Vert \Delta {\stackrel{-}{x}}_{k}\Vert}_{2})$$ | (3) |

Termination criterion (2) is the recommended choice, despite its additional costs per iteration when using the Bi-CGSTAB(ℓ$\ell $) method with ℓ > 1$\ell >1$. Also, if the norm of the initial estimate is much larger than the norm of the solution, that is, if ‖x_{0}‖ ≫ ‖x‖$\Vert {x}_{0}\Vert \gg \Vert x\Vert $, 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 w$w$ of user-specified weights can be used in the computation of the vector norms in termination criterion (2), i.e., ‖v‖_{p}^{(w)} = ‖v^{(w)}‖_{p}${{\Vert v\Vert}_{p}}^{\left(w\right)}={\Vert {v}^{\left(w\right)}\Vert}_{p}$, where (v^{(w)})_{i} = w_{i} v_{i}${\left({v}^{\left(w\right)}\right)}_{\mathit{i}}={w}_{\mathit{i}}{v}_{\mathit{i}}$, for i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$. Note that the use of weights increases the computational costs.

The sequence of calls to the functions comprising the suite is enforced: first, the setup function nag_sparse_complex_gen_basic_setup (f11br) must be called, followed by the solver nag_sparse_complex_gen_basic_solver (f11bs). nag_sparse_complex_gen_basic_diag (f11bt) can be called either when nag_sparse_complex_gen_basic_solver (f11bs) is carrying out a monitoring step or after nag_sparse_complex_gen_basic_solver (f11bs) has completed its tasks. Incorrect sequencing will raise an error condition.

In general, it is not possible to recommend one method in preference to another. RGMRES is often used in the solution of systems arising from PDEs. On the other hand, it can easily stagnate when the size m$m$ of the orthogonal basis is too small, or the preconditioner is not good enough. CGS can be the fastest method, but the computed residuals can exhibit instability which may greatly affect the convergence and quality of the solution. Bi-CGSTAB(ℓ$\ell $) seems robust and reliable, but it can be slower than the other methods: if a preconditioner is used and ℓ > 1$\ell >1$, Bi-CGSTAB(ℓ$\ell $) computes the solution of the preconditioned system x_{k} = Mx_{k}${\stackrel{-}{x}}_{k}=M{x}_{k}$: the preconditioning equations must be solved to obtain the required solution. The algorithm employed limits to 10%$10\%$ or less, when no intermediate monitoring is requested, the number of times the preconditioner has to be thus applied compared with the total number of applications of the preconditioner. TFQMR can be viewed as a more robust variant of the CGS method: it shares the CGS method speed but avoids the CGS fluctuations in the residual, which may give rise to instability. Also, when the termination criterion (2) is used, the CGS, Bi-CGSTAB(ℓ$\ell $) and TFQMR methods will restart the iteration automatically when necessary in order to solve the given problem.

Arnoldi W (1951) The principle of minimized iterations in the solution of the matrix eigenvalue problem *Quart. Appl. Math.* **9** 17–29

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

Freund R W (1993) A transpose-free quasi-minimal residual algorithm for non-Hermitian linear systems *SIAM J. Sci. Comput.* **14** 470–482

Freund R W and Nachtigal N (1991) QMR: a Quasi-Minimal Residual Method for Non-Hermitian Linear Systems *Numer. Math.* **60** 315–339

Golub G H and Van Loan C F (1996) *Matrix Computations* (3rd Edition) Johns Hopkins University Press, Baltimore

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

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 Lanczos-type solver for nonsymmetric linear systems *SIAM J. Sci. Statist. Comput.* **10** 36–52

Van der Vorst H (1989) Bi-CGSTAB, a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems *SIAM J. Sci. Statist. Comput.* **13** 631–644

- 1: method – string
- The iterative method to be used.
- method = 'RGMRES'${\mathbf{method}}=\text{'RGMRES'}$
- Restarted generalized minimum residual method.
- method = 'CGS'${\mathbf{method}}=\text{'CGS'}$
- Conjugate gradient squared method.
- method = 'BICGSTAB'${\mathbf{method}}=\text{'BICGSTAB'}$
- Bi-conjugate gradient stabilized (ℓ$\ell $) method.
- method = 'TFQMR'${\mathbf{method}}=\text{'TFQMR'}$
- Transpose-free quasi-minimal residual method.

*Constraint*: method = 'RGMRES'${\mathbf{method}}=\text{'RGMRES'}$, 'CGS'$\text{'CGS'}$, 'BICGSTAB'$\text{'BICGSTAB'}$ or 'TFQMR'$\text{'TFQMR'}$. - 2: precon – string (length ≥ 1)
- Determines whether preconditioning is used.
- 3: n – int64int32nag_int scalar
- n$n$, the order of the matrix A$A$.
- 4: m – int64int32nag_int scalar
- If method = 'RGMRES'${\mathbf{method}}=\text{'RGMRES'}$, m is the dimension m$m$ of the restart subspace.If method = 'BICGSTAB'${\mathbf{method}}=\text{'BICGSTAB'}$, m is the order ℓ$\ell $ of the polynomial Bi-CGSTAB method.Otherwise, m is not referenced.
*Constraints*:- if method = 'RGMRES'${\mathbf{method}}=\text{'RGMRES'}$, 0 < m ≤ min (n,50)$0<{\mathbf{m}}\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}({\mathbf{n}},50)$;
- if method = 'BICGSTAB'${\mathbf{method}}=\text{'BICGSTAB'}$, 0 < m ≤ min (n,10)$0<{\mathbf{m}}\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}({\mathbf{n}},10)$.

- 5: tol – double scalar
- The tolerance τ$\tau $ for the termination criterion. If tol ≤ 0.0,τ = max (sqrt(ε),sqrt(n)ε)${\mathbf{tol}}\le 0.0,\tau =\mathrm{max}\phantom{\rule{0.125em}{0ex}}(\sqrt{\epsilon},\sqrt{n}\epsilon )$ is used, where ε$\epsilon $ is the machine precision. Otherwise τ = max (tol,10ε,sqrt(n)ε)$\tau =\mathrm{max}\phantom{\rule{0.125em}{0ex}}({\mathbf{tol}},10\epsilon ,\sqrt{n}\epsilon )$ is used.
- 6: maxitn – int64int32nag_int scalar
- The maximum number of iterations.
- 7: anorm – double scalar
- If anorm > 0.0${\mathbf{anorm}}>0.0$, the value of ‖A‖
_{p}${\Vert A\Vert}_{p}$ to be used in the termination criterion (2) (iterm = 1${\mathbf{iterm}}=1$).If anorm ≤ 0.0${\mathbf{anorm}}\le 0.0$, iterm = 1${\mathbf{iterm}}=1$ and norm = '1'${\mathbf{norm}}=\text{'1'}$ or 'I'$\text{'I'}$, then ‖A‖_{1}= ‖A‖_{∞}${\Vert A\Vert}_{1}={\Vert A\Vert}_{\infty}$ is estimated internally by nag_sparse_complex_gen_basic_solver (f11bs). - 8: sigmax – double scalar
- 9: monit – int64int32nag_int scalar
- If monit > 0${\mathbf{monit}}>0$, the frequency at which a monitoring step is executed by nag_sparse_complex_gen_basic_solver (f11bs): if method = 'CGS'${\mathbf{method}}=\text{'CGS'}$ or 'TFQMR'$\text{'TFQMR'}$, a monitoring step is executed every monit iterations; otherwise, a monitoring step is executed every monit super-iterations (groups of up to m$m$ or ℓ$\ell $ iterations for RGMRES or Bi-CGSTAB(ℓ$\ell $), respectively).There are some additional computational costs involved in monitoring the solution and residual vectors when the Bi-CGSTAB(ℓ$\ell $) method is used with ℓ > 1$\ell >1$.
- 10: lwork – int64int32nag_int scalar
- The dimension of the array work as declared in the (sub)program from which nag_sparse_complex_gen_basic_setup (f11br) is called.
**Note:**although the minimum value of lwork ensures the correct functioning of nag_sparse_complex_gen_basic_setup (f11br), a larger value is required by the other functions in the suite, namely nag_sparse_complex_gen_basic_solver (f11bs) and nag_sparse_complex_gen_basic_diag (f11bt). The required value is as follows:

where**Method****Requirements**RGMRES lwork = 120 + n(m + 3) + m(m + 5) + 1${\mathbf{lwork}}=120+n(m+3)+m(m+5)+1$, where m$m$ is the dimension of the basis. CGS lwork = 120 + 7n${\mathbf{lwork}}=120+7n$. Bi-CGSTAB(ℓ$\ell $) lwork = 120 + (2n + ℓ)(ℓ + 2) + p${\mathbf{lwork}}=120+(2n+\ell )(\ell +2)+p$, where ℓ$\ell $ is the order of the method. TFQMR lwork = 120 + 10n${\mathbf{lwork}}=120+10n$, p = 2n$p=2n$ if ℓ > 1$\ell >1$ and iterm = 2${\mathbf{iterm}}=2$ was supplied. p = n$p=n$ if ℓ > 1$\ell >1$ and a preconditioner is used or iterm = 2${\mathbf{iterm}}=2$ was supplied. p = 0$p=0$ otherwise.

- 1: norm_p – string (length ≥ 1)
- Defines the matrix and vector norm to be used in the termination criteria.
*Default*:- if iterm = 1${\mathbf{iterm}}=1$, 'I' $\text{'I'}$;
- otherwise '2' $\text{'2'}$.

- 2: weight – string (length ≥ 1)
- Specifies whether a vector w$w$ of user-supplied weights is to be used in the computation of the vector norms required in termination criterion (2) (iterm = 1${\mathbf{iterm}}=1$): ‖v‖
_{p}^{(w)}= ‖v^{(w)}‖_{p}${{\Vert v\Vert}_{p}}^{\left(w\right)}={\Vert {v}^{\left(w\right)}\Vert}_{p}$, where v_{i}^{(w)}= w_{i}v_{i}${v}_{\mathit{i}}^{\left(w\right)}={w}_{\mathit{i}}{v}_{\mathit{i}}$, for i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$. The suffix p = 1,2,∞$p=1,2,\infty $ denotes the vector norm used, as specified by the parameter norm_p. Note that weights cannot be used when iterm = 2${\mathbf{iterm}}=2$, i.e., when criterion (3) is used.- weight = 'W'${\mathbf{weight}}=\text{'W'}$
- User-supplied weights are to be used and must be supplied on initial entry to nag_sparse_complex_gen_basic_solver (f11bs).
- weight = 'N'${\mathbf{weight}}=\text{'N'}$
- All weights are implicitly set equal to one. Weights do not need to be supplied on initial entry to nag_sparse_complex_gen_basic_solver (f11bs).

*Default*: 'N'$\text{'N'}$ - 3: iterm – int64int32nag_int scalar
- Defines the termination criterion to be used.
*Default*: 1$1$*Constraints*:

None.

- 1: lwreq – int64int32nag_int scalar
- The minimum amount of workspace required by nag_sparse_complex_gen_basic_solver (f11bs). (See also Section [Parameters] in (f11bs).)
- 2: work(lwork) – complex array
- The array work is initialized by nag_sparse_complex_gen_basic_setup (f11br). It must
**not**be modified before calling the next function in the suite, namely nag_sparse_complex_gen_basic_solver (f11bs). - 3: ifail – int64int32nag_int scalar
- ifail = 0${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Errors or warnings detected by the function:

- If ifail = − i${\mathbf{ifail}}=-i$, parameter i$i$ had an illegal value on entry. The parameters are numbered as follows:

- nag_sparse_complex_gen_basic_setup (f11br) has been called out of sequence.

Not applicable.

RGMRES can estimate internally the maximum singular value σ_{1}${\sigma}_{1}$ of the iteration matrix, using σ_{1} ∼ ‖T‖_{1}${\sigma}_{1}\sim {\Vert T\Vert}_{1}$, where T$T$ is the upper triangular matrix obtained by QR$QR$ factorization of the upper Hessenberg matrix generated by the Arnoldi process. The computational costs of this computation are negligible when compared to the overall costs.

Loss of orthogonality in the RGMRES method, or of bi-orthogonality in the Bi-CGSTAB(ℓ$\ell $) method may degrade the solution and speed of convergence. For both methods, the algorithms employed include checks on the basis vectors so that the number of basis vectors used for a given super-iteration may be less than the value specified in the input parameter m. Also, if termination criterion (2) is used, the CGS, Bi-CGSTAB(ℓ$\ell $) and TFQMR methods will restart automatically the computation from the last available iterates, when the stability of the solution process requires it.

Termination criterion (3), when available, involves only the residual (or norm of the residual) produced directly by the iteration process: this may differ from the norm of the true residual r̃_{k} = b − Ax_{k}${\stackrel{~}{r}}_{k}=b-A{x}_{k}$, particularly when the norm of the residual is very small. Also, if the norm of the initial estimate of the solution is much larger than the norm of the exact solution, convergence can be achieved despite very large errors in the solution. On the other hand, termination criterion (3) is cheaper to use and inspects the progress of the actual iteration. Termination criterion (2) should be preferred in most cases, despite its slightly larger costs.

Open in the MATLAB editor: nag_sparse_complex_gen_basic_setup_example

function nag_sparse_complex_gen_basic_setup_examplemethod = 'TFQMR '; precon = 'P'; n = 8; m = int64(1); tol = 1e-08; maxitn = int64(20); anorm = 0; sigmax = 0; monit = int64(2); lwork = int64(6000); nz=24; a=complex(zeros(10000,1)); a(1:24)=[ 2 + 1.i, -1 + 1.i, 1 - 3.i, 4 + 7.i, -3 + 0.i, 2 + 4.i, -7 - 5.i, 2 + 1.i, 3 + 2.i, -4 + 2.i, 0 + 1.i, 5 - 3.i, -1 + 2.i, 8 + 6.i, -3 - 4.i, -6 - 2.i, 5 - 2.i, 2 + 0.i, 0 - 5.i, -1 + 5.i, 6 + 2.i, -1 + 4.i, 2 + 0.i, 3 + 3.i]; irow=zeros(10000,1,'int64'); irow(1:24) = [int64(1),1,1,2,2,2,3,3,4,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8]; icol=zeros(10000,1,'int64'); icol(1:24) = [int64(1),4,8,1,2,5,3,6,1,3,4,7,2,5,7,1,3,6,3,5,7,2,6,8]; ipivp = zeros(n,1,'int64'); ipivq = zeros(n,1,'int64'); irevcm = int64(0); lfill = int64(0); dtol = 0; milu = 'N'; wgt = zeros(8,1); u = complex(zeros(8,1)); v = [ 7 + 11.i; 1 + 24.i; -13 - 18.i; -10 + 3.i; 23 + 14.i; 17 - 7.i; 15 - 3.i; -3 + 20.i]; [a, irow, icol, ipivp, ipivq, istr, idiag, nnzc, npivm, ifail] = ... nag_sparse_complex_gen_precon_ilu(int64(nz), a, irow, icol, lfill, ... dtol, milu, ipivp, ipivq); [lwreq, work, ifail] = ... nag_sparse_complex_gen_basic_setup(method, precon, int64(n), m, tol, ... maxitn, anorm, sigmax, monit, lwork, 'norm_p', '1'); while (irevcm ~= 4) [irevcm, u, v, work, ifail] = ... nag_sparse_complex_gen_basic_solver(irevcm, u, v, wgt, work); if (irevcm == -1) [v, ifail] = ... nag_sparse_complex_gen_matvec('T', a(1:nz), irow(1:nz), icol(1:nz), ... 'N', u); elseif (irevcm == 1) [v, ifail] = ... nag_sparse_complex_gen_matvec('N', a(1:nz), irow(1:nz), icol(1:nz), ... 'N', u); elseif (irevcm == 2) [v, ifail] = ... nag_sparse_complex_gen_precon_ilu_solve('N', a, irow, icol, ipivp, ipivq, istr, idiag, 'N', u); elseif (irevcm == 3) [itn, stplhs, stprhs, anorm, sigmax, work, ifail] = ... nag_sparse_complex_gen_basic_diag(work); fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', ... itn, stplhs); fprintf('\n Solution Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(u(i)), imag(u(i))); end fprintf('\n Residual Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(v(i)), imag(v(i))); end end end % Get information about the computation [itn, stplhs, stprhs, anorm, sigmax, work, ifail] = ... nag_sparse_complex_gen_basic_diag(work); fprintf('\nNumber of iterations for convergence: %d\n', itn); fprintf('Residual norm: %14.4e\n', stplhs); fprintf('Right-hand side of termination criteria: %14.4e\n', stprhs); fprintf('i-norm of matrix a: %14.4e\n', anorm); fprintf('\n Solution Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(u(i)), imag(u(i))); end fprintf('\n Residual Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(v(i)), imag(v(i))); end

Monitoring at iteration number 2 residual norm: 8.2345e+01 Solution Vector +6.9055e-01 + +1.4236e+00I +7.3931e-02 + -1.1880e+00I +1.4778e+00 + +4.7846e-01I +5.6572e+00 + -3.0786e+00I +1.4243e+00 + -1.1246e+00I +1.0374e-01 + +1.9740e+00I +4.4985e-01 + -1.2715e+00I +2.5704e+00 + +1.7578e+00I Residual Vector +1.7772e+00 + +4.6797e+00I +1.0774e+00 + +6.4600e+00I -3.2812e+00 + -1.1314e+01I -3.8698e+00 + -1.6438e+00I +8.9912e+00 + +1.1100e+01I +9.7428e+00 + -4.6218e-01I +3.1668e+00 + +2.8721e+00I -1.0323e+01 + +1.5837e+00I Number of iterations for convergence: 4 Residual norm: 5.6586e-12 Right-hand side of termination criteria: 8.9100e-06 i-norm of matrix a: 2.7000e+01 Solution Vector +1.0000e+00 + +1.0000e+00I +2.0000e+00 + -1.0000e+00I +3.0000e+00 + +1.0000e+00I +4.0000e+00 + -1.0000e+00I +3.0000e+00 + -1.0000e+00I +2.0000e+00 + +1.0000e+00I +1.0000e+00 + -1.0000e+00I -1.8161e-13 + +3.0000e+00I Residual Vector -4.5297e-14 + -3.4817e-13I +4.1744e-13 + -5.4712e-13I +2.5402e-13 + +5.1870e-13I +3.5705e-13 + +1.0658e-14I -5.3291e-13 + -4.6185e-13I -4.6541e-13 + +1.0658e-13I -4.6896e-13 + +4.8317e-13I +6.1284e-13 + -2.8422e-14I

Open in the MATLAB editor: f11br_example

function f11br_examplemethod = 'TFQMR '; precon = 'P'; n = 8; m = int64(1); tol = 1e-08; maxitn = int64(20); anorm = 0; sigmax = 0; monit = int64(2); lwork = int64(6000); nz=24; a=complex(zeros(10000,1)); a(1:24)=[ 2 + 1.i, -1 + 1.i, 1 - 3.i, 4 + 7.i, -3 + 0.i, 2 + 4.i, -7 - 5.i, 2 + 1.i, 3 + 2.i, -4 + 2.i, 0 + 1.i, 5 - 3.i, -1 + 2.i, 8 + 6.i, -3 - 4.i, -6 - 2.i, 5 - 2.i, 2 + 0.i, 0 - 5.i, -1 + 5.i, 6 + 2.i, -1 + 4.i, 2 + 0.i, 3 + 3.i]; irow=zeros(10000,1,'int64'); irow(1:24) = [int64(1),1,1,2,2,2,3,3,4,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8]; icol=zeros(10000,1,'int64'); icol(1:24) = [int64(1),4,8,1,2,5,3,6,1,3,4,7,2,5,7,1,3,6,3,5,7,2,6,8]; ipivp = zeros(n,1,'int64'); ipivq = zeros(n,1,'int64'); irevcm = int64(0); lfill = int64(0); dtol = 0; milu = 'N'; wgt = zeros(8,1); u = complex(zeros(8,1)); v = [ 7 + 11.i; 1 + 24.i; -13 - 18.i; -10 + 3.i; 23 + 14.i; 17 - 7.i; 15 - 3.i; -3 + 20.i]; [a, irow, icol, ipivp, ipivq, istr, idiag, nnzc, npivm, ifail] = ... f11dn(int64(nz), a, irow, icol, lfill, dtol, milu, ipivp, ipivq); [lwreq, work, ifail] = ... f11br(method, precon, int64(n), m, tol, maxitn, anorm, sigmax, monit, lwork, 'norm_p', '1'); while (irevcm ~= 4) [irevcm, u, v, work, ifail] = f11bs(irevcm, u, v, wgt, work); if (irevcm == -1) [v, ifail] = f11xn('T', a(1:nz), irow(1:nz), icol(1:nz), 'N', u); elseif (irevcm == 1) [v, ifail] = f11xn('N', a(1:nz), irow(1:nz), icol(1:nz), 'N', u); elseif (irevcm == 2) [v, ifail] = f11dp('N', a, irow, icol, ipivp, ipivq, istr, idiag, 'N', u); elseif (irevcm == 3) [itn, stplhs, stprhs, anorm, sigmax, work, ifail] = f11bt(work); fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', itn, stplhs); fprintf('\n Solution Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(u(i)), imag(u(i))); end fprintf('\n Residual Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(v(i)), imag(v(i))); end end end % Get information about the computation [itn, stplhs, stprhs, anorm, sigmax, work, ifail] = f11bt(work); fprintf('\nNumber of iterations for convergence: %d\n', itn); fprintf('Residual norm: %14.4e\n', stplhs); fprintf('Right-hand side of termination criteria: %14.4e\n', stprhs); fprintf('i-norm of matrix a: %14.4e\n', anorm); fprintf('\n Solution Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(u(i)), imag(u(i))); end fprintf('\n Residual Vector\n'); for i = 1:n fprintf('%+16.4e + %+16.4eI\n', real(v(i)), imag(v(i))); end

Monitoring at iteration number 2 residual norm: 8.2345e+01 Solution Vector +6.9055e-01 + +1.4236e+00I +7.3931e-02 + -1.1880e+00I +1.4778e+00 + +4.7846e-01I +5.6572e+00 + -3.0786e+00I +1.4243e+00 + -1.1246e+00I +1.0374e-01 + +1.9740e+00I +4.4985e-01 + -1.2715e+00I +2.5704e+00 + +1.7578e+00I Residual Vector +1.7772e+00 + +4.6797e+00I +1.0774e+00 + +6.4600e+00I -3.2812e+00 + -1.1314e+01I -3.8698e+00 + -1.6438e+00I +8.9912e+00 + +1.1100e+01I +9.7428e+00 + -4.6218e-01I +3.1668e+00 + +2.8721e+00I -1.0323e+01 + +1.5837e+00I Number of iterations for convergence: 4 Residual norm: 1.8941e-11 Right-hand side of termination criteria: 8.9100e-06 i-norm of matrix a: 2.7000e+01 Solution Vector +1.0000e+00 + +1.0000e+00I +2.0000e+00 + -1.0000e+00I +3.0000e+00 + +1.0000e+00I +4.0000e+00 + -1.0000e+00I +3.0000e+00 + -1.0000e+00I +2.0000e+00 + +1.0000e+00I +1.0000e+00 + -1.0000e+00I +3.4360e-13 + +3.0000e+00I Residual Vector +5.2047e-13 + -1.0925e-12I +2.8848e-12 + -4.2988e-13I +5.2935e-13 + +1.8510e-12I +1.4637e-12 + +6.0574e-13I -1.4779e-12 + -1.8439e-12I -1.6023e-12 + +7.1054e-14I -1.7515e-12 + +3.0198e-14I +1.5401e-12 + +1.2470e-12I

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013