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_1rhs (f04at)

Purpose

nag_linsys_real_square_solve_1rhs (f04at) calculates the accurate solution of a set of real linear equations with a single right-hand side, using an LU$LU$ factorization with partial pivoting, and iterative refinement.

Syntax

[c, aa, ifail] = f04at(a, b, 'n', n)
[c, aa, ifail] = nag_linsys_real_square_solve_1rhs(a, b, 'n', n)

Description

Given a set of real linear equations, Ax = b$Ax=b$, the function first computes an LU$LU$ factorization of A$A$ with partial pivoting, PA = LU$PA=LU$, where P$P$ is a permutation matrix, L$L$ is lower triangular and U$U$ is unit upper triangular. An approximation to x$x$ is found by forward and backward substitution in Ly = Pb$Ly=Pb$ and Ux = y$Ux=y$. The residual vector r = bAx$r=b-Ax$ is then calculated using additional precision, and a correction d$d$ to x$x$ is found by solving LUd = r$LUd=r$. x$x$ is replaced by x + d$x+d$, and this iterative refinement of the solution is repeated until full machine accuracy is obtained.

References

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

Parameters

Compulsory Input Parameters

1:     a(lda, : $:$) – double array
The first dimension of the array a must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The n$n$ by n$n$ matrix A$A$.
2:     b( : $:$) – double array
Note: the dimension of the array b must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The right-hand side vector b$b$.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the array a The second dimension of the arrays a, b.
n$n$, the order of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.

Input Parameters Omitted from the MATLAB Interface

lda ldaa wks1 wks2

Output Parameters

1:     c(n) – double array
The solution vector x$x$.
2:     aa(ldaa,n${\mathbf{n}}$) – double array
The first dimension of the array aa will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
ldaamax (1,n)$\mathit{ldaa}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The triangular factors L$L$ and U$U$, except that the unit diagonal elements of U$U$ are not stored.
3:     ifail – int64int32nag_int scalar
${\mathrm{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:
ifail = 1${\mathbf{ifail}}=1$
The matrix A$A$ is singular, possibly due to rounding errors.
ifail = 2${\mathbf{ifail}}=2$
Iterative refinement fails to improve the solution, i.e., the matrix A$A$ is too ill-conditioned.
ifail = 3${\mathbf{ifail}}=3$
 On entry, n < 0${\mathbf{n}}<0$, or lda < max (1,n)$\mathit{lda}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or ldaa < max (1,n)$\mathit{ldaa}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.

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_1rhs (f04at) is approximately proportional to n3${n}^{3}$.
The function must not be called with the same name for parameters b and c.

Example

```function nag_linsys_real_square_solve_1rhs_example
a = [33, 16, 72;
-24, -10, -57;
-8, -4, -17];
b = [-359;
281;
85];
[c, aa, ifail] = nag_linsys_real_square_solve_1rhs(a, b)
```
```

c =

1
-2
-5

aa =

-8.0000    0.5000    2.1250
-24.0000    2.0000   -3.0000
33.0000   -0.5000    0.3750

ifail =

0

```
```function f04at_example
a = [33, 16, 72;
-24, -10, -57;
-8, -4, -17];
b = [-359;
281;
85];
[c, aa, ifail] = f04at(a, b)
```
```

c =

1
-2
-5

aa =

-8.0000    0.5000    2.1250
-24.0000    2.0000   -3.0000
33.0000   -0.5000    0.3750

ifail =

0

```