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_real_gen_precon_jacobi (f11dk)

## Purpose

nag_sparse_real_gen_precon_jacobi (f11dk) computes the approximate solution of a real, symmetric or nonsymmetric, sparse system of linear equations applying a number of Jacobi iterations. It is expected that nag_sparse_real_gen_precon_jacobi (f11dk) will be used as a preconditioner for the iterative solution of real sparse systems of equations.

## Syntax

[x, diag, ifail] = f11dk(store, trans, init, niter, a, irow, icol, check, b, diag, 'n', n, 'nnz', nnz)
[x, diag, ifail] = nag_sparse_real_gen_precon_jacobi(store, trans, init, niter, a, irow, icol, check, b, diag, 'n', n, 'nnz', nnz)

## Description

nag_sparse_real_gen_precon_jacobi (f11dk) computes the approximate solution of the real sparse system of linear equations Ax = b$Ax=b$ using niter iterations of the Jacobi algorithm (see also Golub and Van Loan (1996) and Young (1971)):
 xk + 1 = xk + D − 1(b − Axk) $xk+1=xk+D-1(b-Axk)$ (1)
where k = 1,2,,niter$k=1,2,\dots ,{\mathbf{niter}}$ and x0 = 0${x}_{0}=0$.
nag_sparse_real_gen_precon_jacobi (f11dk) can be used both for nonsymmetric and symmetric systems of equations. For symmetric matrices, either all nonzero elements of the matrix A$A$ can be supplied using coordinate storage (CS), or only the nonzero elements of the lower triangle of A$A$, using symmetric coordinate storage (SCS) (see the F11 Chapter Introduction).
It is expected that nag_sparse_real_gen_precon_jacobi (f11dk) will be used as a preconditioner for the iterative solution of real sparse systems of equations, using either the suite comprising the functions nag_sparse_real_symm_basic_setup (f11gd), nag_sparse_real_symm_basic_solver (f11ge) and nag_sparse_real_symm_basic_diag (f11gf), for symmetric systems, or the suite comprising the functions nag_sparse_real_gen_basic_setup (f11bd), nag_sparse_real_gen_basic_solver (f11be) and nag_sparse_real_gen_basic_diag (f11bf), for nonsymmetric 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:     store – string (length ≥ 1)
Specifies whether the matrix A$A$ is stored using symmetric coordinate storage (SCS) (applicable only to a symmetric matrix A$A$) or coordinate storage (CS) (applicable to both symmetric and nonsymmetric matrices).
store = 'N'${\mathbf{store}}=\text{'N'}$
The complete matrix A$A$ is stored in CS format.
store = 'S'${\mathbf{store}}=\text{'S'}$
The lower triangle of the symmetric matrix A$A$ is stored in SCS format.
Constraint: store = 'N'${\mathbf{store}}=\text{'N'}$ or 'S'$\text{'S'}$.
2:     trans – string (length ≥ 1)
If store = 'N'${\mathbf{store}}=\text{'N'}$, specifies whether the approximate solution of Ax = b$Ax=b$ or of ATx = b${A}^{\mathrm{T}}x=b$ is required.
trans = 'N'${\mathbf{trans}}=\text{'N'}$
The approximate solution of Ax = b$Ax=b$ is calculated.
trans = 'T'${\mathbf{trans}}=\text{'T'}$
The approximate solution of ATx = b${A}^{\mathrm{T}}x=b$ is calculated.
Constraint: trans = 'N'${\mathbf{trans}}=\text{'N'}$ or 'T'$\text{'T'}$.
3:     init – string (length ≥ 1)
On first entry, init should be set to 'I', unless the diagonal elements of A$A$ are already stored in the array diag. If diag already contains the diagonal of A$A$, it must be set to 'N'.
init = 'N'${\mathbf{init}}=\text{'N'}$
diag must contain the diagonal of A$A$.
init = 'I'${\mathbf{init}}=\text{'I'}$
diag will store the diagonal of A$A$ on exit.
Constraint: init = 'N'${\mathbf{init}}=\text{'N'}$ or 'I'$\text{'I'}$.
4:     niter – int64int32nag_int scalar
The number of Jacobi iterations requested.
Constraint: niter1${\mathbf{niter}}\ge 1$.
5:     a(nnz) – double array
nnz, the dimension of the array, must satisfy the constraint
• if store = 'N'${\mathbf{store}}=\text{'N'}$, 1nnzn2$1\le {\mathbf{nnz}}\le {{\mathbf{n}}}^{2}$;
• if store = 'S'${\mathbf{store}}=\text{'S'}$, 1nnzn × (n + 1) / 2$1\le {\mathbf{nnz}}\le {\mathbf{n}}×\left({\mathbf{n}}+1\right)/2$.
If store = 'N'${\mathbf{store}}=\text{'N'}$, the nonzero elements in the matrix A$A$ (CS format).
If store = 'S'${\mathbf{store}}=\text{'S'}$, the nonzero elements in the lower triangle of the matrix A$A$ (SCS format).
In both cases, the elements of either A$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_real_gen_sort (f11za) or nag_sparse_real_symm_sort (f11zb) may be used to reorder the elements in this way for CS and SCS storage, respectively.
6:     irow(nnz) – int64int32nag_int array
7:     icol(nnz) – int64int32nag_int array
nnz, the dimension of the array, must satisfy the constraint
• if store = 'N'${\mathbf{store}}=\text{'N'}$, 1nnzn2$1\le {\mathbf{nnz}}\le {{\mathbf{n}}}^{2}$;
• if store = 'S'${\mathbf{store}}=\text{'S'}$, 1nnzn × (n + 1) / 2$1\le {\mathbf{nnz}}\le {\mathbf{n}}×\left({\mathbf{n}}+1\right)/2$.
If store = 'N'${\mathbf{store}}=\text{'N'}$, the row and column indices of the nonzero elements supplied in a.
If store = 'S'${\mathbf{store}}=\text{'S'}$, the row and column indices of the nonzero elements of the lower triangle of the matrix A$A$ supplied in a.
Constraints:
• 1irow(i)n$1\le {\mathbf{irow}}\left(\mathit{i}\right)\le {\mathbf{n}}$, for i = 1,2,,nnz$\mathit{i}=1,2,\dots ,{\mathbf{nnz}}$;
• if store = 'N'${\mathbf{store}}=\text{'N'}$, 1icol(i)n$1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{n}}$, for i = 1,2,,nnz$\mathit{i}=1,2,\dots ,{\mathbf{nnz}}$;
• if store = 'S'${\mathbf{store}}=\text{'S'}$, 1icol(i)irow(i)$1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{irow}}\left(\mathit{i}\right)$, for i = 1,2,,nnz$\mathit{i}=1,2,\dots ,{\mathbf{nnz}}$;
• either irow(i1) < irow(i)${\mathbf{irow}}\left(\mathit{i}-1\right)<{\mathbf{irow}}\left(\mathit{i}\right)$ or both irow(i1) = irow(i)${\mathbf{irow}}\left(\mathit{i}-1\right)={\mathbf{irow}}\left(\mathit{i}\right)$ and icol(i1) < icol(i)${\mathbf{icol}}\left(\mathit{i}-1\right)<{\mathbf{icol}}\left(\mathit{i}\right)$, for i = 2,3,,nnz$\mathit{i}=2,3,\dots ,{\mathbf{nnz}}$.
8:     check – string (length ≥ 1)
Specifies whether or not the CS or SCS representation of the matrix A$A$ should be checked.
check = 'C'${\mathbf{check}}=\text{'C'}$
Checks are carried out on the values of n, nnz, irow, icol; if init = 'N'${\mathbf{init}}=\text{'N'}$, diag is also checked.
check = 'N'${\mathbf{check}}=\text{'N'}$
None of these checks are carried out.
Constraint: check = 'C'${\mathbf{check}}=\text{'C'}$ or 'N'$\text{'N'}$.
9:     b(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The right-hand side vector b$b$.
10:   diag(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
If init = 'N'${\mathbf{init}}=\text{'N'}$, the diagonal elements of A$A$.

### Optional Input Parameters

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

work

### Output Parameters

1:     x(n) – double array
The approximate solution vector ${x}_{{\mathbf{niter}}}$.
2:     diag(n) – double array
If init = 'N'${\mathbf{init}}=\text{'N'}$, unchanged on exit.
If init = 'I'${\mathbf{init}}=\text{'I'}$, the diagonal elements of A$A$.
3:     ifail – int64int32nag_int scalar
${\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:
ifail = 1${\mathbf{ifail}}=1$
 On entry, store ≠ 'N'${\mathbf{store}}\ne \text{'N'}$ or 'S'$\text{'S'}$, or trans ≠ 'N'${\mathbf{trans}}\ne \text{'N'}$ or 'T'$\text{'T'}$, or init ≠ 'N'${\mathbf{init}}\ne \text{'N'}$ or 'I'$\text{'I'}$, or check ≠ 'C'${\mathbf{check}}\ne \text{'C'}$ or 'N'$\text{'N'}$, or niter ≤ 0${\mathbf{niter}}\le 0$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, n < 1${\mathbf{n}}<1$, or nnz < 1${\mathbf{nnz}}<1$, or nnz > n2${\mathbf{nnz}}>{{\mathbf{n}}}^{2}$, if store = 'N'${\mathbf{store}}=\text{'N'}$, or 1 ≤ nnz ≤ [n(n + 1)] / 2$1\le {\mathbf{nnz}}\le \left[{\mathbf{n}}\left({\mathbf{n}}+1\right)\right]/2$, if store = 'S'${\mathbf{store}}=\text{'S'}$.
ifail = 3${\mathbf{ifail}}=3$
On entry, the arrays irow and icol fail to satisfy the following constraints:
• 1irow(i)n$1\le {\mathbf{irow}}\left(i\right)\le {\mathbf{n}}$ and
• if store = 'N'${\mathbf{store}}=\text{'N'}$ then 1icol(i)n$1\le {\mathbf{icol}}\left(i\right)\le {\mathbf{n}}$, or
• if store = 'S'${\mathbf{store}}=\text{'S'}$ then 1icol(i)irow(i)$1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{irow}}\left(\mathit{i}\right)$, for i = 1,2,,nnz$\mathit{i}=1,2,\dots ,{\mathbf{nnz}}$.
• irow(i1) < irow(i)${\mathbf{irow}}\left(\mathit{i}-1\right)<{\mathbf{irow}}\left(\mathit{i}\right)$ or irow(i1) = irow(i)${\mathbf{irow}}\left(\mathit{i}-1\right)={\mathbf{irow}}\left(\mathit{i}\right)$ and icol(i1) < icol(i)${\mathbf{icol}}\left(\mathit{i}-1\right)<{\mathbf{icol}}\left(\mathit{i}\right)$, for i = 2,3,,nnz$\mathit{i}=2,3,\dots ,{\mathbf{nnz}}$.
Therefore a nonzero element has been supplied which does not lie within the matrix A$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 store = 'N'${\mathbf{store}}=\text{'N'}$ or store = 'S'${\mathbf{store}}=\text{'S'}$, respectively.
ifail = 4${\mathbf{ifail}}=4$
On entry, init = 'N'${\mathbf{init}}=\text{'N'}$ and some diagonal elements of A$A$ stored in diag are zero.
ifail = 5${\mathbf{ifail}}=5$
On entry, init = 'I'${\mathbf{init}}=\text{'I'}$ and some diagonal elements of A$A$ are zero.

## 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_real_gen_precon_jacobi (f11dk) is proportional to ${\mathbf{niter}}×{\mathbf{nnz}}$.

### Use of check

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

## Example

```function nag_sparse_real_gen_precon_jacobi_example
store = 'Non symmetric';
trans = 'N';
init = 'I';
niter = int64(4);
n = int64(8);
wgt = zeros(n, 1);
dgnl = zeros(n, 1);
lwork = int64(200);
method = 'bicgstab';
precon = 'p';
m = int64(2);
tol = 1e-6;
maxitn = int64(20);
monit = int64(1);
niter = int64(4);
a = [4;-1;1;4;-5;2;-7;2;2;-1;6;2;-1;8;-2;-2;5;8;-2;-1;7;-1;2;6];
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 = [6; 8; -9; 46; 17; 21; 22; 34];
x = zeros(n, 1);
anorm = 0;
sigmax = 0;

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

% Repeatedly call nag_sparse_real_gen_basic_solver to solve the equations.  On final exit x will
% contain the solution and b the residual vector.

irevcm = int64(0);
init = 'i';
ifail = int64(0);

while irevcm ~= 4
[irevcm, x, b, work, ifail] = nag_sparse_real_gen_basic_solver(irevcm, x, b, wgt, work);

if (irevcm == -1)
[b, ifail1] = nag_sparse_real_gen_matvec('Transpose', a, irow, icol, 'N', x);
elseif (irevcm == 1)
[b, ifail1] = nag_sparse_real_gen_matvec('No Transpose', a, irow, icol, 'N', x);
elseif (irevcm == 2)
[b, dgnl, ifail1] = nag_sparse_real_gen_precon_jacobi(store, trans, init, niter, a, ...
irow, icol, 'Check', x, dgnl);
init = 'n';
elseif (irevcm == 3)
[itn, stplhs, stprhs, anorm, sigmax, ifail1] = nag_sparse_real_gen_basic_diag(work);
fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', itn, stplhs);
end

if (ifail1 ~= 0)
irevcm = 6;
end
end

% Get information about the computation
if ifail == 0
[itn, stplhs, stprhs, anorm, sigmax, ifail] = nag_sparse_real_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  Residual Vector\n');
for i = 1:n
fprintf('%+16.4e %+16.4e\n', x(i), b(i));
end
end
```
```

Number of iterations for convergence: 2
Residual norm:                               1.1177e-04
Right-hand side of termination criteria:     5.4082e-04
i-norm of matrix a:                          1.5000e+01

Solution Vector  Residual Vector
+1.7035e+00      +3.2377e-07
+1.0805e+00      -1.7625e-05
+1.8305e+00      +2.7964e-05
+6.0251e+00      -2.5914e-05
+3.2942e+00      +7.8156e-06
+1.9068e+00      +9.2064e-06
+4.1365e+00      -3.0848e-06
+5.2111e+00      +1.9834e-05

```
```function f11dk_example
store = 'Non symmetric';
trans = 'N';
init = 'I';
niter = int64(4);
n = int64(8);
wgt = zeros(n, 1);
dgnl = zeros(n, 1);
lwork = int64(200);
method = 'bicgstab';
precon = 'p';
m = int64(2);
tol = 1e-6;
maxitn = int64(20);
monit = int64(1);
niter = int64(4);
a = [4;-1;1;4;-5;2;-7;2;2;-1;6;2;-1;8;-2;-2;5;8;-2;-1;7;-1;2;6];
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 = [6; 8; -9; 46; 17; 21; 22; 34];
x = zeros(n, 1);
anorm = 0;
sigmax = 0;

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

% Repeatedly call f11be to solve the equations.  On final exit x will
% contain the solution and b the residual vector.

irevcm = int64(0);
init = 'i';
ifail = int64(0);

while irevcm ~= 4
[irevcm, x, b, work, ifail] = f11be(irevcm, x, b, wgt, work);

if (irevcm == -1)
[b, ifail1] = f11xa('Transpose', a, irow, icol, 'N', x);
elseif (irevcm == 1)
[b, ifail1] = f11xa('No Transpose', a, irow, icol, 'N', x);
elseif (irevcm == 2)
[b, dgnl, ifail1] = f11dk(store, trans, init, niter, a, ...
irow, icol, 'Check', x, dgnl);
init = 'n';
elseif (irevcm == 3)
[itn, stplhs, stprhs, anorm, sigmax, ifail1] = f11bf(work);
fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', itn, stplhs);
end

if (ifail1 ~= 0)
irevcm = 6;
end
end

% Get information about the computation
if ifail == 0
[itn, stplhs, stprhs, anorm, sigmax, ifail] = f11bf(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  Residual Vector\n');
for i = 1:n
fprintf('%+16.4e %+16.4e\n', x(i), b(i));
end
end
```
```

Number of iterations for convergence: 2
Residual norm:                               1.1177e-04
Right-hand side of termination criteria:     5.4082e-04
i-norm of matrix a:                          1.5000e+01

Solution Vector  Residual Vector
+1.7035e+00      +3.2377e-07
+1.0805e+00      -1.7625e-05
+1.8305e+00      +2.7964e-05
+6.0251e+00      -2.5914e-05
+3.2942e+00      +7.8156e-06
+1.9068e+00      +9.2064e-06
+4.1365e+00      -3.0848e-06
+5.2111e+00      +1.9834e-05

```