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_ref (f04ab)

## Purpose

nag_linsys_real_posdef_solve_ref (f04ab) calculates the accurate solution of a set of real symmetric positive definite linear equations with multiple right-hand sides, using a Cholesky factorization and iterative refinement.

## Syntax

[a, c, bb, ifail] = f04ab(a, b, 'n', n, 'm', m)
[a, c, bb, ifail] = nag_linsys_real_posdef_solve_ref(a, b, 'n', n, 'm', m)

## Description

Given a set of real linear equations AX = B$AX=B$, where A$A$ is symmetric positive definite, nag_linsys_real_posdef_solve_ref (f04ab) first computes a Cholesky factorization of A$A$ as A = LLT$A=L{L}^{\mathrm{T}}$, where L$L$ is lower triangular. An approximation to X$X$ is found by forward and backward substitution. The residual matrix R = BAX$R=B-AX$ is then calculated using additional precision, and a correction D$D$ to X$X$ is found by solving LLTD = R$L{L}^{\mathrm{T}}D=R$. X$X$ is replaced by X + D$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:     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 upper triangle of the n$n$ by n$n$ positive definite symmetric matrix A$A$. The elements of the array below the diagonal need not be set.
2:     b(ldb, : $:$) – double array
The first 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 second dimension of the array must be at least max (1,m)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$
The n$n$ by m$m$ right-hand side matrix B$B$.

### Optional Input Parameters

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

### Input Parameters Omitted from the MATLAB Interface

lda ldb ldc wkspce ldbb

### Output Parameters

1:     a(lda, : $:$) – double array
The first dimension of the array a 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)$
ldamax (1,n)$\mathit{lda}\ge \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 A$A$ is unchanged.
2:     c(ldc,m${\mathbf{m}}$) – double array
The first dimension of the array c 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,m)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$
ldcmax (1,n)$\mathit{ldc}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The n$n$ by m$m$ solution matrix X$X$.
3:     bb(ldbb,m${\mathbf{m}}$) – double array
The first dimension of the array bb 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,m)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$
ldbbmax (1,n)$\mathit{ldbb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The final n$n$ by m$m$ residual matrix R = BAX$R=B-AX$.
4:     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 not positive definite, 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 m < 0${\mathbf{m}}<0$, or lda < max (1,n)$\mathit{lda}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or ldb < max (1,n)$\mathit{ldb}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or ldc < max (1,n)$\mathit{ldc}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$, or ldbb < max (1,n)$\mathit{ldbb}<\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 39 of Wilkinson and Reinsch (1971).

The time taken by nag_linsys_real_posdef_solve_ref (f04ab) is approximately proportional to n3${n}^{3}$.
If there is only one right-hand side, it is simpler to use nag_linsys_real_posdef_solve_1rhs (f04as).

## Example

```function nag_linsys_real_posdef_solve_ref_example
a = [5, 7, 6, 5;
7, 10, 8, 7;
6, 8, 10, 9;
5, 7, 9, 10];
b = [23;
32;
33;
31];
[aOut, c, bb, ifail] = nag_linsys_real_posdef_solve_ref(a, b)
```
```

aOut =

5.0000    7.0000    6.0000    5.0000
3.1305   10.0000    8.0000    7.0000
2.6833   -0.8944   10.0000    9.0000
2.2361         0    2.1213   10.0000

c =

1
1
1
1

bb =

0
0
0
0

ifail =

0

```
```function f04ab_example
a = [5, 7, 6, 5;
7, 10, 8, 7;
6, 8, 10, 9;
5, 7, 9, 10];
b = [23;
32;
33;
31];
[aOut, c, bb, ifail] = f04ab(a, b)
```
```

aOut =

5.0000    7.0000    6.0000    5.0000
3.1305   10.0000    8.0000    7.0000
2.6833   -0.8944   10.0000    9.0000
2.2361         0    2.1213   10.0000

c =

1
1
1
1

bb =

0
0
0
0

ifail =

0

```