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
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
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
R
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
