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_linsys_real_square_solve_ref (f04ae)

## Purpose

nag_linsys_real_square_solve_ref (f04ae) calculates the accurate solution of a set of real linear equations with multiple right-hand sides using an $LU$ factorization with partial pivoting, and iterative refinement.

## Syntax

[c, aa, bb, ifail] = f04ae(a, b, 'n', n, 'm', m)
[c, aa, bb, ifail] = nag_linsys_real_square_solve_ref(a, b, 'n', n, 'm', m)

## Description

Given a set of real linear equations $AX=B$, the function first computes an $LU$ factorization of $A$ with partial pivoting, $PA=LU$, where $P$ is a permutation matrix, $L$ is lower triangular and $U$ is unit upper triangular. An approximation to $X$ is found by forward and backward substitution. The residual matrix $R=B-AX$ is then calculated using additional precision, and a correction $D$ to $X$ is found by solving $LUD=PR$. $X$ is replaced by $X+D$ and this iterative refinement of the solution is repeated until full machine accuracy has been obtained.

## References

Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag

## Parameters

### Compulsory Input Parameters

1:     $\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$.
2:     $\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{m}}\right)$.
The $n$ by $m$ right-hand side matrix $B$.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the arrays a, b and the second dimension of the array a.
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
2:     $\mathrm{m}$int64int32nag_int scalar
Default: the second dimension of the array b.
$m$, the number of right-hand sides.
Constraint: ${\mathbf{m}}\ge 0$.

### Output Parameters

1:     $\mathrm{c}\left(\mathit{ldc},:\right)$ – double array
The first dimension of the array c will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array c will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The $n$ by $m$ solution matrix $X$.
2:     $\mathrm{aa}\left(\mathit{ldaa},:\right)$ – double array
The first dimension of the array aa will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array aa will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The triangular factors $L$ and $U$, except that the unit diagonal elements of $U$ are not stored.
3:     $\mathrm{bb}\left(\mathit{ldbb},:\right)$ – double array
The first dimension of the array bb will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The second dimension of the array bb will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The final $n$ by $m$ residual matrix $R=B-AX$.
4:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
The matrix $A$ is singular, possibly due to rounding errors.
${\mathbf{ifail}}=2$
Iterative refinement fails to improve the solution, i.e., the matrix $A$ is too ill-conditioned.
${\mathbf{ifail}}=3$
 On entry, ${\mathbf{n}}<0$, or ${\mathbf{m}}<0$, or $\mathit{lda}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or $\mathit{ldb}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or $\mathit{ldc}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or $\mathit{ldaa}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or $\mathit{ldbb}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

The computed solutions should be correct to full machine accuracy. For a detailed error analysis see page 107 of Wilkinson and Reinsch (1971).

The time taken by nag_linsys_real_square_solve_ref (f04ae) is approximately proportional to ${n}^{3}$.
If there is only one right-hand side, it is simpler to use nag_linsys_real_square_solve_1rhs (f04at).

## Example

This example solves the set of linear equations $AX=B$ where
 $A= 33 16 72 -24 -10 -57 -8 -4 -17 and B= -359 281 85 .$
```function f04ae_example

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

% Solve Ax = b for general matrix A
a = [  33,  16,  72;
-24, -10, -57;
-8,  -4, -17];
b = [-359; 281;  85];

[x, LU, resid, ifail] = f04ae(a, b);

disp('Solution');
disp(x);

```
```f04ae example results

Solution
1
-2
-5

```