Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_linsys_real_sparse_fac_solve (f04ax)

## Purpose

nag_linsys_real_sparse_fac_solve (f04ax) calculates the approximate solution of a set of real sparse linear equations with a single right-hand side, $Ax=b$ or ${A}^{\mathrm{T}}x=b$, where $A$ has been factorized by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs).

## Syntax

[rhs, resid] = f04ax(a, icn, ikeep, rhs, mtype, idisp, 'n', n, 'licn', licn)
[rhs, resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, rhs, mtype, idisp, 'n', n, 'licn', licn)

## Description

To solve a system of real linear equations $Ax=b$ or ${A}^{\mathrm{T}}x=b$, where $A$ is a general sparse matrix, $A$ must first be factorized by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs). nag_linsys_real_sparse_fac_solve (f04ax) then computes $x$ by block forward or backward substitution using simple forward and backward substitution within each diagonal block.
The method is fully described in Duff (1977).
A more recent method is available through solver function nag_sparse_direct_real_gen_solve (f11mf) and factorization function nag_sparse_direct_real_gen_lu (f11me).

## References

Duff I S (1977) MA28 – a set of Fortran subroutines for sparse unsymmetric linear equations AERE Report R8730 HMSO

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{a}\left({\mathbf{licn}}\right)$ – double array
The nonzero elements in the factorization of the matrix $A$, as returned by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs).
2:     $\mathrm{icn}\left({\mathbf{licn}}\right)$int64int32nag_int array
The column indices of the nonzero elements of the factorization, as returned by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs).
3:     $\mathrm{ikeep}\left(5×{\mathbf{n}}\right)$int64int32nag_int array
ikeep provides, on entry, indexing information about the factorization, as returned by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs). Used as internal workspace prior to being restored and hence is unchanged.
4:     $\mathrm{rhs}\left({\mathbf{n}}\right)$ – double array
The right-hand side vector $b$.
5:     $\mathrm{mtype}$int64int32nag_int scalar
mtype specifies the task to be performed.
${\mathbf{mtype}}=1$
Solve $Ax=b$.
${\mathbf{mtype}}\ne 1$
Solve ${A}^{\mathrm{T}}x=b$.
6:     $\mathrm{idisp}\left(2\right)$int64int32nag_int array
The values returned in idisp by nag_matop_real_gen_sparse_lu (f01br).

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the dimension of the array rhs.
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
2:     $\mathrm{licn}$int64int32nag_int scalar
Default: the dimension of the arrays a, icn. (An error is raised if these dimensions are not equal.)
The dimension of the arrays a and icn.

### Output Parameters

1:     $\mathrm{rhs}\left({\mathbf{n}}\right)$ – double array
rhs stores the solution vector $x$.
2:     $\mathrm{resid}$ – double scalar
The value of the maximum residual, $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(\left|{b}_{i}-\sum _{j}^{}{a}_{ij}{x}_{j}\right|\right)$, over all the unsatisfied equations, in case nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs) has been used to factorize a singular or rectangular matrix.

None.

## Accuracy

The accuracy of the computed solution depends on the conditioning of the original matrix. Since nag_linsys_real_sparse_fac_solve (f04ax) is always used with either nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs), you are recommended to set ${\mathbf{grow}}=\mathit{true}$ on entry to these functions and to examine the value of $\mathit{w}\left(1\right)$ on exit (see nag_matop_real_gen_sparse_lu (f01br) and nag_matop_real_gen_sparse_lu_reuse (f01bs)). For a detailed error analysis see page 17 of Duff (1977).
If storage for the original matrix is available then the error can be estimated by calculating the residual
 $r=b-Ax or ​b-ATx$
and calling nag_linsys_real_sparse_fac_solve (f04ax) again to find a correction $\delta$ for $x$ by solving
 $Aδ=r or ​ATδ=r.$

If the factorized form contains $\tau$ nonzeros (${\mathbf{idisp}}\left(2\right)=\tau$) then the time taken is very approximately $2\tau$ times longer than the inner loop of full matrix code. Some advantage is taken of zeros in the right-hand side when solving ${A}^{\mathrm{T}}x=b$ (${\mathbf{mtype}}\ne 1$).

## Example

This example solves the set of linear equations $Ax=b$ where
 $A= 5 0 0 0 0 0 0 2 -1 2 0 0 0 0 3 0 0 0 -2 0 0 1 1 0 -1 0 0 -1 2 -3 -1 -1 0 0 0 6 and b= 15 12 18 3 -6 0 .$
The nonzero elements of $A$ and indexing information are read in by the program, as described in the document for nag_matop_real_gen_sparse_lu (f01br).
```function f04ax_example

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

% Solve sparse system Ax = b
n  = int64(6);
nz = int64(15);
a  = zeros(150,1);
a(1:15) = [ 5;
2; -1; 2;
3;
-2;         1; 1;
-1;        -1; 2; -3;
-1; -1;            6];

irn = zeros(75,1,'int64');
icn = zeros(150,1,'int64');
irn(1:15) = [int64(1);
2; 2; 2;
3;
4;       4; 4;
5;       5; 5; 5;
6; 6;          6];
icn(1:15) = [int64(1);
2; 3; 4;
3;
1;       4; 5;
1;       4; 5; 6;
1; 2;          6];

b   = [15;
12;
18;
3;
-6;
0];

% Factorize A
abort = [true;     true;     false;     true];
[a, irn, icn, ikeep, w, idisp, ifail] = ...
f01br( ...
n, nz, a, irn, icn, abort);

% Solve Ax = b
mtype = int64(1);
[x, resid] = f04ax( ...
a, icn, ikeep, b, mtype, idisp);

disp('Solution');
disp(x);

```
```f04ax example results

Solution
3
3
6
6
3
1

```