This function is specifically designed to precede nag_matop_real_symm_posdef_geneig (f01bv) for the transformation of the symmetric-definite eigenproblem

A x = λ B x

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: $k$ – int64int32nag_int scalar

k

, the change-over point in the decomposition.

Constraint:

m1 - 1 \leq k \leq n

2: $a (lda, n)$ – double array

lda, the first dimension of the array, must satisfy the constraint

lda \geq m1

The upper triangle of the

n

n

symmetric band matrix

A

, with the diagonal of the matrix stored in the

(m + 1)

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

\begin{array}{l} * & * & a_{13} & a_{24} & a_{35} & a_{46} \\ * & a_{12} & a_{23} & a_{34} & a_{45} & a_{56} \\ a_{11} & a_{22} & a_{33} & a_{44} & a_{55} & a_{66} \end{array}

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$ , the order of the matrix $A$ .
2: $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 $m1 ≪ n$ .

Output Parameters

1: $a (lda, n)$ – double array: $A$ stores the corresponding elements of $L$ , $D$ and $U$ .
2: $ifail$ – int64int32nag_int scalar: $ifail = 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$

On entry,

k < m1 - 1

k > n

$ifail = 2$
$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. $ifail = 3$ when the failure occurs in the leading principal sub-matrix of order k and $ifail = 2$ when it occurs in the complement.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$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

U L D L^{T} U^{T} = A + E

where the

2

-norms of

A

and

E

are related by

‖E‖ \leq c {(m + 1)}^{2} ε ‖A‖

where

c

is a constant of order unity and

ε

is the machine precision. In practice, the error is usually appreciably smaller than this.

Further Comments

The time taken by nag_matop_real_symm_posdef_fac (f01bu) is approximately proportional to

n m^{2} + 3 n m

This function is specifically designed for use as the first stage in the solution of the generalized symmetric eigenproblem

A x = λ B x

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^{- T}

is diagonal in its last

n - k

rows and columns,

L

is such that

L^{- 1} U^{- 1} A U^{- T} L^{- T} = D

and

D

is diagonal. To find

U

L

and

D

where

A = U L D L^{T} U^{T}

requires

n m (m + 3) / 2 - m (m + 1) (m + 2) / 3

multiplications and divisions which, is independent of

k

Example

This example finds a

U L D L^{T} U^{T}

decomposition of the real symmetric positive definite matrix

(\begin{array}{r} 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 \end{array}) .

Open in the MATLAB editor: f01bu_example

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

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox