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_ssor_solve (f11dr)

## Purpose

nag_sparse_complex_gen_precon_ssor_solve (f11dr) solves a system of linear equations involving the preconditioning matrix corresponding to SSOR applied to a complex sparse non-Hermitian matrix, represented in coordinate storage format.

## Syntax

[x, ifail] = f11dr(trans, a, irow, icol, rdiag, omega, check, y, 'n', n, 'nz', nz)
[x, ifail] = nag_sparse_complex_gen_precon_ssor_solve(trans, a, irow, icol, rdiag, omega, check, y, 'n', n, 'nz', nz)

## Description

nag_sparse_complex_gen_precon_ssor_solve (f11dr) solves a system of linear equations
 $Mx=y, or MHx=y,$
according to the value of the argument trans, where the matrix
 $M=1ω2-ω D+ω L D-1 D+ω U$
corresponds to symmetric successive-over-relaxation (SSOR) Young (1971) applied to a linear system $Ax=b$, where $A$ is a complex sparse non-Hermitian matrix stored in coordinate storage (CS) format (see Coordinate storage (CS) 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$, $U$ is the strictly upper triangular part of $A$, and $\omega$ is a user-defined relaxation parameter.
It is envisaged that a common use of nag_sparse_complex_gen_precon_ssor_solve (f11dr) will be to carry out the preconditioning step required in the application of nag_sparse_complex_gen_basic_solver (f11bs) to sparse linear systems. For an illustration of this use of nag_sparse_complex_gen_precon_ssor_solve (f11dr) see the example program given in Example. nag_sparse_complex_gen_precon_ssor_solve (f11dr) is also used for this purpose by the Black Box function nag_sparse_complex_gen_solve_jacssor (f11ds).

## References

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

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{trans}$ – string (length ≥ 1)
Specifies whether or not the matrix $M$ is transposed.
${\mathbf{trans}}=\text{'N'}$
$Mx=y$ is solved.
${\mathbf{trans}}=\text{'T'}$
${M}^{\mathrm{H}}x=y$ is solved.
Constraint: ${\mathbf{trans}}=\text{'N'}$ or $\text{'T'}$.
2:     $\mathrm{a}\left({\mathbf{nz}}\right)$ – complex array
The nonzero elements in 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_complex_gen_sort (f11zn) may be used to order the elements in this way.
3:     $\mathrm{irow}\left({\mathbf{nz}}\right)$int64int32nag_int array
4:     $\mathrm{icol}\left({\mathbf{nz}}\right)$int64int32nag_int array
The row and column indices of the nonzero elements supplied in a.
Constraints:
irow and icol must satisfy the following constraints (which may be imposed by a call to nag_sparse_complex_gen_sort (f11zn)):
• $1\le {\mathbf{irow}}\left(\mathit{i}\right)\le {\mathbf{n}}$ and $1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{n}}$, 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}}$.
5:     $\mathrm{rdiag}\left({\mathbf{n}}\right)$ – complex array
The elements of the diagonal matrix ${D}^{-1}$, where $D$ is the diagonal part of $A$.
6:     $\mathrm{omega}$ – double scalar
The relaxation parameter $\omega$.
Constraint: $0.0<{\mathbf{omega}}<2.0$.
7:     $\mathrm{check}$ – string (length ≥ 1)
Specifies whether or not the CS representation of the matrix $M$ should be checked.
${\mathbf{check}}=\text{'C'}$
Checks are carried 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'}$.
8:     $\mathrm{y}\left({\mathbf{n}}\right)$ – complex 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 matrix $A$.
Constraint: $1\le {\mathbf{nz}}\le {{\mathbf{n}}}^{2}$.

### Output Parameters

1:     $\mathrm{x}\left({\mathbf{n}}\right)$ – complex 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{trans}}\ne \text{'N'}$ or $\text{'T'}$, or ${\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}}}^{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(\mathit{i}\right)\le {\mathbf{n}}$ and $1\le {\mathbf{icol}}\left(\mathit{i}\right)\le {\mathbf{n}}$, for $\mathit{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 matrix $A$, is out of order, or has duplicate row and column indices. Call nag_sparse_complex_gen_sort (f11zn) to reorder and sum or remove duplicates.
${\mathbf{ifail}}=4$
On entry, the matrix $A$ has a zero diagonal element. The SSOR preconditioner is not appropriate for this problem.
${\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

If ${\mathbf{trans}}=\text{'N'}$ 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+ωU,$
$c\left(n\right)$ is a modest linear function of $n$, and $\epsilon$ is the machine precision. An equivalent result holds when ${\mathbf{trans}}=\text{'T'}$.

### Timing

The time taken for a call to nag_sparse_complex_gen_precon_ssor_solve (f11dr) is proportional to nz.

### Use of check

It is expected that a common use of nag_sparse_complex_gen_precon_ssor_solve (f11dr) will be to carry out the preconditioning step required in the application of nag_sparse_complex_gen_basic_solver (f11bs) to sparse linear systems. In this situation nag_sparse_complex_gen_precon_ssor_solve (f11dr) 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 ${\mathbf{check}}=\text{'N'}$ for all subsequent calls.

## Example

This example solves a complex sparse linear system of equations
 $Ax=b,$
using RGMRES with SSOR preconditioning.
The RGMRES algorithm itself is implemented by the reverse communication function nag_sparse_complex_gen_basic_solver (f11bs), 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_complex_gen_matvec (f11xn).
• If ${\mathbf{irevcm}}=-1$, a conjugate transposed matrix-vector product $v={A}^{\mathrm{H}}u$ is required in the estimation of the norm of $A$. This is implemented by a call to nag_sparse_complex_gen_matvec (f11xn).
• If ${\mathbf{irevcm}}=2$, a solution of the preconditioning equation $Mv=u$ is required. This is achieved by a call to nag_sparse_complex_gen_precon_ssor_solve (f11dr).
• If ${\mathbf{irevcm}}=4$, nag_sparse_complex_gen_basic_solver (f11bs) 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_complex_gen_basic_solver (f11bs).
```function f11dr_example

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

% Solve sparse system Ax = b using preconditioned CGS

% Sparse A and b
n    = int64(5);
m    = int64(2);
nz   = int64(16);
sparseA = [ 2  3   1    1;
1 -1   1    2;
-1  0   1    4;
0  2   2    2;
-2  1   2    3;
1  0   2    5;
0 -1   3    1;
5  4   3    3;
3 -1   3    4;
1  0   3    5;
-2  2   4    1;
-3  1   4    4;
0  3   4    5;
4 -2   5    2;
-2  0   5    3;
-6  1   5    5];
a(1:nz)    = sparseA(:,1) + i*sparseA(:,2);
irow(1:nz) = int64(sparseA(:,3));
icol(1:nz) = int64(sparseA(:,4));
b = [ -3 +  3i;
-11 +  5i;
23 + 48i;
-41 +  2i;
-28 - 31i];

% Solver setup input argument initialization
method = 'CGS';
precon = 'P';
norm_p = 'I';
tol    = 1e-10;
maxitn = int64(1000);
anorm  = 0;
sigmax = 0;
monit  = int64(0);
lwork  = max([121+n*(3+m)+m*(m+5),120+7*n,120+(2*n+m)*(m+2)+2*n,120+10*n]);

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

rdiag = zeros(n, 1);
% Calculate reciprocal diagonal matrix elements for preconditioner
if precon == 'P' || precon == 'p'
dcount = zeros(n, 1, 'int64');

for j = 1:nz
if irow(j) == icol(j)
dcount(irow(j)) = dcount(irow(j)) + 1;
if a(j) ~= 0
rdiag(irow(j)) = 1/a(j);
else
error('Matrix has a zero diagonal element');
end
end
end

for j = 1:n
if dcount(j)==0
error('Matrix has a missing diagonal element');
elseif dcount(j) > 1
error('Matrix has a multiple diagonal element');
end
end
end
% Other preconditioner inputs
omega = 1.4;
ckdrf = 'C';

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

% matrix-vector input
ckxnf = 'C';

% solve the linear system by revere communication
while (irevcm ~= 4)
[irevcm, x, b, work, ifail] = f11bs( ...
irevcm, x, b, wgt, work);

if (irevcm == -1)
% b = A^T x
[b, ifail] = f11xn( ...
'T', a, irow, icol, ckxnf, x);
ckxnf = 'N';
elseif (irevcm == 1)
% b = Ax
[b, ifail] = f11xn( ...
'N', a, irow, icol, ckxnf, x);
ckxnf = 'N';
elseif (irevcm == 2)
% SSOR preconditioning
[b, ifail] = f11dr( ...
'N', a, irow, icol, rdiag, omega, ckdrf, x);
ckdrf = 'N';
end
end

if ifail == 0
% Successful termination, get solution details from communication array
[itn, stplhs, stprhs, anorm, sigmax, work, ifail] = ...
f11bt(work);

fprintf('Converged in %d iterations\n', itn);
fprintf('Matrix norm         = %16.3e\n', anorm);
fprintf('Final residual norm = %16.4e\n\n', stplhs);
disp('Solution');
disp(x);

end

```
```f11dr example results

Converged in 5 iterations
Matrix norm         =        1.500e+01
Final residual norm =       4.6185e-14

Solution
1.0000 + 2.0000i
2.0000 + 3.0000i
3.0000 + 4.0000i
4.0000 + 5.0000i
5.0000 + 6.0000i

```

Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015