NAG Library Routine Document
F07ABF (DGESVX)
1 Purpose
F07ABF (DGESVX) uses the
$LU$ factorization to compute the solution to a real system of linear equations
where
$A$ is an
$n$ by
$n$ matrix and
$X$ and
$B$ are
$n$ by
$r$ matrices. Error bounds on the solution and a condition estimate are also provided.
2 Specification
SUBROUTINE F07ABF ( 
FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO) 
INTEGER 
N, NRHS, LDA, LDAF, IPIV(*), LDB, LDX, IWORK(N), INFO 
REAL (KIND=nag_wp) 
A(LDA,*), AF(LDAF,*), R(*), C(*), B(LDB,*), X(LDX,*), RCOND, FERR(NRHS), BERR(NRHS), WORK(max(1,4*N)) 
CHARACTER(1) 
FACT, TRANS, EQUED 

The routine may be called by its
LAPACK
name dgesvx.
3 Description
F07ABF (DGESVX) 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
${\mathbf{FACT}}=\text{'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
$AX=B$ and
${A}^{\mathrm{T}}X=B$ are
and
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}^{\mathrm{T}}X=B$ is sought).
 Factorization
The matrix
$A$, or its scaled form, is copied and factored using the
$LU$ decomposition
where
$P$ is a permutation matrix,
$L$ is a unit lower triangular matrix, and
$U$ is upper triangular.
This stage can be bypassed when a factored matrix (with scaled matrices and scaling factors) are supplied; for example, as provided by a previous call to F07ABF (DGESVX) with the same matrix $A$.
 Condition Number Estimation
The $LU$ 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 ${\mathbf{TRANS}}=\text{'N'}$) or ${D}_{R}$ (if ${\mathbf{TRANS}}=\text{'T'}$ or $\text{'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.
 ${\mathbf{FACT}}=\text{'F'}$
 AF and IPIV contain the factorized form of $A$. If ${\mathbf{EQUED}}\ne \text{'N'}$, the matrix $A$ has been equilibrated with scaling factors given by R and C. A, AF and IPIV are not modified.
 ${\mathbf{FACT}}=\text{'N'}$
 The matrix $A$ will be copied to AF and factorized.
 ${\mathbf{FACT}}=\text{'E'}$
 The matrix $A$ will be equilibrated if necessary, then copied to AF and factorized.
Constraint:
${\mathbf{FACT}}=\text{'F'}$, $\text{'N'}$ or $\text{'E'}$.
 2: TRANS – CHARACTER(1)Input
On entry: specifies the form of the system of equations.
 ${\mathbf{TRANS}}=\text{'N'}$
 $AX=B$ (No transpose).
 ${\mathbf{TRANS}}=\text{'T'}$ or $\text{'C'}$
 ${A}^{\mathrm{T}}X=B$ (Transpose).
Constraint:
${\mathbf{TRANS}}=\text{'N'}$, $\text{'T'}$ or $\text{'C'}$.
 3: N – INTEGERInput
On entry: $n$, the number of linear equations, i.e., the order of the matrix $A$.
Constraint:
${\mathbf{N}}\ge 0$.
 4: NRHS – INTEGERInput
On entry: $r$, the number of righthand sides, i.e., the number of columns of the matrix $B$.
Constraint:
${\mathbf{NRHS}}\ge 0$.
 5: A(LDA,$*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array
A
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: the
$n$ by
$n$ matrix
$A$.
If
${\mathbf{FACT}}=\text{'F'}$ and
${\mathbf{EQUED}}\ne \text{'N'}$,
A must have been equilibrated by the scaling factors in
R and/or
C.
On exit: if
${\mathbf{FACT}}=\text{'F'}$ or
$\text{'N'}$, or if
${\mathbf{FACT}}=\text{'E'}$ and
${\mathbf{EQUED}}=\text{'N'}$,
A is not modified.
If
${\mathbf{FACT}}=\text{'E'}$ or
${\mathbf{EQUED}}\ne \text{'N'}$,
$A$ is scaled as follows:
 if ${\mathbf{EQUED}}=\text{'R'}$, $A={D}_{R}A$;
 if ${\mathbf{EQUED}}=\text{'C'}$, $A=A{D}_{C}$;
 if ${\mathbf{EQUED}}=\text{'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 F07ABF (DGESVX) is called.
Constraint:
${\mathbf{LDA}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
 7: AF(LDAF,$*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array
AF
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: if
${\mathbf{FACT}}=\text{'F'}$,
AF contains the factors
$L$ and
$U$ from the factorization
$A=PLU$ as computed by
F07ADF (DGETRF). If
${\mathbf{EQUED}}\ne \text{'N'}$,
AF is the factorized form of the equilibrated matrix
$A$.
If
${\mathbf{FACT}}=\text{'N'}$ or
$\text{'E'}$,
AF need not be set.
On exit: if
${\mathbf{FACT}}=\text{'N'}$,
AF returns the factors
$L$ and
$U$ from the factorization
$A=PLU$ of the original matrix
$A$.
If
${\mathbf{FACT}}=\text{'E'}$,
AF returns the factors
$L$ and
$U$ from the factorization
$A=PLU$ of the equilibrated matrix
$A$ (see the description of
A for the form of the equilibrated matrix).
If
${\mathbf{FACT}}=\text{'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 F07ABF (DGESVX) is called.
Constraint:
${\mathbf{LDAF}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
 9: IPIV($*$) – INTEGER arrayInput/Output

Note: the dimension of the array
IPIV
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: if
${\mathbf{FACT}}=\text{'F'}$,
IPIV contains the pivot indices from the factorization
$A=PLU$ as computed by
F07ADF (DGETRF); at the
$i$th step row
$i$ of the matrix was interchanged with row
${\mathbf{IPIV}}\left(i\right)$.
${\mathbf{IPIV}}\left(i\right)=i$ indicates a row interchange was not required.
If
${\mathbf{FACT}}=\text{'N'}$ or
$\text{'E'}$,
IPIV need not be set.
On exit: if
${\mathbf{FACT}}=\text{'N'}$,
IPIV contains the pivot indices from the factorization
$A=PLU$ of the original matrix
$A$.
If
${\mathbf{FACT}}=\text{'E'}$,
IPIV contains the pivot indices from the factorization
$A=PLU$ of the equilibrated matrix
$A$.
If
${\mathbf{FACT}}=\text{'F'}$,
IPIV is unchanged from entry.
 10: EQUED – CHARACTER(1)Input/Output
On entry: if
${\mathbf{FACT}}=\text{'N'}$ or
$\text{'E'}$,
EQUED need not be set.
If
${\mathbf{FACT}}=\text{'F'}$,
EQUED must specify the form of the equilibration that was performed as follows:
 if ${\mathbf{EQUED}}=\text{'N'}$, no equilibration;
 if ${\mathbf{EQUED}}=\text{'R'}$, row equilibration, i.e., $A$ has been premultiplied by ${D}_{R}$;
 if ${\mathbf{EQUED}}=\text{'C'}$, column equilibration, i.e., $A$ has been postmultiplied by ${D}_{C}$;
 if ${\mathbf{EQUED}}=\text{'B'}$, both row and column equilibration, i.e., $A$ has been replaced by ${D}_{R}A{D}_{C}$.
On exit: if
${\mathbf{FACT}}=\text{'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 ${\mathbf{FACT}}=\text{'F'}$, ${\mathbf{EQUED}}=\text{'N'}$, $\text{'R'}$, $\text{'C'}$ or $\text{'B'}$.
 11: R($*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the dimension of the array
R
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: if
${\mathbf{FACT}}=\text{'N'}$ or
$\text{'E'}$,
R need not be set.
If
${\mathbf{FACT}}=\text{'F'}$ and
${\mathbf{EQUED}}=\text{'R'}$ or
$\text{'B'}$,
R must contain the row scale factors for
$A$,
${D}_{R}$; each element of
R must be positive.
On exit: if
${\mathbf{FACT}}=\text{'F'}$,
R is unchanged from entry.
Otherwise, if no constraints are violated and
${\mathbf{EQUED}}=\text{'R'}$ or
$\text{'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
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
On entry: if
${\mathbf{FACT}}=\text{'N'}$ or
$\text{'E'}$,
C need not be set.
If
${\mathbf{FACT}}=\text{'F'}$ or
${\mathbf{EQUED}}=\text{'C'}$ or
$\text{'B'}$,
C must contain the column scale factors for
$A$,
${D}_{C}$; each element of
C must be positive.
On exit: if
${\mathbf{FACT}}=\text{'F'}$,
C is unchanged from entry.
Otherwise, if no constraints are violated and
${\mathbf{EQUED}}=\text{'C'}$ or
$\text{'B'}$,
C contains the row scale factors for
$A$,
${D}_{C}$; each element of
C is positive.
 13: B(LDB,$*$) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array
B
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{NRHS}}\right)$.
On entry: the $n$ by $r$ righthand side matrix $B$.
On exit: if
${\mathbf{EQUED}}=\text{'N'}$,
B is not modified.
If
${\mathbf{TRANS}}=\text{'N'}$ and
${\mathbf{EQUED}}=\text{'R'}$ or
$\text{'B'}$,
B is overwritten by
${D}_{R}B$.
If
${\mathbf{TRANS}}=\text{'T'}$ or
$\text{'C'}$ and
${\mathbf{EQUED}}=\text{'C'}$ or
$\text{'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 F07ABF (DGESVX) is called.
Constraint:
${\mathbf{LDB}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
 15: X(LDX,$*$) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array
X
must be at least
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{NRHS}}\right)$.
On exit: if ${\mathbf{INFO}}={\mathbf{0}}$ or ${\mathbf{N}+{\mathbf{1}}}$, the $n$ by $r$ solution matrix $X$ to the original system of equations. Note that the arrays $A$ and $B$ are modified on exit if ${\mathbf{EQUED}}\ne \text{'N'}$, and the solution to the equilibrated system is ${D}_{C}^{1}X$ if ${\mathbf{TRANS}}=\text{'N'}$ and ${\mathbf{EQUED}}=\text{'C'}$ or $\text{'B'}$, or ${D}_{R}^{1}X$ if ${\mathbf{TRANS}}=\text{'T'}$ or $\text{'C'}$ and ${\mathbf{EQUED}}=\text{'R'}$ or $\text{'B'}$.
 16: LDX – INTEGERInput
On entry: the first dimension of the array
X as declared in the (sub)program from which F07ABF (DGESVX) is called.
Constraint:
${\mathbf{LDX}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{N}}\right)$.
 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 ${\mathbf{RCOND}}=1.0/\left({\Vert A\Vert}_{1}{\Vert {A}^{1}\Vert}_{1}\right)$.
 18: FERR(NRHS) – REAL (KIND=nag_wp) arrayOutput
On exit: if
${\mathbf{INFO}}={\mathbf{0}}$ or
${\mathbf{N}+{\mathbf{1}}}$, an estimate of the forward error bound for each computed solution vector, such that
${\Vert {\hat{x}}_{j}{x}_{j}\Vert}_{\infty}/{\Vert {x}_{j}\Vert}_{\infty}\le {\mathbf{FERR}}\left(j\right)$ 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 ${\mathbf{INFO}}={\mathbf{0}}$ or ${\mathbf{N}+{\mathbf{1}}}$, an estimate of the componentwise relative backward error of each computed solution vector ${\hat{x}}_{j}$ (i.e., the smallest relative change in any element of $A$ or $B$ that makes ${\hat{x}}_{j}$ an exact solution).
 20: WORK($\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,4\times {\mathbf{N}}\right)$) – REAL (KIND=nag_wp) arrayOutput
On exit:
${\mathbf{WORK}}\left(1\right)$ contains the reciprocal pivot growth factor
$\Vert A\Vert /\Vert U\Vert $. The ‘max absolute element’ norm is used. If
${\mathbf{WORK}}\left(1\right)$ is much less than
$1$, then the stability of the
$LU$ factorization of the (equilibrated) matrix
$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
$\mathbf{INFO}>{\mathbf{0}}\text{and}\mathbf{INFO}\le \mathbf{N}$, then
${\mathbf{WORK}}\left(1\right)$ contains the reciprocal pivot growth factor for the leading
INFO columns of
$A$.
 21: IWORK(N) – INTEGER arrayWorkspace
 22: INFO – INTEGEROutput
On exit:
${\mathbf{INFO}}=0$ unless the routine detects an error (see
Section 6).
6 Error Indicators and Warnings
Errors or warnings detected by the routine:
 ${\mathbf{INFO}}<0$
If ${\mathbf{INFO}}=i$, the $i$th argument had an illegal value. An explanatory message is output, and execution of the program is terminated.
 ${\mathbf{INFO}}>0\text{and}{\mathbf{INFO}}\le {\mathbf{N}}$
If ${\mathbf{INFO}}=i$, ${u}_{ii}$ 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. ${\mathbf{RCOND}}=0.0$ is returned.
 ${\mathbf{INFO}}={\mathbf{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 righthand side vector
$b$, the computed solution
$\hat{x}$ is the exact solution of a perturbed system of equations
$\left(A+E\right)\hat{x}=b$, where
$c\left(n\right)$ is a modest linear function of
$n$, and
$\epsilon $ is the
machine precision. See Section 9.3 of
Higham (2002) for further details.
If
$x$ is the true solution, then the computed solution
$\hat{x}$ satisfies a forward error bound of the form
where
$\mathrm{cond}\left(A,\hat{x},b\right)={\Vert \left{A}^{1}\right\left(\leftA\right\left\hat{x}\right+\leftb\right\right)\Vert}_{\infty}/{\Vert \hat{x}\Vert}_{\infty}\le \mathrm{cond}\left(A\right)={\Vert \left{A}^{1}\right\leftA\right\Vert}_{\infty}\le {\kappa}_{\infty}\left(A\right)$.
If
$\hat{x}$ is the
$j$th column of
$X$, then
${w}_{c}$ is returned in
${\mathbf{BERR}}\left(j\right)$ and a bound on
${\Vert x\hat{x}\Vert}_{\infty}/{\Vert \hat{x}\Vert}_{\infty}$ is returned in
${\mathbf{FERR}}\left(j\right)$. See Section 4.4 of
Anderson et al. (1999) for further details.
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 $Ax=b$ or ${A}^{\mathrm{T}}x=b$; the number is usually $4$ or $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 routine is
F07APF (ZGESVX).
9 Example
This example solves the equations
where
$A$ is the general matrix
and
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.
9.1 Program Text
Program Text (f07abfe.f90)
9.2 Program Data
Program Data (f07abfe.d)
9.3 Program Results
Program Results (f07abfe.r)