PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_sparse_complex_herm_basic_solver (f11gs)
Purpose
nag_sparse_complex_herm_basic_solver (f11gs) is an iterative solver for a complex Hermitian system of simultaneous linear equations;
nag_sparse_complex_herm_basic_solver (f11gs) is the second in a suite of three functions, where the first function,
nag_sparse_complex_herm_basic_setup (f11gr), must be called prior to
nag_sparse_complex_herm_basic_solver (f11gs) to set up the suite, and the third function in the suite,
nag_sparse_complex_herm_basic_diag (f11gt), can be used to return additional information about the computation.
These three functions are suitable for the solution of large sparse complex Hermitian systems of equations.
Syntax
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: lwork has been removed from the interface
.
Description
nag_sparse_complex_herm_basic_solver (f11gs) solves the complex Hermitian system of linear simultaneous equations
Ax = b$Ax=b$ using either 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)) or a preconditioned Lanczos method based upon the algorithm SYMMLQ (see
Paige and Saunders (1975) and
Barrett et al. (1994)).
For a general description of the methods employed you are referred to
Section [Description] in
(f11gr).
nag_sparse_complex_herm_basic_solver (f11gs) can solve the system after the first function in the suite,
nag_sparse_complex_herm_basic_setup (f11gr), has been called to initialize the computation and specify the method of solution. The third function in the suite,
nag_sparse_complex_herm_basic_diag (f11gt), can be used to return additional information generated by the computation during monitoring steps and after
nag_sparse_complex_herm_basic_solver (f11gs) has completed its tasks.
nag_sparse_complex_herm_basic_solver (f11gs) uses
reverse communication, i.e.,
nag_sparse_complex_herm_basic_solver (f11gs) returns repeatedly to the calling program with the parameter
irevcm (see
Section [Parameters]) set to specified values which require the calling program to carry out a specific task: either to compute the matrixvector product
v = Au$v=Au$; to solve the preconditioning equation
Mv = u$Mv=u$; to notify the completion of the computation; or, to 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 function, 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 functions used to perform matrix operations.
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
Parameters
Note: this function 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.
Compulsory Input Parameters
 1:
irevcm – int64int32nag_int scalar
On initial entry:
irevcm = 0${\mathbf{irevcm}}=0$, otherwise an error condition will be raised.
On intermediate reentry:
irevcm must either be unchanged from its previous exit value, or can have one of the following values.
 irevcm = 5${\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. nag_sparse_complex_herm_basic_solver (f11gs) will then return with the termination code irevcm = 4${\mathbf{irevcm}}=4$. Note that before calling nag_sparse_complex_herm_basic_solver (f11gs) with irevcm = 5${\mathbf{irevcm}}=5$ the calling program must have performed the tasks required by the value of irevcm returned by the previous call to nag_sparse_complex_herm_basic_solver (f11gs), otherwise subsequently returned values may be invalid.
 irevcm = 6${\mathbf{irevcm}}=6$
 Immediate termination: nag_sparse_complex_herm_basic_solver (f11gs) will return immediately with termination code irevcm = 4${\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. nag_sparse_complex_herm_basic_solver (f11gs) will then return with the termination code irevcm = 4${\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.
Constraint:
on initial entry,
irevcm = 0${\mathbf{irevcm}}=0$; on reentry, either
irevcm must remain unchanged or be reset to
irevcm = 5${\mathbf{irevcm}}=5$ or
6$6$.
 2:
u( : $:$) – complex array
Note: the dimension of the array
u
must be at least
n$\mathit{n}$.
On initial entry: an initial estimate, x_{0}${x}_{0}$, of the solution of the system of equations Ax = b$Ax=b$.
On intermediate reentry: must remain unchanged.
 3:
v( : $:$) – complex array

Note: the dimension of the array
v
must be at least
n$\mathit{n}$.
On initial entry: the righthand side b$b$ of the system of equations Ax = b$Ax=b$.
On intermediate reentry: the returned value of
irevcm determines the contents of
v in the following way.
If
irevcm = 1${\mathbf{irevcm}}=1$ or
2$2$,
v must store the vector
v$v$, the result of the operation specified by the value of
irevcm returned by the previous call to
nag_sparse_complex_herm_basic_solver (f11gs)
If
irevcm = 3${\mathbf{irevcm}}=3$,
v must remain unchanged.
 4:
wgt( : $:$) – double array

Note: the dimension of the array
wgt
must be at least
max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}(1,\mathit{n})$.
The usersupplied weights, if these are to be used in the computation of the vector norms in the termination criterion (see
Sections [Description] and
[Parameters] in
(f11gr)).
Constraint:
if weights are to be used, at least one element of
wgt must be nonzero.
 5:
work(lwork) – complex array
lwork, the dimension of the array, must satisfy the constraint
lwork ≥ lwreq$\mathit{lwork}\ge {\mathbf{lwreq}}$, where
lwreq is returned by
nag_sparse_complex_herm_basic_setup (f11gr).
On initial entry: the array
work as returned by
nag_sparse_complex_herm_basic_setup (f11gr) (see also
Section [Parameters] in
(f11gr)).
lwork, the dimension of the array, must satisfy the constraint
lwork ≥ lwreq$\mathit{lwork}\ge {\mathbf{lwreq}}$, where
lwreq is returned by
nag_sparse_complex_herm_basic_setup (f11gr).
On intermediate reentry: must remain unchanged.
Optional Input Parameters
None.
Input Parameters Omitted from the MATLAB Interface
 lwork
Output Parameters
 1:
irevcm – int64int32nag_int scalar
On intermediate exit:
has the following meanings.
 irevcm = 1${\mathbf{irevcm}}=1$
 The calling program must compute the matrixvector product v = Au$v=Au$, where u$u$ and v$v$ are stored in u and v, respectively.
 irevcm = 2${\mathbf{irevcm}}=2$
 The calling program must solve the preconditioning equation Mv = u$Mv=u$, where u$u$ and v$v$ are stored in u and v, respectively.
 irevcm = 3${\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 nag_sparse_complex_herm_basic_diag (f11gt) can be called at this step.
On final exit: if
irevcm = 4${\mathbf{irevcm}}=4$,
nag_sparse_complex_herm_basic_solver (f11gs) 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.
 2:
u( : $:$) – complex array
Note: the dimension of the array
u
must be at least
n$\mathit{n}$.
On intermediate exit:
the returned value of
irevcm determines the contents of
u in the following way.
If
irevcm = 1${\mathbf{irevcm}}=1$ or
2$2$,
u holds the vector
u$u$ on which the operation specified by
irevcm is to be carried out.
If
irevcm = 3${\mathbf{irevcm}}=3$,
u holds the current iterate of the solution vector.
On final exit: if
ifail = 3${\mathbf{ifail}}={\mathbf{3}}$ or
− i${{\mathit{i}}}$, the array
u is unchanged from the initial entry to
nag_sparse_complex_herm_basic_solver (f11gs). If
ifail = 1${\mathbf{ifail}}={\mathbf{1}}$, the array
u is unchanged from the last entry to
nag_sparse_complex_herm_basic_solver (f11gs). Otherwise,
u holds the last iterate of the solution of the system of equations, for all returned values of
ifail.
 3:
v( : $:$) – complex array

Note: the dimension of the array
v
must be at least
n$\mathit{n}$.
On intermediate exit:
if
irevcm = 3${\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
ifail = 3${\mathbf{ifail}}={\mathbf{3}}$ or
0${\mathbf{0}}$, the array
v is unchanged from the last entry to
nag_sparse_complex_herm_basic_solver (f11gs). If
ifail = 1${\mathbf{ifail}}={\mathbf{1}}$, the array
v is unchanged from the initial entry to
nag_sparse_complex_herm_basic_solver (f11gs). If
ifail = 0${\mathbf{ifail}}={\mathbf{0}}$ or
2${\mathbf{2}}$, the array
v contains the true residual vector of the system of equations (see also
Section [Error Indicators and Warnings]). Otherwise,
v stores the last iterate of the residual vector unless the Lanczos method (SYMMLQ) was used and
ifail ≥ 5${\mathbf{ifail}}\ge {\mathbf{5}}$, in which case
v is set to
0.0$0.0$.
 4:
work(lwork) – complex array
 5:
ifail – int64int32nag_int scalar
ifail = 0${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see
[Error Indicators and Warnings]).
Error Indicators and Warnings
Errors or warnings detected by the function:
Cases prefixed with W are classified as warnings and
do not generate an error of type NAG:error_n. See nag_issue_warnings.
 ifail = − i${\mathbf{ifail}}=i$
If
ifail = − i${\mathbf{ifail}}=i$, parameter
i$i$ had an illegal value on entry. The parameters are numbered as follows:
1:
irevcm, 2:
u, 3:
v, 4:
wgt, 5:
work, 6:
lwork, 7:
ifail.
It is possible that
ifail refers to a parameter that is omitted from the MATLAB interface. This usually indicates that an error in one of the other input parameters has caused an incorrect value to be inferred.
 W ifail = 1${\mathbf{ifail}}=1$
nag_sparse_complex_herm_basic_solver (f11gs) has been called again after returning the termination code
irevcm = 4${\mathbf{irevcm}}=4$. No further computation has been carried out and all input data and data stored for access by
nag_sparse_complex_herm_basic_diag (f11gt) have remained unchanged.
 W ifail = 2${\mathbf{ifail}}=2$
The required accuracy could not be obtained. However,
nag_sparse_complex_herm_basic_solver (f11gs) has terminated with reasonable accuracy: the last iterate of the residual satisfied the termination criterion but the exact residual
r = b − Ax$r=bAx$, did not. A small number of iterations have been carried out after the iterated residual satisfied the termination criterion, but were unable to 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
nag_sparse_complex_herm_basic_diag (f11gt) to check the values of the left and righthand side of the termination condition.
 ifail = 3${\mathbf{ifail}}=3$
nag_sparse_complex_herm_basic_setup (f11gr) was either not called before calling
nag_sparse_complex_herm_basic_solver (f11gs) or it returned an error. The arrays
u and
v remain unchanged.
 W ifail = 4${\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.
 W ifail = 5${\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.
 ifail = 6${\mathbf{ifail}}=6$
The preconditioner appears not to be positive definite. It is likely that your results are meaningless: both methods require a positive definite preconditioner (see also
Section [Description]). However, the array
u returns the last iterate of the solution, the array
v returns the last iterate of the residual vector, for the conjugate gradient method only.
 W ifail = 7${\mathbf{ifail}}=7$
The matrix of the coefficients appears not to be positive definite (conjugate gradient method only). The arrays
u and
v return the last iterates of the solution and residual vector, respectively. However, you should be warned that the results returned can be be in error.
 W ifail = 8${\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 conjugate gradient method only.
 ifail = 9${\mathbf{ifail}}=9$
The matrix appears to be singular.
 ifail = 10${\mathbf{ifail}}=10$
Usersupplied weights are to be used, but all elements of the array
wgt are zero.
Accuracy
On completion, i.e.,
irevcm = 4${\mathbf{irevcm}}=4$ on exit, the arrays
u and
v will return the solution and residual vectors,
x_{k}${x}_{k}$ and
r_{k} = b − Ax_{k}${r}_{k}=bA{x}_{k}$, respectively, at the
k$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 [Description] in
(f11gr). The computed values of the left and righthand sides of the termination criterion selected can be obtained by a call to
nag_sparse_complex_herm_basic_diag (f11gt).
Further Comments
The number of operations carried out by nag_sparse_complex_herm_basic_solver (f11gs) for each iteration is likely to be principally determined by the computation of the matrixvector products v = Au$v=Au$ and by the solution of the preconditioning equation Mv = u$Mv=u$ in the calling program. Each of these operations is carried out once every iteration.
The number of the remaining operations in nag_sparse_complex_herm_basic_solver (f11gs) for each iteration is approximately proportional to n$\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 A = E^{ − 1}AE^{ − H}$\stackrel{}{A}={E}^{1}A{E}^{\mathrm{H}}$.
Additional matrixvector products are required for the computation of
‖A‖_{1} = ‖A‖_{∞}${\Vert A\Vert}_{1}={\Vert A\Vert}_{\infty}$, when this has not been supplied to
nag_sparse_complex_herm_basic_setup (f11gr) and is required by the termination criterion employed.
The number of operations required to compute
σ_{1}(A)${\sigma}_{1}\left(\stackrel{}{A}\right)$ is negligible for reasonable values of
sigtol and
maxits (see
Sections [Parameters] and
[Further Comments] in
(f11gr)).
If the termination criterion
‖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})$ is used (see
Section [Description] in
(f11gr)) and
‖x_{0}‖ ≫ ‖x_{k}‖$\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:
nag_sparse_complex_herm_basic_solver (f11gs) sets
x_{0} = x_{k}${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}${x}_{0}$ which is as close to the true solution
x̃$\stackrel{~}{x}$ as can be estimated. Otherwise, the iteration should start from
x_{0} = 0${x}_{0}=0$.
Example
Open in the MATLAB editor:
nag_sparse_complex_herm_basic_solver_example
function nag_sparse_complex_herm_basic_solver_example
nz = 23;
a = zeros(10000,1);
a(1:nz) = [ 6 + 0.i;
1 + 1.i;
6 + 0.i;
0 + 1.i;
5 + 0.i;
5 + 0.i;
2  2.i;
4 + 0.i;
1 + 1.i;
2 + 0.i;
6 + 0.i;
4 + 3.i;
0 + 1.i;
1 + 0.i;
6 + 0.i;
1  1.i;
0  1.i;
9 + 0.i;
1 + 3.i;
1 + 2.i;
1 + 0.i;
1 + 4.i;
9 + 0.i];
irow = zeros(10000, 1, 'int64');
irow(1:nz) = [int64(1);2;2;3;3;4;5;5;6;6;6;7;7;7;7;8;8;8;9;9;9;9;9];
icol = zeros(10000, 1, 'int64');
icol(1:nz) = [int64(1);1;2;2;3;4;1;5;3;4;6;2;5;6;7;4;6;8;1;5;6;8;9];
lfill = int64(0);
dtol = 0;
mic = 'N';
dscale = 0;
ipiv = [int64(0);0;0;0;0;0;0;0;0];
method = 'CG ';
precon = 'P';
n = 9;
tol = 1e06;
maxitn = int64(20);
anorm = 0;
sigmax = 0;
maxits = int64(9);
monit = int64(2);
irevcm = int64(0);
u = [complex(0);
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i];
v = [ 8 + 54i;
10  92i;
25 + 27i;
26  28i;
54 + 12i;
26  22i;
47 + 65i;
71  57i;
60 + 70i];
wgt = [0; 0; 0; 0; 0; 0; 0; 0; 0];
[a, irow, icol, ipiv, istr, nnzc, npivm, ifail] = ...
nag_sparse_complex_herm_precon_ilu(int64(nz), a, irow, icol, lfill, ...
dtol, mic, dscale, ipiv);
[lwreq, work, ifail] = ...
nag_sparse_complex_herm_basic_setup(method, precon, int64(n), tol, ...
maxitn, anorm, sigmax, maxits, monit, 'sigcmp', 's', 'norm_p', '1');
while (irevcm ~= 4)
[irevcm, u, v, work, ifail] = ...
nag_sparse_complex_herm_basic_solver(irevcm, u, v, wgt, work);
if (irevcm == 1)
[v, ifail] = ...
nag_sparse_complex_herm_matvec(a(1:nz), irow(1:nz), icol(1:nz), 'N', u);
if (ifail ~= 0)
irecvm = 6;
end
elseif (irevcm == 2)
[v, ifail] = ...
nag_sparse_complex_herm_precon_ilu_solve(a, irow, icol, ipiv, istr, 'N', u);
if (ifail ~= 0)
irecvm = 6;
end
elseif (irevcm == 3)
[itn, stplhs, stprhs, anorm, sigmax, its, sigerr, ifail] = ...
nag_sparse_complex_herm_basic_diag(work);
if (ifail ~= 0)
irecvm = 6;
end
fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', ...
itn, stplhs);
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
[itn, stplhs, stprhs, anorm, sigmax, its, sigerr, ifail] = ...
nag_sparse_complex_herm_basic_diag(work);
fprintf('\nNumber of iterations for convergence: %d\n', itn);
fprintf('Residual norm: %14.4e\n', stplhs);
fprintf('Righthand side of termination criteria: %14.4e\n', stprhs);
fprintf('inorm 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: 1.4937e+01
+2.1423e01 + +4.5333e+00I
1.6589e+00 + 1.2672e+01I
+2.4101e+00 + +7.4551e+00I
+4.4400e+00 + 6.4174e+00I
+9.1135e+00 + +3.7812e+00I
+4.4419e+00 + 4.0382e+00I
+1.4757e+00 + +1.2662e+00I
+8.4872e+00 + 3.5347e+00I
+5.9948e+00 + +9.6851e01I
Residual Vector
1.8370e+00 + +3.6956e+00I
6.5005e01 + +2.5458e01I
1.2616e01 + 1.3625e01I
1.3120e01 + +1.4130e01I
1.1471e+00 + +7.3386e01I
5.5054e01 + 1.0535e+00I
+1.7165e+00 + 1.4614e+00I
3.5829e01 + +2.8764e01I
3.0278e01 + 3.5324e01I
Monitoring at iteration number 4
residual norm: 1.4602e+00
+1.0061e+00 + +8.9847e+00I
+1.9637e+00 + 7.9768e+00I
+3.0067e+00 + +7.0285e+00I
+3.9830e+00 + 5.9636e+00I
+5.0390e+00 + +5.0432e+00I
+6.0488e+00 + 4.0771e+00I
+6.9710e+00 + +3.0168e+00I
+8.0118e+00 + 1.9806e+00I
+9.0074e+00 + +9.6458e01I
Residual Vector
+1.1524e02 + 2.8188e02I
+1.3513e02 + 1.7345e01I
+1.8173e02 + +1.9627e02I
+1.8900e02 + 2.0354e02I
9.0877e02 + 1.0895e01I
2.3890e01 + +3.2440e01I
+1.9031e01 + 1.5499e02I
+5.1611e02 + 4.1435e02I
+4.3615e02 + +5.0884e02I
Number of iterations for convergence: 5
Residual norm: 7.1054e14
Righthand side of termination criteria: 2.7340e03
inorm of matrix a: 2.2000e+01
Solution Vector
+1.0000e+00 + +9.0000e+00I
+2.0000e+00 + 8.0000e+00I
+3.0000e+00 + +7.0000e+00I
+4.0000e+00 + 6.0000e+00I
+5.0000e+00 + +5.0000e+00I
+6.0000e+00 + 4.0000e+00I
+7.0000e+00 + +3.0000e+00I
+8.0000e+00 + 2.0000e+00I
+9.0000e+00 + +1.0000e+00I
Residual Vector
+1.0658e14 + +7.1054e15I
+0.0000e+00 + +0.0000e+00I
+0.0000e+00 + 3.5527e15I
+3.5527e15 + +3.5527e15I
+0.0000e+00 + +1.0658e14I
+0.0000e+00 + 1.0658e14I
+1.4211e14 + +0.0000e+00I
+0.0000e+00 + +0.0000e+00I
7.1054e15 + +0.0000e+00I
Open in the MATLAB editor:
f11gs_example
function f11gs_example
nz = 23;
a = zeros(10000,1);
a(1:nz) = [ 6 + 0.i;
1 + 1.i;
6 + 0.i;
0 + 1.i;
5 + 0.i;
5 + 0.i;
2  2.i;
4 + 0.i;
1 + 1.i;
2 + 0.i;
6 + 0.i;
4 + 3.i;
0 + 1.i;
1 + 0.i;
6 + 0.i;
1  1.i;
0  1.i;
9 + 0.i;
1 + 3.i;
1 + 2.i;
1 + 0.i;
1 + 4.i;
9 + 0.i];
irow = zeros(10000, 1, 'int64');
irow(1:nz) = [int64(1);2;2;3;3;4;5;5;6;6;6;7;7;7;7;8;8;8;9;9;9;9;9];
icol = zeros(10000, 1, 'int64');
icol(1:nz) = [int64(1);1;2;2;3;4;1;5;3;4;6;2;5;6;7;4;6;8;1;5;6;8;9];
lfill = int64(0);
dtol = 0;
mic = 'N';
dscale = 0;
ipiv = [int64(0);0;0;0;0;0;0;0;0];
method = 'CG ';
precon = 'P';
n = 9;
tol = 1e06;
maxitn = int64(20);
anorm = 0;
sigmax = 0;
maxits = int64(9);
monit = int64(2);
irevcm = int64(0);
u = [complex(0);
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i;
0 + 0i];
v = [ 8 + 54i;
10  92i;
25 + 27i;
26  28i;
54 + 12i;
26  22i;
47 + 65i;
71  57i;
60 + 70i];
wgt = [0; 0; 0; 0; 0; 0; 0; 0; 0];
[a, irow, icol, ipiv, istr, nnzc, npivm, ifail] = ...
f11jn(int64(nz), a, irow, icol, lfill, dtol, mic, dscale, ipiv);
[lwreq, work, ifail] = ...
f11gr(method, precon, int64(n), tol, maxitn, anorm, sigmax, maxits, monit, ...
'sigcmp', 's', 'norm_p', '1');
while (irevcm ~= 4)
[irevcm, u, v, work, ifail] = f11gs(irevcm, u, v, wgt, work);
if (irevcm == 1)
[v, ifail] = f11xs(a(1:nz), irow(1:nz), icol(1:nz), 'N', u);
if (ifail ~= 0)
irecvm = 6;
end
elseif (irevcm == 2)
[v, ifail] = f11jp(a, irow, icol, ipiv, istr, 'N', u);
if (ifail ~= 0)
irecvm = 6;
end
elseif (irevcm == 3)
[itn, stplhs, stprhs, anorm, sigmax, its, sigerr, ifail] = f11gt(work);
if (ifail ~= 0)
irecvm = 6;
end
fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', ...
itn, stplhs);
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
[itn, stplhs, stprhs, anorm, sigmax, its, sigerr, ifail] = f11gt(work);
fprintf('\nNumber of iterations for convergence: %d\n', itn);
fprintf('Residual norm: %14.4e\n', stplhs);
fprintf('Righthand side of termination criteria: %14.4e\n', stprhs);
fprintf('inorm 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: 1.4937e+01
+2.1423e01 + +4.5333e+00I
1.6589e+00 + 1.2672e+01I
+2.4101e+00 + +7.4551e+00I
+4.4400e+00 + 6.4174e+00I
+9.1135e+00 + +3.7812e+00I
+4.4419e+00 + 4.0382e+00I
+1.4757e+00 + +1.2662e+00I
+8.4872e+00 + 3.5347e+00I
+5.9948e+00 + +9.6851e01I
Residual Vector
1.8370e+00 + +3.6956e+00I
6.5005e01 + +2.5458e01I
1.2616e01 + 1.3625e01I
1.3120e01 + +1.4130e01I
1.1471e+00 + +7.3386e01I
5.5054e01 + 1.0535e+00I
+1.7165e+00 + 1.4614e+00I
3.5829e01 + +2.8764e01I
3.0278e01 + 3.5324e01I
Monitoring at iteration number 4
residual norm: 1.4602e+00
+1.0061e+00 + +8.9847e+00I
+1.9637e+00 + 7.9768e+00I
+3.0067e+00 + +7.0285e+00I
+3.9830e+00 + 5.9636e+00I
+5.0390e+00 + +5.0432e+00I
+6.0488e+00 + 4.0771e+00I
+6.9710e+00 + +3.0168e+00I
+8.0118e+00 + 1.9806e+00I
+9.0074e+00 + +9.6458e01I
Residual Vector
+1.1524e02 + 2.8188e02I
+1.3513e02 + 1.7345e01I
+1.8173e02 + +1.9627e02I
+1.8900e02 + 2.0354e02I
9.0877e02 + 1.0895e01I
2.3890e01 + +3.2440e01I
+1.9031e01 + 1.5499e02I
+5.1611e02 + 4.1435e02I
+4.3615e02 + +5.0884e02I
Number of iterations for convergence: 5
Residual norm: 7.0166e14
Righthand side of termination criteria: 2.7340e03
inorm of matrix a: 2.2000e+01
Solution Vector
+1.0000e+00 + +9.0000e+00I
+2.0000e+00 + 8.0000e+00I
+3.0000e+00 + +7.0000e+00I
+4.0000e+00 + 6.0000e+00I
+5.0000e+00 + +5.0000e+00I
+6.0000e+00 + 4.0000e+00I
+7.0000e+00 + +3.0000e+00I
+8.0000e+00 + 2.0000e+00I
+9.0000e+00 + +1.0000e+00I
Residual Vector
+1.5099e14 + +1.4211e14I
+1.7764e15 + +0.0000e+00I
+0.0000e+00 + +7.1054e15I
+3.5527e15 + +7.1054e15I
+7.1054e15 + +0.0000e+00I
+0.0000e+00 + 7.1054e15I
+0.0000e+00 + +0.0000e+00I
+0.0000e+00 + 7.1054e15I
+0.0000e+00 + +0.0000e+00I
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013