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_posdef_solve_1rhs (f04as)

## Purpose

nag_linsys_real_posdef_solve_1rhs (f04as) calculates the accurate solution of a set of real symmetric positive definite linear equations with a single right-hand side, $Ax=b$, using a Cholesky factorization and iterative refinement.

## Syntax

[a, c, ifail] = f04as(a, b, 'n', n)
[a, c, ifail] = nag_linsys_real_posdef_solve_1rhs(a, b, 'n', n)

## Description

Given a set of real linear equations $Ax=b$, where $A$ is a symmetric positive definite matrix, nag_linsys_real_posdef_solve_1rhs (f04as) first computes a Cholesky factorization of $A$ as $A=L{L}^{\mathrm{T}}$ where $L$ is lower triangular. An approximation to $x$ is found by forward and backward substitution. The residual vector $r=b-Ax$ is then calculated using additional precision and a correction $d$ to $x$ is found by solving $L{L}^{\mathrm{T}}d=r$. $x$ is then replaced by $x+d$, and this iterative refinement of the solution is repeated until 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:     $\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 upper triangle of the $n$ by $n$ positive definite symmetric matrix $A$. The elements of the array below the diagonal need not be set.
2:     $\mathrm{b}\left(\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)\right)$ – double array
The dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The right-hand side vector $b$.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the array a and the second dimension of the arrays a, b.
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\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)$.
The elements of the array below the diagonal are overwritten; the upper triangle of ${\mathbf{a}}$ is unchanged.
2:     $\mathrm{c}\left(\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)\right)$ – double array
The solution vector $x$.
3:     $\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 not positive definite, 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 $\mathit{lda}<\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 39 of Wilkinson and Reinsch (1971).

The time taken by nag_linsys_real_posdef_solve_1rhs (f04as) is approximately proportional to ${n}^{3}$.
The function must not be called with the same name for arguments b and c.

## Example

This example solves the set of linear equations $Ax=b$ where
 $A= 5 7 6 5 7 10 8 7 6 8 10 9 5 7 9 10 and b= 23 32 33 31 .$
```function f04as_example

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

% Accurate solution to Ax = b, for positive definite A
a = [5,  7,  6,  5;
7, 10,  8,  7;
6,  8, 10,  9;
5,  7,  9, 10];
b = [23;
32;
33;
31];

[afac, x, ifail] = f04as(a, b);

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

```
```f04as example results

Solution
1
1
1
1

```