f11be:: Large Scale Linear Systems (NAG Toolbox)

nag_sparse_real_gen_basic_solver (f11be) solves the general (nonsymmetric) system of linear simultaneous equations

A x = b

of order

n

, where

n

is large and the coefficient matrix

A

is sparse, using one of four available methods: RGMRES, the preconditioned restarted generalized minimum residual method (see Saad and Schultz (1986)); CGS, the preconditioned conjugate gradient squared method (see Sonneveld (1989)); Bi-CGSTAB(

ℓ

), the bi-conjugate gradient stabilized method of order

ℓ

(see Van der Vorst (1989) and Sleijpen and Fokkema (1993)); or TFQMR, the transpose-free quasi-minimal residual method (see Freund and Nachtigal (1991) and Freund (1993)).

nag_sparse_real_gen_basic_solver (f11be) uses reverse communication, i.e., it returns repeatedly to the calling program with the argument irevcm (see Arguments) set to specified values which require the calling program to carry out one of the following tasks:

–	compute the matrix-vector product $v = A u$ or $v = A^{T} u$ (the four methods require the matrix transpose-vector product only if ${‖A‖}_{1}$ or ${‖A‖}_{\infty}$ is estimated internally by Higham's method (see Higham (1988)));
–	solve the preconditioning equation $M v = u$ ;
–	notify the completion of the computation;
–	allow the calling program to monitor the solution.

Reverse communication has the following advantages.

1.	Maximum flexibility in the representation and storage of sparse matrices: all matrix operations are performed outside the solver function, thereby avoiding the need for a complicated interface with enough flexibility to cope with all types of storage schemes and sparsity patterns. This applies also to preconditioners.
2.	Enhanced user interaction: you can closely monitor the progress of the solution and tidy or immediate termination can be requested. This is useful, for example, when alternative termination criteria are to be employed or in case of failure of the external functions used to perform matrix operations.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

Accuracy

Further Comments

If the termination criterion

{‖r_{k}‖}_{p} \leq τ ({‖b‖}_{p} + {‖A‖}_{p} \times {‖x_{k}‖}_{p})

is used (see Description in nag_sparse_real_gen_basic_setup (f11bd)) and

‖x_{0}‖ ≫ ‖x_{k}‖

, then the required accuracy cannot be obtained due to loss of significant digits. The iteration is restarted automatically at some suitable point: nag_sparse_real_gen_basic_solver (f11be) sets

x_{0} = x_{k}

and the computation begins again. For particularly badly scaled problems, more than one restart may be necessary. This does not apply to the RGMRES method which, by its own nature, self-restarts every super-iteration. Naturally, restarting adds to computational costs: it is recommended that the iteration should start from a value

x_{0}

which is as close to the true solution

\tilde{x}

as can be estimated. Otherwise, the iteration should start from

x_{0} = 0

Example

function f11be_example


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

% Solve sparse system Ax = b iteratively using ILU preconditioner

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

% ILU preconditioner
%     Input parameters
lfill = int64(0);
dtol  = 0;
milu  = 'No modification';
ipivp = zeros(n, 1, 'int64');
ipivq = zeros(n, 1, 'int64');
% Compute preconditioner
[a, irow, icol, ipivp, ipivq, istr, idiag, nnzc, npivm, ifail] = ...
    f11da(...
          nz, a, irow, icol, lfill, dtol, milu, ipivp, ipivq);

% Iterative Solution Setup
%      Input parameters
method = 'BICGSTAB';
precon = 'Preconditioned';
lpoly  = int64(1);
tol    = sqrt(x02aj);
maxitn = int64(20);
anorm  = 0;
sigmax = 0;
monit  = int64(1);
lwork  = int64(6000);
% Setup iterative method
[lwreq, work, ifail] = ...
  f11bd(...
        method, precon, n, lpoly, tol, maxitn, anorm, sigmax, ...
        monit, lwork, 'norm_p', '1');

% Reverse communication loop calling iterative solver
irevcm = int64(0);
u      = zeros(8, 1);
v      = b;
wgt    = zeros(8, 1);

while (irevcm ~= 4)

  [irevcm, u, v, work, ifail] = f11be(...
                                      irevcm, u, v, wgt, work);

  if (irevcm == -1)
    % v = A^Tu
    [v, ifail] = f11xa(...
                       'T', a(1:nz), irow(1:nz), icol(1:nz), 'N', u);
  elseif (irevcm == 1)
    % v = Au
    [v, ifail] = f11xa(...
                       'N', a(1:nz), irow(1:nz), icol(1:nz), 'N', u);
  elseif (irevcm == 2)
    % Solve ILU v = u  
    [v, ifail] = f11db('N', a, irow, icol, ipivp, ipivq, istr, idiag, 'N', u);
    if (ifail ~= 0)
      irevcm = 6;
    end
  elseif (irevcm == 3)
    % Monitoring stage
    [itn, stplhs, stprhs, anorm, sigmax, ifail] = ...
      f11bf(work);
    fprintf('\nMonitoring at iteration number %2d\n',itn);
    fprintf('residual norm:              %14.4e\n', stplhs);
    fprintf('\n   Solution Vector  Residual Vector\n');
    for i = 1:n
      fprintf('%16.4f %16.2e\n', u(i), v(i));
    end
  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', u(i), v(i));
end

f11be example results


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

   Solution Vector  Residual Vector
         -4.5858         1.53e+01
          1.0154         2.66e+01
         -2.2234        -8.75e+00
          6.0097         1.86e+01
          1.3827         8.28e+00
         -7.9070         2.04e+01
          0.4427         9.61e+00
          5.9248         3.31e+01

Monitoring at iteration number  2
residual norm:                  3.2742e+01

   Solution Vector  Residual Vector
          4.1642        -2.96e+00
          4.9370        -5.55e+00
          4.8101         8.21e-01
          5.4324        -1.68e+01
          5.8531         5.60e-01
         11.9250        -1.91e+00
          8.4826         1.01e+00
          6.0625        -3.10e+00

Number of iterations for convergence:        3
Residual norm:                               1.0373e-08
Right-hand side of termination criteria:     5.8900e-06
i-norm of matrix a:                          1.1000e+01

   Solution Vector  Residual Vector
          1.0000        -1.36e-09
          2.0000        -2.61e-09
          3.0000         2.25e-10
          4.0000        -3.22e-09
          5.0000         6.30e-10
          6.0000        -5.24e-10
          7.0000         9.58e-10
          8.0000        -8.49e-10

Method	Requirements
RGMRES	$lwork = 100 + n (m + 3) + m (m + 5) + 1$ , where $m$ is the dimension of the basis.
CGS	$lwork = 100 + 7 n$ .
Bi-CGSTAB( $ℓ$ )	$lwork = 100 + (2 n + ℓ) (ℓ + 2) + p$ , where $ℓ$ is the order of the method.
TFQMR	$lwork = 100 + 10 n$ ,

$p = 2 n$	if $ℓ > 1$ and $iterm = 2$ was supplied.
$p = n$	if $ℓ > 1$ and a preconditioner is used or $iterm = 2$ was supplied.
$p = 0$	otherwise.

NAG Toolbox: nag_sparse_real_gen_basic_solver (f11be)

▸▿ Contents

Purpose

Syntax

Description