# NAG Library Routine Document

## 1Purpose

f04axf 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 f01brf or f01bsf.

## 2Specification

Fortran Interface
 Subroutine f04axf ( n, a, licn, icn, rhs, w,
 Integer, Intent (In) :: n, licn, icn(licn), mtype, idisp(2) Integer, Intent (Inout) :: ikeep(5*n) Real (Kind=nag_wp), Intent (In) :: a(licn) Real (Kind=nag_wp), Intent (Inout) :: rhs(n) Real (Kind=nag_wp), Intent (Out) :: w(n), resid
#include nagmk26.h
 void f04axf_ (const Integer *n, const double a[], const Integer *licn, const Integer icn[], Integer ikeep[], double rhs[], double w[], const Integer *mtype, const Integer idisp[], double *resid)

## 3Description

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 f01brf or f01bsf. f04axf 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 routine f11mff and factorization routine f11mef.

## 4References

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

## 5Arguments

1:     $\mathbf{n}$ – IntegerInput
On entry: $n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
2:     $\mathbf{a}\left({\mathbf{licn}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: the nonzero elements in the factorization of the matrix $A$, as returned by f01brf or f01bsf.
3:     $\mathbf{licn}$ – IntegerInput
On entry: the dimension of the arrays a and icn as declared in the (sub)program from which f04axf is called.
4:     $\mathbf{icn}\left({\mathbf{licn}}\right)$ – Integer arrayCommunication Array
On entry: the column indices of the nonzero elements of the factorization, as returned by f01brf or f01bsf.
5:     $\mathbf{ikeep}\left(5×{\mathbf{n}}\right)$ – Integer arrayInput
ikeep provides, on entry, indexing information about the factorization, as returned by f01brf or f01bsf. Used as internal workspace prior to being restored and hence is unchanged.
6:     $\mathbf{rhs}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: the right-hand side vector $b$.
On exit: rhs is overwritten by the solution vector $x$.
7:     $\mathbf{w}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayWorkspace
8:     $\mathbf{mtype}$ – IntegerInput
On entry: mtype specifies the task to be performed.
${\mathbf{mtype}}=1$
Solve $Ax=b$.
${\mathbf{mtype}}\ne 1$
Solve ${A}^{\mathrm{T}}x=b$.
9:     $\mathbf{idisp}\left(2\right)$ – Integer arrayCommunication Array
On entry: the values returned in idisp by f01brf.
10:   $\mathbf{resid}$ – Real (Kind=nag_wp)Output
On exit: 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 f01brf or f01bsf has been used to factorize a singular or rectangular matrix.

## 6Error Indicators and Warnings

If an error is detected in an input argument f04axf will act as if a soft noisy exit has been requested (see Section 3.4.4 in How to Use the NAG Library and its Documentation).

## 7Accuracy

The accuracy of the computed solution depends on the conditioning of the original matrix. Since f04axf is always used with either f01brf or f01bsf, you are recommended to set ${\mathbf{grow}}=\mathrm{.TRUE.}$ on entry to these routines and to examine the value of ${\mathbf{w}}\left(1\right)$ on exit (see f01brf and f01bsf). 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 f04axf again to find a correction $\delta$ for $x$ by solving
 $Aδ=r or ​ATδ=r.$

## 8Parallelism and Performance

f04axf is not threaded in any implementation.

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$).

## 10Example

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 f01brf.

### 10.1Program Text

Program Text (f04axfe.f90)

### 10.2Program Data

Program Data (f04axfe.d)

### 10.3Program Results

Program Results (f04axfe.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017