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_fac (f01bu)

## Purpose

nag_matop_real_symm_posdef_fac (f01bu) performs a $ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}$ decomposition of a real symmetric positive definite band matrix.

## Syntax

[a, ifail] = f01bu(k, a, 'n', n, 'm1', m1)
[a, ifail] = nag_matop_real_symm_posdef_fac(k, a, 'n', n, 'm1', m1)
Note: the interface to this routine has changed since earlier releases of the toolbox:
 At Mark 22: m1 was made optional

## Description

The symmetric positive definite matrix $A$, of order $n$ and bandwidth $2m+1$, is divided into the leading principal sub-matrix of order $k$ and its complement, where $m\le k\le n$. A $UD{U}^{\mathrm{T}}$ decomposition of the latter and an $LD{L}^{\mathrm{T}}$ decomposition of the former are obtained by means of a sequence of elementary transformations, where $U$ is unit upper triangular, $L$ is unit lower triangular and $D$ is diagonal. Thus if $k=n$, an $LD{L}^{\mathrm{T}}$ decomposition of $A$ is obtained.
This function is specifically designed to precede nag_matop_real_symm_posdef_geneig (f01bv) for the transformation of the symmetric-definite eigenproblem $Ax=\lambda Bx$ by the method of Crawford where $A$ and $B$ are of band form. In this context, $k$ is chosen to be close to $n/2$ and the decomposition is applied to the matrix $B$.

## References

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford
Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{k}$int64int32nag_int scalar
$k$, the change-over point in the decomposition.
Constraint: ${\mathbf{m1}}-1\le {\mathbf{k}}\le {\mathbf{n}}$.
2:     $\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{m1}}$.
The upper triangle of the $n$ by $n$ symmetric band matrix $A$, with the diagonal of the matrix stored in the $\left(m+1\right)$th row of the array, and the $m$ superdiagonals within the band stored in the first $m$ rows of the array. Each column of the matrix is stored in the corresponding column of the array. For example, if $n=6$ and $m=2$, the storage scheme is
 $* * a13 a24 a35 a46 * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66$
Elements in the top left corner of the array are not used. The following code assigns the matrix elements within the band to the correct elements of the array:
``` for j=1:n
for i=max(1,j-m1+1):j
a(i-j+m1,j) = matrix (i,j);
end
end ```

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the second dimension of the array a.
$n$, the order of the matrix $A$.
2:     $\mathrm{m1}$int64int32nag_int scalar
Default: the first dimension of the array a.
$m+1$, where $m$ is the number of nonzero superdiagonals in $A$. Normally ${\mathbf{m1}}\ll {\mathbf{n}}$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},{\mathbf{n}}\right)$ – double array
$A$ stores the corresponding elements of $L$, $D$ and $U$.
2:     $\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$
 On entry, ${\mathbf{k}}<{\mathbf{m1}}-1$ or ${\mathbf{k}}>{\mathbf{n}}$.
${\mathbf{ifail}}=2$
${\mathbf{ifail}}=3$
The matrix $A$ is not positive definite, perhaps as a result of rounding errors, giving an element of $D$ which is zero or negative. ${\mathbf{ifail}}={\mathbf{3}}$ when the failure occurs in the leading principal sub-matrix of order k and ${\mathbf{ifail}}={\mathbf{2}}$ when it occurs in the complement.
${\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 Cholesky decomposition of a positive definite matrix is known for its remarkable numerical stability (see Wilkinson (1965)). The computed $U$, $L$ and $D$ satisfy the relation $ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}=A+E$ where the $2$-norms of $A$ and $E$ are related by $‖E‖\le c{\left(m+1\right)}^{2}\epsilon ‖A‖$ where $c$ is a constant of order unity and $\epsilon$ is the machine precision. In practice, the error is usually appreciably smaller than this.

The time taken by nag_matop_real_symm_posdef_fac (f01bu) is approximately proportional to $n{m}^{2}+3nm$.
This function is specifically designed for use as the first stage in the solution of the generalized symmetric eigenproblem $Ax=\lambda Bx$ by Crawford's method which preserves band form in the transformation to a similar standard problem. In this context, for maximum efficiency, $k$ should be chosen as the multiple of $m$ nearest to $n/2$.
The matrix $U$ is such that ${U}^{-1}A{U}^{-\mathrm{T}}$ is diagonal in its last $n-k$ rows and columns, $L$ is such that ${L}^{-1}{U}^{-1}A{U}^{-\mathrm{T}}{L}^{-\mathrm{T}}=D$ and $D$ is diagonal. To find $U$, $L$ and $D$ where $A=ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}$ requires $nm\left(m+3\right)/2-m\left(m+1\right)\left(m+2\right)/3$ multiplications and divisions which, is independent of $k$.

## Example

This example finds a $ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}$ decomposition of the real symmetric positive definite matrix
 $3 -9 6 -9 31 -2 -4 6 -2 123 -66 15 -4 -66 145 -24 4 15 -24 61 -74 -18 4 -74 98 24 -18 24 6 .$
```function f01bu_example

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

% A Banded matrix A in banded storage format
m1 = int64(3); n = int64(7);
a = [0,  0,   6,  -4,  15,   4, -18;
0, -9,  -2, -66, -24, -74,  24;
3, 31, 123, 145,  61,  98,   6];

k = int64(4);
[a, ifail] = f01bu(k, a);

ptitle = 'Factorized form of the matrix';

[ifail] = x04ce( ...
n, n, int64(0), m1-1, a, ptitle);

```
```f01bu example results

Factorized form of the matrix
1          2          3          4          5          6          7
1      3.0000    -3.0000     2.0000
2                 4.0000     4.0000    -1.0000
3                            2.0000     5.0000     3.0000
4                                       3.0000    -4.0000     2.0000
5                                                  5.0000    -1.0000    -3.0000
6                                                             2.0000     4.0000
7                                                                        6.0000
```