NAG Library Routine Document

f04axf (real_sparse_fac_solve)

1
Purpose

f04axf calculates the approximate solution of a set of real sparse linear equations with a single right-hand side, Ax=b or ATx=b, where A has been factorized by f01brf or f01bsf.

2
Specification

Fortran Interface
Subroutine f04axf ( n, a, licn, icn, ikeep, rhs, w, mtype, idisp, resid)
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 Ax=b or ATx=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: n0.
2:     alicn – 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:     icnlicn – Integer arrayCommunication Array
On entry: the column indices of the nonzero elements of the factorization, as returned by f01brf or f01bsf.
5:     ikeep5×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:     rhsn – 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:     wn – Real (Kind=nag_wp) arrayWorkspace
8:     mtype – IntegerInput
On entry: mtype specifies the task to be performed.
mtype=1
Solve Ax=b.
mtype1
Solve ATx=b.
9:     idisp2 – 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, maxbi-jaijxj, 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 w1 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-Axor ​b-ATx  
and calling f04axf again to find a correction δ for x by solving
Aδ=ror ​ATδ=r.  

8
Parallelism and Performance

f04axf is not threaded in any implementation.

9
Further Comments

If the factorized form contains τ nonzeros (idisp2=τ) 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 ATx=b (mtype1).

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

10.1
Program Text

Program Text (f04axfe.f90)

10.2
Program Data

Program Data (f04axfe.d)

10.3
Program Results

Program Results (f04axfe.r)