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$Ax=b$ or ATx = b${A}^{\mathrm{T}}x=b$, where A$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$Ax=b$ or ATx = b${A}^{\mathrm{T}}x=b$, where A$A$ is a general sparse matrix, A$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$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:     a(licn) – double array
The nonzero elements in the factorization of the matrix A$A$, as returned by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs).
2:     icn(licn) – 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:     ikeep(5 × n$5×{\mathbf{n}}$) – 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:     rhs(n) – double array
n, the dimension of the array, must satisfy the constraint n0${\mathbf{n}}\ge 0$.
The right-hand side vector b$b$.
5:     mtype – int64int32nag_int scalar
mtype specifies the task to be performed.
mtype = 1${\mathbf{mtype}}=1$
Solve Ax = b$Ax=b$.
mtype1${\mathbf{mtype}}\ne 1$
Solve ATx = b${A}^{\mathrm{T}}x=b$.
6:     idisp(2$2$) – int64int32nag_int array
The values returned in idisp by nag_matop_real_gen_sparse_lu (f01br).

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the array rhs.
n$n$, the order of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.
2:     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 as declared in the (sub)program from which nag_linsys_real_sparse_fac_solve (f04ax) is called.

w

### Output Parameters

1:     rhs(n) – double array
rhs stores the solution vector x$x$.
2:     resid – double scalar
The value of the maximum residual, max (|bijaijxj|)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(|{b}_{i}-\sum _{j}^{}{a}_{ij}{x}_{j}|\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.

## Error Indicators and Warnings

If an error is detected in an input parameter nag_linsys_real_sparse_fac_solve (f04ax) will act as if a soft noisy exit has been requested (see Section [Soft Fail Option] in the (essin)).

## 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 grow = true${\mathbf{grow}}=\mathbf{true}$ on entry to these functions and to examine the value of w(1)$\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) $r=b-Ax (or ​b-ATx)$
and calling nag_linsys_real_sparse_fac_solve (f04ax) again to find a correction δ$\delta$ for x$x$ by solving
 Aδ = r (or ​ATδ = r). $Aδ=r (or ​ATδ=r).$

If the factorized form contains τ$\tau$ nonzeros (idisp(2) = τ${\mathbf{idisp}}\left(2\right)=\tau$) then the time taken is very approximately 2τ$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 ATx = b${A}^{\mathrm{T}}x=b$ (mtype1${\mathbf{mtype}}\ne 1$).

## Example

```function nag_linsys_real_sparse_fac_solve_example
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');
irn(1:15) = [int64(1);2;2;2;3;4; ...
4;4;5;5;5;5; ...
6;6;6];
icn = zeros(150,1,'int64');
icn(1:15) = [int64(1);2;3;4;3;1; ...
4;5;1;4;5;6; ...
1;2;6];
abort = [true;
true;
false;
true];

rhs = [15;
12;
18;
3;
-6;
0];
mtype = int64(1);
[a, irn, icn, ikeep, w, idisp, ifail] = nag_matop_real_gen_sparse_lu(n, nz, a, irn, icn, abort);
[rhsOut, resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, rhs, mtype, idisp)
```
```

rhsOut =

3
3
6
6
3
1

resid =

0

```
```function f04ax_example
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');
irn(1:15) = [int64(1);2;2;2;3;4; ...
4;4;5;5;5;5; ...
6;6;6];
icn = zeros(150,1,'int64');
icn(1:15) = [int64(1);2;3;4;3;1; ...
4;5;1;4;5;6; ...
1;2;6];
abort = [true;
true;
false;
true];

rhs = [15;
12;
18;
3;
-6;
0];
mtype = int64(1);
[a, irn, icn, ikeep, w, idisp, ifail] = f01br(n, nz, a, irn, icn, abort);
[rhsOut, resid] = f04ax(a, icn, ikeep, rhs, mtype, idisp)
```
```

rhsOut =

3
3
6
6
3
1

resid =

0

```