NAG Library Routine Document

F07APF (ZGESVX)

INTEGER	N, NRHS, LDA, LDAF, IPIV(*), LDB, LDX, INFO
REAL (KIND=nag_wp)	R(), C(), RCOND, FERR(NRHS), BERR(NRHS), RWORK(max(1,2*N))
COMPLEX (KIND=nag_wp)	A(LDA,), AF(LDAF,), B(LDB,), X(LDX,), WORK(2*N)
CHARACTER(1)	FACT, TRANS, EQUED

The routine may be called by its LAPACK name zgesvx.

3 Description

F07APF (ZGESVX) performs the following steps:

Equilibration
The linear system to be solved may be badly scaled. However, the system can be equilibrated as a first stage by setting $FACT ='E'$ . In this case, real scaling factors are computed and these factors then determine whether the system is to be equilibrated. Equilibrated forms of the systems $A X = B$ , $A^{T} X = B$ and $A^{H} X = B$ are
$(D_{R} A D_{C}) (D_{C}^{- 1} X) = D_{R} B,$
${(D_{R} A D_{C})}^{T} (D_{R}^{- 1} X) = D_{C} B,$
and
${(D_{R} A D_{C})}^{H} (D_{R}^{- 1} X) = D_{C} B,$
respectively, where $D_{R}$ and $D_{C}$ are diagonal matrices, with positive diagonal elements, formed from the computed scaling factors.

When equilibration is used, $A$ will be overwritten by $D_{R} A D_{C}$ and $B$ will be overwritten by $D_{R} B$ (or $D_{C} B$ when the solution of $A^{T} X = B$ or $A^{H} X = B$ is sought).
Factorization
The matrix $A$ , or its scaled form, is copied and factored using the $L U$ decomposition
$A = P L U,$
where $P$ is a permutation matrix, $L$ is a unit lower triangular matrix, and $U$ is upper triangular.

This stage can be by-passed when a factored matrix (with scaled matrices and scaling factors) are supplied; for example, as provided by a previous call to F07APF (ZGESVX) with the same matrix $A$ .
Condition Number Estimation
The $L U$ factorization of $A$ determines whether a solution to the linear system exists. If some diagonal element of $U$ is zero, then $U$ is exactly singular, no solution exists and the routine returns with a failure. Otherwise the factorized 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 then a warning code is returned on final exit.
Solution
The (equilibrated) system is solved for $X$ ( $D_{C}^{- 1} X$ or $D_{R}^{- 1} X$ ) using the factored form of $A$ ( $D_{R} A D_{C}$ ).
Iterative Refinement
Iterative refinement is applied to improve the computed solution matrix and to calculate error bounds and backward error estimates for the computed solution.
Construct Solution Matrix $X$
If equilibration was used, the matrix $X$ is premultiplied by $D_{C}$ (if $TRANS ='N'$ ) or $D_{R}$ (if $TRANS ='T'$ or $'C'$ ) 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'$: AF and IPIV contain the factorized form of $A$ . If $EQUED \neq'N'$ , the matrix $A$ has been equilibrated with scaling factors given by R and C. A, AF and IPIV are not modified.
$FACT ='N'$: The matrix $A$ will be copied to AF and factorized.
$FACT ='E'$: The matrix $A$ will be equilibrated if necessary, then copied to AF and factorized.

Constraint:

FACT ='F'

'N'

'E'

2: TRANS – CHARACTER(1)Input

On entry: specifies the form of the system of equations.

$TRANS ='N'$: $A X = B$ (No transpose).
$TRANS ='T'$: $A^{T} X = B$ (Transpose).
$TRANS ='C'$: $A^{H} X = B$ (Conjugate transpose).

Constraint:

TRANS ='N'

'T'

'C'

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: A(LDA, $*$ ) – COMPLEX (KIND=nag_wp) arrayInput/Output

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

\max (1, N)

On entry: the

n

n

matrix

A

FACT ='F'

and

EQUED \neq'N'

, A must have been equilibrated by the scaling factors in R and/or C.

On exit: if

FACT ='F'

'N'

, or if

FACT ='E'

and

EQUED ='N'

, A is not modified.

FACT ='E'

EQUED \neq'N'

A

is scaled as follows:

if $EQUED ='R'$ , $A = D_{R} A$ ;
if $EQUED ='C'$ , $A = A D_{C}$ ;
if $EQUED ='B'$ , $A = D_{R} A D_{C}$ .

6: LDA – INTEGERInput

On entry: the first dimension of the array A as declared in the (sub)program from which F07APF (ZGESVX) is called.

Constraint:

LDA \geq \max (1, N)

7: AF(LDAF, $*$ ) – COMPLEX (KIND=nag_wp) arrayInput/Output

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

\max (1, N)

On entry: if

FACT ='F'

, AF contains the factors

L

and

U

from the factorization

A = P L U

as computed by F07ARF (ZGETRF). If

EQUED \neq'N'

, AF is the factorized form of the equilibrated matrix

A

FACT ='N'

'E'

, AF need not be set.

On exit: if

FACT ='N'

, AF returns the factors

L

and

U

from the factorization

A = P L U

of the original matrix

A

FACT ='E'

, AF returns the factors

L

and

U

from the factorization

A = P L U

of the equilibrated matrix

A

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

FACT ='F'

, AF is unchanged from entry.

8: LDAF – INTEGERInput

On entry: the first dimension of the array AF as declared in the (sub)program from which F07APF (ZGESVX) is called.

Constraint:

LDAF \geq \max (1, N)

9: IPIV( $*$ ) – INTEGER arrayInput/Output

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

\max (1, N)

On entry: if

FACT ='F'

, IPIV contains the pivot indices from the factorization

A = P L U

as computed by F07ARF (ZGETRF); at the

i

th step row

i

of the matrix was interchanged with row

IPIV (i)

IPIV (i) = i

indicates a row interchange was not required.

FACT ='N'

'E'

, IPIV need not be set.

On exit: if

FACT ='N'

, IPIV contains the pivot indices from the factorization

A = P L U

of the original matrix

A

FACT ='E'

, IPIV contains the pivot indices from the factorization

A = P L U

of the equilibrated matrix

A

FACT ='F'

, IPIV is unchanged from entry.

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 ='R'$ , row equilibration, i.e., $A$ has been premultiplied by $D_{R}$ ;
if $EQUED ='C'$ , column equilibration, i.e., $A$ has been postmultiplied by $D_{C}$ ;
if $EQUED ='B'$ , both row and column equilibration, i.e., $A$ has been replaced by $D_{R} A D_{C}$ .

On exit: if

FACT ='F'

, EQUED is unchanged from entry.

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

Constraint: if

FACT ='F'

EQUED ='N'

'R'

'C'

'B'

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

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

\max (1, N)

On entry: if

FACT ='N'

'E'

, R need not be set.

FACT ='F'

and

EQUED ='R'

'B'

, R must contain the row scale factors for

A

D_{R}

; each element of R must be positive.

On exit: if

FACT ='F'

, R is unchanged from entry.

Otherwise, if no constraints are violated and

EQUED ='R'

'B'

, R contains the row scale factors for

A

D_{R}

, such that

A

is multiplied on the left by

D_{R}

; each element of R is positive.

12: C( $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

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

\max (1, N)

On entry: if

FACT ='N'

'E'

, C need not be set.

FACT ='F'

EQUED ='C'

'B'

, C must contain the column scale factors for

A

D_{C}

; each element of C must be positive.

On exit: if

FACT ='F'

, C is unchanged from entry.

Otherwise, if no constraints are violated and

EQUED ='C'

'B'

, C contains the row scale factors for

A

D_{C}

; each element of C is positive.

13: B(LDB, $*$ ) – COMPLEX (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.

TRANS ='N'

and

EQUED ='R'

'B'

, B is overwritten by

D_{R} B

TRANS ='T'

'C'

and

EQUED ='C'

'B'

, B is overwritten by

D_{C} B

14: LDB – INTEGERInput

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

Constraint:

LDB \geq \max (1, N)

15: X(LDX, $*$ ) – COMPLEX (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 \neq'N'

, and the solution to the equilibrated system is

D_{C}^{- 1} X

TRANS ='N'

and

EQUED ='C'

'B'

, or

D_{R}^{- 1} X

TRANS ='T'

'C'

and

EQUED ='R'

'B'

16: LDX – INTEGERInput

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

Constraint:

LDX \geq \max (1, N)

17: 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})

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

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

20: WORK( $2 \times N$ ) – COMPLEX (KIND=nag_wp) arrayWorkspace

21: RWORK( $\max (1, 2 \times N)$ ) – REAL (KIND=nag_wp) arrayOutput

On exit:

RWORK (1)

contains the reciprocal pivot growth factor

‖A‖ / ‖U‖

. The ‘max absolute element’ norm is used. If

RWORK (1)

is much less than

1

, then the stability of the

L U

factorization of the (equilibrated) matrix

A

could be poor. This also means that the solution X, condition estimator RCOND, and forward error bound FERR could be unreliable. If factorization fails with

INFO > 0 and INFO \leq N

, then

RWORK (1)

contains the reciprocal pivot growth factor for the leading INFO columns of

A

22: 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$ , $u_{i i}$ is exactly zero. The factorization has been completed, but the factor $U$ is exactly singular, so the solution and error bounds could not be computed. $RCOND = 0.0$ is returned.

$INFO = N + 1$: The triangular matrix $U$ 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| \leq c (n) ε P |L| |U|,

c (n)

is a modest linear function of

n

, and

ε

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

x

is the true solution, then the computed solution

\hat{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{8}{3} n^{3}

floating point operations.

Estimating the forward error involves solving a number of systems of linear equations of the form

A x = b

A^{T} x = b

; the number is usually

4

5

and never more than

11

. Each solution involves approximately

8 n^{2}

operations.

In practice the condition number estimator is very reliable, but it can underestimate the true condition number; see Section 15.3 of Higham (2002) for further details.

The real analogue of this routine is F07ABF (DGESVX).

9 Example

This example solves the equations

A X = B,

where

A

is the general matrix

A = (\begin{array}{r} - 1.34 + 2.55 i & 0.28 + 3.17 i & - 6.39 - 2.20 i & 0.72 - 0.92 i \\ - 1.70 - 14.10 i & 33.10 - 1.50 i & - 1.50 + 13.40 i & 12.90 + 13.80 i \\ - 3.29 - 2.39 i & - 1.91 + 4.42 i & - 0.14 - 1.35 i & 1.72 + 1.35 i \\ 2.41 + 0.39 i & - 0.56 + 1.47 i & - 0.83 - 0.69 i & - 1.96 + 0.67 i \end{array})

and

B = (\begin{array}{r} 26.26 + 51.78 i & 31.32 - 6.70 i \\ 64.30 - 86.80 i & 158.60 - 14.20 i \\ - 5.75 + 25.31 i & - 2.15 + 30.19 i \\ 1.16 + 2.57 i & - 2.56 + 7.55 i \end{array}) .

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

A

and an estimate of the reciprocal of the pivot growth factor for the factorization of

A

are also output.

NAG Library Routine DocumentF07APF (ZGESVX)

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

F07APF (ZGESVX)