      SUBROUTINE HOKDRD(NR,N,X,F,G,A,UDIAG,P,XPLS,FPLS,FCN,SX,STEPMX,
     +     STEPTL,DLT,IRETCD,MXTAKE,AMU,DLTP,PHI,PHIP0,
     +     SC,XPLSP,WRK0,EPSM,ITNCNT,IPR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C PURPOSE
C -------
C FIND A NEXT NEWTON ITERATE (XPLS) BY THE MORE-HEBDON METHOD
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C X(N)         --> OLD ITERATE X[K-1]
C F            --> FUNCTION VALUE AT OLD ITERATE, F(X)
C G(N)         --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE
C A(N,N)       --> CHOLESKY DECOMPOSITION OF HESSIAN IN LOWER
C                  TRIANGULAR PART AND DIAGONAL.
C                  HESSIAN IN UPPER TRIANGULAR PART AND UDIAG.
C UDIAG(N)     --> DIAGONAL OF HESSIAN IN A(.,.)
C P(N)         --> NEWTON STEP
C XPLS(N)     <--  NEW ITERATE X[K]
C FPLS        <--  FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL       --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C                  CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C DLT         <--> TRUST REGION RADIUS
C IRETCD      <--  RETURN CODE
C                    =0 SATISFACTORY XPLS FOUND
C                    =1 FAILED TO FIND SATISFACTORY XPLS SUFFICIENTLY
C                       DISTINCT FROM X
C MXTAKE      <--  BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C AMU         <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C DLTP        <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C PHI         <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C PHIP0       <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C SC(N)        --> WORKSPACE
C XPLSP(N)     --> WORKSPACE
C WRK0(N)      --> WORKSPACE
C EPSM         --> MACHINE EPSILON
C ITNCNT       --> ITERATION COUNT
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION X(N),G(N),P(N),XPLS(N),SX(N)
      DIMENSION A(NR,1),UDIAG(N)
      DIMENSION SC(N),XPLSP(N),WRK0(N)
      LOGICAL MXTAKE,NWTAKE
      LOGICAL FSTIME
      EXTERNAL FCN
C
      IRETCD=4
      FSTIME=.TRUE.
      TMP=0.D0
      DO 5 I=1,N
        TMP=TMP+SX(I)*SX(I)*P(I)*P(I)
    5 CONTINUE
      RNWTLN=SQRT(TMP)
C$    WRITE(IPR,954) RNWTLN
C
      IF(ITNCNT.GT.1) GO TO 100
C     IF(ITNCNT.EQ.1)
C     THEN
        AMU=0.D0
C
C       IF FIRST ITERATION AND TRUST REGION NOT PROVIDED BY USER,
C       COMPUTE INITIAL TRUST REGION.
C
        IF(DLT.NE. (-1.D0)) GO TO 100
C       IF(DLT.EQ. (-1.))
C       THEN
          ALPHA=0.D0
          DO 10 I=1,N
            ALPHA=ALPHA+(G(I)*G(I))/(SX(I)*SX(I))
   10     CONTINUE
          BETA=0.0D0
          DO 30 I=1,N
            TMP=0.D0
            DO 20 J=I,N
              TMP=TMP + (A(J,I)*G(J))/(SX(J)*SX(J))
   20       CONTINUE
            BETA=BETA+TMP*TMP
   30     CONTINUE
          DLT=ALPHA*SQRT(ALPHA)/BETA
          DLT = MIN(DLT, STEPMX)
C$        WRITE(IPR,950)
C$        WRITE(IPR,951) ALPHA,BETA,DLT
C       ENDIF
C     ENDIF
C
  100 CONTINUE
C
C FIND NEW STEP BY MORE-HEBDON ALGORITHM
      CALL HOKSTD(NR,N,G,A,UDIAG,P,SX,RNWTLN,DLT,AMU,
     +     DLTP,PHI,PHIP0,FSTIME,SC,NWTAKE,WRK0,EPSM,IPR)
      DLTP=DLT
C
C CHECK NEW POINT AND UPDATE TRUST REGION
      CALL TRGUPD(NR,N,X,F,G,A,FCN,SC,SX,NWTAKE,STEPMX,STEPTL,
     +         DLT,IRETCD,XPLSP,FPLSP,XPLS,FPLS,MXTAKE,IPR,3,UDIAG)
      IF(IRETCD.LE.1) RETURN
      GO TO 100
C
  950 FORMAT(43H HOKDRD    INITIAL TRUST REGION NOT GIVEN. ,
     +       21H COMPUTE CAUCHY STEP.)
  951 FORMAT(18H HOKDRD    ALPHA =,E20.13/
     +       18H HOKDRD    BETA  =,E20.13/
     +       18H HOKDRD    DLT   =,E20.13)
  952 FORMAT(28H HOKDRD    CURRENT STEP (SC))
  954 FORMAT(18H0HOKDRD    RNWTLN=,E20.13)
  955 FORMAT(14H HOKDRD       ,5(E20.13,3X))
      END
