      SUBROUTINE HOKSTD(NR,N,G,A,UDIAG,P,SX,RNWTLN,DLT,AMU,
     +     DLTP,PHI,PHIP0,FSTIME,SC,NWTAKE,WRK0,EPSM,IPR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C PURPOSE
C -------
C FIND NEW STEP BY MORE-HEBDON ALGORITHM
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C G(N)         --> GRADIENT AT CURRENT ITERATE, G(X)
C A(N,N)       --> CHOLESKY DECOMPOSITION OF HESSIAN IN
C                  LOWER TRIANGULAR PART AND DIAGONAL.
C                  HESSIAN OR APPROX IN UPPER TRIANGULAR PART
C UDIAG(N)     --> DIAGONAL OF HESSIAN IN A(.,.)
C P(N)         --> NEWTON STEP
C SX(N)        --> DIAGONAL SCALING MATRIX FOR N
C RNWTLN       --> NEWTON STEP LENGTH
C DLT         <--> TRUST REGION RADIUS
C AMU         <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C DLTP         --> TRUST REGION RADIUS AT LAST EXIT FROM THIS ROUTINE
C PHI         <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C PHIP0       <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C FSTIME      <--> BOOLEAN. =.TRUE. IF FIRST ENTRY TO THIS ROUTINE
C                  DURING K-TH ITERATION
C SC(N)       <--  CURRENT STEP
C NWTAKE      <--  BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
C WRK0         --> WORKSPACE
C EPSM         --> MACHINE EPSILON
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION G(N),P(N),SX(N),SC(N),WRK0(N)
      DIMENSION A(NR,1),UDIAG(N)
      LOGICAL NWTAKE,DONE
      LOGICAL FSTIME
C
C HI AND ALO ARE CONSTANTS USED IN THIS ROUTINE.
C CHANGE HERE IF OTHER VALUES ARE TO BE SUBSTITUTED.
      IPR=IPR
      HI=1.5D0
      ALO=.75D0
C -----
      IF(RNWTLN.GT.HI*DLT) GO TO 15
C     IF(RNWTLN.LE.HI*DLT)
C     THEN
C
C       TAKE NEWTON STEP
C
        NWTAKE=.TRUE.
        DO 10 I=1,N
          SC(I)=P(I)
   10   CONTINUE
        DLT=MIN(DLT,RNWTLN)
        AMU=0.D0
C$      WRITE(IPR,951)
        RETURN
C     ELSE
C
C       NEWTON STEP NOT TAKEN
C
   15   CONTINUE
C$      WRITE(IPR,952)
        NWTAKE=.FALSE.
        IF(AMU.LE.0.D0) GO TO 20
C       IF(AMU.GT.0.)
C       THEN
          AMU=AMU- (PHI+DLTP) *((DLTP-DLT)+PHI)/(DLT*PHIP)
C$        WRITE(IPR,956) AMU
C       ENDIF
   20   CONTINUE
        PHI=RNWTLN-DLT
        IF(.NOT.FSTIME) GO TO 28
C       IF(FSTIME)
C       THEN
          DO 25 I=1,N
            WRK0(I)=SX(I)*SX(I)*P(I)
   25     CONTINUE
C
C         SOLVE L*Y = (SX**2)*P
C
          CALL FORSLD(NR,N,A,WRK0,WRK0)
          PHIP0=-DNRM2(N,WRK0,1)**2/RNWTLN
          FSTIME=.FALSE.
C       ENDIF
   28   PHIP=PHIP0
        AMULO=-PHI/PHIP
        AMUUP=0.0D0
        DO 30 I=1,N
          AMUUP=AMUUP+(G(I)*G(I))/(SX(I)*SX(I))
   30   CONTINUE
        AMUUP=SQRT(AMUUP)/DLT
        DONE=.FALSE.
C$      WRITE(IPR,956) AMU
C$      WRITE(IPR,959) PHI
C$      WRITE(IPR,960) PHIP
C$      WRITE(IPR,957) AMULO
C$      WRITE(IPR,958) AMUUP
C
C       TEST VALUE OF AMU; GENERATE NEXT AMU IF NECESSARY
C
  100   CONTINUE
        IF(DONE) RETURN
C$      WRITE(IPR,962)
        IF(AMU.GE.AMULO .AND. AMU.LE.AMUUP) GO TO 110
C       IF(AMU.LT.AMULO .OR.  AMU.GT.AMUUP)
C       THEN
          AMU=MAX(SQRT(AMULO*AMUUP),AMUUP*1.0D-3)
C$        WRITE(IPR,956) AMU
C       ENDIF
  110   CONTINUE
C
C       COPY (H,UDIAG) TO L
C       WHERE H <-- H+AMU*(SX**2) [DO NOT ACTUALLY CHANGE (H,UDIAG)]
        DO 130 J=1,N
          A(J,J)=UDIAG(J) + AMU*SX(J)*SX(J)
          IF(J.EQ.N) GO TO 130
          JP1=J+1
          DO 120 I=JP1,N
            A(I,J)=A(J,I)
  120     CONTINUE
  130   CONTINUE
C
C       FACTOR H=L(L+)
C
        CALL CHLDCD(NR,N,A,0.0D0,SQRT(EPSM),ADDMAX)
C
C       SOLVE H*P = L(L+)*SC = -G
C
        DO 140 I=1,N
          WRK0(I)=-G(I)
  140   CONTINUE
        CALL LLTSLD(NR,N,A,SC,WRK0)
C$      WRITE(IPR,955)
C$      WRITE(IPR,963) (SC(I),I=1,N)
C
C       RESET H.  NOTE SINCE UDIAG HAS NOT BEEN DESTROYED WE NEED DO
C       NOTHING HERE.  H IS IN THE UPPER PART AND IN UDIAG, STILL INTACT
C
        STEPLN=0.D0
        DO 150 I=1,N
          STEPLN=STEPLN + SX(I)*SX(I)*SC(I)*SC(I)
  150   CONTINUE
        STEPLN=SQRT(STEPLN)
        PHI=STEPLN-DLT
        DO 160 I=1,N
          WRK0(I)=SX(I)*SX(I)*SC(I)
  160   CONTINUE
        CALL FORSLD(NR,N,A,WRK0,WRK0)
        PHIP=-DNRM2(N,WRK0,1)**2/STEPLN
C$      WRITE(IPR,961) DLT,STEPLN
C$      WRITE(IPR,959) PHI
C$      WRITE(IPR,960) PHIP
        IF((ALO*DLT.GT.STEPLN .OR. STEPLN.GT.HI*DLT) .AND.
     +       (AMUUP-AMULO.GT.0.D0)) GO TO 170
C       IF((ALO*DLT.LE.STEPLN .AND. STEPLN.LE.HI*DLT) .OR.
C            (AMUUP-AMULO.LE.0.))
C       THEN
C
C         SC IS ACCEPTABLE HOKSTDEP
C
C$        WRITE(IPR,954)
          DONE=.TRUE.
          GO TO 100
C       ELSE
C
C         SC NOT ACCEPTABLE HOKSTDEP.  SELECT NEW AMU
C
  170     CONTINUE
C$        WRITE(IPR,953)
          AMULO=MAX(AMULO,AMU-(PHI/PHIP))
          IF(PHI.LT.0.D0) AMUUP=MIN(AMUUP,AMU)
          AMU=AMU-(STEPLN*PHI)/(DLT*PHIP)
C$        WRITE(IPR,956) AMU
C$        WRITE(IPR,957) AMULO
C$        WRITE(IPR,958) AMUUP
          GO TO 100
C       ENDIF
C     ENDIF
C
  951 FORMAT(27H0HOKSTD    TAKE NEWTON STEP)
  952 FORMAT(32H0HOKSTD    NEWTON STEP NOT TAKEN)
  953 FORMAT(31H HOKSTD    SC IS NOT ACCEPTABLE)
  954 FORMAT(27H HOKSTD    SC IS ACCEPTABLE)
  955 FORMAT(28H HOKSTD    CURRENT STEP (SC))
  956 FORMAT(18H HOKSTD    AMU   =,E20.13)
  957 FORMAT(18H HOKSTD    AMULO =,E20.13)
  958 FORMAT(18H HOKSTD    AMUUP =,E20.13)
  959 FORMAT(18H HOKSTD    PHI   =,E20.13)
  960 FORMAT(18H HOKSTD    PHIP  =,E20.13)
  961 FORMAT(18H HOKSTD    DLT   =,E20.13/
     +       18H HOKSTD    STEPLN=,E20.13)
  962 FORMAT(23H0HOKSTD    FIND NEW AMU)
  963 FORMAT(14H HOKSTD       ,5(E20.13,3X))
      END
