NAG Library Routine Document

F07HBF (DPBSVX)

INTEGER	N, KD, NRHS, LDAB, LDAFB, LDB, LDX, IWORK(N), INFO
REAL (KIND=nag_wp)	AB(LDAB,), AFB(LDAFB,), S(), B(LDB,), X(LDX,), RCOND, FERR(NRHS), BERR(NRHS), WORK(3N)
CHARACTER(1)	FACT, UPLO, EQUED

The routine may be called by its LAPACK name dpbsvx.

3 Description

F07HBF (DPBSVX) performs the following steps:

If $FACT ='E'$ , real diagonal scaling factors, $D_{S}$ , are computed to equilibrate the system:
$(D_{S} A D_{S}) (D_{S}^{- 1} X) = D_{S} B .$
Whether or not the system will be equilibrated depends on the scaling of the matrix $A$ , but if equilibration is used, $A$ is overwritten by $D_{S} A D_{S}$ and $B$ by $D_{S} B$ .
If $FACT ='N'$ or $'E'$ , the Cholesky decomposition is used to factor the matrix $A$ (after equilibration if $FACT ='E'$ ) as $A = U^{T} U$ if $UPLO ='U'$ or $A = L L^{T}$ if $UPLO ='L'$ , where $U$ is an upper triangular matrix and $L$ is a lower triangular matrix.
If the leading $i$ by $i$ principal minor of $A$ is not positive definite, then the routine returns with $INFO = i$ . Otherwise, the factored form of $A$ is used to estimate the condition number of the matrix $A$ . If the reciprocal of the condition number is less than machine precision, $INFO = N + 1$ is returned as a warning, but the routine still goes on to solve for $X$ and compute error bounds as described below.
The system of equations is solved for $X$ using the factored form of $A$ .
Iterative refinement is applied to improve the computed solution matrix and to calculate error bounds and backward error estimates for it.
If equilibration was used, the matrix $X$ is premultiplied by $D_{S}$ so that it solves the original system before equilibration.

4 References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Higham N J (2002) Accuracy and Stability of Numerical Algorithms (2nd Edition) SIAM, Philadelphia

5 Parameters

1: FACT – CHARACTER(1)Input

On entry: specifies whether or not the factorized form of the matrix

A

is supplied on entry, and if not, whether the matrix

A

should be equilibrated before it is factorized.

$FACT ='F'$: AFB contains the factorized form of $A$ . If $EQUED ='Y'$ , the matrix $A$ has been equilibrated with scaling factors given by S. AB and AFB will not be modified.
$FACT ='N'$: The matrix $A$ will be copied to AFB and factorized.
$FACT ='E'$: The matrix $A$ will be equilibrated if necessary, then copied to AFB and factorized.

Constraint:

FACT ='F'

'N'

'E'

2: UPLO – CHARACTER(1)Input

On entry: if

UPLO ='U'

, the upper triangle of

A

is stored.

UPLO ='L'

, the lower triangle of

A

is stored.

Constraint:

UPLO ='U'

'L'

3: N – INTEGERInput

On entry:

n

, the number of linear equations, i.e., the order of the matrix

A

Constraint:

N \geq 0

4: KD – INTEGERInput

On entry:

k_{d}

, the number of superdiagonals of the matrix

A

UPLO ='U'

, or the number of subdiagonals if

UPLO ='L'

Constraint:

KD \geq 0

5: NRHS – INTEGERInput

On entry:

r

, the number of right-hand sides, i.e., the number of columns of the matrix

B

Constraint:

NRHS \geq 0

6: AB(LDAB, $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array AB must be at least

\max (1, N)

On entry: the upper or lower triangle of the symmetric band matrix

A

, except if

FACT ='F'

and

EQUED ='Y'

, in which case AB must contain the equilibrated matrix

D_{S} A D_{S}

The matrix is stored in rows

1

k_{d} + 1

, more precisely,

if $UPLO ='U'$ , the elements of the upper triangle of $A$ within the band must be stored with element $A_{i j}$ in $AB (k_{d} + 1 + i - j, j) for \max (1, j - k_{d}) \leq i \leq j$ ;
if $UPLO ='L'$ , the elements of the lower triangle of $A$ within the band must be stored with element $A_{i j}$ in $AB (1 + i - j, j) for j \leq i \leq \min (n, j + k_{d}) .$

On exit: if

FACT ='E'

and

EQUED ='Y'

, AB is overwritten by

D_{S} A D_{S}

7: LDAB – INTEGERInput

On entry: the first dimension of the array AB as declared in the (sub)program from which F07HBF (DPBSVX) is called.

Constraint:

LDAB \geq KD + 1

8: AFB(LDAFB, $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array AFB must be at least

\max (1, N)

On entry: if

FACT ='F'

, AFB contains the triangular factor

U

L

from the Cholesky factorization

A = U^{T} U

A = L L^{T}

of the band matrix

A

, in the same storage format as

A

. If

EQUED ='Y'

, AFB is the factorized form of the equilibrated matrix

A

On exit: if

FACT ='N'

, AFB returns the triangular factor

U

L

from the Cholesky factorization

A = U^{T} U

A = L L^{T}

FACT ='E'

, AFB returns the triangular factor

U

L

from the Cholesky factorization

A = U^{T} U

A = L L^{T}

of the equilibrated matrix

A

(see the description of AB for the form of the equilibrated matrix).

9: LDAFB – INTEGERInput

On entry: the first dimension of the array AFB as declared in the (sub)program from which F07HBF (DPBSVX) is called.

Constraint:

LDAFB \geq KD + 1

10: EQUED – CHARACTER(1)Input/Output

On entry: if

FACT ='N'

'E'

, EQUED need not be set.

FACT ='F'

, EQUED must specify the form of the equilibration that was performed as follows:

if $EQUED ='N'$ , no equilibration;
if $EQUED ='Y'$ , equilibration was performed, i.e., $A$ has been replaced by $D_{S} A D_{S}$ .

On exit: if

FACT ='F'

, EQUED is unchanged from entry.

Otherwise, if no constraints are violated, EQUED specifies the form of the equilibration that was performed as specified above.

Constraint: if

FACT ='F'

EQUED ='N'

'Y'

11: S( $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

Note: the dimension of the array S must be at least

\max (1, N)

On entry: if

FACT ='N'

'E'

, S need not be set.

FACT ='F'

and

EQUED ='Y'

, S must contain the scale factors,

D_{S}

, for

A

; each element of S must be positive.

On exit: if

FACT ='F'

, S is unchanged from entry.

Otherwise, if no constraints are violated and

EQUED ='Y'

, S contains the scale factors,

D_{S}

, for

A

; each element of S is positive.

12: B(LDB, $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array B must be at least

\max (1, NRHS)

On entry: the

n

r

right-hand side matrix

B

On exit: if

EQUED ='N'

, B is not modified.

EQUED ='Y'

, B is overwritten by

D_{S} B

13: LDB – INTEGERInput

On entry: the first dimension of the array B as declared in the (sub)program from which F07HBF (DPBSVX) is called.

Constraint:

LDB \geq \max (1, N)

14: X(LDX, $*$ ) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array X must be at least

\max (1, NRHS)

On exit: if

INFO = 0

N + 1

, the

n

r

solution matrix

X

to the original system of equations. Note that the arrays

A

and

B

are modified on exit if

EQUED ='Y'

, and the solution to the equilibrated system is

D_{S}^{- 1} X

15: LDX – INTEGERInput

On entry: the first dimension of the array X as declared in the (sub)program from which F07HBF (DPBSVX) is called.

Constraint:

LDX \geq \max (1, N)

16: RCOND – REAL (KIND=nag_wp)Output

On exit: if no constraints are violated, an estimate of the reciprocal condition number of the matrix

A

(after equilibration if that is performed), computed as

RCOND = 1.0 / ({‖A‖}_{1} {‖A^{- 1}‖}_{1})

17: FERR(NRHS) – REAL (KIND=nag_wp) arrayOutput

On exit: if

INFO = 0

N + 1

, an estimate of the forward error bound for each computed solution vector, such that

{‖{\hat{x}}_{j} - x_{j}‖}_{\infty} / {‖x_{j}‖}_{\infty} \leq FERR (j)

where

{\hat{x}}_{j}

is the

j

th column of the computed solution returned in the array X and

x_{j}

is the corresponding column of the exact solution

X

. The estimate is as reliable as the estimate for RCOND, and is almost always a slight overestimate of the true error.

18: BERR(NRHS) – REAL (KIND=nag_wp) arrayOutput

On exit: if

INFO = 0

N + 1

, an estimate of the component-wise relative backward error of each computed solution vector

{\hat{x}}_{j}

(i.e., the smallest relative change in any element of

A

B

that makes

{\hat{x}}_{j}

an exact solution).

19: WORK( $3 \times N$ ) – REAL (KIND=nag_wp) arrayWorkspace

20: IWORK(N) – INTEGER arrayWorkspace

21: INFO – INTEGEROutput

On exit:

INFO = 0

unless the routine detects an error (see Section 6).

6 Error Indicators and Warnings

Errors or warnings detected by the routine:

$INFO < 0$: If $INFO = - i$ , the $i$ th argument had an illegal value. An explanatory message is output, and execution of the program is terminated.

$INFO > 0 and INFO \leq N$: If $INFO = i$ and $i \leq N$ , the leading minor of order $i$ of $A$ is not positive definite, so the factorization could not be completed, and the solution has not been computed. $RCOND = 0.0$ is returned.

$INFO = N + 1$: The triangular matrix $U$ (or $L$ ) is nonsingular, but RCOND is less than machine precision, meaning that the matrix is singular to working precision. Nevertheless, the solution and error bounds are computed because there are a number of situations where the computed solution can be more accurate than the value of RCOND would suggest.

7 Accuracy

For each right-hand side vector

b

, the computed solution

x

is the exact solution of a perturbed system of equations

(A + E) x = b

, where

if $UPLO ='U'$ , $|E| \leq c (n) ε |U^{T}| |U|$ ;
if $UPLO ='L'$ , $|E| \leq c (n) ε |L| |L^{T}|$ ,

c (n)

is a modest linear function of

n

, and

ε

is the machine precision. See Section 10.1 of Higham (2002) for further details.

\hat{x}

is the true solution, then the computed solution

x

satisfies a forward error bound of the form

\frac{{‖x - \hat{x}‖}_{\infty}}{{‖\hat{x}‖}_{\infty}} \leq w_{c} cond (A, \hat{x}, b)

where

cond (A, \hat{x}, b) = {‖|A^{- 1}| (|A| |\hat{x}| + |b|)‖}_{\infty} / {‖\hat{x}‖}_{\infty} \leq cond (A) = {‖|A^{- 1}| |A|‖}_{\infty} \leq κ_{\infty} (A)

. If

\hat{x}

is the

j

th column of

X

, then

w_{c}

is returned in

BERR (j)

and a bound on

{‖x - \hat{x}‖}_{\infty} / {‖\hat{x}‖}_{\infty}

is returned in

FERR (j)

. See Section 4.4 of Anderson et al. (1999) for further details.

8 Further Comments

When

n ≫ k

, the factorization of

A

requires approximately

n {(k + 1)}^{2}

floating point operations, where

k

is the number of superdiagonals.

For each right-hand side, computation of the backward error involves a minimum of

8 n k

floating point operations. Each step of iterative refinement involves an additional

12 n k

operations. At most five steps of iterative refinement are performed, but usually only one or two steps are required. Estimating the forward error involves solving a number of systems of equations of the form

A x = b

; the number is usually

4

5

and never more than

11

. Each solution involves approximately

4 n k

operations.

The complex analogue of this routine is F07HPF (ZPBSVX).

9 Example

This example solves the equations

A X = B,

where

A

is the symmetric positive definite band matrix

A = (\begin{array}{r} 5.49 & 2.68 & 0 & 0 \\ 2.68 & 5.63 & - 2.39 & 0 \\ 0 & - 2.39 & 2.60 & - 2.22 \\ 0 & 0 & - 2.22 & 5.17 \end{array})

and

B = (\begin{matrix} 22.09 & 5.10 \\ 9.31 & 30.81 \\ - 5.24 & - 25.82 \\ 11.83 & 22.90 \end{matrix}) .

Error estimates for the solutions, information on equilibration and an estimate of the reciprocal of the condition number of the scaled matrix

A

are also output.

NAG Library Routine DocumentF07HBF (DPBSVX)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

F07HBF (DPBSVX)