NAG Toolbox: nag_sparse_real_symm_precon_ssor_solve (f11jd)

Purpose

Syntax

Description

References

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.

See also Use of check.

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}\times \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

   $\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}\times \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$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Timing

Use of check

Example

Open in the MATLAB editor: f11jd_example

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