CHFAC

Computes an upper triangular factorization of a real symmetric nonnegative definite matrix.

Required Arguments

A  N by N symmetric nonnegative definite matrix for which an upper triangular factorization is desired. (Input)
Only elements in the upper triangle of A are referenced.

IRANK  Rank of A. (Output)
N  IRANK is the number of effective zero pivots.

R  N by N upper triangular matrix containing the R matrix from a Cholesky decomposition RTR of A. (Output)
The elements of the appropriate rows of R are set to 0.0 if linear dependence of the columns of A is declared. (There are N  IRANK rows of R whose elements are set to 0.0.) If A is not needed, then R and A can share the same storage locations.

Optional Arguments

N  Order of the matrix. (Input)
Default: N = size (A,2).

LDA  Leading dimension of A exactly as specified in the dimension statement in the calling program. (Input)
Default: LDA = size (A,1).

TOL  Tolerance used in determining linear dependence. (Input)
TOL = 100 * AMACH(4) is a common choice. See documentation for routine AMACH.
Default: TOL = 1.19e –5 for single precision and 2.2d –14 for double precision.

LDR  Leading dimension of R exactly as specified in the dimension statement in the calling program. (Input)
Default: LDR = size (R,1).

FORTRAN 90 Interface

Generic: CALL CHFAC (A, IRANK, R [])

Specific: The specific interface names are S_CHFAC and D_CHFAC.

FORTRAN 77 Interface

Single: CALL CHFAC (N, A, LDA, TOL, IRANK, R, LDR)

Double: The double precision name is DCHFAC.

Description

Routine CHFAC computes a Cholesky factorization RTR = A of an n × n symmetric nonnegative definite matrix A. The matrix R is taken to be an upper triangular matrix. The diagonal elements of R are taken to be nonnegative. If A is singular and has rank r, n  r rows of R have all their elements taken to be zero.

The algorithm is based on the work of Healy (1968). The algorithm proceeds sequentially by columns. The i‑th column is declared to be linearly dependent on the first i  1 columns if

 

where ɛ (stored in TOL) is the input tolerance. When a linear dependence is declared, all elements in the i‑th row of R are set to zero.

Modifications due to Farebrother and Berry (1974) and Barrett and Healy (1978) for checking for matrices that are not nonnegative definite are also incorporated. Routine CHFAC declares A to not be nonnegative definite and issues an error message with an error code of 1 if either of the
following conditions is satisfied:

 

Healy’s (1968) algorithm and CHFAC permit the matrices A and R to occupy the same storage locations. Barrett and Healy (1978) in their remark neglect this fact. Routine CHFAC uses

 

for aii in the Condition 2 above to remedy this problem.

Comments

1. Informational error

 

Type

Code

Description

3

1

The input matrix is not nonnegative definite within the tolerance defined by TOL.

2. Elements of row i of R are set to 0.0 if a linear dependence is declared. Linear dependence is declared if

 

 

Example

A Cholesky factorization of a 5 × 5 symmetric nonnegative definite matrix is computed. Maindonald (1984, pages 8586) discusses in detail the computations for this problem.

 

! SPECIFICATIONS FOR PARAMETERS

USE IMSL_LIBRARIES

 

IMPLICIT NONE

INTEGER LDA, LDR, N, J

PARAMETER (N=5, LDA=N, LDR=N)

!

INTEGER IRANK, NOUT

REAL A(LDA,N), R(LDR,N), TOL

!

DATA (A(1,J),J=1,N)/36.0, 12.0, 30.0, 6.0, 18.0/

DATA (A(2,J),J=1,N)/12.0, 20.0, 2.0, 10.0, 22.0/

DATA (A(3,J),J=1,N)/30.0, 2.0, 29.0, 1.0, 7.0/

DATA (A(4,J),J=1,N)/ 6.0, 10.0, 1.0, 14.0, 20.0/

DATA (A(5,J),J=1,N)/ 8.0, 22.0, 7.0, 20.0, 40.0/

!

CALL CHFAC (A, IRANK, R)

CALL UMACH (2, NOUT)

WRITE (NOUT,*) 'IRANK = ', IRANK

CALL WRRRN ('R', R)

END

Output

 

IRANK = 4

 

R

1 2 3 4 5

1 6.000 2.000 5.000 1.000 3.000

2 0.000 4.000 -2.000 2.000 4.000

3 0.000 0.000 0.000 0.000 0.000

4 0.000 0.000 0.000 3.000 3.000

5 0.000 0.000 0.000 0.000 2.449