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

## Purpose

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

## Syntax

[x, ifail] = f11dd(trans, a, irow, icol, rdiag, omega, check, y, 'n', n, 'nnz', nnz)
[x, ifail] = nag_sparse_real_gen_precon_ssor_solve(trans, a, irow, icol, rdiag, omega, check, y, 'n', n, 'nnz', nnz)

## Description

nag_sparse_real_gen_precon_ssor_solve (f11dd) solves a system of linear equations
 Mx = y, or MTx = y,
$Mx=y, or MTx=y,$
according to the value of the parameter trans, where the matrix
 M = 1/(ω(2 − ω))(D + ωL) D − 1 (D + ωU) $M=1ω(2-ω) (D+ω L) D-1 (D+ω U)$
corresponds to symmetric successive-over-relaxation (SSOR) (see Young (1971)) applied to a linear system Ax = b$Ax=b$, where A$A$ is a real sparse nonsymmetric matrix stored in coordinate storage (CS) format (see Section [Coordinate storage (CS) format] in the F11 Chapter Introduction).
In the definition of M$M$ given above D$D$ is the diagonal part of A$A$, L$L$ is the strictly lower triangular part of A$A$, U$U$ is the strictly upper triangular part of A$A$, and ω$\omega$ is a user-defined relaxation parameter.
It is envisaged that a common use of nag_sparse_real_gen_precon_ssor_solve (f11dd) will be to carry out the preconditioning step required in the application of nag_sparse_real_gen_basic_solver (f11be) to sparse linear systems. For an illustration of this use of nag_sparse_real_gen_precon_ssor_solve (f11dd) see the example program given in Section [Example]. nag_sparse_real_gen_precon_ssor_solve (f11dd) is also used for this purpose by the Black Box function nag_sparse_real_gen_solve_jacssor (f11de).

## References

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

## Parameters

### Compulsory Input Parameters

1:     trans – string (length ≥ 1)
Specifies whether or not the matrix M$M$ is transposed.
trans = 'N'${\mathbf{trans}}=\text{'N'}$
Mx = y$Mx=y$ is solved.
trans = 'T'${\mathbf{trans}}=\text{'T'}$
MTx = y${M}^{\mathrm{T}}x=y$ is solved.
Constraint: trans = 'N'${\mathbf{trans}}=\text{'N'}$ or 'T'$\text{'T'}$.
2:     a(nnz) – double array
nnz, the dimension of the array, must satisfy the constraint 1nnzn2$1\le {\mathbf{nnz}}\le {{\mathbf{n}}}^{2}$.
The nonzero elements in the matrix A$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_gen_sort (f11za) may be used to order the elements in this way.
3:     irow(nnz) – int64int32nag_int array
4:     icol(nnz) – int64int32nag_int array
nnz, the dimension of the array, must satisfy the constraint 1nnzn2$1\le {\mathbf{nnz}}\le {{\mathbf{n}}}^{2}$.
The row and column indices of the nonzero elements supplied in array a.
Constraints:
irow and icol must satisfy the following constraints (which may be imposed by a call to nag_sparse_real_gen_sort (f11za)):
• 1irow(i)n$1\le {\mathbf{irow}}\left(\mathit{i}\right)\le {\mathbf{n}}$ and 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}}$;
• 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}}$.
5:     rdiag(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The elements of the diagonal matrix D1${D}^{-1}$, where D$D$ is the diagonal part of A$A$.
6:     omega – double scalar
The relaxation parameter ω$\omega$.
Constraint: 0.0 < omega < 2.0$0.0<{\mathbf{omega}}<2.0$.
7:     check – string (length ≥ 1)
Specifies whether or not the CS representation of the matrix M$M$ should be checked.
check = 'C'${\mathbf{check}}=\text{'C'}$
Checks are carried on the values of n, nnz, irow, icol and omega.
check = 'N'${\mathbf{check}}=\text{'N'}$
None of these checks are carried out.
Constraint: check = 'C'${\mathbf{check}}=\text{'C'}$ or 'N'$\text{'N'}$.
8:     y(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The right-hand side vector y$y$.

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays rdiag, y. (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.)
The number of nonzero elements in the matrix A$A$.
Constraint: 1nnzn2$1\le {\mathbf{nnz}}\le {{\mathbf{n}}}^{2}$.

iwork

### Output Parameters

1:     x(n) – double array
The solution vector x$x$.
2:     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, trans ≠ 'N'${\mathbf{trans}}\ne \text{'N'}$ or 'T'$\text{'T'}$, or check ≠ 'C'${\mathbf{check}}\ne \text{'C'}$ or 'N'$\text{'N'}$.
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}$, or omega lies outside the interval (0.0,2.0)$\left(0.0,2.0\right)$,
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(\mathit{i}\right)\le {\mathbf{n}}$ and 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}}$;
• irow(i1) < irow(i)${\mathbf{irow}}\left(i-1\right)<{\mathbf{irow}}\left(i\right)$ or irow(i1) = irow(i)${\mathbf{irow}}\left(i-1\right)={\mathbf{irow}}\left(i\right)$ and icol(i1) < icol(i)${\mathbf{icol}}\left(i-1\right)<{\mathbf{icol}}\left(i\right)$, for i = 2,3,,nnz$i=2,3,\dots ,{\mathbf{nnz}}$.
Therefore a nonzero element has been supplied which does not lie in the matrix A$A$, is out of order, or has duplicate row and column indices. Call nag_sparse_real_gen_sort (f11za) to reorder and sum or remove duplicates.
ifail = 4${\mathbf{ifail}}=4$
On entry, the matrix A$A$ has a zero diagonal element. The SSOR preconditioner is not appropriate for this problem.

## Accuracy

If trans = 'N'${\mathbf{trans}}=\text{'N'}$ the computed solution x$x$ is the exact solution of a perturbed system of equations (M + δM)x = y$\left(M+\delta M\right)x=y$, where
 |δM| ≤ c(n)ε|D + ωL||D − 1||D + ωU|, $|δM|≤c(n)ε|D+ωL||D-1||D+ωU|,$
c(n)$c\left(n\right)$ is a modest linear function of n$n$, and ε$\epsilon$ is the machine precision. An equivalent result holds when trans = 'T'${\mathbf{trans}}=\text{'T'}$.

### Timing

The time taken for a call to nag_sparse_real_gen_precon_ssor_solve (f11dd) is proportional to nnz.

### Use of check

It is expected that a common use of nag_sparse_real_gen_precon_ssor_solve (f11dd) will be to carry out the preconditioning step required in the application of nag_sparse_real_gen_basic_solver (f11be) to sparse linear systems. In this situation nag_sparse_real_gen_precon_ssor_solve (f11dd) is likely to be called many times with the same matrix M$M$. 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 for all subsequent calls set check = 'N'${\mathbf{check}}=\text{'N'}$.

## Example

function nag_sparse_real_gen_precon_ssor_solve_example
n = int64(5);
m = int64(2);
nz = int64(16);
method = 'rgmres';
precon = 'P';
iterm = int64(1);
tol = 1e-10;
maxitn = int64(1000);
anorm = 0;
sigmax = 0;
trans = 'N';
omega = 1.1;
check = 'C';
weight = 'N';
monit = int64(0);
lwork = max([n*(m+3)+m*(m+5)+101,7*n+100,(2*n+m)*(m+2)+n+100,10*n+100]);
a = zeros(3*nz, 1);
a(1:nz) = [2; 1; -1; -3; -2; 1; 1; 5; 3; 1; -2; -3; -1; 4; -2; -6];
irow = zeros(3*nz, 1, 'int64');
irow(1:nz) = [1; 1; 1; 2; 2; 2; 3; 3; 3; 3; 4; 4; 4; 5; 5; 5];
icol = zeros(3*nz, 1, 'int64');
icol(1:nz) = [1; 2; 4; 2; 3; 5; 1; 3; 4; 5; 1; 4; 5; 2; 3; 5];
b = [0; -7; 33; -19; -28];
x = zeros(n, 1);
iwork = zeros(2*n + 1, 1, 'int64');
rdiag = zeros(n, 1);
wgt = zeros(n, 1);
irevcm = int64(0);
ckxaf = 'C';
ckddf = 'C';

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

% Calculate reciprocal diagonal matrix elements if necessary
if strcmpi(precon, 'P')

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 element');
elseif iwork(i) >= 2
error('Matrix has a multiple diagonal element');
end
end

end

% Solve the linear system
while irevcm ~= 4
[irevcm, x, b, work, ifail] = ...
nag_sparse_real_gen_basic_solver(irevcm, x, b, wgt, work);

if (irevcm == -1)
% Compute transposed matrix-vector product
[b, ifail] = nag_sparse_real_gen_matvec('T', a(1:nz), irow(1:nz), ...
icol(1:nz), ckxaf, x);
ckxaf = 'N';
elseif (irevcm == 1)
% Compute matrix-vector product
[b, ifail] = nag_sparse_real_gen_matvec('N', a(1:nz), irow(1:nz), ...
icol(1:nz), ckxaf, x);
ckxaf = 'N';
elseif (irevcm == 2)
% SSOR preconditioning
[b, ifail] = nag_sparse_real_gen_precon_ssor_solve(trans, a(1:nz), ...
irow(1:nz), icol(1:nz), rdiag, omega, ckddf, x);
ckddf = 'N';
end
end

% Get information about the computation
[itn, stplhs, stprhs, anorm, sigmax, ifail] = ...
nag_sparse_real_gen_basic_diag(work);
fprintf('\nConverged in %d iterations\n', itn);
fprintf('Matrix norm = %16.3e\n', anorm);
fprintf('Final residual norm = %16.3e\n', stplhs);
fprintf('      X\n');
disp(x);

Converged in 12 iterations
Matrix norm =        1.200e+01
Final residual norm =        3.841e-09
X
1.0000
2.0000
3.0000
4.0000
5.0000

function f11dd_example
n = int64(5);
m = int64(2);
nz = int64(16);
method = 'rgmres';
precon = 'P';
iterm = int64(1);
tol = 1e-10;
maxitn = int64(1000);
anorm = 0;
sigmax = 0;
trans = 'N';
omega = 1.1;
check = 'C';
weight = 'N';
monit = int64(0);
lwork = max([n*(m+3)+m*(m+5)+101,7*n+100,(2*n+m)*(m+2)+n+100,10*n+100]);
a = zeros(3*nz, 1);
a(1:nz) = [2; 1; -1; -3; -2; 1; 1; 5; 3; 1; -2; -3; -1; 4; -2; -6];
irow = zeros(3*nz, 1, 'int64');
irow(1:nz) = [1; 1; 1; 2; 2; 2; 3; 3; 3; 3; 4; 4; 4; 5; 5; 5];
icol = zeros(3*nz, 1, 'int64');
icol(1:nz) = [1; 2; 4; 2; 3; 5; 1; 3; 4; 5; 1; 4; 5; 2; 3; 5];
b = [0; -7; 33; -19; -28];
x = zeros(n, 1);
iwork = zeros(2*n + 1, 1, 'int64');
rdiag = zeros(n, 1);
wgt = zeros(n, 1);
irevcm = int64(0);
ckxaf = 'C';
ckddf = 'C';

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

% Calculate reciprocal diagonal matrix elements if necessary
if strcmpi(precon, 'P')

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 element');
elseif iwork(i) >= 2
error('Matrix has a multiple diagonal element');
end
end

end

% Solve the linear system
while irevcm ~= 4
[irevcm, x, b, work, ifail] = f11be(irevcm, x, b, wgt, work);

if (irevcm == -1)
% Compute transposed matrix-vector product
[b, ifail] = f11xa('T', a(1:nz), irow(1:nz), icol(1:nz), ckxaf, x);
ckxaf = 'N';
elseif (irevcm == 1)
% Compute matrix-vector product
[b, ifail] = f11xa('N', a(1:nz), irow(1:nz), icol(1:nz), ckxaf, x);
ckxaf = 'N';
elseif (irevcm == 2)
% SSOR preconditioning
[b, ifail] = f11dd(trans, a(1:nz), irow(1:nz), icol(1:nz), rdiag,...
omega, ckddf, x);
ckddf = 'N';
end
end

% Get information about the computation
[itn, stplhs, stprhs, anorm, sigmax, ifail] = f11bf(work);
fprintf('\nConverged in %d iterations\n', itn);
fprintf('Matrix norm = %16.3e\n', anorm);
fprintf('Final residual norm = %16.3e\n', stplhs);
fprintf('      X\n');
disp(x);

Converged in 12 iterations
Matrix norm =        1.200e+01
Final residual norm =        3.841e-09
X
1.0000
2.0000
3.0000
4.0000
5.0000