Integer, Intent (In)	::	n, lal
Integer, Intent (Inout)	::	nrow(n), ifail
Real (Kind=nag_wp), Intent (In)	::	a(lal)
Real (Kind=nag_wp), Intent (Out)	::	al(lal), d(n)

C Header Interface

#include nagmk26.h

void	f01mcf_ ( const Integer n, const double a[], const Integer lal, Integer nrow[], double al[], double d[], Integer *ifail)

3

Description

f01mcf determines the unit lower triangular matrix

L

and the diagonal matrix

D

in the Cholesky factorization

A = L D L^{T}

of a symmetric positive definite variable-bandwidth matrix

A

of order

n

. (Such a matrix is sometimes called a ‘sky-line’ matrix.)

The matrix

A

is represented by the elements lying within the envelope of its lower triangular part, that is, between the first nonzero of each row and the diagonal (see Section 10 for an example). The width

nrow (i)

of the

i

th row is the number of elements between the first nonzero element and the element on the diagonal, inclusive. Although, of course, any matrix possesses an envelope as defined, this routine is primarily intended for the factorization of symmetric positive definite matrices with an average bandwidth which is small compared with

n

(also see Section 9).

The method is based on the property that during Cholesky factorization there is no fill-in outside the envelope.

The determination of

L

and

D

is normally the first of two steps in the solution of the system of equations

A x = b

. The remaining step, viz. the solution of

L D L^{T} x = b

, may be carried out using f04mcf.

4

References

Jennings A (1966) A compact storage scheme for the solution of symmetric linear simultaneous equations Comput. J. 9 281–285

Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag

5

Arguments

1: $n$ – IntegerInput

On entry:

n

, the order of the matrix

A

Constraint:

n \geq 1

2: $a (lal)$ – Real (Kind=nag_wp) arrayInput

On entry: the elements within the envelope of the lower triangle of the positive definite symmetric matrix

A

, taken in row by row order. The following code assigns the matrix elements within the envelope to the correct elements of the array:

   k = 0
   Do 20 i = 1, n
      Do 10 j = i-nrow(i)+1, i
         k = k + 1
         a(k) = matrix (i,j)
10    Continue
20 Continue

6

Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$

On entry,	$n < 1$ ,
or	for some $i$ , $nrow (i) < 1$ or $nrow (i) > i$ ,
or	$lal < nrow (1) + nrow (2) + \dots + nrow (n)$ .

$ifail = 2$: $A$ is not positive definite, or this property has been destroyed by rounding errors. The factorization has not been completed.

$ifail = 3$: $A$ is not positive definite, or this property has been destroyed by rounding errors. The factorization has not been completed.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7

Accuracy

ifail = 0

on exit, then the computed

L

and

D

satisfy the relation

L D L^{T} = A + F

, where

{‖F‖}_{2} \leq k m^{2} ε \times \max_{i} a_{i i}

and

{‖F‖}_{2} \leq k m^{2} ε \times {‖A‖}_{2},

where

k

is a constant of order unity,

m

is the largest value of

nrow (i)

, and

ε

is the machine precision. See pages 25–27 and 54–55 of Wilkinson and Reinsch (1971).

8

Parallelism and Performance

f01mcf is not threaded in any implementation.

9

Further Comments

The time taken by f01mcf is approximately proportional to the sum of squares of the values of

nrow (i)

The distribution of row widths may be very non-uniform without undue loss of efficiency. Moreover, the routine has been designed to be as competitive as possible in speed with routines designed for full or uniformly banded matrices, when applied to such matrices.

Unless otherwise stated in the Users' Note for your implementation, the routine may be called with the same actual array supplied for arguments a and al, in which case

L

overwrites the lower triangle of a. However this is not standard Fortran and may not work in all implementations.

10

Example

This example obtains the Cholesky factorization of the symmetric matrix, whose lower triangle is:

(\begin{array}{r} 1 \\ 2 & 5 \\ 0 & 3 & 13 \\ 0 & 0 & 0 & 16 \\ 5 & 14 & 18 & 8 & 55 \\ 0 & 0 & 0 & 24 & 17 & 77 \end{array}) .

For this matrix, the elements of nrow must be set to

1

2

2

1

5

3

, and the elements within the envelope must be supplied in row order as:

1, 2, 5, 3, 13, 16, 5, 14, 18, 8, 55, 24, 17, 77 .

NAG Library Routine Document

f01mcf (real_vband_posdef_fac)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification