hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_sparse_real_gen_precon_jacobi (f11dk)

Purpose

nag_sparse_real_gen_precon_jacobi (f11dk) computes the approximate solution of a real, symmetric or nonsymmetric, sparse system of linear equations applying a number of Jacobi iterations. It is expected that nag_sparse_real_gen_precon_jacobi (f11dk) will be used as a preconditioner for the iterative solution of real sparse systems of equations.

Syntax

[x, diag, ifail] = f11dk(store, trans, init, niter, a, irow, icol, check, b, diag, 'n', n, 'nnz', nnz)
[x, diag, ifail] = nag_sparse_real_gen_precon_jacobi(store, trans, init, niter, a, irow, icol, check, b, diag, 'n', n, 'nnz', nnz)

Description

nag_sparse_real_gen_precon_jacobi (f11dk) computes the approximate solution of the real sparse system of linear equations Ax = bAx=b using niter iterations of the Jacobi algorithm (see also Golub and Van Loan (1996) and Young (1971)):
xk + 1 = xk + D1(bAxk)
xk+1=xk+D-1(b-Axk)
(1)
where k = 1,2,,niterk=1,2,,niter and x0 = 0x0=0.
nag_sparse_real_gen_precon_jacobi (f11dk) can be used both for nonsymmetric and symmetric systems of equations. For symmetric matrices, either all nonzero elements of the matrix AA can be supplied using coordinate storage (CS), or only the nonzero elements of the lower triangle of AA, using symmetric coordinate storage (SCS) (see the F11 Chapter Introduction).
It is expected that nag_sparse_real_gen_precon_jacobi (f11dk) will be used as a preconditioner for the iterative solution of real sparse systems of equations, using either the suite comprising the functions nag_sparse_real_symm_basic_setup (f11gd), nag_sparse_real_symm_basic_solver (f11ge) and nag_sparse_real_symm_basic_diag (f11gf), for symmetric systems, or the suite comprising the functions nag_sparse_real_gen_basic_setup (f11bd), nag_sparse_real_gen_basic_solver (f11be) and nag_sparse_real_gen_basic_diag (f11bf), for nonsymmetric systems of equations.

References

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Young D (1971) Iterative Solution of Large Linear Systems Academic Press, New York

Parameters

Compulsory Input Parameters

1:     store – string (length ≥ 1)
Specifies whether the matrix AA is stored using symmetric coordinate storage (SCS) (applicable only to a symmetric matrix AA) or coordinate storage (CS) (applicable to both symmetric and nonsymmetric matrices).
store = 'N'store='N'
The complete matrix AA is stored in CS format.
store = 'S'store='S'
The lower triangle of the symmetric matrix AA is stored in SCS format.
Constraint: store = 'N'store='N' or 'S''S'.
2:     trans – string (length ≥ 1)
If store = 'N'store='N', specifies whether the approximate solution of Ax = bAx=b or of ATx = bATx=b is required.
trans = 'N'trans='N'
The approximate solution of Ax = bAx=b is calculated.
trans = 'T'trans='T'
The approximate solution of ATx = bATx=b is calculated.
Constraint: trans = 'N'trans='N' or 'T''T'.
3:     init – string (length ≥ 1)
On first entry, init should be set to 'I', unless the diagonal elements of AA are already stored in the array diag. If diag already contains the diagonal of AA, it must be set to 'N'.
init = 'N'init='N'
diag must contain the diagonal of AA.
init = 'I'init='I'
diag will store the diagonal of AA on exit.
Constraint: init = 'N'init='N' or 'I''I'.
4:     niter – int64int32nag_int scalar
The number of Jacobi iterations requested.
Constraint: niter1niter1.
5:     a(nnz) – double array
nnz, the dimension of the array, must satisfy the constraint
  • if store = 'N'store='N', 1nnzn21nnzn2;
  • if store = 'S'store='S', 1nnzn × (n + 1) / 21nnzn×(n+1)/2.
If store = 'N'store='N', the nonzero elements in the matrix AA (CS format).
If store = 'S'store='S', the nonzero elements in the lower triangle of the matrix AA (SCS format).
In both cases, the elements of either AA 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:     irow(nnz) – int64int32nag_int array
7:     icol(nnz) – int64int32nag_int array
nnz, the dimension of the array, must satisfy the constraint
  • if store = 'N'store='N', 1nnzn21nnzn2;
  • if store = 'S'store='S', 1nnzn × (n + 1) / 21nnzn×(n+1)/2.
If store = 'N'store='N', the row and column indices of the nonzero elements supplied in a.
If store = 'S'store='S', the row and column indices of the nonzero elements of the lower triangle of the matrix AA supplied in a.
Constraints:
  • 1irow(i)n1irowin, for i = 1,2,,nnzi=1,2,,nnz;
  • if store = 'N'store='N', 1icol(i)n1icolin, for i = 1,2,,nnzi=1,2,,nnz;
  • if store = 'S'store='S', 1icol(i)irow(i)1icoliirowi, for i = 1,2,,nnzi=1,2,,nnz;
  • either irow(i1) < irow(i)irowi-1<irowi or both irow(i1) = irow(i)irowi-1=irowi and icol(i1) < icol(i)icoli-1<icoli, for i = 2,3,,nnzi=2,3,,nnz.
8:     check – string (length ≥ 1)
Specifies whether or not the CS or SCS representation of the matrix AA should be checked.
check = 'C'check='C'
Checks are carried out on the values of n, nnz, irow, icol; if init = 'N'init='N', diag is also checked.
check = 'N'check='N'
None of these checks are carried out.
Constraint: check = 'C'check='C' or 'N''N'.
9:     b(n) – double array
n, the dimension of the array, must satisfy the constraint n1n1.
The right-hand side vector bb.
10:   diag(n) – double array
n, the dimension of the array, must satisfy the constraint n1n1.
If init = 'N'init='N', the diagonal elements of AA.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays b, diag. (An error is raised if these dimensions are not equal.)
nn, the order of the matrix AA.
Constraint: n1n1.
2:     nnz – int64int32nag_int scalar
Default: The dimension of the arrays a, irow, icol. (An error is raised if these dimensions are not equal.)
If store = 'N'store='N', the number of nonzero elements in the matrix AA.
If store = 'S'store='S', the number of nonzero elements in the lower triangle of the matrix AA.
Constraints:
  • if store = 'N'store='N', 1nnzn21nnzn2;
  • if store = 'S'store='S', 1nnzn × (n + 1) / 21nnzn×(n+1)/2.

Input Parameters Omitted from the MATLAB Interface

work

Output Parameters

1:     x(n) – double array
The approximate solution vector xniterxniter.
2:     diag(n) – double array
If init = 'N'init='N', unchanged on exit.
If init = 'I'init='I', the diagonal elements of AA.
3:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,store'N'store'N' or 'S''S',
ortrans'N'trans'N' or 'T''T',
orinit'N'init'N' or 'I''I',
orcheck'C'check'C' or 'N''N',
orniter0niter0.
  ifail = 2ifail=2
On entry,n < 1n<1,
ornnz < 1nnz<1,
ornnz > n2nnz>n2, if store = 'N'store='N',
or1nnz[n(n + 1)] / 21nnz[nn+1]/2, if store = 'S'store='S'.
  ifail = 3ifail=3
On entry, the arrays irow and icol fail to satisfy the following constraints:
  • 1irow(i)n1irowin and
    • if store = 'N'store='N' then 1icol(i)n1icolin, or
    • if store = 'S'store='S' then 1icol(i)irow(i)1icoliirowi, for i = 1,2,,nnzi=1,2,,nnz.
  • irow(i1) < irow(i)irowi-1<irowi or irow(i1) = irow(i)irowi-1=irowi and icol(i1) < icol(i)icoli-1<icoli, for i = 2,3,,nnzi=2,3,,nnz.
Therefore a nonzero element has been supplied which does not lie within the matrix AA, 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 store = 'N'store='N' or store = 'S'store='S', respectively.
  ifail = 4ifail=4
On entry, init = 'N'init='N' and some diagonal elements of AA stored in diag are zero.
  ifail = 5ifail=5
On entry, init = 'I'init='I' and some diagonal elements of AA are zero.

Accuracy

In general, the Jacobi method cannot be used on its own to solve systems of linear equations. The rate of convergence is bound by its spectral properties (see, for example, Golub and Van Loan (1996)) and as a solver, the Jacobi method can only be applied to a limited set of matrices. One condition that guarantees convergence is strict diagonal dominance.
However, the Jacobi method can be used successfully as a preconditioner to a wider class of systems of equations. The Jacobi method has good vector/parallel properties, hence it can be applied very efficiently. Unfortunately, it is not possible to provide criteria which define the applicability of the Jacobi method as a preconditioner, and its usefulness must be judged for each case.

Further Comments

Timing

The time taken for a call to nag_sparse_real_gen_precon_jacobi (f11dk) is proportional to niter × nnzniter×nnz.

Use of check

It is expected that a common use of nag_sparse_real_gen_precon_jacobi (f11dk) will be as preconditioner for the iterative solution of real, symmetric or nonsymmetric, linear systems. In this situation, nag_sparse_real_gen_precon_jacobi (f11dk) is likely to be called many times. In the interests of both reliability and efficiency, you are recommended to set check = 'C'check='C' for the first of such calls, and to set check = 'N'check='N' for all subsequent calls.

Example

function nag_sparse_real_gen_precon_jacobi_example
store = 'Non symmetric';
trans = 'N';
init = 'I';
niter = int64(4);
n = int64(8);
wgt = zeros(n, 1);
dgnl = zeros(n, 1);
lwork = int64(200);
method = 'bicgstab';
precon = 'p';
m = int64(2);
tol = 1e-6;
maxitn = int64(20);
monit = int64(1);
niter = int64(4);
a = [4;-1;1;4;-5;2;-7;2;2;-1;6;2;-1;8;-2;-2;5;8;-2;-1;7;-1;2;6];
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 = [6; 8; -9; 46; 17; 21; 22; 34];
x = zeros(n, 1);
anorm = 0;
sigmax = 0;


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

% Repeatedly call nag_sparse_real_gen_basic_solver to solve the equations.  On final exit x will
% contain the solution and b the residual vector.

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

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

  if (irevcm == -1)
    [b, ifail1] = nag_sparse_real_gen_matvec('Transpose', a, irow, icol, 'N', x);
  elseif (irevcm == 1)
    [b, ifail1] = nag_sparse_real_gen_matvec('No Transpose', a, irow, icol, 'N', x);
  elseif (irevcm == 2)
    [b, dgnl, ifail1] = nag_sparse_real_gen_precon_jacobi(store, trans, init, niter, a, ...
                             irow, icol, 'Check', x, dgnl);
    init = 'n';
  elseif (irevcm == 3)
    [itn, stplhs, stprhs, anorm, sigmax, ifail1] = nag_sparse_real_gen_basic_diag(work);
    fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', itn, stplhs);
  end

  if (ifail1 ~= 0)
    irevcm = 6;
  end
end

% Get information about the computation
if ifail == 0
  [itn, stplhs, stprhs, anorm, sigmax, ifail] = nag_sparse_real_gen_basic_diag(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  Residual Vector\n');
  for i = 1:n
    fprintf('%+16.4e %+16.4e\n', x(i), b(i));
  end
end
 

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

   Solution Vector  Residual Vector
     +1.7035e+00      +3.2377e-07
     +1.0805e+00      -1.7625e-05
     +1.8305e+00      +2.7964e-05
     +6.0251e+00      -2.5914e-05
     +3.2942e+00      +7.8156e-06
     +1.9068e+00      +9.2064e-06
     +4.1365e+00      -3.0848e-06
     +5.2111e+00      +1.9834e-05

function f11dk_example
store = 'Non symmetric';
trans = 'N';
init = 'I';
niter = int64(4);
n = int64(8);
wgt = zeros(n, 1);
dgnl = zeros(n, 1);
lwork = int64(200);
method = 'bicgstab';
precon = 'p';
m = int64(2);
tol = 1e-6;
maxitn = int64(20);
monit = int64(1);
niter = int64(4);
a = [4;-1;1;4;-5;2;-7;2;2;-1;6;2;-1;8;-2;-2;5;8;-2;-1;7;-1;2;6];
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 = [6; 8; -9; 46; 17; 21; 22; 34];
x = zeros(n, 1);
anorm = 0;
sigmax = 0;


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

% Repeatedly call f11be to solve the equations.  On final exit x will
% contain the solution and b the residual vector.

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, ifail1] = f11xa('Transpose', a, irow, icol, 'N', x);
  elseif (irevcm == 1)
    [b, ifail1] = f11xa('No Transpose', a, irow, icol, 'N', x);
  elseif (irevcm == 2)
    [b, dgnl, ifail1] = f11dk(store, trans, init, niter, a, ...
                             irow, icol, 'Check', x, dgnl);
    init = 'n';
  elseif (irevcm == 3)
    [itn, stplhs, stprhs, anorm, sigmax, ifail1] = f11bf(work);
    fprintf('\nMonitoring at iteration number %d\nresidual norm: %14.4e\n', itn, stplhs);
  end

  if (ifail1 ~= 0)
    irevcm = 6;
  end
end

% Get information about the computation
if ifail == 0
  [itn, stplhs, stprhs, anorm, sigmax, ifail] = f11bf(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  Residual Vector\n');
  for i = 1:n
    fprintf('%+16.4e %+16.4e\n', x(i), b(i));
  end
end
 

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

   Solution Vector  Residual Vector
     +1.7035e+00      +3.2377e-07
     +1.0805e+00      -1.7625e-05
     +1.8305e+00      +2.7964e-05
     +6.0251e+00      -2.5914e-05
     +3.2942e+00      +7.8156e-06
     +1.9068e+00      +9.2064e-06
     +4.1365e+00      -3.0848e-06
     +5.2111e+00      +1.9834e-05


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

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