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, 'nz', nz)
[x, diag, ifail] = nag_sparse_real_gen_precon_jacobi(store, trans, init, niter, a, irow, icol, check, b, diag, 'n', n, 'nz', nz)

## Description

nag_sparse_real_gen_precon_jacobi (f11dk) computes the approximate solution of the real 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_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$ 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_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:     $\mathrm{store}$ – string (length ≥ 1)
Specifies whether the matrix $A$ is stored using symmetric coordinate storage (SCS) (applicable only to a symmetric matrix $A$) or coordinate storage (CS) (applicable to both symmetric and non-symmetric matrices).
${\mathbf{store}}=\text{'N'}$
The complete matrix $A$ is stored in CS format.
${\mathbf{store}}=\text{'S'}$
The lower triangle of the symmetric 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 symmetric 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{T}}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{T}}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)$ – double 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_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:     $\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)$ – double array
The right-hand side vector $b$.
10:   $\mathrm{diag}\left({\mathbf{n}}\right)$ – double 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)$ – double array
The approximate solution vector ${x}_{{\mathbf{niter}}}$.
2:     $\mathrm{diag}\left({\mathbf{n}}\right)$ – double 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_real_gen_precon_jacobi (f11dk) is proportional to ${\mathbf{niter}}×{\mathbf{nz}}$.

### 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 ${\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 real sparse nonsymmetric system of equations $Ax=b$ iteratively using nag_sparse_real_gen_precon_jacobi (f11dk) as a preconditioner.
```function f11dk_example

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

% Solve sparse system Ax = b iteratively using Jacobi preconditioner

% Define sparse matrix A and RHS B
n    = int64(8);
nz   = int64(24);
irow(1:nz) = 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(1:nz) = 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]);
a(1:nz)    = [4;-1; 1; 4;-5; 2;-7; 2; 2;-1; 6; 2;
-1; 8;-2;-2; 5; 8;-2;-1; 7;-1; 2; 6];
b          = [6; 8;-9;46;17;21;22;34];

% Iterative Solution Setup
%      Input parameters
method = 'BICGSTAB';
precon = 'Preconditioned';
lpoly  = int64(2);
tol    = sqrt(x02aj)^(3/8);
maxitn = int64(20);
anorm  = 0;
sigmax = 0;
monit  = int64(1);
lwork  = int64(200);

% Setup iterative method
[lwreq, work, ifail] = ...
f11bd(method, precon, n, lpoly, tol, maxitn, anorm, sigmax, ...
monit, lwork, 'norm_p', '1');

% Preconditioner parameters
store = 'Non symmetric';
trans = 'N';
init  = 'I';
niter = int64(4);
wgt   = zeros(n, 1);
dgnl  = zeros(n, 1);
x     = zeros(n, 1);

% Repeatedly call f11be to solve the equations.
% On final exit x will contain the solution and b the residual vector.
% Jacobi preconditioner is via f11dk when irevcm==2.

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 = A^Tx
[b, ifail1] = f11xa(...
'Transpose', a, irow, icol, 'N', x);
elseif (irevcm == 1)
% b = Ax
[b, ifail1] = f11xa(...
'No Transpose', a, irow, icol, 'N', x);
elseif (irevcm == 2)
% Jacobi preconditioning Jb = x
[b, dgnl, ifail1] = f11dk(...
store, trans, init, niter, a(1:nz), irow, icol, 'Check', x, dgnl);
init = 'n';
elseif (irevcm == 3)
% Monitoring
[itn, stplhs, stprhs, anorm, sigmax, ifail1] = ...
f11bf(work);
fprintf('\nMonitoring at iteration number %2d\n',itn);
fprintf('residual norm:              %14.4e\n', stplhs);
end

end

% Get information about the computation
[itn, stplhs, stprhs, anorm, sigmax, ifail] = ...
f11bf(work);

fprintf('\nNumber of iterations for convergence:     %4d\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.4f %16.2e\n', x(i), b(i));
end

```
```f11dk example results

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

Solution Vector  Residual Vector
1.7035         3.24e-07
1.0805        -1.76e-05
1.8305         2.80e-05
6.0251        -2.59e-05
3.2942         7.82e-06
1.9068         9.21e-06
4.1365        -3.08e-06
5.2111         1.98e-05
```