Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_lapack_dgesvx (f07ab)

## Purpose

nag_lapack_dgesvx (f07ab) uses the $LU$ factorization to compute the solution to a real system of linear equations
 $AX=B or ATX=B ,$
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.

## Syntax

[a, af, ipiv, equed, r, c, b, x, rcond, ferr, berr, work, info] = f07ab(fact, trans, a, af, ipiv, equed, r, c, b, 'n', n, 'nrhs_p', nrhs_p)
[a, af, ipiv, equed, r, c, b, x, rcond, ferr, berr, work, info] = nag_lapack_dgesvx(fact, trans, a, af, ipiv, equed, r, c, b, 'n', n, 'nrhs_p', nrhs_p)

## Description

nag_lapack_dgesvx (f07ab) performs the following steps:
1. 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
 $DR A DC DC-1X = DR B$
and
 $DR A DC T DR-1 X = DC 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}^{\mathrm{T}}X=B$ is sought).
2. Factorization
The matrix $A$, or its scaled form, is copied and factored using the $LU$ decomposition
 $A=PLU ,$
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_lapack_dgesvx (f07ab) with the same matrix $A$.
3. 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 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.
4. 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}$).
5. 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.
6. 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.

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

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{fact}$ – string (length ≥ 1)
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:     $\mathrm{trans}$ – string (length ≥ 1)
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:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array
The first dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
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.
4:     $\mathrm{af}\left(\mathit{ldaf},:\right)$ – double array
The first dimension of the array af must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array af must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If ${\mathbf{fact}}=\text{'F'}$, af contains the factors $L$ and $U$ from the factorization $A=PLU$ as computed by nag_lapack_dgetrf (f07ad). 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.
5:     $\mathrm{ipiv}\left(:\right)$int64int32nag_int array
The dimension of the array ipiv must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
If ${\mathbf{fact}}=\text{'F'}$, ipiv contains the pivot indices from the factorization $A=PLU$ as computed by nag_lapack_dgetrf (f07ad); 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.
6:     $\mathrm{equed}$ – string (length ≥ 1)
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}$.
Constraint: if ${\mathbf{fact}}=\text{'F'}$, ${\mathbf{equed}}=\text{'N'}$, $\text{'R'}$, $\text{'C'}$ or $\text{'B'}$.
7:     $\mathrm{r}\left(:\right)$ – double array
The dimension of the array r must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
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.
8:     $\mathrm{c}\left(:\right)$ – double array
The dimension of the array c must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
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.
9:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – double array
The first dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs_p}}\right)$.
The $n$ by $r$ right-hand side matrix $B$.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the arrays a, af, b and the second dimension of the arrays a, af, ipiv, r, c.
$n$, the number of linear equations, i.e., the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
2:     $\mathrm{nrhs_p}$int64int32nag_int scalar
Default: the second dimension of the array b.
$r$, the number of right-hand sides, i.e., the number of columns of the matrix $B$.
Constraint: ${\mathbf{nrhs_p}}\ge 0$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double array
The first dimension of the array a will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array a will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
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}$.
2:     $\mathrm{af}\left(\mathit{ldaf},:\right)$ – double array
The first dimension of the array af will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array af will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
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.
3:     $\mathrm{ipiv}\left(:\right)$int64int32nag_int array
The dimension of the array ipiv will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
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.
4:     $\mathrm{equed}$ – string (length ≥ 1)
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.
5:     $\mathrm{r}\left(:\right)$ – double array
The dimension of the array r will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
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.
6:     $\mathrm{c}\left(:\right)$ – double array
The dimension of the array c will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
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.
7:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – double array
The first dimension of the array b will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array b will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs_p}}\right)$.
If ${\mathbf{equed}}=\text{'N'}$, b is not modified.
If ${\mathbf{trans}}=\text{'N'}$ and ${\mathbf{equed}}=\text{'R'}$ or $\text{'B'}$, b stores ${D}_{R}B$.
If ${\mathbf{trans}}=\text{'T'}$ or $\text{'C'}$ and ${\mathbf{equed}}=\text{'C'}$ or $\text{'B'}$, b stores ${D}_{C}B$.
8:     $\mathrm{x}\left(\mathit{ldx},:\right)$ – double array
The first dimension of the array x will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array x will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs_p}}\right)$.
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'}$.
9:     $\mathrm{rcond}$ – double scalar
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({‖A‖}_{1}{‖{A}^{-1}‖}_{1}\right)$.
10:   $\mathrm{ferr}\left({\mathbf{nrhs_p}}\right)$ – double array
If ${\mathbf{info}}={\mathbf{0}}$ or $\mathbf{n}+{\mathbf{1}}$, an estimate of the forward error bound for each computed solution vector, such that ${‖{\stackrel{^}{x}}_{j}-{x}_{j}‖}_{\infty }/{‖{x}_{j}‖}_{\infty }\le {\mathbf{ferr}}\left(j\right)$ where ${\stackrel{^}{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.
11:   $\mathrm{berr}\left({\mathbf{nrhs_p}}\right)$ – double array
If ${\mathbf{info}}={\mathbf{0}}$ or $\mathbf{n}+{\mathbf{1}}$, an estimate of the component-wise relative backward error of each computed solution vector ${\stackrel{^}{x}}_{j}$ (i.e., the smallest relative change in any element of $A$ or $B$ that makes ${\stackrel{^}{x}}_{j}$ an exact solution).
12:   $\mathrm{work}\left(\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,4×{\mathbf{n}}\right)\right)$ – double array
${\mathbf{work}}\left(1\right)$ contains the reciprocal pivot growth factor $‖A‖/‖U‖$. 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$.
13:   $\mathrm{info}$int64int32nag_int scalar
${\mathbf{info}}=0$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{info}}<0$
If ${\mathbf{info}}=-i$, argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.
W  ${\mathbf{info}}>0 \text{and} {\mathbf{info}}\le {\mathbf{n}}$
Element $_$ of the diagonal 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.
W  ${\mathbf{info}}={\mathbf{n}}+1$
$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.

## Accuracy

For each right-hand side vector $b$, the computed solution $\stackrel{^}{x}$ is the exact solution of a perturbed system of equations $\left(A+E\right)\stackrel{^}{x}=b$, where
 $E≤cnεPLU ,$
$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 $\stackrel{^}{x}$ satisfies a forward error bound of the form
 $x-x^∞ x^∞ ≤ wc condA,x^,b$
where $\mathrm{cond}\left(A,\stackrel{^}{x},b\right)={‖\left|{A}^{-1}\right|\left(\left|A\right|\left|\stackrel{^}{x}\right|+\left|b\right|\right)‖}_{\infty }/{‖\stackrel{^}{x}‖}_{\infty }\le \mathrm{cond}\left(A\right)={‖\left|{A}^{-1}\right|\left|A\right|‖}_{\infty }\le {\kappa }_{\infty }\left(A\right)$. If $\stackrel{^}{x}$ is the $j$th column of $X$, then ${w}_{c}$ is returned in ${\mathbf{berr}}\left(j\right)$ and a bound on ${‖x-\stackrel{^}{x}‖}_{\infty }/{‖\stackrel{^}{x}‖}_{\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 function is nag_lapack_zgesvx (f07ap).

## Example

This example solves the equations
 $AX=B ,$
where $A$ is the general matrix
 $A= 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$
and
 $B = 9.52 18.47 2435.00 225.00 0.77 -13.28 -6.22 -6.21 .$
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.
```function f07ab_example

fprintf('f07ab example results\n\n');

% Linear Problem
n    = 4;
nrhs = 2;
a = [  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];
b = [   9.52,  18.47;
2435.00, 225.00;
0.77, -13.28;
-6.22,  -6.21];

% Input parameter initialization
fact  = 'Equilibration';
trans = 'No transpose';
equed = 'N';
af   = zeros(n, n);
ipiv = zeros(n, 1, 'int64');
r    = zeros(n, 1);
c    = zeros(n, 1);

% Solve
[a, af, ipiv, equed, r, c, b, x, rcond, ferr, berr, work, info] = ...
f07ab(...
fact, trans, a, af, ipiv, equed, r, c, b);

fprintf('Solution is x:\n');
disp(x);
fprintf('\nApproximate condition number = %8.3f\n',1/rcond);
fprintf('Approximate forward  errors  :\n');
fprintf('                               %11.3e\n',ferr);
fprintf('Approximate backward errors  :\n')
fprintf('                               %11.3e\n',berr);

```
```f07ab example results

Solution is x:
1.0000    3.0000
-1.0000    2.0000
3.0000    4.0000
-5.0000    1.0000

Approximate condition number =   54.967
Approximate forward  errors  :
2.384e-14
3.301e-14
Approximate backward errors  :
6.800e-17
8.040e-17
```