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_symm_precon_ssor_solve (f11jd)

## Purpose

nag_sparse_real_symm_precon_ssor_solve (f11jd) solves a system of linear equations involving the preconditioning matrix corresponding to SSOR applied to a real sparse symmetric matrix, represented in symmetric coordinate storage format.

## Syntax

[x, ifail] = f11jd(a, irow, icol, rdiag, omega, check, y, 'n', n, 'nz', nz)
[x, ifail] = nag_sparse_real_symm_precon_ssor_solve(a, irow, icol, rdiag, omega, check, y, 'n', n, 'nz', nz)

## Description

nag_sparse_real_symm_precon_ssor_solve (f11jd) solves a system of equations
 $Mx=y$
involving the preconditioning matrix
 $M=1ω2-ω D+ω L D-1 D+ω LT$
corresponding to symmetric successive-over-relaxation (SSOR) (see Young (1971)) on a linear system $Ax=b$, where $A$ is a sparse symmetric matrix stored in symmetric coordinate storage (SCS) format (see Symmetric coordinate storage (SCS) format in the F11 Chapter Introduction).
In the definition of $M$ given above $D$ is the diagonal part of $A$, $L$ is the strictly lower triangular part of $A$, and $\omega$ is a user-defined relaxation parameter.
It is envisaged that a common use of nag_sparse_real_symm_precon_ssor_solve (f11jd) will be to carry out the preconditioning step required in the application of nag_sparse_real_symm_basic_solver (f11ge) to sparse linear systems. For an illustration of this use of nag_sparse_real_symm_precon_ssor_solve (f11jd) see the example program given in Example. nag_sparse_real_symm_precon_ssor_solve (f11jd) is also used for this purpose by the Black Box function nag_sparse_real_symm_solve_jacssor (f11je).

## References

Young D (1971) Iterative Solution of Large Linear Systems Academic Press, New York

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{a}\left({\mathbf{nz}}\right)$ – double array
The nonzero elements in the lower triangular part of the matrix $A$, ordered by increasing row index, and by increasing column index within each row. Multiple entries for the same row and column indices are not permitted. The function nag_sparse_real_symm_sort (f11zb) may be used to order the elements in this way.
2:     $\mathrm{irow}\left({\mathbf{nz}}\right)$int64int32nag_int array
3:     $\mathrm{icol}\left({\mathbf{nz}}\right)$int64int32nag_int array
The row and column indices of the nonzero elements supplied in array a.
Constraints:
irow and icol must satisfy these constraints (which may be imposed by a call to nag_sparse_real_symm_sort (f11zb)):
• $1\le {\mathbf{irow}}\left(\mathit{i}\right)\le {\mathbf{n}}$ and $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}}$.
4:     $\mathrm{rdiag}\left({\mathbf{n}}\right)$ – double array
The elements of the diagonal matrix ${D}^{-1}$, where $D$ is the diagonal part of $A$.
5:     $\mathrm{omega}$ – double scalar
The relaxation parameter $\omega$.
Constraint: $0.0<{\mathbf{omega}}<2.0$.
6:     $\mathrm{check}$ – string (length ≥ 1)
Specifies whether or not the input data should be checked.
${\mathbf{check}}=\text{'C'}$
Checks are carried out on the values of n, nz, irow, icol and omega.
${\mathbf{check}}=\text{'N'}$
None of these checks are carried out.
Constraint: ${\mathbf{check}}=\text{'C'}$ or $\text{'N'}$.
7:     $\mathrm{y}\left({\mathbf{n}}\right)$ – double array
The right-hand side vector $y$.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the dimension of the arrays rdiag, y. (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.)
The number of nonzero elements in the lower triangular part of $A$.
Constraint: $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 solution vector $x$.
2:     $\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{check}}\ne \text{'C'}$ or $\text{'N'}$.
${\mathbf{ifail}}=2$
 On entry, ${\mathbf{n}}<1$, or ${\mathbf{nz}}<1$, or ${\mathbf{nz}}>{\mathbf{n}}×\left({\mathbf{n}}+1\right)/2$, or omega lies outside the interval $\left(0.0,2.0\right)$,
${\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 $1\le {\mathbf{icol}}\left(i\right)\le {\mathbf{irow}}\left(i\right)$, for $i=1,2,\dots ,{\mathbf{nz}}$;
• ${\mathbf{irow}}\left(i-1\right)<{\mathbf{irow}}\left(i\right)$ or ${\mathbf{irow}}\left(i-1\right)={\mathbf{irow}}\left(i\right)$ and ${\mathbf{icol}}\left(i-1\right)<{\mathbf{icol}}\left(i\right)$, for $i=2,3,\dots ,{\mathbf{nz}}$.
Therefore a nonzero element has been supplied which does not lie in the lower triangular part of $A$, is out of order, or has duplicate row and column indices. Call nag_sparse_real_symm_sort (f11zb) to reorder and sum or remove duplicates.
${\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

The computed solution $x$ is the exact solution of a perturbed system of equations $\left(M+\delta M\right)x=y$, where
 $δM≤cnεD+ωLD-1D+ωLT,$
$c\left(n\right)$ is a modest linear function of $n$, and $\epsilon$ is the machine precision.

### Timing

The time taken for a call to nag_sparse_real_symm_precon_ssor_solve (f11jd) is proportional to nz.

### Use of check

It is expected that a common use of nag_sparse_real_symm_precon_ssor_solve (f11jd) will be to carry out the preconditioning step required in the application of nag_sparse_real_symm_basic_solver (f11ge) to sparse symmetric linear systems. In this situation nag_sparse_real_symm_precon_ssor_solve (f11jd) is likely to be called many times with the same matrix $M$. 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 a sparse symmetric linear system of equations
 $Ax=b,$
using the conjugate-gradient (CG) method with SSOR preconditioning.
The CG algorithm itself is implemented by the reverse communication function nag_sparse_real_symm_basic_solver (f11ge), which returns repeatedly to the calling program with various values of the argument irevcm. This argument indicates the action to be taken by the calling program.
• If ${\mathbf{irevcm}}=1$, a matrix-vector product $v=Au$ is required. This is implemented by a call to nag_sparse_real_symm_matvec (f11xe).
• If ${\mathbf{irevcm}}=2$, a solution of the preconditioning equation $Mv=u$ is required. This is achieved by a call to nag_sparse_real_symm_precon_ssor_solve (f11jd).
• If ${\mathbf{irevcm}}=4$, nag_sparse_real_symm_basic_solver (f11ge) has completed its tasks. Either the iteration has terminated, or an error condition has arisen.
For further details see the function document for nag_sparse_real_symm_basic_solver (f11ge).
```function f11jd_example

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

% Solve Ax = b by CG with SSOR preconditioning

% Matrix A and RHS b

n    = int64(7);
nz   = int64(16);
a = zeros(nz, 1);
irow = zeros(nz, 1, 'int64');
icol = irow;
a          = [4; 1; 5; 2; 2; 3;-1; 1; 4; 1;-2; 3; 2;-1;-2; 5];
irow(1:nz) = [1; 2; 2; 3; 4; 4; 5; 5; 5; 6; 6; 6; 7; 7; 7; 7];
icol(1:nz) = [1; 1; 2; 3; 2; 4; 1; 4; 5; 2; 5; 6; 1; 2; 3; 7];

b = [15; 18; -8; 21; 11; 10; 29];

% Initialise solver
method = 'CG';
precon = 'P';
tol    = 1e-6;
maxitn = int64(100);
anorm  = 0;
sigmax = 0;
maxits = int64(10);
monit = int64(0);
[lwreq, work, ifail] = ...
f11gd( ...
method, precon, n, tol, maxitn, anorm, sigmax, ...
maxits, monit, 'sigcmp', 'N', 'norm_p', 'I');

% SSOR preconditioner
% Calculate reciprocal diagonal matrix elements.
iwork = zeros(n+1, 1, 'int64');
rdiag = zeros(n, 1);
for i=1:nz
if irow(i) == icol(i)
iwork(irow(i)) = iwork(irow(i)) + 1;
if a(i) ~= 0
rdiag(irow(i)) = 1/a(i);
else
error('Matrix has a zero diagonal element');
end
end
end
for i=1:n
if iwork(i) == 0
error('Matrix has a missing diagonal elemen');
elseif iwork(i) > 2
error('Matrix has a multiple diagonal element');
end
end

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

% Preconditioner inputs
omega  = 1;
check  = 'C';

% Solve the linear system
while irevcm ~= 4
[irevcm, x, b, work, ifail] = f11ge( ...
irevcm, x, b, wgt, work);
if (irevcm == 1)
% Compute matrix vector product
[b, ifail] = f11xe( ...
a, irow, icol, 'N', x);
elseif (irevcm == 2)
% SSOR preconditioning
[b, ifail] = f11jd( ...
a, irow, icol, rdiag, omega, check, x);
end
end

% Get information about the computation
[itn, stplhs, stprhs, anorm, sigmax, its, sigerr, ifail] = ...
f11gf(work);
fprintf('Converged in %d iterations\n', itn);
fprintf('Final residual norm = %14.4e\n\n', stplhs);
disp('Solution');
disp(x);

```
```f11jd example results

Converged in 6 iterations
Final residual norm =     7.1054e-15

Solution
1.0000
2.0000
3.0000
4.0000
5.0000
6.0000
7.0000

```