NAG Library Routine Document

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

C Header Interface

#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)

3

Description

To solve a system of real linear equations

A x = b

A^{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.

4

References

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

5

Arguments

1: $n$ – IntegerInput

On entry:

n

, the order of the matrix

A

Constraint:

n \geq 0

2: $a (licn)$ – Real (Kind=nag_wp) arrayInput

On entry: the nonzero elements in the factorization of the matrix

A

, as returned by f01brf or f01bsf.

3: $licn$ – IntegerInput

On entry: the dimension of the arrays a and icn as declared in the (sub)program from which f04axf is called.

4: $icn (licn)$ – Integer arrayCommunication Array

On entry: the column indices of the nonzero elements of the factorization, as returned by f01brf or f01bsf.

5: $ikeep (5 \times n)$ – 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: $rhs (n)$ – 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: $w (n)$ – Real (Kind=nag_wp) arrayWorkspace

8: $mtype$ – IntegerInput

On entry: mtype specifies the task to be performed.

$mtype = 1$: Solve $A x = b$ .
$mtype \neq 1$: Solve $A^{T} x = b$ .

9: $idisp (2)$ – Integer arrayCommunication Array

On entry: the values returned in idisp by f01brf.

10: $resid$ – Real (Kind=nag_wp)Output

On exit: the value of the maximum residual,

\max (|b_{i} - \sum_{j}^{} a_{i j} x_{j}|)

, over all the unsatisfied equations, in case f01brf or f01bsf has been used to factorize a singular or rectangular matrix.

6

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

7

Accuracy

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

grow = .TRUE.

on entry to these routines and to examine the value of

w (1)

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 - A x (or ​ b - A^{T} x)

and calling f04axf again to find a correction

δ

for

x

by solving

A δ = r (or ​ A^{T} δ = r) .

8

Parallelism and Performance

f04axf is not threaded in any implementation.

9

Further Comments

If the factorized form contains

τ

nonzeros (

idisp (2) = τ

) then the time taken is very approximately

2 τ

times longer than the inner loop of full matrix code. Some advantage is taken of zeros in the right-hand side when solving

A^{T} x = b

(

mtype \neq 1

10

Example

This example solves the set of linear equations

A x = b

where

A = (\begin{array}{r} 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 \end{array}) and b = (\begin{array}{r} 15 \\ 12 \\ 18 \\ 3 \\ - 6 \\ 0 \end{array}) .

The nonzero elements of

A

and indexing information are read in by the program, as described in the document for f01brf.

NAG Library Routine Document

f04axf (real_sparse_fac_solve)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification