NAG Toolbox: nag_sparse_complex_gen_precon_jacobi (f11dx)

Purpose

Syntax

Description

References

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 Hermitian matrix $A$) or coordinate storage (CS) (applicable to both Hermitian and non-Hermitian matrices).

$\mathbf{store}=\text{'N'}$

The complete matrix $A$ is stored in CS format.

$\mathbf{store}=\text{'S'}$

The lower triangle of the Hermitian 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 Hermitian 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{H}}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{H}}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)$ – complex 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_complex_gen_sort (f11zn) or nag_sparse_complex_herm_sort (f11zp) 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.

See also Use of check.

Constraint: $\mathbf{check}=\text{'C'}$ or $\text{'N'}$.

9:     $\mathrm{b}\left(\mathbf{n}\right)$ – complex array

The right-hand side vector $b$.

10:   $\mathrm{diag}\left(\mathbf{n}\right)$ – complex 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}\times \left(\mathbf{n}+1\right)/2$.

Output Parameters

1:     $\mathrm{x}\left(\mathbf{n}\right)$ – complex array

The approximate solution vector $x_{\mathbf{niter}}$.

2:     $\mathrm{diag}\left(\mathbf{n}\right)$ – complex 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

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

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: f11dx_example

function f11dx_example


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

% Solve sparse system Ax = b iteratively using Jacobi preconditioner

% Define sparse matrix A and RHS B
nz   = int64(24);
n    = int64(8);
a = [ 2 + 1i, -1 + 1i,  1 - 3i,  4 + 7i, -3 + 0i,  2 + 4i, ...
     -7 - 5i,  2 + 1i,  3 + 2i, -4 + 2i,  0 + 1i,  5 - 3i, ...
     -1 + 2i,  8 + 6i, -3 - 4i, -6 - 2i,  5 - 2i,  2 + 0i, ...
      0 - 5i, -1 + 5i,  6 + 2i, -1 + 4i,  2 + 0i,  3 + 3i];
b = [ 7 + 11i;
      1 + 24i;
    -13 - 18i;
    -10 +  3i;
     23 + 14i;
     17 -  7i;
     15 -  3i;
     -3 + 20i];
irow = 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 = 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]);

b = [  7 + 11i;
       1 + 24i;
     -13 - 18i;
     -10 +  3i;
      23 + 14i;
      17 -  7i;
      15 -  3i;
      -3 + 20i];


% Iterative Solution Setup
%      Input parameters
method = 'TFQMR';
precon = 'p';
m      = int64(2);
tol    = 1e-6;
maxitn = int64(20);
anorm  = 0;
sigmax = 0;
monit  = int64(1);
lwork = int64(300);

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

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

% Preconditioner inputs
store = 'Non Hermitian';
trans = 'N';
init  = 'I';
niter = int64(2);
check = 'Check';
diag  = complex(zeros(n, 1));

% Solve by reverse communication
while (irevcm ~= 4)
  [irevcm, x, b, work, ifail] = f11bs( ...
                                       irevcm, x, b, wgt, work);

  if (irevcm == -1)
    % b = A^Tx
    [b, ifail] = f11xn('T', a, irow, icol, 'N', x);
  elseif (irevcm == 1)
    % b = Ax
    [b, ifail] = f11xn('N', a, irow, icol, 'N', x);
  elseif (irevcm == 2)
    % Jacobi preconditioner
    [b, diag, ifail] = f11dx( ...
                              store, trans, init, niter, ...
                              a, irow, icol, check, x, diag);
  elseif (irevcm == 3)
    [itn, stplhs, stprhs, anorm, sigmax, work, ifail] = ...
    f11bt(work);
    fprintf('\nMonitoring at iteration number %d\n', itn);
    fprintf('residual norm: %14.4e\n', stplhs);
  end
end

% Obtain information about the computation
[itn, stplhs, stprhs, anorm, sigmax, work, ifail] = ...
f11bt(work);
fprintf('\nNumber of iterations for convergence: %d\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\n');
disp(x);
fprintf('\n   Residual Vector\n');
disp(b);

f11dx example results


Monitoring at iteration number 1
residual norm:     1.5062e+02

Monitoring at iteration number 2
residual norm:     1.5704e+02

Monitoring at iteration number 3
residual norm:     1.4803e+02

Monitoring at iteration number 4
residual norm:     8.5215e+01

Monitoring at iteration number 5
residual norm:     4.2951e+01

Monitoring at iteration number 6
residual norm:     2.5055e+01

Monitoring at iteration number 7
residual norm:     1.9090e-01

Number of iterations for convergence: 8
Residual norm:                               9.5485e-08
Right-hand side of termination criteria:     8.9100e-04
i-norm of matrix a:                          2.7000e+01

   Solution Vector
   1.0000 + 1.0000i
   2.0000 - 1.0000i
   3.0000 + 1.0000i
   4.0000 - 1.0000i
   3.0000 - 1.0000i
   2.0000 + 1.0000i
   1.0000 - 1.0000i
   0.0000 + 3.0000i


   Residual Vector
   1.0e-07 *

   0.0471 - 0.0394i
   0.0877 - 0.0981i
  -0.0287 + 0.0416i
   0.0358 - 0.1112i
  -0.0692 - 0.0869i
  -0.0225 + 0.0849i
   0.0052 - 0.0334i
  -0.0558 - 0.1073i