hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_matop_real_symm_posdef_fac (f01bu)


    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example


nag_matop_real_symm_posdef_fac (f01bu) performs a ULDLTUT decomposition of a real symmetric positive definite band matrix.


[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


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 mkn. A UDUT decomposition of the latter and an LDLT 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 LDLT 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=λ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.


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


Compulsory Input Parameters

1:     k int64int32nag_int scalar
k, the change-over point in the decomposition.
Constraint: m1-1kn.
2:     aldan – double array
lda, the first dimension of the array, must satisfy the constraint ldam1.
The upper triangle of the n by n symmetric band matrix A, with the diagonal of the matrix stored in the m+1th 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);

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 m1n.

Output Parameters

1:     aldan – 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:
On entry,k<m1-1 or k>n.
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.
An unexpected error has been triggered by this routine. Please contact NAG.
Your licence key may have expired or may not have been installed correctly.
Dynamic memory allocation failed.


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 ULDLTUT=A+E where the 2-norms of A and E are related by Ec 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 nm2+3nm.
This function is specifically designed for use as the first stage in the solution of the generalized symmetric eigenproblem Ax=λ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-1AU-T is diagonal in its last n-k rows and columns, L is such that L-1U-1AU-TL-T=D and D is diagonal. To find U, L and D where A=ULDLTUT requires nmm+3/2-mm+1m+2/3 multiplications and divisions which, is independent of k.


This example finds a ULDLTUT 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

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015