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 ULDLTUT$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:
Mark 22: m1 has been made optional
.

## Description

The symmetric positive definite matrix A$A$, of order n$n$ and bandwidth 2m + 1$2m+1$, is divided into the leading principal sub-matrix of order k$k$ and its complement, where mkn$m\le k\le n$. A UDUT$UD{U}^{\mathrm{T}}$ decomposition of the latter and an LDLT$LD{L}^{\mathrm{T}}$ decomposition of the former are obtained by means of a sequence of elementary transformations, where U$U$ is unit upper triangular, L$L$ is unit lower triangular and D$D$ is diagonal. Thus if k = n$k=n$, an LDLT$LD{L}^{\mathrm{T}}$ decomposition of A$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 = λBx$Ax=\lambda Bx$ by the method of Crawford where A$A$ and B$B$ are of band form. In this context, k$k$ is chosen to be close to n / 2$n/2$ and the decomposition is applied to the matrix B$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:     k – int64int32nag_int scalar
k$k$, the change-over point in the decomposition.
Constraint: m11kn${\mathbf{m1}}-1\le {\mathbf{k}}\le {\mathbf{n}}$.
2:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldam1$\mathit{lda}\ge {\mathbf{m1}}$.
The upper triangle of the n$n$ by n$n$ symmetric band matrix A$A$, with the diagonal of the matrix stored in the (m + 1)$\left(m+1\right)$th row of the array, and the m$m$ superdiagonals within the band stored in the first m$m$ rows of the array. Each column of the matrix is stored in the corresponding column of the array. For example, if n = 6$n=6$ and m = 2$m=2$, the storage scheme is
 * * a13 a24 a35 a46 * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66
$* * 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:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
n$n$, the order of the matrix A$A$.
2:     m1 – int64int32nag_int scalar
Default: The first dimension of the array a.
m + 1$m+1$, where m$m$ is the number of nonzero superdiagonals in A$A$. Normally m1n${\mathbf{m1}}\ll {\mathbf{n}}$.

lda w

### Output Parameters

1:     a(lda,n) – double array
ldam1$\mathit{lda}\ge {\mathbf{m1}}$.
A$A$ stores the corresponding elements of L$L$, D$D$ and U$U$.
2:     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$
 On entry, k < m1${\mathbf{k}}<{\mathbf{m1}}$ or k > n${\mathbf{k}}>{\mathbf{n}}$.
ifail = 2${\mathbf{ifail}}=2$
ifail = 3${\mathbf{ifail}}=3$
The matrix A$A$ is not positive definite, perhaps as a result of rounding errors, giving an element of D$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.

## Accuracy

The Cholesky decomposition of a positive definite matrix is known for its remarkable numerical stability (see Wilkinson (1965)). The computed U$U$, L$L$ and D$D$ satisfy the relation ULDLTUT = A + E$ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}=A+E$ where the 2$2$-norms of A$A$ and E$E$ are related by Ec(m + 1)2εA$‖E‖\le c{\left(m+1\right)}^{2}\epsilon ‖A‖$ where c$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 nm2 + 3nm$n{m}^{2}+3nm$.
This function is specifically designed for use as the first stage in the solution of the generalized symmetric eigenproblem Ax = λBx$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$k$ should be chosen as the multiple of m$m$ nearest to n / 2$n/2$.
The matrix U$U$ is such that U1AUT${U}^{-1}A{U}^{-\mathrm{T}}$ is diagonal in its last nk$n-k$ rows and columns, L$L$ is such that L1U1AUTLT = D${L}^{-1}{U}^{-1}A{U}^{-\mathrm{T}}{L}^{-\mathrm{T}}=D$ and D$D$ is diagonal. To find U$U$, L$L$ and D$D$ where A = ULDLTUT$A=ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}$ requires nm(m + 3) / 2m(m + 1)(m + 2) / 3$nm\left(m+3\right)/2-m\left(m+1\right)\left(m+2\right)/3$ multiplications and divisions which, is independent of k$k$.

## Example

```function nag_matop_real_symm_posdef_fac_example
k = int64(4);
a = [0, 0, 6, -4, 15, 4, -18;
0, -9, -2, -66, -24, -74, 24;
3, 31, 123, 145, 61, 98, 6];
[aOut, ifail] = nag_matop_real_symm_posdef_fac(k, a)
```
```

aOut =

0     0     2    -1     3     2    -3
0    -3     4     5    -4    -1     4
3     4     2     3     5     2     6

ifail =

0

```
```function f01bu_example
k = int64(4);
a = [0, 0, 6, -4, 15, 4, -18;
0, -9, -2, -66, -24, -74, 24;
3, 31, 123, 145, 61, 98, 6];
[aOut, ifail] = f01bu(k, a)
```
```

aOut =

0     0     2    -1     3     2    -3
0    -3     4     5    -4    -1     4
3     4     2     3     5     2     6

ifail =

0

```