NAG Library Routine Document
F01MCF
1 Purpose
F01MCF computes the Cholesky factorization of a real symmetric positive definite variable-bandwidth matrix.
2 Specification
INTEGER |
N, LAL, NROW(N), IFAIL |
REAL (KIND=nag_wp) |
A(LAL), AL(LAL), D(N) |
|
3 Description
F01MCF determines the unit lower triangular matrix and the diagonal matrix in the Cholesky factorization of a symmetric positive definite variable-bandwidth matrix of order . (Such a matrix is sometimes called a ‘sky-line’ matrix.)
The matrix
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 9 for an example). The
width
of the
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
(also see
Section 8).
The method is based on the property that during Cholesky factorization there is no fill-in outside the envelope.
The determination of
and
is normally the first of two steps in the solution of the system of equations
. The remaining step, viz. the solution of
, 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 Parameters
- 1: N – INTEGERInput
On entry: , the order of the matrix .
Constraint:
.
- 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
, 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
- 3: LAL – INTEGERInput
On entry: the dimension of the arrays
A and
AL as declared in the (sub)program from which F01MCF is called.
Constraint:
.
- 4: NROW(N) – INTEGER arrayInput
On entry: must contain the width of row of the matrix , i.e., the number of elements between the first (leftmost) nonzero element and the element on the diagonal, inclusive.
Constraint:
, for .
- 5: AL(LAL) – REAL (KIND=nag_wp) arrayOutput
On exit: the elements within the envelope of the lower triangular matrix
, taken in row by row order. The envelope of
is identical to that of the lower triangle of
. The unit diagonal elements of
are stored explicitly. See also
Section 8.
- 6: D(N) – REAL (KIND=nag_wp) arrayOutput
On exit: the diagonal elements of the diagonal matrix . Note that the determinant of is equal to the product of these diagonal elements. If the value of the determinant is required it should not be determined by forming the product explicitly, because of the possibility of overflow or underflow. The logarithm of the determinant may safely be formed from the sum of the logarithms of the diagonal elements.
- 7: IFAIL – INTEGERInput/Output
-
On entry:
IFAIL must be set to
,
. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
.
When the value is used it is essential to test the value of IFAIL on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
On entry, | , |
or | for some , or , |
or | . |
is not positive definite, or this property has been destroyed by rounding errors. The factorization has not been completed.
is not positive definite, or this property has been destroyed by rounding errors. The factorization has not been completed.
7 Accuracy
If
on exit, then the
computed
and
satisfy the relation
, where
and
where
is a constant of order unity,
is the largest value of
, and
is the
machine precision. See pages 25–27 and 54–55 of
Wilkinson and Reinsch (1971).
The time taken by F01MCF is approximately proportional to the sum of squares of the values of .
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 parameters
A and
AL, in which case
overwrites the lower triangle of
A. However this is not standard Fortran and may not work in all implementations.
9 Example
This example obtains the Cholesky factorization of the symmetric matrix, whose lower triangle is:
For this matrix, the elements of
NROW must be set to
,
,
,
,
,
, and the elements within the envelope must be supplied in row order as:
9.1 Program Text
Program Text (f01mcfe.f90)
9.2 Program Data
Program Data (f01mcfe.d)
9.3 Program Results
Program Results (f01mcfe.r)