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)

Purpose

nag_matop_real_symm_posdef_fac (f01bu) performs a ULDLTUTULDLTUT 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 AA, of order nn and bandwidth 2m + 12m+1, is divided into the leading principal sub-matrix of order kk and its complement, where mknmkn. A UDUTUDUT decomposition of the latter and an LDLTLDLT decomposition of the former are obtained by means of a sequence of elementary transformations, where UU is unit upper triangular, LL is unit lower triangular and DD is diagonal. Thus if k = nk=n, an LDLTLDLT decomposition of AA 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 = λBxAx=λBx by the method of Crawford where AA and BB are of band form. In this context, kk is chosen to be close to n / 2n/2 and the decomposition is applied to the matrix BB.

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
kk, the change-over point in the decomposition.
Constraint: m11knm1-1kn.
2:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldam1ldam1.
The upper triangle of the nn by nn symmetric band matrix AA, with the diagonal of the matrix stored in the (m + 1)(m+1)th row of the array, and the mm superdiagonals within the band stored in the first mm rows of the array. Each column of the matrix is stored in the corresponding column of the array. For example, if n = 6n=6 and m = 2m=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.
nn, the order of the matrix AA.
2:     m1 – int64int32nag_int scalar
Default: The first dimension of the array a.
m + 1m+1, where mm is the number of nonzero superdiagonals in AA. Normally m1nm1n.

Input Parameters Omitted from the MATLAB Interface

lda w

Output Parameters

1:     a(lda,n) – double array
ldam1ldam1.
AA stores the corresponding elements of LL, DD and UU.
2:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,k < m1k<m1 or k > nk>n.
  ifail = 2ifail=2
  ifail = 3ifail=3
The matrix AA is not positive definite, perhaps as a result of rounding errors, giving an element of DD which is zero or negative. ifail = 3ifail=3 when the failure occurs in the leading principal sub-matrix of order k and ifail = 2ifail=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 UU, LL and DD satisfy the relation ULDLTUT = A + EULDLTUT=A+E where the 22-norms of AA and EE are related by Ec(m + 1)2εAEc (m+1) 2εA where cc 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 + 3nmnm2+3nm.
This function is specifically designed for use as the first stage in the solution of the generalized symmetric eigenproblem Ax = λBxAx=λBx by Crawford's method which preserves band form in the transformation to a similar standard problem. In this context, for maximum efficiency, kk should be chosen as the multiple of mm nearest to n / 2n/2.
The matrix UU is such that U1AUTU-1AU-T is diagonal in its last nkn-k rows and columns, LL is such that L1U1AUTLT = DL-1U-1AU-TL-T=D and DD is diagonal. To find UU, LL and DD where A = ULDLTUTA=ULDLTUT requires nm(m + 3) / 2m(m + 1)(m + 2) / 3nm(m+3)/2-m(m+1)(m+2)/3 multiplications and divisions which, is independent of kk.

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



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–2013