      SUBROUTINE OPTDRV(NR,N,X,FCN,D1FCN,D2FCN,TYPSIZ,FSCALE,
     +     METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,
     +     DLT,GRADTL,STEPMX,STEPTL,
     +     XPLS,FPLS,GPLS,ITRMCD,
     +     A,UDIAG,G,P,SX,WRK0,WRK1,WRK2,WRK3)
C
C PURPOSE
C -------
C DRIVER FOR NON-LINEAR OPTIMIZATION PROBLEM
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C X(N)         --> ON ENTRY: ESTIMATE TO A ROOT OF FCN
C FCN          --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION
C                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C                            FCN: R(N) --> R(1)
C D1FCN        --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE GRADIENT
C                  OF FCN.  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C D2FCN        --> (OPTIONAL) NAME OF SUBROUTINE TO EVALUATE HESSIAN OF
C                  OF FCN.  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C TYPSIZ(N)    --> TYPICAL SIZE FOR EACH COMPONENT OF X
C FSCALE       --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION
C METHOD       --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C                    =1 LINE SEARCH
C                    =2 DOUBLE DOGLEG
C                    =3 MORE-HEBDON
C IEXP         --> =1 IF OPTIMIZATION FUNCTION FCN IS EXPENSIVE TO
C                  EVALUATE, =0 OTHERWISE.  IF SET THEN HESSIAN WILL
C                  BE EVALUATED BY SECANT UPDATE INSTEAD OF
C                  ANALYTICALLY OR BY FINITE DIFFERENCES
C MSG         <--> ON INPUT:  (.GT.0) MESSAGE TO INHIBIT CERTAIN
C                    AUTOMATIC CHECKS
C                  ON OUTPUT: (.LT.0) ERROR CODE; =0 NO ERROR
C NDIGIT       --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN
C ITNLIM       --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C IAGFLG       --> =1 IF ANALYTIC GRADIENT SUPPLIED
C IAHFLG       --> =1 IF ANALYTIC HESSIAN SUPPLIED
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C DLT          --> TRUST REGION RADIUS
C GRADTL       --> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE
C                  ENOUGH TO ZERO TO TERMINATE ALGORITHM
C STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL       --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C                  CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C XPLS(N)     <--> ON EXIT:  XPLS IS LOCAL MINIMUM
C FPLS        <--> ON EXIT:  FUNCTION VALUE AT SOLUTION, XPLS
C GPLS(N)     <--> ON EXIT:  GRADIENT AT SOLUTION XPLS
C ITRMCD      <--  TERMINATION CODE
C A(N,N)       --> WORKSPACE FOR HESSIAN (OR ESTIMATE)
C                  AND ITS CHOLESKY DECOMPOSITION
C UDIAG(N)     --> WORKSPACE [FOR DIAGONAL OF HESSIAN]
C G(N)         --> WORKSPACE (FOR GRADIENT AT CURRENT ITERATE)
C P(N)         --> WORKSPACE FOR STEP
C SX(N)        --> WORKSPACE (FOR DIAGONAL SCALING MATRIX)
C WRK0(N)      --> WORKSPACE
C WRK1(N)      --> WORKSPACE
C WRK2(N)      --> WORKSPACE
C WRK3(N)      --> WORKSPACE
C
C
C INTERNAL VARIABLES
C ------------------
C ANALTL           TOLERANCE FOR COMPARISON OF ESTIMATED AND
C                  ANALYTICAL GRADIENTS AND HESSIANS
C EPSM             MACHINE EPSILON
C F                FUNCTION VALUE: FCN(X)
C ITNCNT           CURRENT ITERATION, K
C RNF              RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN.
C                       NOISE=10.**(-NDIGIT)
C
      DIMENSION X(N),XPLS(N),G(N),GPLS(N),P(N)
      DIMENSION TYPSIZ(N),SX(N)
      DIMENSION A(NR,1),UDIAG(N)
      DIMENSION WRK0(N),WRK1(N),WRK2(N),WRK3(N)
      LOGICAL MXTAKE,NOUPDT
      EXTERNAL FCN,D1FCN,D2FCN
C
C INITIALIZATION
C --------------
      DO 10 I=1,N
        P(I)=0.
   10 CONTINUE
      ITNCNT=0
      IRETCD=-1
      EPSM=R1MACH(4)
      CALL OPTCHK(N,X,TYPSIZ,SX,FSCALE,GRADTL,ITNLIM,NDIGIT,EPSM,
     +     DLT,METHOD,IEXP,IAGFLG,IAHFLG,STEPMX,MSG,IPR)
      IF(MSG.LT.0) RETURN
C***
      RNF=MAX(10.0**(-NDIGIT),EPSM)
C***
      ANALTL=MAX(1.0E-2,SQRT(RNF))
C
      IF(MOD(MSG/8,2).EQ.1) GO TO 15
      WRITE(IPR,901)
      WRITE(IPR,900) (TYPSIZ(I),I=1,N)
      WRITE(IPR,902)
      WRITE(IPR,900) (SX(I),I=1,N)
      WRITE(IPR,903) FSCALE
      WRITE(IPR,904) NDIGIT,IAGFLG,IAHFLG,IEXP,METHOD,ITNLIM,EPSM
      WRITE(IPR,905) STEPMX,STEPTL,GRADTL,DLT,RNF,ANALTL
   15 CONTINUE
C
C EVALUATE FCN(X)
C
      CALL FCN(N,X,F)
C
C EVALUATE ANALYTIC OR FINITE DIFFERENCE GRADIENT AND CHECK ANALYTIC
C GRADIENT, IF REQUESTED.
C
      IF (IAGFLG .EQ. 1) GO TO 20
C     IF (IAGFLG .EQ. 0)
C     THEN
        CALL FSTOFD (1, 1, N, X, FCN, F, G, SX, RNF, WRK, 1)
        GO TO 25
C
   20 CALL D1FCN (N, X, G)
      IF (MOD(MSG/2,2) .EQ. 1) GO TO 25
C     IF (MOD(MSG/2,2).EQ.0)
C     THEN
        CALL GRDCHK (N, X, FCN, F, G, TYPSIZ, SX, FSCALE,
     1    RNF, ANALTL, WRK1, MSG, IPR)
        IF (MSG .LT. 0) RETURN
   25 CONTINUE
C
      CALL OPTSTP(N,X,F,G,WRK1,ITNCNT,ICSCMX,
     +            ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,
     +            IPR,MSG)
      IF(ITRMCD.NE.0) GO TO 700
C
      IF(IEXP.NE.1) GO TO 80
C
C IF OPTIMIZATION FUNCTION EXPENSIVE TO EVALUATE (IEXP=1), THEN
C HESSIAN WILL BE OBTAINED BY SECANT UPDATES.  GET INITIAL HESSIAN.
C
      CALL HSNINT(NR,N,A,SX,METHOD)
      GO TO 90
   80 CONTINUE
C
C EVALUATE ANALYTIC OR FINITE DIFFERENCE HESSIAN AND CHECK ANALYTIC
C HESSIAN IF REQUESTED (ONLY IF USER-SUPPLIED ANALYTIC HESSIAN
C ROUTINE D2FCN FILLS ONLY LOWER TRIANGULAR PART AND DIAGONAL OF A).
C
      IF (IAHFLG .EQ. 1) GO TO 82
C     IF (IAHFLG .EQ. 0)
C     THEN
         IF (IAGFLG .EQ. 1) CALL FSTOFD (NR, N, N, X, D1FCN, G, A, SX,
     1      RNF, WRK1, 3)
         IF (IAGFLG .NE. 1) CALL SNDOFD (NR, N, X, FCN, F, A, SX, RNF,
     1      WRK1, WRK2)
         GO TO 88
C
C     ELSE
   82    IF (MOD(MSG/4,2).EQ.0) GO TO 85
C        IF (MOD(MSG/4, 2) .EQ. 1)
C        THEN
            CALL D2FCN (NR, N, X, A)
            GO TO 88
C
C        ELSE
   85       CALL HESCHK (NR, N, X, FCN, D1FCN, D2FCN, F, G, A, TYPSIZ,
     1         SX, RNF, ANALTL, IAGFLG, UDIAG, WRK1, WRK2, MSG, IPR)
C
C           HESCHK EVALUATES D2FCN AND CHECKS IT AGAINST THE FINITE
C           DIFFERENCE HESSIAN WHICH IT CALCULATES BY CALLING FSTOFD
C           (IF IAGFLG .EQ. 1) OR SNDOFD (OTHERWISE).
C
            IF (MSG .LT. 0) RETURN
   88 CONTINUE
C
   90 IF(MOD(MSG/8,2).EQ.0)
     +     CALL RESULT(NR,N,X,F,G,A,P,ITNCNT,1,IPR)
C
C
C ITERATION
C ---------
  100 ITNCNT=ITNCNT+1
C
C FIND PERTURBED LOCAL MODEL HESSIAN AND ITS LL+ DECOMPOSITION
C (SKIP THIS STEP IF LINE SEARCH OR DOGSTEP TECHNIQUES BEING USED WITH
C SECANT UPDATES.  CHOLESKY DECOMPOSITION L ALREADY OBTAINED FROM
C SECFAC.)
C
      IF(IEXP.EQ.1 .AND. METHOD.NE.3) GO TO 105
  103   CALL CHLHSN(NR,N,A,EPSM,SX,UDIAG)
  105 CONTINUE
C
C SOLVE FOR NEWTON STEP:  AP=-G
C
      DO 110 I=1,N
        WRK1(I)=-G(I)
  110 CONTINUE
      CALL LLTSLV(NR,N,A,P,WRK1)
C
C DECIDE WHETHER TO ACCEPT NEWTON STEP  XPLS=X + P
C OR TO CHOOSE XPLS BY A GLOBAL STRATEGY.
C
      IF (IAGFLG .NE. 0 .OR. METHOD .EQ. 1) GO TO 111
      DLTSAV = DLT
      IF (METHOD .EQ. 2) GO TO 111
      AMUSAV = AMU
      DLPSAV = DLTP
      PHISAV = PHI
      PHPSAV = PHIP0
  111 IF(METHOD.EQ.1)
     +     CALL LNSRCH(N,X,F,G,P,XPLS,FPLS,FCN,MXTAKE,IRETCD,
     +     STEPMX,STEPTL,SX,IPR)
      IF(METHOD.EQ.2)
     +     CALL DOGDRV(NR,N,X,F,G,A,P,XPLS,FPLS,FCN,SX,STEPMX,
     +     STEPTL,DLT,IRETCD,MXTAKE,WRK0,WRK1,WRK2,WRK3,IPR)
      IF(METHOD.EQ.3)
     +     CALL HOOKDR(NR,N,X,F,G,A,UDIAG,P,XPLS,FPLS,FCN,SX,STEPMX,
     +     STEPTL,DLT,IRETCD,MXTAKE,AMU,DLTP,PHI,PHIP0,WRK0,
     +     WRK1,WRK2,EPSM,ITNCNT,IPR)
C
C IF COULD NOT FIND SATISFACTORY STEP AND FORWARD DIFFERENCE
C GRADIENT WAS USED, RETRY USING CENTRAL DIFFERENCE GRADIENT.
C
      IF (IRETCD .NE. 1 .OR. IAGFLG .NE. 0) GO TO 112
C     IF (IRETCD .EQ. 1 .AND. IAGFLG .EQ. 0)
C     THEN
C
C        SET IAGFLG FOR CENTRAL DIFFERENCES
C
         IAGFLG = -1
         WRITE(IPR,906) ITNCNT
C
         CALL FSTOCD (N, X, FCN, SX, RNF, G)
         IF (METHOD .EQ. 1) GO TO 105
         DLT = DLTSAV
         IF (METHOD .EQ. 2) GO TO 105
         AMU = AMUSAV
         DLTP = DLPSAV
         PHI = PHISAV
         PHIP0 = PHPSAV
         GO TO 103
C     ENDIF
C
C CALCULATE STEP FOR OUTPUT
C
  112 CONTINUE
      DO 114 I = 1, N
         P(I) = XPLS(I) - X(I)
  114 CONTINUE
C
C CALCULATE GRADIENT AT XPLS
C
      IF (IAGFLG .EQ. (-1)) GO TO 116
      IF (IAGFLG .EQ. 0) GO TO 118
C
C ANALYTIC GRADIENT
      CALL D1FCN (N, XPLS, GPLS)
      GO TO 120
C
C CENTRAL DIFFERENCE GRADIENT
  116 CALL FSTOCD (N, XPLS, FCN, SX, RNF, GPLS)
      GO TO 120
C
C FORWARD DIFFERENCE GRADIENT
  118 CALL FSTOFD (1, 1, N, XPLS, FCN, FPLS, GPLS, SX, RNF, WRK, 1)
  120 CONTINUE
C
C CHECK WHETHER STOPPING CRITERIA SATISFIED
C
      CALL OPTSTP(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX,
     +            ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,
     +            IPR,MSG)
      IF(ITRMCD.NE.0) GO TO 690
C
C EVALUATE HESSIAN AT XPLS
C
      IF(IEXP.EQ.0) GO TO 130
      IF(METHOD.EQ.3)
     +     CALL SECUNF(NR,N,X,G,A,UDIAG,XPLS,GPLS,EPSM,ITNCNT,RNF,
     +     IAGFLG,NOUPDT,WRK1,WRK2,WRK3)
      IF(METHOD.NE.3)
     +     CALL SECFAC(NR,N,X,G,A,XPLS,GPLS,EPSM,ITNCNT,RNF,IAGFLG,
     +     NOUPDT,WRK0,WRK1,WRK2,WRK3)
      GO TO 150
  130 IF(IAHFLG.EQ.1) GO TO 140
      IF(IAGFLG.EQ.1)
     +     CALL FSTOFD(NR,N,N,XPLS,D1FCN,GPLS,A,SX,RNF,WRK1,3)
      IF(IAGFLG.NE.1) CALL SNDOFD(NR,N,XPLS,FCN,FPLS,A,SX,RNF,WRK1,WRK2)
      GO TO 150
  140 CALL D2FCN(NR,N,XPLS,A)
  150 CONTINUE
      IF(MOD(MSG/16,2).EQ.1)
     +     CALL RESULT(NR,N,XPLS,FPLS,GPLS,A,P,ITNCNT,1,IPR)
C
C X <-- XPLS  AND  G <-- GPLS  AND  F <-- FPLS
C
      F=FPLS
      DO 160 I=1,N
        X(I)=XPLS(I)
        G(I)=GPLS(I)
  160 CONTINUE
      GO TO 100
C
C TERMINATION
C -----------
C RESET XPLS,FPLS,GPLS,  IF PREVIOUS ITERATE SOLUTION
C
  690 IF(ITRMCD.NE.3) GO TO 710
  700 CONTINUE
      FPLS=F
      DO 705 I=1,N
        XPLS(I)=X(I)
        GPLS(I)=G(I)
  705 CONTINUE
C
C PRINT RESULTS
C
  710 CONTINUE
      IF(MOD(MSG/8,2).EQ.0)
     +     CALL RESULT(NR,N,XPLS,FPLS,GPLS,A,P,ITNCNT,0,IPR)
      MSG=0
      RETURN
C
  900 FORMAT(14H OPTDRV       ,5(E20.13,3X))
  901 FORMAT(20H0OPTDRV    TYPICAL X)
  902 FORMAT(40H OPTDRV    DIAGONAL SCALING MATRIX FOR X)
  903 FORMAT(22H OPTDRV    TYPICAL F =,E20.13)
  904 FORMAT(40H0OPTDRV    NUMBER OF GOOD DIGITS IN FCN=,I5/
     +       27H OPTDRV    GRADIENT FLAG  =,I5,18H   (=1 IF ANALYTIC,
     +       19H GRADIENT SUPPLIED)/
     +       27H OPTDRV    HESSIAN FLAG   =,I5,18H   (=1 IF ANALYTIC,
     +       18H HESSIAN SUPPLIED)/
     +       27H OPTDRV    EXPENSE FLAG   =,I5, 9H   (=1 IF,
     +       45H MINIMIZATION FUNCTION EXPENSIVE TO EVALUATE)/
     +       27H OPTDRV    METHOD TO USE  =,I5,19H   (=1,2,3 FOR LINE,
     +       49H SEARCH, DOUBLE DOGLEG, MORE-HEBDON RESPECTIVELY)/
     +       27H OPTDRV    ITERATION LIMIT=,I5/
     +       27H OPTDRV    MACHINE EPSILON=,E20.13)
  905 FORMAT(30H0OPTDRV    MAXIMUM STEP SIZE =,E20.13/
     +       30H OPTDRV    STEP TOLERANCE    =,E20.13/
     +       30H OPTDRV    GRADIENT TOLERANCE=,E20.13/
     +       30H OPTDRV    TRUST REG RADIUS  =,E20.13/
     +       30H OPTDRV    REL NOISE IN FCN  =,E20.13/
     +       30H OPTDRV    ANAL-FD TOLERANCE =,E20.13)
  906 FORMAT(52H OPTDRV    SHIFT FROM FORWARD TO CENTRAL DIFFERENCES,
     1   14H IN ITERATION , I5)
      END
