      SUBROUTINE CHLHSD(NR,N,A,EPSM,SX,UDIAG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C PURPOSE
C -------
C FIND THE L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION OF THE PERTURBED
C MODEL HESSIAN MATRIX A+MU*I(WHERE MU\0 AND I IS THE IDENTITY MATRIX)
C WHICH IS SAFELY POSITIVE DEFINITE.  IF A IS SAFELY POSITIVE DEFINITE
C UPON ENTRY, THEN MU=0.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)      <--> ON ENTRY; "A" IS MODEL HESSIAN (ONLY LOWER
C                  TRIANGULAR PART AND DIAGONAL STORED)
C                  ON EXIT:  A CONTAINS L OF LL+ DECOMPOSITION OF
C                  PERTURBED MODEL HESSIAN IN LOWER TRIANGULAR
C                  PART AND DIAGONAL AND CONTAINS HESSIAN IN UPPER
C                  TRIANGULAR PART AND UDIAG
C EPSM         --> MACHINE EPSILON
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C UDIAG(N)    <--  ON EXIT: CONTAINS DIAGONAL OF HESSIAN
C
C INTERNAL VARIABLES
C ------------------
C TOL              TOLERANCE
C DIAGMN           MINIMUM ELEMENT ON DIAGONAL OF A
C DIAGMX           MAXIMUM ELEMENT ON DIAGONAL OF A
C OFFMAX           MAXIMUM OFF-DIAGONAL ELEMENT OF A
C OFFROW           SUM OF OFF-DIAGONAL ELEMENTS IN A ROW OF A
C EVMIN            MINIMUM EIGENVALUE OF A
C EVMAX            MAXIMUM EIGENVALUE OF A
C
C DESCRIPTION
C -----------
C 1. IF "A" HAS ANY NEGATIVE DIAGONAL ELEMENTS, THEN CHOOSE MU>0
C SUCH THAT THE DIAGONAL OF A:=A+MU*I IS ALL POSITIVE
C WITH THE RATIO OF ITS SMALLEST TO LARGEST ELEMENT ON THE
C ORDER OF SQRT(EPSM).
C
C 2. "A" UNDERGOES A PERTURBED CHOLESKY DECOMPOSITION WHICH
C RESULTS IN AN LL+ DECOMPOSITION OF A+D, WHERE D IS A
C NON-NEGATIVE DIAGONAL MATRIX WHICH IS IMPLICITLY ADDED TO
C "A" DURING THE DECOMPOSITION IF "A" IS NOT POSITIVE DEFINITE.
C "A" IS RETAINED AND NOT CHANGED DURING THIS PROCESS BY
C COPYING L INTO THE UPPER TRIANGULAR PART OF "A" AND THE
C DIAGONAL INTO UDIAG.  THEN THE CHOLESKY DECOMPOSITION ROUTINE
C IS CALLED.  ON RETURN, ADDMAX CONTAINS MAXIMUM ELEMENT OF D.
C
C 3. IF ADDMAX=0, "A" WAS POSITIVE DEFINITE GOING INTO STEP 2
C AND RETURN IS MADE TO CALLING PROGRAM.  OTHERWISE,
C THE MINIMUM NUMBER SDD WHICH MUST BE ADDED TO THE
C DIAGONAL OF A TO MAKE IT SAFELY STRICTLY DIAGONALLY DOMINANT
C IS CALCULATED.  SINCE A+ADDMAX*I AND A+SDD*I ARE SAFELY
C POSITIVE DEFINITE, CHOOSE MU=MIN(ADDMAX,SDD) AND DECOMPOSE
C A+MU*I TO OBTAIN L.
C
      DIMENSION A(NR,1),SX(N),UDIAG(N)
C
C SCALE HESSIAN
C PRE- AND POST- MULTIPLY "A" BY INV(SX)
C
      DO 20 J=1,N
        DO 10 I=J,N
          A(I,J)=A(I,J)/(SX(I)*SX(J))
   10   CONTINUE
   20 CONTINUE
C
C STEP1
C -----
C NOTE:  IF A DIFFERENT TOLERANCE IS DESIRED THROUGHOUT THIS
C ALGORITHM, CHANGE TOLERANCE HERE:
      TOL=SQRT(EPSM)
C
      DIAGMX=A(1,1)
      DIAGMN=A(1,1)
      IF(N.EQ.1) GO TO 35
      DO 30 I=2,N
        IF(A(I,I).LT.DIAGMN) DIAGMN=A(I,I)
        IF(A(I,I).GT.DIAGMX) DIAGMX=A(I,I)
   30 CONTINUE
   35 POSMAX=MAX(DIAGMX,0.D0)
C
C DIAGMN .LE. 0
C
      IF(DIAGMN.GT.POSMAX*TOL) GO TO 100
C     IF(DIAGMN.LE.POSMAX*TOL)
C     THEN
        AMU=TOL*(POSMAX-DIAGMN)-DIAGMN
        IF(AMU.NE.0.D0) GO TO 60
C       IF(AMU.EQ.0.)
C       THEN
C
C FIND LARGEST OFF-DIAGONAL ELEMENT OF A
          OFFMAX=0.D0
          IF(N.EQ.1) GO TO 50
          DO 45 I=2,N
            IM1=I-1
            DO 40 J=1,IM1
              IF(ABS(A(I,J)).GT.OFFMAX) OFFMAX=ABS(A(I,J))
   40       CONTINUE
   45     CONTINUE
   50     AMU=OFFMAX
          IF(AMU.NE.0.D0) GO TO 55
C         IF(AMU.EQ.0.)
C         THEN
            AMU=1.0D0
            GO TO 60
C         ELSE
   55       AMU=AMU*(1.0D0+TOL)
C         ENDIF
C       ENDIF
C
C A=A + MU*I
C
   60   DO 65 I=1,N
          A(I,I)=A(I,I)+AMU
   65   CONTINUE
        DIAGMX=DIAGMX+AMU
C     ENDIF
C
C STEP2
C -----
C COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART
C AND DIAGONAL OF "A" TO UDIAG
C
  100 CONTINUE
      DO 110 J=1,N
        UDIAG(J)=A(J,J)
        IF(J.EQ.N) GO TO 110
        JP1=J+1
        DO 105 I=JP1,N
          A(J,I)=A(I,J)
  105   CONTINUE
  110 CONTINUE
C
      CALL CHLDCD(NR,N,A,DIAGMX,TOL,ADDMAX)
C
C
C STEP3
C -----
C IF ADDMAX=0, "A" WAS POSITIVE DEFINITE GOING INTO STEP 2,
C THE LL+ DECOMPOSITION HAS BEEN DONE, AND WE RETURN.
C OTHERWISE, ADDMAX>0.  PERTURB "A" SO THAT IT IS SAFELY
C DIAGONALLY DOMINANT AND FIND LL+ DECOMPOSITION
C
      IF(ADDMAX.LE.0.D0) GO TO 170
C     IF(ADDMAX.GT.0.)
C     THEN
C
C RESTORE ORIGINAL "A" (LOWER TRIANGULAR PART AND DIAGONAL)
C
        DO 120 J=1,N
          A(J,J)=UDIAG(J)
          IF(J.EQ.N) GO TO 120
          JP1=J+1
          DO 115 I=JP1,N
            A(I,J)=A(J,I)
  115     CONTINUE
  120   CONTINUE
C
C FIND SDD SUCH THAT A+SDD*I IS SAFELY POSITIVE DEFINITE
C NOTE:  EVMIN<0 SINCE A IS NOT POSITIVE DEFINITE;
C
        EVMIN=0.D0
        EVMAX=A(1,1)
        DO 150 I=1,N
          OFFROW=0.D0
          IF(I.EQ.1) GO TO 135
          IM1=I-1
          DO 130 J=1,IM1
            OFFROW=OFFROW+ABS(A(I,J))
  130     CONTINUE
  135     IF(I.EQ.N) GO TO 145
          IP1=I+1
          DO 140 J=IP1,N
            OFFROW=OFFROW+ABS(A(J,I))
  140     CONTINUE
  145     EVMIN=MIN(EVMIN,A(I,I)-OFFROW)
          EVMAX=MAX(EVMAX,A(I,I)+OFFROW)
  150   CONTINUE
        SDD=TOL*(EVMAX-EVMIN)-EVMIN
C
C PERTURB "A" AND DECOMPOSE AGAIN
C
        AMU=MIN(SDD,ADDMAX)
        DO 160 I=1,N
          A(I,I)=A(I,I)+AMU
          UDIAG(I)=A(I,I)
  160   CONTINUE
C
C "A" NOW GUARANTEED SAFELY POSITIVE DEFINITE
C
        CALL CHLDCD(NR,N,A,0.0D0,TOL,ADDMAX)
C     ENDIF
C
C UNSCALE HESSIAN AND CHOLESKY DECOMPOSITION MATRIX
C
  170 DO 190 J=1,N
        DO 175 I=J,N
          A(I,J)=SX(I)*A(I,J)
  175   CONTINUE
        IF(J.EQ.1) GO TO 185
        JM1=J-1
        DO 180 I=1,JM1
          A(I,J)=SX(I)*SX(J)*A(I,J)
  180   CONTINUE
  185   UDIAG(J)=UDIAG(J)*SX(J)*SX(J)
  190 CONTINUE
      RETURN
      END
