NAG Library Function Document

nag_dgesvx (f07abc)

nag_dgesvx (Nag_OrderType order, Nag_FactoredFormType fact, Nag_TransType trans, Integer n, Integer nrhs, double a[], Integer pda, double af[], Integer pdaf, Integer ipiv[], Nag_EquilibrationType *equed, double r[], double c[], double b[], Integer pdb, double x[], Integer pdx, double *rcond, double ferr[], double berr[], double *recip_growth_factor, NagError *fail)

3 Description

nag_dgesvx (f07abc) 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 = Nag_EquilibrateAndFactor$ . 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$ and $A^{T} X = B$ are
$(D_{R} A D_{C}) (D_{C}^{- 1} X) = D_{R} B$
and
${(D_{R} A D_{C})}^{T} (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$ 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 nag_dgesvx (f07abc) 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 function 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 = Nag_NoTrans$ ) or $D_{R}$ (if $trans = Nag_Trans$ or $Nag_ConjTrans$ ) 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 Arguments

1: order – Nag_OrderTypeInput

On entry: the order argument specifies the two-dimensional storage scheme being used, i.e., row-major ordering or column-major ordering. C language defined storage is specified by

order = Nag_RowMajor

. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.

Constraint:

order = Nag_RowMajor

or Nag_ColMajor.

2: fact – Nag_FactoredFormTypeInput

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 = Nag_Factored$: af and ipiv contain the factorized form of $A$ . If $equed \neq Nag_NoEquilibration$ , the matrix $A$ has been equilibrated with scaling factors given by r and c. a, af and ipiv are not modified.
$fact = Nag_NotFactored$: The matrix $A$ will be copied to af and factorized.
$fact = Nag_EquilibrateAndFactor$: The matrix $A$ will be equilibrated if necessary, then copied to af and factorized.

Constraint:

fact = Nag_Factored

Nag_NotFactored

Nag_EquilibrateAndFactor

3: trans – Nag_TransTypeInput

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

$trans = Nag_NoTrans$: $A X = B$ (No transpose).
$trans = Nag_Trans$ or $Nag_ConjTrans$: $A^{T} X = B$ (Transpose).

Constraint:

trans = Nag_NoTrans

Nag_Trans

Nag_ConjTrans

4: n – IntegerInput

On entry:

n

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

A

Constraint:

n \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: a[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array a must be at least

\max (1, pda \times n)

The

(i, j)

th element of the matrix

A

is stored in

$a [(j - 1) \times pda + i - 1]$ when $order = Nag_ColMajor$ ;
$a [(i - 1) \times pda + j - 1]$ when $order = Nag_RowMajor$ .

On entry: the

n

n

matrix

A

fact = Nag_Factored

and

equed \neq Nag_NoEquilibration

, a must have been equilibrated by the scaling factors in r and/or c.

On exit: if

fact = Nag_Factored

Nag_NotFactored

, or if

fact = Nag_EquilibrateAndFactor

and

equed = Nag_NoEquilibration

, a is not modified.

fact = Nag_EquilibrateAndFactor

equed \neq Nag_NoEquilibration

A

is scaled as follows:

if $equed = Nag_RowEquilibration$ , $A = D_{R} A$ ;
if $equed = Nag_ColumnEquilibration$ , $A = A D_{C}$ ;
if $equed = Nag_RowAndColumnEquilibration$ , $A = D_{R} A D_{C}$ .

7: pda – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array a.

Constraint:

pda \geq \max (1, n)

8: af[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array af must be at least

\max (1, pdaf \times n)

The

(i, j)

th element of the matrix is stored in

$af [(j - 1) \times pdaf + i - 1]$ when $order = Nag_ColMajor$ ;
$af [(i - 1) \times pdaf + j - 1]$ when $order = Nag_RowMajor$ .

On entry: if

fact = Nag_Factored

, af contains the factors

L

and

U

from the factorization

A = P L U

as computed by nag_dgetrf (f07adc). If

equed \neq Nag_NoEquilibration

, af is the factorized form of the equilibrated matrix

A

fact = Nag_NotFactored

Nag_EquilibrateAndFactor

, af need not be set.

On exit: if

fact = Nag_NotFactored

, af returns the factors

L

and

U

from the factorization

A = P L U

of the original matrix

A

fact = Nag_EquilibrateAndFactor

, 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 = Nag_Factored

, af is unchanged from entry.

9: pdaf – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array af.

Constraint:

pdaf \geq \max (1, n)

10: ipiv[ $\dim$ ] – IntegerInput/Output

Note: the dimension, dim, of the array ipiv must be at least

\max (1, n)

On entry: if

fact = Nag_Factored

, ipiv contains the pivot indices from the factorization

A = P L U

as computed by nag_dgetrf (f07adc); at the

i

th step row

i

of the matrix was interchanged with row

ipiv [i - 1]

ipiv [i - 1] = i

indicates a row interchange was not required.

fact = Nag_NotFactored

Nag_EquilibrateAndFactor

, ipiv need not be set.

On exit: if

fact = Nag_NotFactored

, ipiv contains the pivot indices from the factorization

A = P L U

of the original matrix

A

fact = Nag_EquilibrateAndFactor

, ipiv contains the pivot indices from the factorization

A = P L U

of the equilibrated matrix

A

fact = Nag_Factored

, ipiv is unchanged from entry.

11: equed – Nag_EquilibrationType*Input/Output

On entry: if

fact = Nag_NotFactored

Nag_EquilibrateAndFactor

, equed need not be set.

fact = Nag_Factored

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

if $equed = Nag_NoEquilibration$ , no equilibration;
if $equed = Nag_RowEquilibration$ , row equilibration, i.e., $A$ has been premultiplied by $D_{R}$ ;
if $equed = Nag_ColumnEquilibration$ , column equilibration, i.e., $A$ has been postmultiplied by $D_{C}$ ;
if $equed = Nag_RowAndColumnEquilibration$ , both row and column equilibration, i.e., $A$ has been replaced by $D_{R} A D_{C}$ .

On exit: if

fact = Nag_Factored

, 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 = Nag_Factored

equed = Nag_NoEquilibration

Nag_RowEquilibration

Nag_ColumnEquilibration

Nag_RowAndColumnEquilibration

12: r[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array r must be at least

\max (1, n)

On entry: if

fact = Nag_NotFactored

Nag_EquilibrateAndFactor

, r need not be set.

fact = Nag_Factored

and

equed = Nag_RowEquilibration

Nag_RowAndColumnEquilibration

, r must contain the row scale factors for

A

D_{R}

; each element of r must be positive.

On exit: if

fact = Nag_Factored

, r is unchanged from entry.

Otherwise, if no constraints are violated and

equed = Nag_RowEquilibration

Nag_RowAndColumnEquilibration

, 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.

13: c[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array c must be at least

\max (1, n)

On entry: if

fact = Nag_NotFactored

Nag_EquilibrateAndFactor

, c need not be set.

fact = Nag_Factored

equed = Nag_ColumnEquilibration

Nag_RowAndColumnEquilibration

, c must contain the column scale factors for

A

D_{C}

; each element of c must be positive.

On exit: if

fact = Nag_Factored

, c is unchanged from entry.

Otherwise, if no constraints are violated and

equed = Nag_ColumnEquilibration

Nag_RowAndColumnEquilibration

, c contains the row scale factors for

A

D_{C}

; each element of c is positive.

14: b[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array b must be at least

$\max (1, pdb \times nrhs)$ when $order = Nag_ColMajor$ ;
$\max (1, n \times pdb)$ when $order = Nag_RowMajor$ .

The

(i, j)

th element of the matrix

B

is stored in

$b [(j - 1) \times pdb + i - 1]$ when $order = Nag_ColMajor$ ;
$b [(i - 1) \times pdb + j - 1]$ when $order = Nag_RowMajor$ .

On entry: the

n

r

right-hand side matrix

B

On exit: if

equed = Nag_NoEquilibration

, b is not modified.

trans = Nag_NoTrans

and

equed = Nag_RowEquilibration

Nag_RowAndColumnEquilibration

, b is overwritten by

D_{R} B

trans = Nag_Trans

Nag_ConjTrans

and

equed = Nag_ColumnEquilibration

Nag_RowAndColumnEquilibration

, b is overwritten by

D_{C} B

15: pdb – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array b.

Constraints:

if $order = Nag_ColMajor$ , $pdb \geq \max (1, n)$ ;
if $order = Nag_RowMajor$ , $pdb \geq \max (1, nrhs)$ .

16: x[ $\dim$ ] – doubleOutput

Note: the dimension, dim, of the array x must be at least

$\max (1, pdx \times nrhs)$ when $order = Nag_ColMajor$ ;
$\max (1, n \times pdx)$ when $order = Nag_RowMajor$ .

The

(i, j)

th element of the matrix

X

is stored in

$x [(j - 1) \times pdx + i - 1]$ when $order = Nag_ColMajor$ ;
$x [(i - 1) \times pdx + j - 1]$ when $order = Nag_RowMajor$ .

On exit: if

fail . code =

NE_NOERROR or NE_SINGULAR_WP, 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 Nag_NoEquilibration

, and the solution to the equilibrated system is

D_{C}^{- 1} X

trans = Nag_NoTrans

and

equed = Nag_ColumnEquilibration

Nag_RowAndColumnEquilibration

, or

D_{R}^{- 1} X

trans = Nag_Trans

Nag_ConjTrans

and

equed = Nag_RowEquilibration

Nag_RowAndColumnEquilibration

17: pdx – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array x.

Constraints:

if $order = Nag_ColMajor$ , $pdx \geq \max (1, n)$ ;
if $order = Nag_RowMajor$ , $pdx \geq \max (1, nrhs)$ .

18: rcond – double *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})

19: ferr[nrhs] – doubleOutput

On exit: if

fail . code =

NE_NOERROR or NE_SINGULAR_WP, 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 - 1]

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.

20: berr[nrhs] – doubleOutput

On exit: if

fail . code =

NE_NOERROR or NE_SINGULAR_WP, 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).

21: recip_growth_factor – double *Output

On exit: if

fail . code =

NE_NOERROR, the reciprocal pivot growth factor

‖A‖ / ‖U‖

, where

‖.‖

denotes the maximum absolute element norm. If

recip_growth_factor ≪ 1

, the stability of the

L U

factorization of (equilibrated)

A

could be poor. This also means that the solution x, condition estimate rcond, and forward error bound ferr could be unreliable. If the factorization fails with

fail . code =

NE_SINGULAR, then

recip_growth_factor

contains the reciprocal pivot growth factor for the leading

fail . errnum

columns of

A

22: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $〈value〉$ had an illegal value.
NE_INT: On entry, $n = 〈value〉$ .
Constraint: $n \geq 0$ .; On entry, $nrhs = 〈value〉$ .
Constraint: $nrhs \geq 0$ .; On entry, $pda = 〈value〉$ .
Constraint: $pda > 0$ .; On entry, $pdaf = 〈value〉$ .
Constraint: $pdaf > 0$ .; On entry, $pdb = 〈value〉$ .
Constraint: $pdb > 0$ .; On entry, $pdx = 〈value〉$ .
Constraint: $pdx > 0$ .
NE_INT_2: On entry, $pda = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pda \geq \max (1, n)$ .; On entry, $pdaf = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pdaf \geq \max (1, n)$ .; On entry, $pdb = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pdb \geq \max (1, n)$ .; On entry, $pdb = 〈value〉$ and $nrhs = 〈value〉$ .
Constraint: $pdb \geq \max (1, nrhs)$ .; On entry, $pdx = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pdx \geq \max (1, n)$ .; On entry, $pdx = 〈value〉$ and $nrhs = 〈value〉$ .
Constraint: $pdx \geq \max (1, nrhs)$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
NE_SINGULAR: $U (〈value〉, 〈value〉)$ 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.
NE_SINGULAR_WP: $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 - 1]

and a bound on

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

is returned in

ferr [j - 1]

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

8 Further Comments

The factorization of

A

requires approximately

\frac{2}{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

2 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 complex analogue of this function is nag_zgesvx (f07apc).

9 Example

This example solves the equations

A X = B,

where

A

is the general matrix

A = (\begin{matrix} 1.80 & 2.88 & 2.05 & - 0.89 \\ 525.00 & - 295.00 & - 95.00 & - 380.00 \\ 1.58 & - 2.69 & - 2.90 & - 1.04 \\ - 1.11 & - 0.66 & - 0.59 & - 0.80 \end{matrix})

and

B = (\begin{matrix} 9.52 & 18.47 \\ 2435.00 & 225.00 \\ 0.77 & - 13.28 \\ - 6.22 & - 6.21 \end{matrix}) .

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 Function Documentnag_dgesvx (f07abc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

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 Function Document

nag_dgesvx (f07abc)