NAG Toolbox: nag_sparse_real_gen_precon_jacobi (f11dk)

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

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

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

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

The lower triangle of the symmetric 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 symmetric 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{T}}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{T}}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)$ – double 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_real_gen_sort (f11za) or nag_sparse_real_symm_sort (f11zb) 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)$ – double array

The right-hand side vector $b$.

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

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

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

function f11dk_example


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

% Solve sparse system Ax = b iteratively using Jacobi preconditioner

% Define sparse matrix A and RHS B
n    = int64(8);
nz   = int64(24);
irow(1:nz) = 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(1:nz) = 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]);
a(1:nz)    = [4;-1; 1; 4;-5; 2;-7; 2; 2;-1; 6; 2;
             -1; 8;-2;-2; 5; 8;-2;-1; 7;-1; 2; 6];
b          = [6; 8;-9;46;17;21;22;34];

% Iterative Solution Setup
%      Input parameters
method = 'BICGSTAB';
precon = 'Preconditioned';
lpoly  = int64(2);
tol    = sqrt(x02aj)^(3/8);
maxitn = int64(20);
anorm  = 0;
sigmax = 0;
monit  = int64(1);
lwork  = int64(200);

% Setup iterative method
[lwreq, work, ifail] = ...
  f11bd(method, precon, n, lpoly, tol, maxitn, anorm, sigmax, ...
        monit, lwork, 'norm_p', '1');

% Preconditioner parameters
store = 'Non symmetric';
trans = 'N';
init  = 'I';
niter = int64(4);
wgt   = zeros(n, 1);
dgnl  = zeros(n, 1);
x     = zeros(n, 1);

% Repeatedly call f11be to solve the equations.
% On final exit x will contain the solution and b the residual vector.
% Jacobi preconditioner is via f11dk when irevcm==2.

irevcm = int64(0);
init = 'i';
ifail = int64(0);

while irevcm ~= 4
  [irevcm, x, b, work, ifail] = ...
  f11be(...
        irevcm, x, b, wgt, work);

  if (irevcm == -1)
    % b = A^Tx
    [b, ifail1] = f11xa(...
        'Transpose', a, irow, icol, 'N', x);
  elseif (irevcm == 1)
    % b = Ax
    [b, ifail1] = f11xa(...
        'No Transpose', a, irow, icol, 'N', x);
  elseif (irevcm == 2)
    % Jacobi preconditioning Jb = x
    [b, dgnl, ifail1] = f11dk(...
        store, trans, init, niter, a(1:nz), irow, icol, 'Check', x, dgnl);
    init = 'n';
  elseif (irevcm == 3)
    % Monitoring
    [itn, stplhs, stprhs, anorm, sigmax, ifail1] = ...
        f11bf(work);
    fprintf('\nMonitoring at iteration number %2d\n',itn);
    fprintf('residual norm:              %14.4e\n', stplhs);
  end

end

% Get information about the computation
[itn, stplhs, stprhs, anorm, sigmax, ifail] = ...
    f11bf(work);

fprintf('\nNumber of iterations for convergence:     %4d\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  Residual Vector\n');
for i = 1:n
  fprintf('%16.4f %16.2e\n', x(i), b(i));
end

f11dk example results


Number of iterations for convergence:        2
Residual norm:                               1.1177e-04
Right-hand side of termination criteria:     5.5153e-01
i-norm of matrix a:                          1.5000e+01

   Solution Vector  Residual Vector
          1.7035         3.24e-07
          1.0805        -1.76e-05
          1.8305         2.80e-05
          6.0251        -2.59e-05
          3.2942         7.82e-06
          1.9068         9.21e-06
          4.1365        -3.08e-06
          5.2111         1.98e-05