MCHOL

Computes an upper triangular factorization of a real symmetric matrix A plus a diagonal matrix D, where D is determined sequentially during the Cholesky factorization in order to make A + D nonnegative definite.

Required Arguments

A ‑ N by N symmetric matrix for which a Cholesky factorization is attempted. (Input)
Only elements in the upper triangle and diagonal of A are referenced.

IRANK ‑ Rank of A + D. (Output)

R ‑ N by N upper triangular matrix containing the R matrix from a Cholesky decomposition RTR of A + D. (Output)
The lower triangle of R is not referenced. If A is not needed, then R and A can share the same storage locations.

DMAX ‑ Largest diagonal element of D. (Output)
If DMAX equals 0.0, then A is nonnegative definite, and R is a Cholesky factorization of A. If DMAX is positive, then A is indefinite, i.e., A has at least one negative eigenvalue. In this case, DMAX is an upper bound on the absolute value of the minimum eigenvalue.

IND ‑ Index for subsequent computation of a descent direction in the case of a saddle point. (Output)
If IND = 0, then A is nonnegative definite. For positive IND, let e be a unit vector with IND‑th element 1 and remaining elements 0. The solution s of Rs = e is a direction of negative curvature, i.e., sT As is negative.

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.e‑5 for single precision and 2.d –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 MCHOL (A, IRANK, R, DMAX, IND [, …])

Specific: The specific interface names are S_MCHOL and D_MCHOL.

FORTRAN 77 Interface

Single: CALL MCHOL (N, A, LDA, TOL, IRANK, R, LDR, DMAX, IND)

Double: The double precision name is DMCHOL.

Description

Routine MCHOL computes a Cholesky factorization, RTR, of A + D where A is symmetric and D is a diagonal matrix with sufficiently large diagonal elements such that A + D is nonnegative definite. The routine is similar to one described by Gill, Murray, and Wright (1981, pages 108‑111). Here, though, we allow A + D to be singular.

The algorithm proceeds sequentially by rows. If A + D is singular, the Cholesky factor R is taken to have some rows that are entirely zero. The i‑th row of A + D is declared to be linearly dependent on the first i ‑ 1 rows if the following two conditions are satisfied:

where ɛ is the input argument TOL.

The routine MCHOL is often used to find a descent direction in a minimization problem. Let A and g be the current Hessian and gradient, respectively, associated with the minimization problem. The solution s of As = ‑g may not give a descent direction if A is not nonnegative definite. Instead, in order to guarantee a descent direction, a solution s of (A + D)s = ‑g can be found where A + D is nonnegative definite. Routine MCHOL is useful for computing the upper triangular Cholesky factor R of A + D so that routine GIRTS can be invoked to compute the descent direction s by solving successively the two triangular linear systems RTx = ‑g and Rs = x for x and then s. Also if g = 0 and A is not nonnegative definite, i.e., the current solution is a saddle point, GIRTS can be used to compute a descent direction s from the linear system Rs = e where e is a unit vector with

Examples

Example 1

A Cholesky factorization of a 5 × 5 symmetric nonnegative definite matrix is computed. Maindonald (1984, pages 85‑86) discusses the example.

! SPECIFICATIONS FOR PARAMETERS

USE MCHOL_INT

USE UMACH_INT

USE WRRRN_INT

IMPLICIT NONE

INTEGER LDA, LDR, N, J

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

INTEGER IND, IRANK, NOUT

REAL A(LDA,N), DMAX, 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/

TOL = 0.00001

CALL MCHOL (A, IRANK, R, DMAX, IND, TOL=TOL)

CALL UMACH (2, NOUT)

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

WRITE (NOUT,99999) ' DMAX = ', DMAX

WRITE (NOUT,99998) ' IND = ', IND

99998 FORMAT (A, I3)

99999 FORMAT (A, 1PE10.3)

CALL WRRRN ('R', R)

END

Output

IRANK = 4

DMAX = 0.000E+00

IND = 0

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

Example 2

A modified Cholesky factorization of a 3 × 3 symmetric indefinite matrix A is computed. A solution of Rs = e3 is also obtained using routine GIRTS. Note that sT As is negative as verified by using routine BLINF (IMSL MATH/LIBRARY). Gill, Murray, and Wright (1981, page 111) discuss the example.

! SPECIFICATIONS FOR PARAMETERS

USE IMSL_LIBRARIES

IMPLICIT NONE

INTEGER LDA, LDR, N, J

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

INTEGER IND, IRANK, NOUT

REAL A(LDA,N), DMAX, E(N,1), R(LDR,N), S(N, 1), SPAS, TOL

DATA (A(1,J),J=1,N)/1, 1, 2/

DATA (A(2,J),J=1,N)/1, 1, 3/

DATA (A(3,J),J=1,N)/2, 3, 1/

TOL = 0.00001

CALL MCHOL (A, IRANK, R, DMAX, IND, TOL=TOL)

CALL UMACH (2, NOUT)

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

WRITE (NOUT,99999) ' DMAX = ', DMAX

WRITE (NOUT,99998) ' IND = ', IND

CALL WRRRN ('R', R)

IF (IND .GT. 0) THEN

E= 0.0

E(IND, 1) = 1.0

CALL GIRTS (R, E, IRANK, S)

SPAS = BLINF(A, S(:,1), S(:,1))

WRITE (NOUT,*) ' '

WRITE (NOUT,99999) ' trans(s)*A*s = ', SPAS

END IF

99998 FORMAT (A, I3)

99999 FORMAT (A, F10.3)

END

Output

IRANK = 3

DMAX = 5.016

IND = 3

1 2 3

1 1.942 0.515 1.030

2 0.000 2.398 1.030

3 0.000 0.000 1.059

trans(s)*A*s = -2.254