NAG Library Routine Document

F07PBF (DSPSVX)

INTEGER	N, NRHS, IPIV(N), LDB, LDX, IWORK(N), INFO
REAL (KIND=nag_wp)	AP(), AFP(), B(LDB,), X(LDX,), RCOND, FERR(NRHS), BERR(NRHS), WORK(3*N)
CHARACTER(1)	FACT, UPLO

The routine may be called by its LAPACK name dspsvx.

3 Description

F07PBF (DSPSVX) performs the following steps:

If $FACT ='N'$ , the diagonal pivoting method is used to factor $A$ as $A = U D U^{T}$ if $UPLO ='U'$ or $A = L D L^{T}$ if $UPLO ='L'$ , where $U$ (or $L$ ) is a product of permutation and unit upper (lower) triangular matrices and $D$ is symmetric and block diagonal with $1$ by $1$ and $2$ by $2$ diagonal blocks.
If some $d_{i i} = 0$ , so that $D$ is exactly singular, 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.

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

has been supplied.

$FACT ='F'$: AFP and IPIV contain the factorized form of the matrix $A$ . AFP and IPIV will not be modified.
$FACT ='N'$: The matrix $A$ will be copied to AFP and factorized.

Constraint:

FACT ='F'

'N'

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: 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

5: AP( $*$ ) – REAL (KIND=nag_wp) arrayInput

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

\max (1, N \times (N + 1) / 2)

On entry: the

n

n

symmetric matrix

A

, packed by columns.

More precisely,

if $UPLO ='U'$ , the upper triangle of $A$ must be stored with element $A_{i j}$ in $AP (i + j (j - 1) / 2)$ for $i \leq j$ ;
if $UPLO ='L'$ , the lower triangle of $A$ must be stored with element $A_{i j}$ in $AP (i + (2 n - j) (j - 1) / 2)$ for $i \geq j$ .

6: AFP( $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

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

\max (1, N \times (N + 1) / 2)

On entry: if

FACT ='F'

, AFP contains the block diagonal matrix

D

and the multipliers used to obtain the factor

U

L

from the factorization

A = U D U^{T}

A = L D L^{T}

as computed by F07PDF (DSPTRF), stored as a packed triangular matrix in the same storage format as

A

On exit: if

FACT ='N'

, AFP contains the block diagonal matrix

D

and the multipliers used to obtain the factor

U

L

from the factorization

A = U D U^{T}

A = L D L^{T}

as computed by F07PDF (DSPTRF), stored as a packed triangular matrix in the same storage format as

A

7: IPIV(N) – INTEGER arrayInput/Output

On entry: if

FACT ='F'

, IPIV contains details of the interchanges and the block structure of

D

, as determined by F07PDF (DSPTRF).

if $IPIV (i) = k > 0$ , $d_{i i}$ is a $1$ by $1$ pivot block and the $i$ th row and column of $A$ were interchanged with the $k$ th row and column;
if $UPLO ='U'$ and $IPIV (i - 1) = IPIV (i) = - l < 0$ , $(\begin{matrix} d_{i - 1, i - 1} & {\bar{d}}_{i, i - 1} \\ {\bar{d}}_{i, i - 1} & d_{i i} \end{matrix})$ is a $2$ by $2$ pivot block and the $(i - 1)$ th row and column of $A$ were interchanged with the $l$ th row and column;
if $UPLO ='L'$ and $IPIV (i) = IPIV (i + 1) = - m < 0$ , $(\begin{matrix} d_{i i} & d_{i + 1, i} \\ d_{i + 1, i} & d_{i + 1, i + 1} \end{matrix})$ is a $2$ by $2$ pivot block and the $(i + 1)$ th row and column of $A$ were interchanged with the $m$ th row and column.

On exit: if

FACT ='N'

, IPIV contains details of the interchanges and the block structure of

D

, as determined by F07PDF (DSPTRF), as described above.

8: B(LDB, $*$ ) – REAL (KIND=nag_wp) arrayInput

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

9: LDB – INTEGERInput

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

Constraint:

LDB \geq \max (1, N)

10: 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

11: LDX – INTEGERInput

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

Constraint:

LDX \geq \max (1, N)

12: RCOND – REAL (KIND=nag_wp)Output

On exit: the estimate of the reciprocal condition number of the matrix

A

. If

RCOND = 0.0

, the matrix may be exactly singular. This condition is indicated by

INFO > 0 and INFO \leq N

. Otherwise, if RCOND is less than the machine precision, the matrix is singular to working precision. This condition is indicated by

INFO = N + 1

13: 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.

14: 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).

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

16: IWORK(N) – INTEGER arrayWorkspace

17: 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 \leq N$ , $d (i, i)$ is exactly zero. The factorization has been completed, but the factor $D$ is exactly singular, so the solution and error bounds could not be computed. $RCOND = 0.0$ is returned.

$INFO = N + 1$: $D$ 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

\hat{x}

is the exact solution of a perturbed system of equations

(A + E) \hat{x} = b

, where

{‖E‖}_{1} = O (ε) {‖A‖}_{1},

where

ε

is the machine precision. See Chapter 11 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

The factorization of

A

requires approximately

\frac{1}{3} n^{3}

floating point operations.

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

4 n^{2}

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

6 n^{2}

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

2 n^{2}

operations.

The complex analogues of this routine are F07PPF (ZHPSVX) for Hermitian matrices, and F07QPF (ZSPSVX) for symmetric matrices.

9 Example

This example solves the equations

A X = B,

where

A

is the symmetric matrix

A = (\begin{matrix} - 1.81 & 2.06 & 0.63 & - 1.15 \\ 2.06 & 1.15 & 1.87 & 4.20 \\ 0.63 & 1.87 & - 0.21 & 3.87 \\ - 1.15 & 4.20 & 3.87 & 2.07 \end{matrix}) and B = (\begin{matrix} 0.96 & 3.93 \\ 6.07 & 19.25 \\ 8.38 & 9.90 \\ 9.50 & 27.85 \end{matrix}) .

Error estimates for the solutions, and an estimate of the reciprocal of the condition number of the matrix

A

are also output.

NAG Library Routine DocumentF07PBF (DSPSVX)

+− 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

F07PBF (DSPSVX)