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

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_sparse_complex_gen_precon_jacobi (f11dx)

## Purpose

nag_sparse_complex_gen_precon_jacobi (f11dx) computes the approximate solution of a complex, Hermitian or non-Hermitian, sparse system of linear equations applying a number of Jacobi iterations. It is expected that nag_sparse_complex_gen_precon_jacobi (f11dx) will be used as a preconditioner for the iterative solution of complex sparse systems of equations.

## Syntax

[x, diag, ifail] = f11dx(store, trans, init, niter, a, irow, icol, check, b, diag, 'n', n, 'nz', nz)
[x, diag, ifail] = nag_sparse_complex_gen_precon_jacobi(store, trans, init, niter, a, irow, icol, check, b, diag, 'n', n, 'nz', nz)

## Description

nag_sparse_complex_gen_precon_jacobi (f11dx) computes the approximate solution of the complex sparse system of linear equations $Ax=b$ using niter iterations of the Jacobi algorithm (see also Golub and Van Loan (1996) and Young (1971)):
 $xk+1=xk+D-1b-Axk$ (1)
where $k=1,2,\dots ,{\mathbf{niter}}$ and ${x}_{0}=0$.
nag_sparse_complex_gen_precon_jacobi (f11dx) can be used both for non-Hermitian and Hermitian systems of equations. For Hermitian matrices, either all nonzero elements of the matrix $A$ can be supplied using coordinate storage (CS), or only the nonzero elements of the lower triangle of $A$, using symmetric coordinate storage (SCS) (see the F11 Chapter Introduction).
It is expected that nag_sparse_complex_gen_precon_jacobi (f11dx) will be used as a preconditioner for the iterative solution of complex sparse systems of equations, using either the suite comprising the functions nag_sparse_complex_herm_basic_setup (f11gr), nag_sparse_complex_herm_basic_solver (f11gs) and nag_sparse_complex_herm_basic_diag (f11gt), for Hermitian systems, or the suite comprising the functions nag_sparse_complex_gen_basic_setup (f11br), nag_sparse_complex_gen_basic_solver (f11bs) and nag_sparse_complex_gen_basic_diag (f11bt), for non-Hermitian systems of equations.

## References

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Young D (1971) Iterative Solution of Large Linear Systems Academic Press, New York

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{store}$ – string (length ≥ 1)
Specifies whether the matrix $A$ is stored using symmetric coordinate storage (SCS) (applicable only to a Hermitian matrix $A$) or coordinate storage (CS) (applicable to both Hermitian and non-Hermitian matrices).
${\mathbf{store}}=\text{'N'}$
The complete matrix $A$ is stored in CS format.
${\mathbf{store}}=\text{'S'}$
The lower triangle of the Hermitian matrix $A$ is stored in SCS format.
Constraint: ${\mathbf{store}}=\text{'N'}$ or $\text{'S'}$.
2:     $\mathrm{trans}$ – string (length ≥ 1)
Suggested value: if the matrix $A$ is Hermitian and stored in CS format, it is recommended that ${\mathbf{trans}}=\text{'N'}$ for reasons of efficiency.
If ${\mathbf{store}}=\text{'N'}$, specifies whether the approximate solution of $Ax=b$ or of ${A}^{\mathrm{H}}x=b$ is required.
${\mathbf{trans}}=\text{'N'}$
The approximate solution of $Ax=b$ is calculated.
${\mathbf{trans}}=\text{'T'}$
The approximate solution of ${A}^{\mathrm{H}}x=b$ is calculated.
Constraint: ${\mathbf{trans}}=\text{'N'}$ or $\text{'T'}$.
3:     $\mathrm{init}$ – string (length ≥ 1)
Suggested value: ${\mathbf{init}}=\text{'I'}$ on first entry; ${\mathbf{init}}=\text{'N'}$, subsequently, unless diag has been overwritten.
On first entry, init should be set to 'I', unless the diagonal elements of $A$ are already stored in the array diag. If diag already contains the diagonal of $A$, it must be set to 'N'.
${\mathbf{init}}=\text{'N'}$
diag must contain the diagonal of $A$.
${\mathbf{init}}=\text{'I'}$
diag will store the diagonal of $A$ on exit.
Constraint: ${\mathbf{init}}=\text{'N'}$ or $\text{'I'}$.
4:     $\mathrm{niter}$int64int32nag_int scalar
The number of Jacobi iterations requested.
Constraint: ${\mathbf{niter}}\ge 1$.
5:     $\mathrm{a}\left({\mathbf{nz}}\right)$ – complex array
If ${\mathbf{store}}=\text{'N'}$, the nonzero elements in the matrix $A$ (CS format).
If ${\mathbf{store}}=\text{'S'}$, the nonzero elements in the lower triangle of the matrix $A$ (SCS format).
In both cases, the elements of either $A$ or of its lower triangle must be ordered by increasing row index and by increasing column index within each row. Multiple entries for the same row and columns indices are not permitted. The function nag_sparse_complex_gen_sort (f11zn) or nag_sparse_complex_herm_sort (f11zp) may be used to reorder the elements in this way for CS and SCS storage, respectively.
6:     $\mathrm{irow}\left({\mathbf{nz}}\right)$int64int32nag_int array
7:     $\mathrm{icol}\left({\mathbf{nz}}\right)$int64int32nag_int array
If ${\mathbf{store}}=\text{'N'}$, the row and column indices of the nonzero elements supplied in a.
If ${\mathbf{store}}=\text{'S'}$, the row and column indices of the nonzero elements of the lower triangle of the matrix $A$ supplied in a.
Constraints:
• $1\le {\mathbf{irow}}\left(\mathit{i}\right)\le {\mathbf{n}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{nz}}$;
• if ${\mathbf{store}}=\text{'N'}$, $1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{n}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{nz}}$;
• if ${\mathbf{store}}=\text{'S'}$, $1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{irow}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{nz}}$;
• either ${\mathbf{irow}}\left(\mathit{i}-1\right)<{\mathbf{irow}}\left(\mathit{i}\right)$ or both ${\mathbf{irow}}\left(\mathit{i}-1\right)={\mathbf{irow}}\left(\mathit{i}\right)$ and ${\mathbf{icol}}\left(\mathit{i}-1\right)<{\mathbf{icol}}\left(\mathit{i}\right)$, for $\mathit{i}=2,3,\dots ,{\mathbf{nz}}$.
8:     $\mathrm{check}$ – string (length ≥ 1)
Specifies whether or not the CS or SCS representation of the matrix $A$ should be checked.
${\mathbf{check}}=\text{'C'}$
Checks are carried out on the values of n, nz, irow, icol; if ${\mathbf{init}}=\text{'N'}$, diag is also checked.
${\mathbf{check}}=\text{'N'}$
None of these checks are carried out.
Constraint: ${\mathbf{check}}=\text{'C'}$ or $\text{'N'}$.
9:     $\mathrm{b}\left({\mathbf{n}}\right)$ – complex array
The right-hand side vector $b$.
10:   $\mathrm{diag}\left({\mathbf{n}}\right)$ – complex array
If ${\mathbf{init}}=\text{'N'}$, the diagonal elements of $A$.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the dimension of the arrays b, diag. (An error is raised if these dimensions are not equal.)
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 1$.
2:     $\mathrm{nz}$int64int32nag_int scalar
Default: the dimension of the arrays a, irow, icol. (An error is raised if these dimensions are not equal.)
If ${\mathbf{store}}=\text{'N'}$, the number of nonzero elements in the matrix $A$.
If ${\mathbf{store}}=\text{'S'}$, the number of nonzero elements in the lower triangle of the matrix $A$.
Constraints:
• if ${\mathbf{store}}=\text{'N'}$, $1\le {\mathbf{nz}}\le {{\mathbf{n}}}^{2}$;
• if ${\mathbf{store}}=\text{'S'}$, $1\le {\mathbf{nz}}\le {\mathbf{n}}×\left({\mathbf{n}}+1\right)/2$.

### Output Parameters

1:     $\mathrm{x}\left({\mathbf{n}}\right)$ – complex array
The approximate solution vector ${x}_{{\mathbf{niter}}}$.
2:     $\mathrm{diag}\left({\mathbf{n}}\right)$ – complex array
If ${\mathbf{init}}=\text{'N'}$, unchanged on exit.
If ${\mathbf{init}}=\text{'I'}$, the diagonal elements of $A$.
3:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{store}}\ne \text{'N'}$ or $\text{'S'}$, or ${\mathbf{trans}}\ne \text{'N'}$ or $\text{'T'}$, or ${\mathbf{init}}\ne \text{'N'}$ or $\text{'I'}$, or ${\mathbf{check}}\ne \text{'C'}$ or $\text{'N'}$, or ${\mathbf{niter}}\le 0$.
${\mathbf{ifail}}=2$
 On entry, ${\mathbf{n}}<1$, or ${\mathbf{nz}}<1$, or ${\mathbf{nz}}>{{\mathbf{n}}}^{2}$, if ${\mathbf{store}}=\text{'N'}$, or $1\le {\mathbf{nz}}\le \left[{\mathbf{n}}\left({\mathbf{n}}+1\right)\right]/2$, if ${\mathbf{store}}=\text{'S'}$.
${\mathbf{ifail}}=3$
On entry, the arrays irow and icol fail to satisfy the following constraints:
• $1\le {\mathbf{irow}}\left(i\right)\le {\mathbf{n}}$ and
• if ${\mathbf{store}}=\text{'N'}$ then $1\le {\mathbf{icol}}\left(i\right)\le {\mathbf{n}}$, or
• if ${\mathbf{store}}=\text{'S'}$ then $1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{irow}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{nz}}$.
• ${\mathbf{irow}}\left(\mathit{i}-1\right)<{\mathbf{irow}}\left(\mathit{i}\right)$ or ${\mathbf{irow}}\left(\mathit{i}-1\right)={\mathbf{irow}}\left(\mathit{i}\right)$ and ${\mathbf{icol}}\left(\mathit{i}-1\right)<{\mathbf{icol}}\left(\mathit{i}\right)$, for $\mathit{i}=2,3,\dots ,{\mathbf{nz}}$.
Therefore a nonzero element has been supplied which does not lie within the matrix $A$, is out of order, or has duplicate row and column indices. Call either nag_sparse_real_gen_sort (f11za) or nag_sparse_real_symm_sort (f11zb) to reorder and sum or remove duplicates when ${\mathbf{store}}=\text{'N'}$ or ${\mathbf{store}}=\text{'S'}$, respectively.
${\mathbf{ifail}}=4$
On entry, ${\mathbf{init}}=\text{'N'}$ and some diagonal elements of $A$ stored in diag are zero.
${\mathbf{ifail}}=5$
On entry, ${\mathbf{init}}=\text{'I'}$ and some diagonal elements of $A$ are zero.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

In general, the Jacobi method cannot be used on its own to solve systems of linear equations. The rate of convergence is bound by its spectral properties (see, for example, Golub and Van Loan (1996)) and as a solver, the Jacobi method can only be applied to a limited set of matrices. One condition that guarantees convergence is strict diagonal dominance.
However, the Jacobi method can be used successfully as a preconditioner to a wider class of systems of equations. The Jacobi method has good vector/parallel properties, hence it can be applied very efficiently. Unfortunately, it is not possible to provide criteria which define the applicability of the Jacobi method as a preconditioner, and its usefulness must be judged for each case.

### Timing

The time taken for a call to nag_sparse_complex_gen_precon_jacobi (f11dx) is proportional to ${\mathbf{niter}}×{\mathbf{nz}}$.

### Use of check

It is expected that a common use of nag_sparse_complex_gen_precon_jacobi (f11dx) will be as preconditioner for the iterative solution of complex, Hermitian or non-Hermitian, linear systems. In this situation, nag_sparse_complex_gen_precon_jacobi (f11dx) is likely to be called many times. In the interests of both reliability and efficiency, you are recommended to set ${\mathbf{check}}=\text{'C'}$ for the first of such calls, and to set ${\mathbf{check}}=\text{'N'}$ for all subsequent calls.

## Example

This example solves the complex sparse non-Hermitian system of equations $Ax=b$ iteratively using nag_sparse_complex_gen_precon_jacobi (f11dx) as a preconditioner.
```function f11dx_example

fprintf('f11dx example results\n\n');

% Solve sparse system Ax = b iteratively using Jacobi preconditioner

% Define sparse matrix A and RHS B
nz   = int64(24);
n    = int64(8);
a = [ 2 + 1i, -1 + 1i,  1 - 3i,  4 + 7i, -3 + 0i,  2 + 4i, ...
-7 - 5i,  2 + 1i,  3 + 2i, -4 + 2i,  0 + 1i,  5 - 3i, ...
-1 + 2i,  8 + 6i, -3 - 4i, -6 - 2i,  5 - 2i,  2 + 0i, ...
0 - 5i, -1 + 5i,  6 + 2i, -1 + 4i,  2 + 0i,  3 + 3i];
b = [ 7 + 11i;
1 + 24i;
-13 - 18i;
-10 +  3i;
23 + 14i;
17 -  7i;
15 -  3i;
-3 + 20i];
irow = 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 = 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]);

b = [  7 + 11i;
1 + 24i;
-13 - 18i;
-10 +  3i;
23 + 14i;
17 -  7i;
15 -  3i;
-3 + 20i];

% Iterative Solution Setup
%      Input parameters
method = 'TFQMR';
precon = 'p';
m      = int64(2);
tol    = 1e-6;
maxitn = int64(20);
anorm  = 0;
sigmax = 0;
monit  = int64(1);
lwork = int64(300);

% Initialise Solver
[lwreq, work, ifail] = f11br( ...
method, precon, n, m, tol, maxitn, anorm, ...
sigmax, monit, lwork, 'norm_p', '1');

% Solver inputs
irevcm = int64(0);
x = complex(zeros(n, 1));
wgt= zeros(n, 1);

% Preconditioner inputs
store = 'Non Hermitian';
trans = 'N';
init  = 'I';
niter = int64(2);
check = 'Check';
diag  = complex(zeros(n, 1));

% Solve by reverse communication
while (irevcm ~= 4)
[irevcm, x, b, work, ifail] = f11bs( ...
irevcm, x, b, wgt, work);

if (irevcm == -1)
% b = A^Tx
[b, ifail] = f11xn('T', a, irow, icol, 'N', x);
elseif (irevcm == 1)
% b = Ax
[b, ifail] = f11xn('N', a, irow, icol, 'N', x);
elseif (irevcm == 2)
% Jacobi preconditioner
[b, diag, ifail] = f11dx( ...
store, trans, init, niter, ...
a, irow, icol, check, x, diag);
elseif (irevcm == 3)
[itn, stplhs, stprhs, anorm, sigmax, work, ifail] = ...
f11bt(work);
fprintf('\nMonitoring at iteration number %d\n', itn);
fprintf('residual norm: %14.4e\n', stplhs);
end
end

% Obtain 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');
disp(x);
fprintf('\n   Residual Vector\n');
disp(b);

```
```f11dx example results

Monitoring at iteration number 1
residual norm:     1.5062e+02

Monitoring at iteration number 2
residual norm:     1.5704e+02

Monitoring at iteration number 3
residual norm:     1.4803e+02

Monitoring at iteration number 4
residual norm:     8.5215e+01

Monitoring at iteration number 5
residual norm:     4.2951e+01

Monitoring at iteration number 6
residual norm:     2.5055e+01

Monitoring at iteration number 7
residual norm:     1.9090e-01

Number of iterations for convergence: 8
Residual norm:                               9.5485e-08
Right-hand side of termination criteria:     8.9100e-04
i-norm of matrix a:                          2.7000e+01

Solution Vector
1.0000 + 1.0000i
2.0000 - 1.0000i
3.0000 + 1.0000i
4.0000 - 1.0000i
3.0000 - 1.0000i
2.0000 + 1.0000i
1.0000 - 1.0000i
0.0000 + 3.0000i

Residual Vector
1.0e-07 *

0.0471 - 0.0394i
0.0877 - 0.0981i
-0.0287 + 0.0416i
0.0358 - 0.1112i
-0.0692 - 0.0869i
-0.0225 + 0.0849i
0.0052 - 0.0334i
-0.0558 - 0.1073i

```