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$, where $A$ is symmetric positive definite, nag_linsys_real_posdef_solve_ref (f04ab) 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 matrix $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 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 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(\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{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 $A$ is unchanged.
2:     $\mathrm{c}\left(\mathit{ldc},{\mathbf{m}}\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$.
3:     $\mathrm{bb}\left(\mathit{ldbb},{\mathbf{m}}\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 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 ${\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{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 39 of Wilkinson and Reinsch (1971).

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

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

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

% Solve 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, resid, ifail] = f04ab(a, b);

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

f04ab example results

Solution
1
1
1
1