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_matop_real_symm_posdef_inv (f01ab)

## Purpose

nag_matop_real_symm_posdef_inv (f01ab) calculates the accurate inverse of a real symmetric positive definite matrix, using a Cholesky factorization and iterative refinement.

## Syntax

[a, b, ifail] = f01ab(a, 'n', n)
[a, b, ifail] = nag_matop_real_symm_posdef_inv(a, 'n', n)

## Description

To compute the inverse $X$ of a real symmetric positive definite matrix $A$, nag_matop_real_symm_posdef_inv (f01ab) 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 computing ${L}^{-1}$ and then the product ${L}^{-\mathrm{T}}{L}^{-1}$. The residual matrix $R=I-AX$ is 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 inverse 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},{\mathbf{n}}\right)$ – double array
lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge {\mathbf{n}}+1$.
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.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the second dimension of the array a.
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 1$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},{\mathbf{n}}\right)$ – double array
The lower triangle of the inverse matrix $X$ is stored in the elements of the array below the diagonal, in rows $2$ to $n+1$; ${x}_{ij}$ is stored in ${\mathbf{a}}\left(i+1,j\right)$ for $i\ge j$. The upper triangle of the original matrix is unchanged.
2:     $\mathrm{b}\left(\mathit{ldb},{\mathbf{n}}\right)$ – double array
The lower triangle of the inverse matrix $X$, with ${x}_{ij}$ stored in ${\mathbf{b}}\left(i,j\right)$, for $i\ge j$.
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$
The refinement process fails to converge, i.e., the matrix $A$ is ill-conditioned.
${\mathbf{ifail}}=3$
${\mathbf{n}}<1$, or $\mathit{lda}<{\mathbf{n}}+1$, or $\mathit{ldb}<{\mathbf{n}}$.
${\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 inverse should be correct to full machine accuracy. For a detailed error analysis see page 40 of Wilkinson and Reinsch (1971).

The time taken by nag_matop_real_symm_posdef_inv (f01ab) is approximately proportional to ${n}^{3}$.

## Example

This example finds the inverse of the $4$ by $4$ matrix:
 $5 7 6 5 7 10 8 7 6 8 10 9 5 7 9 10 .$
```function f01ab_example

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

a = [  5,  7,  6,  5;
7, 10,  8,  7;
6,  8, 10,  9;
5,  7,  9, 10];

a = [a; 0 0 0 0];

[X, L, ifail] = f01ab(a);

matrix = 'Lower';
diag   = 'Non-unit';
xtitl  = 'Lower triangle of inverse:';
[ifail] = x04ca( ...
matrix, diag, L, xtitl);

```
```f01ab example results

Lower triangle of inverse:
1          2          3          4
1     68.0000
2    -41.0000    25.0000
3    -17.0000    10.0000     5.0000
4     10.0000    -6.0000    -3.0000     2.0000
```