      SUBROUTINE UNCMIN (N,X0,FCN,X,F,INFO,W,LW)
C***BEGIN PROLOGUE  UNCMIN
C***DATE WRITTEN   870923    (YYMMDD)
C***REVISION DATE  871222    (YYMMDD)
C***CATEGORY NO.  G1B1A1
C***KEYWORDS  UNCONSTRAINED MINIMIZATION
C***AUTHOR  NASH, S.G., (GEORGE MASON UNIVERSITY)
C***PURPOSE  UNCMIN minimizes a smooth nonlinear function of n variables.
C     A subroutine that computes the function value at any point
C     must be supplied.  Derivative values are not required.
C     This subroutine provides the simplest interface to the UNCMIN
C     minimization package.  user has no control over options.
C
C***DESCRIPTION
C     This routine uses a quasi-Newton algorithm with line search
C     to minimize the function represented by the subroutine FCN.
C     At each iteration, the nonlinear function is approximated
C     by a quadratic function derived from a Taylor series.
C     The quadratic function is minimized to obtain a search direction,
C     and an approximate minimum of the nonlinear function along
C     the search direction is found using a line search.  The
C     algorithm computes an approximation to the second derivative
C     matrix of the nonlinear function using quasi-Newton techniques.
C
C     The UNCMIN package is quite general, and provides many options
C     for the user.  However, this subroutine is designed to be
C     easy to use, with few choices allowed.  For example:
C
C     1.  Only function values need be computed.  First derivative
C     values are obtained by finite-differencing.  This can be
C     very costly when the number of variables is large.
C
C     2.  It is assumed that the function values can be obtained
C     accurately (to an accuracy comparable to the precision of
C     the computer arithmetic).
C
C     3.  At most 150 iterations are allowed.
C
C     4.  It is assumed that the function values are well-scaled,
C     that is, that the optimal function value is not pathologically
C     large or small.
C
C     For more information, see the reference listed below.
C
C PARAMETERS
C ----------
C N            --> INTEGER
C                  Dimension of problem
C X0(N)        --> REAL
C                  Initial estimate of minimum
C FCN          --> Name of routine to evaluate minimization function.
C                  Must be declared EXTERNAL in calling routine, and
C                  have calling sequence
C                      SUBROUTINE FCN(N, X, F)
C                  with N and X as here, F the computed function value.
C X(N)        <--  REAL
C                  Local minimum
C F           <--  REAL
C                  Function value at local minimum X
C INFO        <--  INTEGER
C                  Termination code
C                      INFO =  0:  Optimal solution found
C                      INFO =  1:  Terminated with gradient small,
C                                  X is probably optimal
C                      INFO =  2:  Terminated with stepsize small,
C                                  X is probably optimal
C                      INFO =  3:  Lower point cannot be found,
C                                  X is probably optimal
C                      INFO =  4:  Iteration limit (150) exceeded
C                      INFO =  5:  Too many large steps,
C                                  function may be unbounded
C                      INFO = -1:  Insufficient workspace
C W(LW)        --> REAL
C                  Workspace
C LW           --> INTEGER
C                  Size of workspace, at least N*(N+10)
C
C***REFERENCES  R.B. SCHNABEL, J.E. KOONTZ, AND BE.E. WEISS, A MODULAR
C                 SYSTEM OF ALGORITHMS FOR UNCONSTRAINED MINIMIZATION,
C                 REPORT CU-CS-240-82, COMP. SCI. DEPT., UNIV. OF
C                 COLORADO AT BOULDER, 1982.
C***ROUTINES CALLED  OPTDRV, XERROR
C***END PROLOGUE  UNCMIN
      DIMENSION X0(N),X(N),W(LW)
      CHARACTER ERRMSG*68
      EXTERNAL  FCN, D1FCN, D2FCN
C----------------------------------------------------------------
C SUBDIVIDE WORKSPACE
C----------------------------------------------------------------
C***FIRST EXECUTABLE STATEMENT  UNCMIN
      IG  = 1
      IT  = IG  + N
      IW1 = IT  + N
      IW2 = IW1 + N
      IW3 = IW2 + N
      IW4 = IW3 + N
      IW5 = IW4 + N
      IW6 = IW5 + N
      IW7 = IW6 + N
      IW8 = IW7 + N
      IA  = IW8 + N
      LWMIN = IA + N*N-1
      IF (LWMIN .GT. LW) THEN
          INFO = -1
          WRITE(ERRMSG, '(
     *      ''UNCMIN ERROR (INFO=-1) -- INSUFFICIENT WORKSPACE'',
     *      '', LW = '', I5 )' ) LW
          CALL XERROR(ERRMSG(1:60), 60, -1, 0)
          RETURN
      END IF
C----------------------------------------------------------------
C SET UP PARAMETERS FOR OPTDRV
C----------------------------------------------------------------
C PARAMETERS THAT SHOULD NOT BE CHANGED WHEN USING CONDENSED CODE
C
C NR     = PARAMETER USED TO DIVIDE WORKSPACE
C METHOD = 1 (LINE SEARCH) -- DO NOT CHANGE
C MSG    = 9 => NO PRINTING, N=1 ALLOWED
C IAHFLG = 1 => ANALYTIC HESSIAN  SUPPLIED (0 OTHERWISE)
C IPR    = DEVICE FOR OUTPUT (IRRELEVANT IN CURRENT VERSION)
C DLT    = (IRRELEVANT PARAMETER FOR METHOD = 1)
C EPSM   = MACHINE EPSILON
C
      NR = N
      METHOD = 1
      MSG = 9
      IAHFLG = 0
      IPR = I1MACH(2)
      DLT = -1.0
      EPSM = R1MACH(4)
C
C PARAMETERS THAT MAY BE CHANGED:
C
C IEXP   = 1 => FUNCTION EXPENSIVE TO EVALUATE (IEXP = 0 => CHEAP)
C IAGFLG = 1 => ANALYTIC GRADIENT SUPPLIED (0 OTHERWISE)
C NDIGIT = -1 => OPTDRV ASSUMES F IS FULLY ACCURATE
C ITNLIM = 150 = MAXIMUM NUMBER OF ITERATIONS ALLOWED
C GRADTL = ZERO TOLERANCE FOR GRADIENT, FOR CONVERGENCE TESTS
C STEPMX = MAXIMUM ALLOWABLE STEP SIZE
C STEPTL = ZERO TOLERANCE FOR STEP, FOR CONVERGENCE TESTS
C FSCALE = TYPICAL ORDER-OF-MAGNITUDE SIZE OF FUNCTION
C TYPSIZ = TYPICAL ORDER-OF-MAGNITUDE SIZE OF X (STORED IN W(LT))
C
      IEXP = 1
      IAGFLG = 0
      NDIGIT = -1
      ITNLIM = 150
      GRADTL = EPSM**(1.0/3.0)
      STEPMX = 0.0
      STEPTL = SQRT(EPSM)
      FSCALE = 1.0
      DO 10 LT = IT,IT+N-1
          W(LT) = 1.0
   10 CONTINUE
C
C MINIMIZE FUNCTION
C
      CALL OPTDRV (NR, N, X0, FCN, D1FCN, D2FCN, W(IT), FSCALE,
     +             METHOD, IEXP, MSG, NDIGIT, ITNLIM, IAGFLG, IAHFLG,
     +             IPR, DLT, GRADTL, STEPMX, STEPTL,
     +             X, F, W(IG), INFO, W(IA),
     +             W(IW1), W(IW2), W(IW3), W(IW4),
     +             W(IW5), W(IW6), W(IW7), W(IW8))
C
      IF (INFO .EQ. 1) THEN
          WRITE(ERRMSG, '(
     *      ''UNCMIN WARNING -- INFO = 1'',
     *      '': PROBABLY CONVERGED, GRADIENT SMALL'')' )
          CALL XERROR(ERRMSG(1:62), 62, INFO, 0)
      END IF
      IF (INFO .EQ. 2) THEN
          WRITE(ERRMSG, '(
     *      ''UNCMIN WARNING -- INFO = 2'',
     *      '': PROBABLY CONVERGED, STEPSIZE SMALL'')' )
          CALL XERROR(ERRMSG(1:62), 62, INFO, 0)
      END IF
      IF (INFO .EQ. 3) THEN
          WRITE(ERRMSG, '(
     *      ''UNCMIN WARNING -- INFO = 3'',
     *      '': CANNOT FIND LOWER POINT'')' )
          CALL XERROR(ERRMSG(1:51), 51, INFO, 0)
      END IF
      IF (INFO .EQ. 4) THEN
          WRITE(ERRMSG, '(
     *      ''UNCMIN WARNING -- INFO = 4'',
     *      '': TOO MANY ITERATIONS'')' )
          CALL XERROR(ERRMSG(1:47), 47, INFO, 0)
      END IF
      IF (INFO .EQ. 5) THEN
          WRITE(ERRMSG, '(
     *      ''UNCMIN WARNING -- INFO = 5'',
     *      '': TOO MANY LARGE STEPS, POSSIBLY UNBOUNDED'')' )
          CALL XERROR(ERRMSG(1:68), 68, INFO, 0)
      END IF
C
      RETURN
      END
C SINGLE PRECISION VERSION OF UNCMIN
C  CHANGES MADE:
C     REMOVAL OF: IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C             AT BEGINNING OF ALL ROUTINES
C     CHANGE OF D1MACH TO R1MACH
C     CHANGE OF DOUBLE-PRECISION CONSTANTS TO SINGLE-PRECISION ONES
C             AT LINES AFTER C***
C     CHANGE OF DOUBLE-PRECISION BLAS TO SINGLE-PRECISION ONES
C             AT LINES AFTER C+++
C     NEW DRIVER:  UNCMIN
C     PROBLEMS OF SIZE N=1 ALLOWED (SEE DFAULT "C---")
C
      SUBROUTINE OPTIF0(NR,N,X,FCN,XPLS,FPLS,GPLS,ITRMCD,A,WRK)
C
C PURPOSE
C -------
C PROVIDE SIMPLEST INTERFACE TO MINIMIZATION PACKAGE.
C USER HAS NO CONTROL OVER OPTIONS.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C X(N)         --> INITIAL ESTIMATE OF MINIMUM
C FCN          --> NAME OF ROUTINE TO EVALUATE MINIMIZATION FUNCTION.
C                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE.
C XPLS(N)     <--  LOCAL MINIMUM
C FPLS        <--  FUNCTION VALUE AT LOCAL MINIMUM XPLS
C GPLS(N)     <--  GRADIENT AT LOCAL MINIMUM XPLS
C ITRMCD      <--  TERMINATION CODE
C A(N,N)       --> WORKSPACE
C WRK(N,9)     --> WORKSPACE
C
      DIMENSION X(N),XPLS(N),GPLS(N)
      DIMENSION A(NR,1),WRK(NR,1)
      EXTERNAL FCN,D1FCN,D2FCN
C
C EQUIVALENCE WRK(N,1) = UDIAG(N)
C             WRK(N,2) = G(N)
C             WRK(N,3) = P(N)
C             WRK(N,4) = TYPSIZ(N)
C             WRK(N,5) = SX(N)
C             WRK(N,6) = WRK0(N)
C             WRK(N,7) = WRK1(N)
C             WRK(N,8) = WRK2(N)
C             WRK(N,9) = WRK3(N)
C
      CALL DFAULT(N,X,WRK(1,4),FSCALE,METHOD,IEXP,MSG,NDIGIT,
     +     ITNLIM,IAGFLG,IAHFLG,IPR,DLT,GRADTL,STEPMX,STEPTL)
      CALL OPTDRV(NR,N,X,FCN,D1FCN,D2FCN,WRK(1,4),FSCALE,
     +     METHOD,IEXP,MSG,NDIGIT,ITNLIM,IAGFLG,IAHFLG,IPR,
     +     DLT,GRADTL,STEPMX,STEPTL,
     +     XPLS,FPLS,GPLS,ITRMCD,
     +     A,WRK(1,1),WRK(1,2),WRK(1,3),WRK(1,5),WRK(1,6),
     +     WRK(1,7),WRK(1,8),WRK(1,9))
      RETURN
      END
      SUBROUTINE BAKSLV(NR,N,A,X,B)
C
C PURPOSE
C -------
C SOLVE  AX=B  WHERE A IS UPPER TRIANGULAR MATRIX.
C NOTE THAT A IS INPUT AS A LOWER TRIANGULAR MATRIX AND
C THAT THIS ROUTINE TAKES ITS TRANSPOSE IMPLICITLY.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)       --> LOWER TRIANGULAR MATRIX (PRESERVED)
C X(N)        <--  SOLUTION VECTOR
C B(N)         --> RIGHT-HAND SIDE VECTOR
C
C NOTE
C ----
C IF B IS NO LONGER REQUIRED BY CALLING ROUTINE,
C THEN VECTORS B AND X MAY SHARE THE SAME STORAGE.
C
      DIMENSION A(NR,1),X(N),B(N)
C
C SOLVE (L-TRANSPOSE)X=B. (BACK SOLVE)
C
      I=N
      X(I)=B(I)/A(I,I)
      IF(N.EQ.1) RETURN
   30 IP1=I
      I=I-1
      SUM=0.
      DO 40 J=IP1,N
        SUM=SUM+A(J,I)*X(J)
   40 CONTINUE
      X(I)=(B(I)-SUM)/A(I,I)
      IF(I.GT.1) GO TO 30
      RETURN
      END
      SUBROUTINE CHLHSN(NR,N,A,EPSM,SX,UDIAG)
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
C***
   35 POSMAX=MAX(DIAGMX,0.0)
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.) GO TO 60
C       IF(AMU.EQ.0.)
C       THEN
C
C FIND LARGEST OFF-DIAGONAL ELEMENT OF A
          OFFMAX=0.
          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.) GO TO 55
C         IF(AMU.EQ.0.)
C         THEN
            AMU=1.0
            GO TO 60
C         ELSE
   55       AMU=AMU*(1.0+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 CHOLDC(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.) 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.
        EVMAX=A(1,1)
        DO 150 I=1,N
          OFFROW=0.
          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
C***
        CALL CHOLDC(NR,N,A,0.0,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
      SUBROUTINE CHOLDC(NR,N,A,DIAGMX,TOL,ADDMAX)
C
C PURPOSE
C -------
C FIND THE PERTURBED L(L-TRANSPOSE) [WRITTEN LL+] DECOMPOSITION
C OF A+D, WHERE D IS A NON-NEGATIVE DIAGONAL MATRIX ADDED TO A IF
C NECESSARY TO ALLOW THE CHOLESKY DECOMPOSITION TO CONTINUE.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)      <--> ON ENTRY: MATRIX FOR WHICH TO FIND PERTURBED
C                       CHOLESKY DECOMPOSITION
C                  ON EXIT:  CONTAINS L OF LL+ DECOMPOSITION
C                  IN LOWER TRIANGULAR PART AND DIAGONAL OF "A"
C DIAGMX       --> MAXIMUM DIAGONAL ELEMENT OF "A"
C TOL          --> TOLERANCE
C ADDMAX      <--  MAXIMUM AMOUNT IMPLICITLY ADDED TO DIAGONAL OF "A"
C                  IN FORMING THE CHOLESKY DECOMPOSITION OF A+D
C INTERNAL VARIABLES
C ------------------
C AMINL    SMALLEST ELEMENT ALLOWED ON DIAGONAL OF L
C AMNLSQ   =AMINL**2
C OFFMAX   MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN OF A
C
C
C DESCRIPTION
C -----------
C THE NORMAL CHOLESKY DECOMPOSITION IS PERFORMED.  HOWEVER, IF AT ANY
C POINT THE ALGORITHM WOULD ATTEMPT TO SET L(I,I)=SQRT(TEMP)
C WITH TEMP < TOL*DIAGMX, THEN L(I,I) IS SET TO SQRT(TOL*DIAGMX)
C INSTEAD.  THIS IS EQUIVALENT TO ADDING TOL*DIAGMX-TEMP TO A(I,I)
C
C
      DIMENSION A(NR,1)
C
      ADDMAX=0.
      AMINL=SQRT(DIAGMX*TOL)
      AMNLSQ=AMINL*AMINL
C
C FORM COLUMN J OF L
C
      DO 100 J=1,N
C FIND DIAGONAL ELEMENTS OF L
        SUM=0.
        IF(J.EQ.1) GO TO 20
        JM1=J-1
        DO 10 K=1,JM1
          SUM=SUM + A(J,K)*A(J,K)
   10   CONTINUE
   20   TEMP=A(J,J)-SUM
        IF(TEMP.LT.AMNLSQ) GO TO 30
C       IF(TEMP.GE.AMINL**2)
C       THEN
          A(J,J)=SQRT(TEMP)
          GO TO 40
C       ELSE
C
C FIND MAXIMUM OFF-DIAGONAL ELEMENT IN COLUMN
   30     OFFMAX=0.
          IF(J.EQ.N) GO TO 37
          JP1=J+1
          DO 35 I=JP1,N
            IF(ABS(A(I,J)).GT.OFFMAX) OFFMAX=ABS(A(I,J))
   35     CONTINUE
   37     IF(OFFMAX.LE.AMNLSQ) OFFMAX=AMNLSQ
C
C ADD TO DIAGONAL ELEMENT  TO ALLOW CHOLESKY DECOMPOSITION TO CONTINUE
          A(J,J)=SQRT(OFFMAX)
          ADDMAX=MAX(ADDMAX,OFFMAX-TEMP)
C       ENDIF
C
C FIND I,J ELEMENT OF LOWER TRIANGULAR MATRIX
   40   IF(J.EQ.N) GO TO 100
        JP1=J+1
        DO 70 I=JP1,N
          SUM=0.0
          IF(J.EQ.1) GO TO 60
          JM1=J-1
          DO 50 K=1,JM1
            SUM=SUM+A(I,K)*A(J,K)
   50     CONTINUE
   60     A(I,J)=(A(I,J)-SUM)/A(J,J)
   70   CONTINUE
  100 CONTINUE
      RETURN
      END
      SUBROUTINE D1FCN(N,X,G)
C
C PURPOSE
C -------
C DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC
C WHEN SPECIFIC ANALYTIC GRADIENT FUNCTION NOT SUPPLIED.
C
      DIMENSION X(N),G(N)
      G(N)=G(N)
      X(N)=X(N)
      STOP
      END
      SUBROUTINE D2FCN(NR,N,X,H)
C
C PURPOSE
C -------
C DUMMY ROUTINE TO PREVENT UNSATISFIED EXTERNAL DIAGNOSTIC
C WHEN SPECIFIC ANALYTIC HESSIAN FUNCTION NOT SUPPLIED.
C
      DIMENSION X(N),H(NR,1)
      H(NR,1)=H(NR,1)
      X(N)=X(N)
      STOP
      END
      SUBROUTINE DFAULT(N,X,TYPSIZ,FSCALE,METHOD,IEXP,MSG,NDIGIT,
     +     ITNLIM,IAGFLG,IAHFLG,IPR,DLT,GRADTL,STEPMX,STEPTL)
C
C PURPOSE
C -------
C SET DEFAULT VALUES FOR EACH INPUT VARIABLE TO
C MINIMIZATION ALGORITHM.
C
C PARAMETERS
C ----------
C N            --> DIMENSION OF PROBLEM
C X(N)         --> INITIAL GUESS TO SOLUTION (TO COMPUTE MAX STEP SIZE)
C TYPSIZ(N)   <--  TYPICAL SIZE FOR EACH COMPONENT OF X
C FSCALE      <--  ESTIMATE OF SCALE OF MINIMIZATION FUNCTION
C METHOD      <--  ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C IEXP        <--  =0 IF MINIMIZATION FUNCTION NOT EXPENSIVE TO EVALUATE
C MSG         <--  MESSAGE TO INHIBIT CERTAIN AUTOMATIC CHECKS + OUTPUT
C NDIGIT      <--  NUMBER OF GOOD DIGITS IN MINIMIZATION FUNCTION
C ITNLIM      <--  MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C IAGFLG      <--  =0 IF ANALYTIC GRADIENT NOT SUPPLIED
C IAHFLG      <--  =0 IF ANALYTIC HESSIAN NOT SUPPLIED
C IPR         <--  DEVICE TO WHICH TO SEND OUTPUT
C DLT         <--  TRUST REGION RADIUS
C GRADTL      <--  TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE ENOUGH
C                  TO ZERO TO TERMINATE ALGORITHM
C STEPMX      <--  VALUE OF ZERO TO TRIP DEFAULT MAXIMUM IN OPTCHK
C STEPTL      <--  TOLERANCE AT WHICH SUCCESSIVE ITERATES CONSIDERED
C                  CLOSE ENOUGH TO TERMINATE ALGORITHM
C
      DIMENSION TYPSIZ(N),X(N)
      X(N)=X(N)
C
C SET TYPICAL SIZE OF X AND MINIMIZATION FUNCTION
      DO 10 I=1,N
        TYPSIZ(I)=1.0
   10 CONTINUE
      FSCALE=1.0
C
C SET TOLERANCES
      DLT=-1.0
      EPSM=R1MACH(4)
      GRADTL=EPSM**(1.0/3.0)
      STEPMX=0.0
      STEPTL=SQRT(EPSM)
C
C SET FLAGS
      METHOD=1
      IEXP=1
C---
C     MSG=0
      MSG=9
      NDIGIT=-1
      ITNLIM=150
      IAGFLG=0
      IAHFLG=0
      IPR=I1MACH(2)
C
      RETURN
      END
      SUBROUTINE DOGDRV(NR,N,X,F,G,A,P,XPLS,FPLS,FCN,SX,STEPMX,
     +     STEPTL,DLT,IRETCD,MXTAKE,SC,WRK1,WRK2,WRK3,IPR)
C
C PURPOSE
C -------
C FIND A NEXT NEWTON ITERATE (XPLS) BY THE DOUBLE DOGLEG 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
C                  IN LOWER TRIANGULAR PART AND DIAGONAL
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                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
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 SC(N)        --> WORKSPACE [CURRENT STEP]
C WRK1(N)      --> WORKSPACE (AND PLACE HOLDING ARGUMENT TO TREGUP)
C WRK2(N)      --> WORKSPACE
C WRK3(N)      --> WORKSPACE
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION X(N),XPLS(N),G(N),P(N)
      DIMENSION SX(N)
      DIMENSION SC(N),WRK1(N),WRK2(N),WRK3(N)
      DIMENSION A(NR,1)
      LOGICAL FSTDOG,NWTAKE,MXTAKE
      EXTERNAL FCN
C
      IRETCD=4
      FSTDOG=.TRUE.
      TMP=0.
      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
  100 CONTINUE
C
C FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM
      CALL DOGSTP(NR,N,G,A,P,SX,RNWTLN,DLT,NWTAKE,FSTDOG,
     +     WRK1,WRK2,CLN,ETA,SC,IPR,STEPMX)
C
C CHECK NEW POINT AND UPDATE TRUST REGION
      CALL TREGUP(NR,N,X,F,G,A,FCN,SC,SX,NWTAKE,STEPMX,STEPTL,DLT,
     +     IRETCD,WRK3,FPLSP,XPLS,FPLS,MXTAKE,IPR,2,WRK1)
      IF(IRETCD.LE.1) RETURN
      GO TO 100
  950 FORMAT(42H DOGDRV    INITIAL TRUST REGION NOT GIVEN.,
     +       22H  COMPUTE CAUCHY STEP.)
  951 FORMAT(18H DOGDRV    ALPHA =,E20.13/
     +       18H DOGDRV    BETA  =,E20.13/
     +       18H DOGDRV    DLT   =,E20.13/
     +       18H DOGDRV    NWTAKE=,L1    )
  952 FORMAT(28H DOGDRV    CURRENT STEP (SC))
  954 FORMAT(18H0DOGDRV    RNWTLN=,E20.13)
  955 FORMAT(14H DOGDRV       ,5(E20.13,3X))
      END
      SUBROUTINE DOGSTP(NR,N,G,A,P,SX,RNWTLN,DLT,NWTAKE,FSTDOG,
     +     SSD,V,CLN,ETA,SC,IPR,STEPMX)
C
C PURPOSE
C -------
C FIND NEW STEP BY DOUBLE DOGLEG 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 PART AND DIAGONAL
C P(N)         --> NEWTON STEP
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C RNWTLN       --> NEWTON STEP LENGTH
C DLT         <--> TRUST REGION RADIUS
C NWTAKE      <--> BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
C FSTDOG      <--> BOOLEAN, =.TRUE. IF ON FIRST LEG OF DOGLEG
C SSD(N)      <--> WORKSPACE [CAUCHY STEP TO THE MINIMUM OF THE
C                  QUADRATIC MODEL IN THE SCALED STEEPEST DESCENT
C                  DIRECTION] [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C V(N)        <--> WORKSPACE  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C CLN         <--> CAUCHY LENGTH
C                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C ETA              [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C SC(N)       <--  CURRENT STEP
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE
C
C INTERNAL VARIABLES
C ------------------
C CLN              LENGTH OF CAUCHY STEP
C
      DIMENSION G(N),P(N)
      DIMENSION SX(N)
      DIMENSION SC(N),SSD(N),V(N)
      DIMENSION A(NR,1)
      LOGICAL NWTAKE,FSTDOG
      IPR=IPR
C
C CAN WE TAKE NEWTON STEP
C
      IF(RNWTLN.GT.DLT) GO TO 100
C     IF(RNWTLN.LE.DLT)
C     THEN
        NWTAKE=.TRUE.
        DO 10 I=1,N
          SC(I)=P(I)
   10   CONTINUE
        DLT=RNWTLN
C$      WRITE(IPR,951)
        GO TO 700
C     ELSE
C
C NEWTON STEP TOO LONG
C CAUCHY STEP IS ON DOUBLE DOGLEG CURVE
C
  100   NWTAKE=.FALSE.
        IF(.NOT.FSTDOG) GO TO 200
C       IF(FSTDOG)
C       THEN
C
C         CALCULATE DOUBLE DOGLEG CURVE (SSD)
          FSTDOG=.FALSE.
          ALPHA=0.
          DO 110 I=1,N
            ALPHA=ALPHA + (G(I)*G(I))/(SX(I)*SX(I))
  110     CONTINUE
          BETA=0.
          DO 130 I=1,N
            TMP=0.
            DO 120 J=I,N
              TMP=TMP + (A(J,I)*G(J))/(SX(J)*SX(J))
  120       CONTINUE
            BETA=BETA+TMP*TMP
  130     CONTINUE
          DO 140 I=1,N
            SSD(I)=-(ALPHA/BETA)*G(I)/SX(I)
  140     CONTINUE
          CLN=ALPHA*SQRT(ALPHA)/BETA
C+++
          ETA=.2 + (.8*ALPHA*ALPHA)/(-BETA*SDOT(N,G,1,P,1))
          DO 150 I=1,N
            V(I)=ETA*SX(I)*P(I) - SSD(I)
  150     CONTINUE
          IF (DLT .EQ. (-1.0)) DLT = MIN(CLN, STEPMX)
C$        WRITE(IPR,954) ALPHA,BETA,CLN,ETA
C$        WRITE(IPR,955)
C$        WRITE(IPR,960) (SSD(I),I=1,N)
C$        WRITE(IPR,956)
C$        WRITE(IPR,960) (V(I),I=1,N)
C       ENDIF
  200   IF(ETA*RNWTLN.GT.DLT) GO TO 220
C       IF(ETA*RNWTLN .LE. DLT)
C       THEN
C
C         TAKE PARTIAL STEP IN NEWTON DIRECTION
C
          DO 210 I=1,N
            SC(I)=(DLT/RNWTLN)*P(I)
  210     CONTINUE
C$        WRITE(IPR,957)
          GO TO 700
C       ELSE
  220     IF(CLN.LT.DLT) GO TO 240
C         IF(CLN.GE.DLT)
C         THEN
C           TAKE STEP IN STEEPEST DESCENT DIRECTION
C
            DO 230 I=1,N
              SC(I)=(DLT/CLN)*SSD(I)/SX(I)
  230       CONTINUE
C$          WRITE(IPR,958)
            GO TO 700
C         ELSE
C
C           CALCULATE CONVEX COMBINATION OF SSD AND ETA*P
C           WHICH HAS SCALED LENGTH DLT
C
C+++
  240       DOT1=SDOT(N,V,1,SSD,1)
C+++
            DOT2=SDOT(N,V,1,V,1)
            ALAM=(-DOT1+SQRT((DOT1*DOT1)-DOT2*(CLN*CLN-DLT*DLT)))/DOT2
            DO 250 I=1,N
              SC(I)=(SSD(I) + ALAM*V(I))/SX(I)
  250       CONTINUE
C$          WRITE(IPR,959)
C         ENDIF
C       ENDIF
C     ENDIF
  700 CONTINUE
C$    WRITE(IPR,952) FSTDOG,NWTAKE,RNWTLN,DLT
C$    WRITE(IPR,953)
C$    WRITE(IPR,960) (SC(I),I=1,N)
      RETURN
C
  951 FORMAT(27H0DOGSTP    TAKE NEWTON STEP)
  952 FORMAT(18H DOGSTP    FSTDOG=,L1/
     +       18H DOGSTP    NWTAKE=,L1/
     +       18H DOGSTP    RNWTLN=,E20.13/
     +       18H DOGSTP    DLT   =,E20.13)
  953 FORMAT(28H DOGSTP    CURRENT STEP (SC))
  954 FORMAT(18H DOGSTP    ALPHA =,E20.13/
     +       18H DOGSTP    BETA  =,E20.13/
     +       18H DOGSTP    CLN   =,E20.13/
     +       18H DOGSTP    ETA   =,E20.13)
  955 FORMAT(28H DOGSTP    CAUCHY STEP (SSD))
  956 FORMAT(12H DOGSTP    V)
  957 FORMAT(48H0DOGSTP    TAKE PARTIAL STEP IN NEWTON DIRECTION)
  958 FORMAT(50H0DOGSTP    TAKE STEP IN STEEPEST DESCENT DIRECTION)
  959 FORMAT(39H0DOGSTP    TAKE CONVEX COMBINATION STEP)
  960 FORMAT(14H DOGSTP       ,5(E20.13,3X))
      END
      SUBROUTINE FORSLV(NR,N,A,X,B)
C
C PURPOSE
C -------
C SOLVE  AX=B  WHERE A IS LOWER TRIANGULAR MATRIX
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)       --> LOWER TRIANGULAR MATRIX (PRESERVED)
C X(N)        <--  SOLUTION VECTOR
C B(N)         --> RIGHT-HAND SIDE VECTOR
C
C NOTE
C ----
C IF B IS NO LONGER REQUIRED BY CALLING ROUTINE,
C THEN VECTORS B AND X MAY SHARE THE SAME STORAGE.
C
      DIMENSION A(NR,1),X(N),B(N)
C
C SOLVE LX=B. (FOREWARD SOLVE)
C
      X(1)=B(1)/A(1,1)
      IF(N.EQ.1) RETURN
      DO 20 I=2,N
        SUM=0.0
        IM1=I-1
        DO 10 J=1,IM1
          SUM=SUM+A(I,J)*X(J)
   10   CONTINUE
        X(I)=(B(I)-SUM)/A(I,I)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE FSTOCD (N, X, FCN, SX, RNOISE, G)
C PURPOSE
C -------
C FIND CENTRAL DIFFERENCE APPROXIMATION G TO THE FIRST DERIVATIVE
C (GRADIENT) OF THE FUNCTION DEFINED BY FCN AT THE POINT X.
C
C PARAMETERS
C ----------
C N            --> DIMENSION OF PROBLEM
C X            --> POINT AT WHICH GRADIENT IS TO BE APPROXIMATED.
C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION.
C SX           --> DIAGONAL SCALING MATRIX FOR X.
C RNOISE       --> RELATIVE NOISE IN FCN [F(X)].
C G           <--  CENTRAL DIFFERENCE APPROXIMATION TO GRADIENT.
C
C
      DIMENSION X(N)
      DIMENSION SX(N)
      DIMENSION G(N)
      EXTERNAL FCN
C
C FIND I TH  STEPSIZE, EVALUATE TWO NEIGHBORS IN DIRECTION OF I TH
C UNIT VECTOR, AND EVALUATE I TH  COMPONENT OF GRADIENT.
C
      THIRD = 1.0/3.0
      DO 10 I = 1, N
         STEPI = RNOISE**THIRD * MAX(ABS(X(I)), 1.0/SX(I))
         XTEMPI = X(I)
         X(I) = XTEMPI + STEPI
         CALL FCN (N, X, FPLUS)
         X(I) = XTEMPI - STEPI
         CALL FCN (N, X, FMINUS)
         X(I) = XTEMPI
         G(I) = (FPLUS - FMINUS)/(2.0*STEPI)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE FSTOFD(NR,M,N,XPLS,FCN,FPLS,A,SX,RNOISE,FHAT,ICASE)
C PURPOSE
C -------
C FIND FIRST ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A" TO THE
C FIRST DERIVATIVE OF THE FUNCTION DEFINED BY THE SUBPROGRAM "FNAME"
C EVALUATED AT THE NEW ITERATE "XPLS".
C
C
C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE:
C 1) THE FIRST DERIVATIVE (GRADIENT) OF THE OPTIMIZATION FUNCTION "FCN
C    ANALYTIC USER ROUTINE HAS BEEN SUPPLIED;
C 2) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION
C    IF NO ANALYTIC USER ROUTINE HAS BEEN SUPPLIED FOR THE HESSIAN BUT
C    ONE HAS BEEN SUPPLIED FOR THE GRADIENT ("FCN") AND IF THE
C    OPTIMIZATION FUNCTION IS INEXPENSIVE TO EVALUATE
C
C NOTE
C ----
C _M=1 (OPTIMIZATION) ALGORITHM ESTIMATES THE GRADIENT OF THE FUNCTION
C      (FCN).   FCN(X) # F: R(N)-->R(1)
C _M=N (SYSTEMS) ALGORITHM ESTIMATES THE JACOBIAN OF THE FUNCTION
C      FCN(X) # F: R(N)-->R(N).
C _M=N (OPTIMIZATION) ALGORITHM ESTIMATES THE HESSIAN OF THE OPTIMIZATIO
C      FUNCTION, WHERE THE HESSIAN IS THE FIRST DERIVATIVE OF "FCN"
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C M            --> NUMBER OF ROWS IN A
C N            --> NUMBER OF COLUMNS IN A; DIMENSION OF PROBLEM
C XPLS(N)      --> NEW ITERATE:  X[K]
C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C FPLS(M)      --> _M=1 (OPTIMIZATION) FUNCTION VALUE AT NEW ITERATE:
C                       FCN(XPLS)
C                  _M=N (OPTIMIZATION) VALUE OF FIRST DERIVATIVE
C                       (GRADIENT) GIVEN BY USER FUNCTION FCN
C                  _M=N (SYSTEMS)  FUNCTION VALUE OF ASSOCIATED
C                       MINIMIZATION FUNCTION
C A(NR,N)     <--  FINITE DIFFERENCE APPROXIMATION (SEE NOTE).  ONLY
C                  LOWER TRIANGULAR MATRIX AND DIAGONAL ARE RETURNED
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C RNOISE       --> RELATIVE NOISE IN FCN [F(X)]
C FHAT(M)      --> WORKSPACE
C ICASE        --> =1 OPTIMIZATION (GRADIENT)
C                  =2 SYSTEMS
C                  =3 OPTIMIZATION (HESSIAN)
C
C INTERNAL VARIABLES
C ------------------
C STEPSZ - STEPSIZE IN THE J-TH VARIABLE DIRECTION
C
      DIMENSION XPLS(N),FPLS(M)
      DIMENSION FHAT(M)
      DIMENSION SX(N)
      DIMENSION A(NR,1)
C
C FIND J-TH COLUMN OF A
C EACH COLUMN IS DERIVATIVE OF F(FCN) WITH RESPECT TO XPLS(J)
C
      DO 30 J=1,N
        STEPSZ=SQRT(RNOISE)*MAX(ABS(XPLS(J)),1./SX(J))
        XTMPJ=XPLS(J)
        XPLS(J)=XTMPJ+STEPSZ
        CALL FCN(N,XPLS,FHAT)
        XPLS(J)=XTMPJ
        DO 20 I=1,M
          A(I,J)=(FHAT(I)-FPLS(I))/STEPSZ
   20   CONTINUE
   30 CONTINUE
      IF(ICASE.NE.3) RETURN
C
C IF COMPUTING HESSIAN, A MUST BE SYMMETRIC
C
      IF(N.EQ.1) RETURN
      NM1=N-1
      DO 50 J=1,NM1
        JP1=J+1
        DO 40 I=JP1,M
          A(I,J)=(A(I,J)+A(J,I))/2.0
   40   CONTINUE
   50 CONTINUE
      RETURN
      END
      SUBROUTINE GRDCHK(N,X,FCN,F,G,TYPSIZ,SX,FSCALE,RNF,
     +     ANALTL,WRK1,MSG,IPR)
C
C PURPOSE
C -------
C CHECK ANALYTIC GRADIENT AGAINST ESTIMATED GRADIENT
C
C PARAMETERS
C ----------
C N            --> DIMENSION OF PROBLEM
C X(N)         --> 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 F            --> FUNCTION VALUE:  FCN(X)
C G(N)         --> GRADIENT:  G(X)
C TYPSIZ(N)    --> TYPICAL SIZE FOR EACH COMPONENT OF X
C SX(N)        --> DIAGONAL SCALING MATRIX:  SX(I)=1./TYPSIZ(I)
C FSCALE       --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN
C RNF          --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C ANALTL       --> TOLERANCE FOR COMPARISON OF ESTIMATED AND
C                  ANALYTICAL GRADIENTS
C WRK1(N)      --> WORKSPACE
C MSG         <--  MESSAGE OR ERROR CODE
C                    ON OUTPUT: =-21, PROBABLE CODING ERROR OF GRADIENT
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION X(N),G(N)
      DIMENSION SX(N),TYPSIZ(N)
      DIMENSION WRK1(N)
      EXTERNAL FCN
C
C COMPUTE FIRST ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO
C ANALYTIC GRADIENT.
C
      CALL FSTOFD(1,1,N,X,FCN,F,WRK1,SX,RNF,WRK,1)
      KER=0
      DO 5 I=1,N
        GS=MAX(ABS(F),FSCALE)/MAX(ABS(X(I)),TYPSIZ(I))
        IF(ABS(G(I)-WRK1(I)).GT.MAX(ABS(G(I)),GS)*ANALTL) KER=1
    5 CONTINUE
      IF(KER.EQ.0) GO TO 20
        WRITE(IPR,901)
        WRITE(IPR,902) (I,G(I),WRK1(I),I=1,N)
        MSG=-21
   20 CONTINUE
      RETURN
  901 FORMAT(47H0GRDCHK    PROBABLE ERROR IN CODING OF ANALYTIC,
     +       19H GRADIENT FUNCTION./
     +       16H GRDCHK     COMP,12X,8HANALYTIC,12X,8HESTIMATE)
  902 FORMAT(11H GRDCHK    ,I5,3X,E20.13,3X,E20.13)
      END
      SUBROUTINE HESCHK(NR,N,X,FCN,D1FCN,D2FCN,F,G,A,TYPSIZ,SX,RNF,
     +     ANALTL,IAGFLG,UDIAG,WRK1,WRK2,MSG,IPR)
C
C PURPOSE
C -------
C CHECK ANALYTIC HESSIAN AGAINST ESTIMATED HESSIAN
C  (THIS MAY BE DONE ONLY IF THE USER SUPPLIED ANALYTIC HESSIAN
C   D2FCN FILLS ONLY THE LOWER TRIANGULAR PART AND DIAGONAL OF A)
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C X(N)         --> 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        --> NAME OF SUBROUTINE TO EVALUATE GRADIENT OF FCN.
C                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C D2FCN        --> NAME OF SUBROUTINE TO EVALUATE HESSIAN OF FCN.
C                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C F            --> FUNCTION VALUE:  FCN(X)
C G(N)        <--  GRADIENT:  G(X)
C A(N,N)      <--  ON EXIT:  HESSIAN IN LOWER TRIANGULAR PART AND DIAG
C TYPSIZ(N)    --> TYPICAL SIZE FOR EACH COMPONENT OF X
C SX(N)        --> DIAGONAL SCALING MATRIX:  SX(I)=1./TYPSIZ(I)
C RNF          --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C ANALTL       --> TOLERANCE FOR COMPARISON OF ESTIMATED AND
C                  ANALYTICAL GRADIENTS
C IAGFLG       --> =1 IF ANALYTIC GRADIENT SUPPLIED
C UDIAG(N)     --> WORKSPACE
C WRK1(N)      --> WORKSPACE
C WRK2(N)      --> WORKSPACE
C MSG         <--> MESSAGE OR ERROR CODE
C                    ON INPUT : IF =1XX DO NOT COMPARE ANAL + EST HESS
C                    ON OUTPUT: =-22, PROBABLE CODING ERROR OF HESSIAN
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION X(N),G(N),A(NR,1)
      DIMENSION TYPSIZ(N),SX(N)
      DIMENSION UDIAG(N),WRK1(N),WRK2(N)
      EXTERNAL FCN,D1FCN
C
C COMPUTE FINITE DIFFERENCE APPROXIMATION A TO THE HESSIAN.
C
      IF(IAGFLG.EQ.1) CALL FSTOFD(NR,N,N,X,D1FCN,G,A,SX,RNF,WRK1,3)
      IF(IAGFLG.NE.1) CALL SNDOFD(NR,N,X,FCN,F,A,SX,RNF,WRK1,WRK2)
C
      KER=0
C
C COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART
C AND DIAGONAL OF "A" TO UDIAG
C
      DO 30 J=1,N
        UDIAG(J)=A(J,J)
        IF(J.EQ.N) GO TO 30
        JP1=J+1
        DO 25 I=JP1,N
          A(J,I)=A(I,J)
   25   CONTINUE
   30 CONTINUE
C
C COMPUTE ANALYTIC HESSIAN AND COMPARE TO FINITE DIFFERENCE
C APPROXIMATION.
C
      CALL D2FCN(NR,N,X,A)
      DO 40 J=1,N
C***
        HS=MAX(ABS(G(J)),1.0)/MAX(ABS(X(J)),TYPSIZ(J))
        IF(ABS(A(J,J)-UDIAG(J)).GT.MAX(ABS(UDIAG(J)),HS)*ANALTL)
     +       KER=1
        IF(J.EQ.N) GO TO 40
        JP1=J+1
        DO 35 I=JP1,N
          IF(ABS(A(I,J)-A(J,I)).GT.MAX(ABS(A(I,J)),HS)*ANALTL) KER=1
   35   CONTINUE
   40 CONTINUE
C
      IF(KER.EQ.0) GO TO 90
        WRITE(IPR,901)
        DO 50 I=1,N
          IF(I.EQ.1) GO TO 45
          IM1=I-1
          DO 43 J=1,IM1
            WRITE(IPR,902) I,J,A(I,J),A(J,I)
   43     CONTINUE
   45     WRITE(IPR,902) I,I,A(I,I),UDIAG(I)
   50   CONTINUE
        MSG=-22
C     ENDIF
   90 CONTINUE
      RETURN
  901 FORMAT(47H HESCHK    PROBABLE ERROR IN CODING OF ANALYTIC,
     +       18H HESSIAN FUNCTION./
     +       21H HESCHK      ROW  COL,14X,8HANALYTIC,14X,10H(ESTIMATE))
  902 FORMAT(11H HESCHK    ,2I5,2X,E20.13,2X,1H(,E20.13,1H))
      END
      SUBROUTINE HOOKDR(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)
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.
      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.
C
C       IF FIRST ITERATION AND TRUST REGION NOT PROVIDED BY USER,
C       COMPUTE INITIAL TRUST REGION.
C
        IF(DLT.NE. (-1.)) GO TO 100
C       IF(DLT.EQ. (-1.))
C       THEN
          ALPHA=0.
          DO 10 I=1,N
            ALPHA=ALPHA+(G(I)*G(I))/(SX(I)*SX(I))
   10     CONTINUE
          BETA=0.0
          DO 30 I=1,N
            TMP=0.
            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 HOOKST(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 TREGUP(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 HOOKDR    INITIAL TRUST REGION NOT GIVEN. ,
     +       21H COMPUTE CAUCHY STEP.)
  951 FORMAT(18H HOOKDR    ALPHA =,E20.13/
     +       18H HOOKDR    BETA  =,E20.13/
     +       18H HOOKDR    DLT   =,E20.13)
  952 FORMAT(28H HOOKDR    CURRENT STEP (SC))
  954 FORMAT(18H0HOOKDR    RNWTLN=,E20.13)
  955 FORMAT(14H HOOKDR       ,5(E20.13,3X))
      END
      SUBROUTINE HOOKST(NR,N,G,A,UDIAG,P,SX,RNWTLN,DLT,AMU,
     +     DLTP,PHI,PHIP0,FSTIME,SC,NWTAKE,WRK0,EPSM,IPR)
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.5
      ALO=.75
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.
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.) 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 FORSLV(NR,N,A,WRK0,WRK0)
C+++
          PHIP0=-SNRM2(N,WRK0,1)**2/RNWTLN
          FSTIME=.FALSE.
C       ENDIF
   28   PHIP=PHIP0
        AMULO=-PHI/PHIP
        AMUUP=0.0
        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.0E-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
C***
        CALL CHOLDC(NR,N,A,0.0,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 LLTSLV(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.
        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 FORSLV(NR,N,A,WRK0,WRK0)
C+++
        PHIP=-SNRM2(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.)) 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 HOOKSTEP
C
C$        WRITE(IPR,954)
          DONE=.TRUE.
          GO TO 100
C       ELSE
C
C         SC NOT ACCEPTABLE HOOKSTEP.  SELECT NEW AMU
C
  170     CONTINUE
C$        WRITE(IPR,953)
          AMULO=MAX(AMULO,AMU-(PHI/PHIP))
          IF(PHI.LT.0.) 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(27H0HOOKST    TAKE NEWTON STEP)
  952 FORMAT(32H0HOOKST    NEWTON STEP NOT TAKEN)
  953 FORMAT(31H HOOKST    SC IS NOT ACCEPTABLE)
  954 FORMAT(27H HOOKST    SC IS ACCEPTABLE)
  955 FORMAT(28H HOOKST    CURRENT STEP (SC))
  956 FORMAT(18H HOOKST    AMU   =,E20.13)
  957 FORMAT(18H HOOKST    AMULO =,E20.13)
  958 FORMAT(18H HOOKST    AMUUP =,E20.13)
  959 FORMAT(18H HOOKST    PHI   =,E20.13)
  960 FORMAT(18H HOOKST    PHIP  =,E20.13)
  961 FORMAT(18H HOOKST    DLT   =,E20.13/
     +       18H HOOKST    STEPLN=,E20.13)
  962 FORMAT(23H0HOOKST    FIND NEW AMU)
  963 FORMAT(14H HOOKST       ,5(E20.13,3X))
      END
      SUBROUTINE HSNINT(NR,N,A,SX,METHOD)
C
C PURPOSE
C -------
C PROVIDE INITIAL HESSIAN WHEN USING SECANT UPDATES
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)      <--  INITIAL HESSIAN (LOWER TRIANGULAR MATRIX)
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C METHOD       --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C                    =1,2 FACTORED SECANT METHOD USED
C                    =3   UNFACTORED SECANT METHOD USED
C
      DIMENSION A(NR,1),SX(N)
C
      DO 100 J=1,N
        IF(METHOD.EQ.3) A(J,J)=SX(J)*SX(J)
        IF(METHOD.NE.3) A(J,J)=SX(J)
        IF(J.EQ.N) GO TO 100
        JP1=J+1
        DO 90 I=JP1,N
          A(I,J)=0.
   90   CONTINUE
  100 CONTINUE
      RETURN
      END
      SUBROUTINE LLTSLV(NR,N,A,X,B)
C
C PURPOSE
C -------
C SOLVE AX=B WHERE A HAS THE FORM L(L-TRANSPOSE)
C BUT ONLY THE LOWER TRIANGULAR PART, L, IS STORED.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)       --> MATRIX OF FORM L(L-TRANSPOSE).
C                  ON RETURN A IS UNCHANGED.
C X(N)        <--  SOLUTION VECTOR
C B(N)         --> RIGHT-HAND SIDE VECTOR
C
C NOTE
C ----
C IF B IS NOT REQUIRED BY CALLING PROGRAM, THEN
C B AND X MAY SHARE THE SAME STORAGE.
C
      DIMENSION A(NR,1),X(N),B(N)
C
C FORWARD SOLVE, RESULT IN X
C
      CALL FORSLV(NR,N,A,X,B)
C
C BACK SOLVE, RESULT IN X
C
      CALL BAKSLV(NR,N,A,X,X)
      RETURN
      END
      SUBROUTINE LNSRCH(N,X,F,G,P,XPLS,FPLS,FCN,MXTAKE,
     +   IRETCD,STEPMX,STEPTL,SX,IPR)
C PURPOSE
C -------
C FIND A NEXT NEWTON ITERATE BY LINE SEARCH.
C
C PARAMETERS
C ----------
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 P(N)         --> NON-ZERO 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 IRETCD      <--  RETURN CODE
C MXTAKE      <--  BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE
C STEPTL       --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C                  CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
C INTERNAL VARIABLES
C ------------------
C SLN              NEWTON LENGTH
C RLN              RELATIVE LENGTH OF NEWTON STEP
C
      INTEGER N,IRETCD
      DIMENSION SX(N)
      DIMENSION X(N),G(N),P(N)
      DIMENSION XPLS(N)
      LOGICAL MXTAKE
C
      IPR=IPR
      MXTAKE=.FALSE.
      IRETCD=2
C$    WRITE(IPR,954)
C$    WRITE(IPR,955) (P(I),I=1,N)
      TMP=0.0
      DO 5 I=1,N
        TMP=TMP+SX(I)*SX(I)*P(I)*P(I)
    5 CONTINUE
      SLN=SQRT(TMP)
      IF(SLN.LE.STEPMX) GO TO 10
C
C NEWTON STEP LONGER THAN MAXIMUM ALLOWED
        SCL=STEPMX/SLN
        CALL SCLMUL(N,SCL,P,P)
        SLN=STEPMX
C$      WRITE(IPR,954)
C$      WRITE(IPR,955) (P(I),I=1,N)
   10 CONTINUE
C+++
      SLP=SDOT(N,G,1,P,1)
      RLN=0.
      DO 15 I=1,N
        RLN=MAX(RLN,ABS(P(I))/MAX(ABS(X(I)),1./SX(I)))
   15 CONTINUE
      RMNLMB=STEPTL/RLN
      ALMBDA=1.0
C$    WRITE(IPR,952) SLN,SLP,RMNLMB,STEPMX,STEPTL
C
C LOOP
C CHECK IF NEW ITERATE SATISFACTORY.  GENERATE NEW LAMBDA IF NECESSARY.
C
  100 CONTINUE
      IF(IRETCD.LT.2) RETURN
      DO 105 I=1,N
        XPLS(I)=X(I) + ALMBDA*P(I)
  105 CONTINUE
      CALL FCN(N,XPLS,FPLS)
C$    WRITE(IPR,950) ALMBDA
C$    WRITE(IPR,951)
C$    WRITE(IPR,955) (XPLS(I),I=1,N)
C$    WRITE(IPR,953) FPLS
      IF(FPLS.GT. F+SLP*1.E-4*ALMBDA) GO TO 130
C     IF(FPLS.LE. F+SLP*1.E-4*ALMBDA)
C     THEN
C
C SOLUTION FOUND
C
        IRETCD=0
        IF(ALMBDA.EQ.1.0 .AND. SLN.GT. .99*STEPMX) MXTAKE=.TRUE.
        GO TO 100
C
C SOLUTION NOT (YET) FOUND
C
C     ELSE
  130   IF(ALMBDA .GE. RMNLMB) GO TO 140
C       IF(ALMBDA .LT. RMNLMB)
C       THEN
C
C NO SATISFACTORY XPLS FOUND SUFFICIENTLY DISTINCT FROM X
C
          IRETCD=1
          GO TO 100
C       ELSE
C
C CALCULATE NEW LAMBDA
C
  140     IF(ALMBDA.NE.1.0) GO TO 150
C         IF(ALMBDA.EQ.1.0)
C         THEN
C
C FIRST BACKTRACK: QUADRATIC FIT
C
            TLMBDA=-SLP/(2.*(FPLS-F-SLP))
            GO TO 170
C         ELSE
C
C ALL SUBSEQUENT BACKTRACKS: CUBIC FIT
C
  150       T1=FPLS-F-ALMBDA*SLP
            T2=PFPLS-F-PLMBDA*SLP
            T3=1.0/(ALMBDA-PLMBDA)
            A=T3*(T1/(ALMBDA*ALMBDA) - T2/(PLMBDA*PLMBDA))
            B=T3*(T2*ALMBDA/(PLMBDA*PLMBDA)
     +           - T1*PLMBDA/(ALMBDA*ALMBDA) )
            DISC=B*B-3.0*A*SLP
            IF(DISC.LE. B*B) GO TO 160
C           IF(DISC.GT. B*B)
C           THEN
C
C ONLY ONE POSITIVE CRITICAL POINT, MUST BE MINIMUM
C
C***
              TLMBDA=(-B+SIGN(1.0,A)*SQRT(DISC))/(3.0*A)
              GO TO 165
C           ELSE
C
C BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM
C
C***
  160         TLMBDA=(-B-SIGN(1.0,A)*SQRT(DISC))/(3.0*A)
C           ENDIF
  165       IF(TLMBDA.GT. .5*ALMBDA) TLMBDA=.5*ALMBDA
C         ENDIF
  170     PLMBDA=ALMBDA
          PFPLS=FPLS
          IF(TLMBDA.GE. ALMBDA*.1) GO TO 180
C         IF(TLMBDA.LT.ALMBDA/10.)
C         THEN
            ALMBDA=ALMBDA*.1
            GO TO 190
C         ELSE
  180       ALMBDA=TLMBDA
C         ENDIF
C       ENDIF
C     ENDIF
  190 GO TO 100
  950 FORMAT(18H LNSRCH    ALMBDA=,E20.13)
  951 FORMAT(29H LNSRCH    NEW ITERATE (XPLS))
  952 FORMAT(18H LNSRCH    SLN   =,E20.13/
     +       18H LNSRCH    SLP   =,E20.13/
     +       18H LNSRCH    RMNLMB=,E20.13/
     +       18H LNSRCH    STEPMX=,E20.13/
     +       18H LNSRCH    STEPTL=,E20.13)
  953 FORMAT(19H LNSRCH    F(XPLS)=,E20.13)
  954 FORMAT(26H0LNSRCH    NEWTON STEP (P))
  955 FORMAT(14H LNSRCH       ,5(E20.13,3X))
      END
      SUBROUTINE MVMLTL(NR,N,A,X,Y)
C
C PURPOSE
C -------
C COMPUTE Y=LX
C WHERE L IS A LOWER TRIANGULAR MATRIX STORED IN A
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)       --> LOWER TRIANGULAR (N*N) MATRIX
C X(N)         --> OPERAND VECTOR
C Y(N)        <--  RESULT VECTOR
C
C NOTE
C ----
C X AND Y CANNOT SHARE STORAGE
C
      DIMENSION A(NR,1),X(N),Y(N)
      DO 30 I=1,N
        SUM=0.
        DO 10 J=1,I
          SUM=SUM+A(I,J)*X(J)
   10   CONTINUE
        Y(I)=SUM
   30 CONTINUE
      RETURN
      END
      SUBROUTINE MVMLTS(NR,N,A,X,Y)
C
C PURPOSE
C -------
C COMPUTE Y=AX
C WHERE "A" IS A SYMMETRIC (N*N) MATRIX STORED IN ITS LOWER
C TRIANGULAR PART AND X,Y ARE N-VECTORS
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)       --> SYMMETRIC (N*N) MATRIX STORED IN
C                  LOWER TRIANGULAR PART AND DIAGONAL
C X(N)         --> OPERAND VECTOR
C Y(N)        <--  RESULT VECTOR
C
C NOTE
C ----
C X AND Y CANNOT SHARE STORAGE.
C
      DIMENSION A(NR,1),X(N),Y(N)
      DO 30 I=1,N
        SUM=0.
        DO 10 J=1,I
          SUM=SUM+A(I,J)*X(J)
   10   CONTINUE
        IF(I.EQ.N) GO TO 25
        IP1=I+1
        DO 20 J=IP1,N
          SUM=SUM+A(J,I)*X(J)
   20   CONTINUE
   25   Y(I)=SUM
   30 CONTINUE
      RETURN
      END
      SUBROUTINE MVMLTU(NR,N,A,X,Y)
C
C PURPOSE
C -------
C COMPUTE Y=(L+)X
C WHERE L IS A LOWER TRIANGULAR MATRIX STORED IN A
C (L-TRANSPOSE (L+) IS TAKEN IMPLICITLY)
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(NR,1)       --> LOWER TRIANGULAR (N*N) MATRIX
C X(N)         --> OPERAND VECTOR
C Y(N)        <--  RESULT VECTOR
C
C NOTE
C ----
C X AND Y CANNOT SHARE STORAGE
C
      DIMENSION A(NR,1),X(N),Y(N)
      DO 30 I=1,N
        SUM=0.
        DO 10 J=I,N
          SUM=SUM+A(J,I)*X(J)
   10   CONTINUE
        Y(I)=SUM
   30 CONTINUE
      RETURN
      END
      SUBROUTINE OPTCHK(N,X,TYPSIZ,SX,FSCALE,GRADTL,ITNLIM,NDIGIT,EPSM,
     +     DLT,METHOD,IEXP,IAGFLG,IAHFLG,STEPMX,MSG,IPR)
C
C PURPOSE
C -------
C CHECK INPUT FOR REASONABLENESS
C
C PARAMETERS
C ----------
C N            --> DIMENSION OF PROBLEM
C X(N)         --> ON ENTRY, ESTIMATE TO ROOT OF FCN
C TYPSIZ(N)   <--> TYPICAL SIZE OF EACH COMPONENT OF X
C SX(N)       <--  DIAGONAL SCALING MATRIX FOR X
C FSCALE      <--> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN
C GRADTL       --> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE
C                  ENOUGH TO ZERO TO TERMINATE ALGORITHM
C ITNLIM      <--> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C NDIGIT      <--> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN
C EPSM         --> MACHINE EPSILON
C DLT         <--> TRUST REGION RADIUS
C METHOD      <--> ALGORITHM INDICATOR
C IEXP        <--> EXPENSE FLAG
C IAGFLG      <--> =1 IF ANALYTIC GRADIENT SUPPLIED
C IAHFLG      <--> =1 IF ANALYTIC HESSIAN SUPPLIED
C STEPMX      <--> MAXIMUM STEP SIZE
C MSG         <--> MESSAGE AND ERROR CODE
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION X(N),TYPSIZ(N),SX(N)
C
C CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES.
C IF NOT, SET THEM TO DEFAULT VALUES.
      IF(METHOD.LT.1 .OR. METHOD.GT.3) METHOD=1
      IF(IAGFLG.NE.1) IAGFLG=0
      IF(IAHFLG.NE.1) IAHFLG=0
      IF(IEXP.NE.0) IEXP=1
      IF(MOD(MSG/2,2).EQ.1 .AND. IAGFLG.EQ.0) GO TO 830
      IF(MOD(MSG/4,2).EQ.1 .AND. IAHFLG.EQ.0) GO TO 835
C
C CHECK DIMENSION OF PROBLEM
C
      IF(N.LE.0) GO TO 805
      IF(N.EQ.1 .AND. MOD(MSG,2).EQ.0) GO TO 810
C
C COMPUTE SCALE MATRIX
C
      DO 10 I=1,N
        IF(TYPSIZ(I).EQ.0.) TYPSIZ(I)=1.0
        IF(TYPSIZ(I).LT.0.) TYPSIZ(I)=-TYPSIZ(I)
        SX(I)=1.0/TYPSIZ(I)
   10 CONTINUE
C
C CHECK MAXIMUM STEP SIZE
C
      IF (STEPMX .GT. 0.0) GO TO 20
      STPSIZ = 0.0
      DO 15 I = 1, N
         STPSIZ = STPSIZ + X(I)*X(I)*SX(I)*SX(I)
   15 CONTINUE
      STPSIZ = SQRT(STPSIZ)
C***
      STEPMX = MAX(1.0E3*STPSIZ, 1.0E3)
   20 CONTINUE
C CHECK FUNCTION SCALE
      IF(FSCALE.EQ.0.) FSCALE=1.0
      IF(FSCALE.LT.0.) FSCALE=-FSCALE
C
C CHECK GRADIENT TOLERANCE
      IF(GRADTL.LT.0.) GO TO 815
C
C CHECK ITERATION LIMIT
      IF(ITNLIM.LE.0) GO TO 820
C
C CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN
      IF(NDIGIT.EQ.0) GO TO 825
      IF(NDIGIT.LT.0) NDIGIT=-LOG10(EPSM)
C
C CHECK TRUST REGION RADIUS
      IF(DLT.LE.0.) DLT=-1.0
      IF (DLT .GT. STEPMX) DLT = STEPMX
      RETURN
C
C ERROR EXITS
C
  805 WRITE(IPR,901) N
      MSG=-1
      GO TO 895
  810 WRITE(IPR,902)
      MSG=-2
      GO TO 895
  815 WRITE(IPR,903) GRADTL
      MSG=-3
      GO TO 895
  820 WRITE(IPR,904) ITNLIM
      MSG=-4
      GO TO 895
  825 WRITE(IPR,905) NDIGIT
      MSG=-5
      GO TO 895
  830 WRITE(IPR,906) MSG,IAGFLG
      MSG=-6
      GO TO 895
  835 WRITE(IPR,907) MSG,IAHFLG
      MSG=-7
  895 RETURN
  901 FORMAT(32H0OPTCHK    ILLEGAL DIMENSION, N=,I5)
  902 FORMAT(55H0OPTCHK    +++ WARNING +++  THIS PACKAGE IS INEFFICIENT,
     +       26H FOR PROBLEMS OF SIZE N=1./
     +       48H OPTCHK    CHECK INSTALLATION LIBRARIES FOR MORE,
     +       22H APPROPRIATE ROUTINES./
     +       41H OPTCHK    IF NONE, SET MSG AND RESUBMIT.)
  903 FORMAT(38H0OPTCHK    ILLEGAL TOLERANCE.  GRADTL=,E20.13)
  904 FORMAT(44H0OPTCHK    ILLEGAL ITERATION LIMIT.  ITNLIM=,I5)
  905 FORMAT(52H0OPTCHK    MINIMIZATION FUNCTION HAS NO GOOD DIGITS.,
     +        9H  NDIGIT=,I5)
  906 FORMAT(50H0OPTCHK    USER REQUESTS THAT ANALYTIC GRADIENT BE,
     +       33H ACCEPTED AS PROPERLY CODED (MSG=,I5, 2H),/
     +       45H OPTCHK    BUT ANALYTIC GRADIENT NOT SUPPLIED,
     +        9H (IAGFLG=,I5, 2H).)
  907 FORMAT(49H0OPTCHK    USER REQUESTS THAT ANALYTIC HESSIAN BE,
     +       33H ACCEPTED AS PROPERLY CODED (MSG=,I5, 2H),/
     +       44H OPTCHK    BUT ANALYTIC HESSIAN NOT SUPPLIED,
     +        9H (IAHFLG=,I5, 2H).)
      END
      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
      SUBROUTINE OPTSTP(N,XPLS,FPLS,GPLS,X,ITNCNT,ICSCMX,
     +      ITRMCD,GRADTL,STEPTL,SX,FSCALE,ITNLIM,IRETCD,MXTAKE,IPR,MSG)
C
C UNCONSTRAINED MINIMIZATION STOPPING CRITERIA
C --------------------------------------------
C FIND WHETHER THE ALGORITHM SHOULD TERMINATE, DUE TO ANY
C OF THE FOLLOWING:
C 1) PROBLEM SOLVED WITHIN USER TOLERANCE
C 2) CONVERGENCE WITHIN USER TOLERANCE
C 3) ITERATION LIMIT REACHED
C 4) DIVERGENCE OR TOO RESTRICTIVE MAXIMUM STEP (STEPMX) SUSPECTED
C
C
C PARAMETERS
C ----------
C N            --> DIMENSION OF PROBLEM
C XPLS(N)      --> NEW ITERATE X[K]
C FPLS         --> FUNCTION VALUE AT NEW ITERATE F(XPLS)
C GPLS(N)      --> GRADIENT AT NEW ITERATE, G(XPLS), OR APPROXIMATE
C X(N)         --> OLD ITERATE X[K-1]
C ITNCNT       --> CURRENT ITERATION K
C ICSCMX      <--> NUMBER CONSECUTIVE STEPS .GE. STEPMX
C                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C ITRMCD      <--  TERMINATION CODE
C GRADTL       --> TOLERANCE AT WHICH RELATIVE GRADIENT CONSIDERED CLOSE
C                  ENOUGH TO ZERO TO TERMINATE ALGORITHM
C STEPTL       --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES
C                  CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C FSCALE       --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION
C ITNLIM       --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS
C IRETCD       --> RETURN CODE
C MXTAKE       --> BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C MSG          --> IF MSG INCLUDES A TERM 8, SUPPRESS OUTPUT
C
C
      INTEGER N,ITNCNT,ICSCMX,ITRMCD,ITNLIM
      DIMENSION SX(N)
      DIMENSION XPLS(N),GPLS(N),X(N)
      LOGICAL MXTAKE
C
      ITRMCD=0
C
C LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X
      IF(IRETCD.NE.1) GO TO 50
C     IF(IRETCD.EQ.1)
C     THEN
        JTRMCD=3
        GO TO 600
C     ENDIF
   50 CONTINUE
C
C FIND DIRECTION IN WHICH RELATIVE GRADIENT MAXIMUM.
C CHECK WHETHER WITHIN TOLERANCE
C
      D=MAX(ABS(FPLS),FSCALE)
      RGX=0.0
      DO 100 I=1,N
        RELGRD=ABS(GPLS(I))*MAX(ABS(XPLS(I)),1./SX(I))/D
        RGX=MAX(RGX,RELGRD)
  100 CONTINUE
      JTRMCD=1
      IF(RGX.LE.GRADTL) GO TO 600
C
      IF(ITNCNT.EQ.0) RETURN
C
C FIND DIRECTION IN WHICH RELATIVE STEPSIZE MAXIMUM
C CHECK WHETHER WITHIN TOLERANCE.
C
      RSX=0.0
      DO 120 I=1,N
        RELSTP=ABS(XPLS(I)-X(I))/MAX(ABS(XPLS(I)),1./SX(I))
        RSX=MAX(RSX,RELSTP)
  120 CONTINUE
      JTRMCD=2
      IF(RSX.LE.STEPTL) GO TO 600
C
C CHECK ITERATION LIMIT
C
      JTRMCD=4
      IF(ITNCNT.GE.ITNLIM) GO TO 600
C
C CHECK NUMBER OF CONSECUTIVE STEPS \ STEPMX
C
      IF(MXTAKE) GO TO 140
C     IF(.NOT.MXTAKE)
C     THEN
        ICSCMX=0
        RETURN
C     ELSE
  140   CONTINUE
        IF (MOD(MSG/8,2) .EQ. 0) WRITE(IPR,900)
        ICSCMX=ICSCMX+1
        IF(ICSCMX.LT.5) RETURN
        JTRMCD=5
C     ENDIF
C
C
C PRINT TERMINATION CODE
C
  600 ITRMCD=JTRMCD
      IF (MOD(MSG/8,2) .EQ. 0) GO TO(601,602,603,604,605), ITRMCD
      GO TO 700
  601 WRITE(IPR,901)
      GO TO 700
  602 WRITE(IPR,902)
      GO TO 700
  603 WRITE(IPR,903)
      GO TO 700
  604 WRITE(IPR,904)
      GO TO 700
  605 WRITE(IPR,905)
C
  700 RETURN
C
  900 FORMAT(48H0OPTSTP    STEP OF MAXIMUM LENGTH (STEPMX) TAKEN)
  901 FORMAT(43H0OPTSTP    RELATIVE GRADIENT CLOSE TO ZERO./
     +       48H OPTSTP    CURRENT ITERATE IS PROBABLY SOLUTION.)
  902 FORMAT(48H0OPTSTP    SUCCESSIVE ITERATES WITHIN TOLERANCE./
     +       48H OPTSTP    CURRENT ITERATE IS PROBABLY SOLUTION.)
  903 FORMAT(52H0OPTSTP    LAST GLOBAL STEP FAILED TO LOCATE A POINT,
     +       14H LOWER THAN X./
     +       51H OPTSTP    EITHER X IS AN APPROXIMATE LOCAL MINIMUM,
     +       17H OF THE FUNCTION,/
     +       50H OPTSTP    THE FUNCTION IS TOO NON-LINEAR FOR THIS,
     +       11H ALGORITHM,/
     +       34H OPTSTP    OR STEPTL IS TOO LARGE.)
  904 FORMAT(36H0OPTSTP    ITERATION LIMIT EXCEEDED./
     +       28H OPTSTP    ALGORITHM FAILED.)
  905 FORMAT(39H0OPTSTP    MAXIMUM STEP SIZE EXCEEDED 5,
     +       19H CONSECUTIVE TIMES./
     +       50H OPTSTP    EITHER THE FUNCTION IS UNBOUNDED BELOW,/
     +       47H OPTSTP    BECOMES ASYMPTOTIC TO A FINITE VALUE,
     +       30H FROM ABOVE IN SOME DIRECTION,/
     +       33H OPTSTP    OR STEPMX IS TOO SMALL)
      END
      SUBROUTINE QRAUX1(NR,N,R,I)
C
C PURPOSE
C -------
C INTERCHANGE ROWS I,I+1 OF THE UPPER HESSENBERG MATRIX R,
C COLUMNS I TO N
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF MATRIX
C R(N,N)      <--> UPPER HESSENBERG MATRIX
C I            --> INDEX OF ROW TO INTERCHANGE (I.LT.N)
C
      DIMENSION R(NR,1)
      DO 10 J=I,N
        TMP=R(I,J)
        R(I,J)=R(I+1,J)
        R(I+1,J)=TMP
   10 CONTINUE
      RETURN
      END
      SUBROUTINE QRAUX2(NR,N,R,I,A,B)
C
C PURPOSE
C -------
C PRE-MULTIPLY R BY THE JACOBI ROTATION J(I,I+1,A,B)
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF MATRIX
C R(N,N)      <--> UPPER HESSENBERG MATRIX
C I            --> INDEX OF ROW
C A            --> SCALAR
C B            --> SCALAR
C
      DIMENSION R(NR,1)
      DEN=SQRT(A*A + B*B)
      C=A/DEN
      S=B/DEN
      DO 10 J=I,N
        Y=R(I,J)
        Z=R(I+1,J)
        R(I,J)=C*Y - S*Z
        R(I+1,J)=S*Y + C*Z
   10 CONTINUE
      RETURN
      END
      SUBROUTINE QRUPDT(NR,N,A,U,V)
C
C PURPOSE
C -------
C FIND AN ORTHOGONAL (N*N) MATRIX (Q*) AND AN UPPER TRIANGULAR (N*N)
C MATRIX (R*) SUCH THAT (Q*)(R*)=R+U(V+)
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C A(N,N)      <--> ON INPUT:  CONTAINS R
C                  ON OUTPUT: CONTAINS (R*)
C U(N)         --> VECTOR
C V(N)         --> VECTOR
C
      DIMENSION A(NR,1)
      DIMENSION U(N),V(N)
C
C DETERMINE LAST NON-ZERO IN U(.)
C
      K=N
   10 IF(U(K).NE.0. .OR. K.EQ.1) GO TO 20
C     IF(U(K).EQ.0. .AND. K.GT.1)
C     THEN
        K=K-1
        GO TO 10
C     ENDIF
C
C (K-1) JACOBI ROTATIONS TRANSFORM
C     R + U(V+) --> (R*) + (U(1)*E1)(V+)
C WHICH IS UPPER HESSENBERG
C
   20 IF(K.LE.1) GO TO 40
        KM1=K-1
        DO 30 II=1,KM1
          I=KM1-II+1
          IF(U(I).NE.0.) GO TO 25
C         IF(U(I).EQ.0.)
C         THEN
            CALL QRAUX1(NR,N,A,I)
            U(I)=U(I+1)
            GO TO 30
C         ELSE
   25       CALL QRAUX2(NR,N,A,I,U(I),-U(I+1))
            U(I)=SQRT(U(I)*U(I) + U(I+1)*U(I+1))
C         ENDIF
   30   CONTINUE
C     ENDIF
C
C R <-- R + (U(1)*E1)(V+)
C
   40 DO 50 J=1,N
        A(1,J)=A(1,J) +U(1)*V(J)
   50 CONTINUE
C
C (K-1) JACOBI ROTATIONS TRANSFORM UPPER HESSENBERG R
C TO UPPER TRIANGULAR (R*)
C
      IF(K.LE.1) GO TO 100
        KM1=K-1
        DO 80 I=1,KM1
          IF(A(I,I).NE.0.) GO TO 70
C         IF(A(I,I).EQ.0.)
C         THEN
            CALL QRAUX1(NR,N,A,I)
            GO TO 80
C         ELSE
   70       T1=A(I,I)
            T2=-A(I+1,I)
            CALL QRAUX2(NR,N,A,I,T1,T2)
C         ENDIF
   80   CONTINUE
C     ENDIF
  100 RETURN
      END
      SUBROUTINE RESULT(NR,N,X,F,G,A,P,ITNCNT,IFLG,IPR)
C
C PURPOSE
C -------
C PRINT INFORMATION
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C X(N)         --> ITERATE X[K]
C F            --> FUNCTION VALUE AT X[K]
C G(N)         --> GRADIENT AT X[K]
C A(N,N)       --> HESSIAN AT X[K]
C P(N)         --> STEP TAKEN
C ITNCNT       --> ITERATION NUMBER K
C IFLG         --> FLAG CONTROLLING INFO TO PRINT
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION X(N),G(N),P(N),A(NR,1)
C PRINT ITERATION NUMBER
      WRITE(IPR,903) ITNCNT
      IF(IFLG.EQ.0) GO TO 120
C
C PRINT STEP
      WRITE(IPR,907)
      WRITE(IPR,905) (P(I),I=1,N)
C
C PRINT CURRENT ITERATE
  120 CONTINUE
      WRITE(IPR,904)
      WRITE(IPR,905) (X(I),I=1,N)
C
C PRINT FUNCTION VALUE
      WRITE(IPR,906)
      WRITE(IPR,905) F
C
C PRINT GRADIENT
      WRITE(IPR,908)
      WRITE(IPR,905) (G(I),I=1,N)
C
C PRINT HESSIAN FROM ITERATION K
      IF(IFLG.EQ.0) GO TO 140
      WRITE(IPR,901)
      DO 130 I=1,N
        WRITE(IPR,900) I
        WRITE(IPR,902) (A(I,J),J=1,I)
  130 CONTINUE
C
  140 RETURN
  900 FORMAT(15H RESULT     ROW,I5)
  901 FORMAT(29H RESULT       HESSIAN AT X(K))
  902 FORMAT(14H RESULT       ,5(2X,E20.13))
  903 FORMAT(/21H0RESULT    ITERATE K=,I5)
  904 FORMAT(18H RESULT       X(K))
  905 FORMAT(22H RESULT               ,5(2X,E20.13) )
  906 FORMAT(30H RESULT       FUNCTION AT X(K))
  907 FORMAT(18H RESULT       STEP)
  908 FORMAT(30H RESULT       GRADIENT AT X(K))
      END
      SUBROUTINE SCLMUL(N,S,V,Z)
C
C PURPOSE
C -------
C MULTIPLY VECTOR BY SCALAR
C RESULT VECTOR MAY BE OPERAND VECTOR
C
C PARAMETERS
C ----------
C N            --> DIMENSION OF VECTORS
C S            --> SCALAR
C V(N)         --> OPERAND VECTOR
C Z(N)        <--  RESULT VECTOR
      DIMENSION V(N),Z(N)
      DO 100 I=1,N
        Z(I)=S*V(I)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE SECFAC(NR,N,X,G,A,XPLS,GPLS,EPSM,ITNCNT,RNF,
     +     IAGFLG,NOUPDT,S,Y,U,W)
C
C PURPOSE
C -------
C UPDATE HESSIAN BY THE BFGS FACTORED 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 G(N)         --> GRADIENT OR APPROXIMATE AT OLD ITERATE
C A(N,N)      <--> ON ENTRY: CHOLESKY DECOMPOSITION OF HESSIAN IN
C                    LOWER PART AND DIAGONAL.
C                  ON EXIT:  UPDATED CHOLESKY DECOMPOSITION OF HESSIAN
C                    IN LOWER TRIANGULAR PART AND DIAGONAL
C XPLS(N)      --> NEW ITERATE, X[K]
C GPLS(N)      --> GRADIENT OR APPROXIMATE AT NEW ITERATE
C EPSM         --> MACHINE EPSILON
C ITNCNT       --> ITERATION COUNT
C RNF          --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C IAGFLG       --> =1 IF ANALYTIC GRADIENT SUPPLIED, =0 ITHERWISE
C NOUPDT      <--> BOOLEAN: NO UPDATE YET
C                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C S(N)         --> WORKSPACE
C Y(N)         --> WORKSPACE
C U(N)         --> WORKSPACE
C W(N)         --> WORKSPACE
C
      DIMENSION X(N),XPLS(N),G(N),GPLS(N)
      DIMENSION A(NR,1)
      DIMENSION S(N),Y(N),U(N),W(N)
      LOGICAL NOUPDT,SKPUPD
C
      IF(ITNCNT.EQ.1) NOUPDT=.TRUE.
      DO 10 I=1,N
        S(I)=XPLS(I)-X(I)
        Y(I)=GPLS(I)-G(I)
   10 CONTINUE
C+++
      DEN1=SDOT(N,S,1,Y,1)
C+++
      SNORM2=SNRM2(N,S,1)
C+++
      YNRM2=SNRM2(N,Y,1)
      IF(DEN1.LT.SQRT(EPSM)*SNORM2*YNRM2) GO TO 110
C     IF(DEN1.GE.SQRT(EPSM)*SNORM2*YNRM2)
C     THEN
        CALL MVMLTU(NR,N,A,S,U)
C+++
        DEN2=SDOT(N,U,1,U,1)
C
C       L <-- SQRT(DEN1/DEN2)*L
C
        ALP=SQRT(DEN1/DEN2)
        IF(.NOT.NOUPDT) GO TO 50
C       IF(NOUPDT)
C       THEN
          DO 30 J=1,N
            U(J)=ALP*U(J)
            DO 20 I=J,N
              A(I,J)=ALP*A(I,J)
   20       CONTINUE
   30     CONTINUE
          NOUPDT=.FALSE.
          DEN2=DEN1
          ALP=1.0
C       ENDIF
   50   SKPUPD=.TRUE.
C
C       W = L(L+)S = HS
C
        CALL MVMLTL(NR,N,A,U,W)
        I=1
        IF(IAGFLG.NE.0) GO TO 55
C       IF(IAGFLG.EQ.0)
C       THEN
          RELTOL=SQRT(RNF)
          GO TO 60
C       ELSE
   55     RELTOL=RNF
C       ENDIF
   60   IF(I.GT.N .OR. .NOT.SKPUPD) GO TO 70
C       IF(I.LE.N .AND. SKPUPD)
C       THEN
          IF(ABS(Y(I)-W(I)) .LT. RELTOL*MAX(ABS(G(I)),ABS(GPLS(I))))
     +         GO TO 65
C         IF(ABS(Y(I)-W(I)) .GE. RELTOL*AMAX1(ABS(G(I)),ABS(GPLS(I))))
C         THEN
            SKPUPD=.FALSE.
            GO TO 60
C         ELSE
   65       I=I+1
            GO TO 60
C         ENDIF
C       ENDIF
   70   IF(SKPUPD) GO TO 110
C       IF(.NOT.SKPUPD)
C       THEN
C
C         W=Y-ALP*L(L+)S
C
          DO 75 I=1,N
            W(I)=Y(I)-ALP*W(I)
   75     CONTINUE
C
C         ALP=1/SQRT(DEN1*DEN2)
C
          ALP=ALP/DEN1
C
C         U=(L+)/SQRT(DEN1*DEN2) = (L+)S/SQRT((Y+)S * (S+)L(L+)S)
C
          DO 80 I=1,N
            U(I)=ALP*U(I)
   80     CONTINUE
C
C         COPY L INTO UPPER TRIANGULAR PART.  ZERO L.
C
          IF(N.EQ.1) GO TO 93
          DO 90 I=2,N
            IM1=I-1
            DO 85 J=1,IM1
              A(J,I)=A(I,J)
              A(I,J)=0.
   85       CONTINUE
   90     CONTINUE
C
C         FIND Q, (L+) SUCH THAT  Q(L+) = (L+) + U(W+)
C
   93     CALL QRUPDT(NR,N,A,U,W)
C
C         UPPER TRIANGULAR PART AND DIAGONAL OF A NOW CONTAIN UPDATED
C         CHOLESKY DECOMPOSITION OF HESSIAN.  COPY BACK TO LOWER
C         TRIANGULAR PART.
C
          IF(N.EQ.1) GO TO 110
          DO 100 I=2,N
            IM1=I-1
            DO 95 J=1,IM1
              A(I,J)=A(J,I)
   95       CONTINUE
  100     CONTINUE
C       ENDIF
C     ENDIF
  110 RETURN
      END
      SUBROUTINE SECUNF(NR,N,X,G,A,UDIAG,XPLS,GPLS,EPSM,ITNCNT,
     +     RNF,IAGFLG,NOUPDT,S,Y,T)
C
C PURPOSE
C -------
C UPDATE HESSIAN BY THE BFGS UNFACTORED 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 G(N)         --> GRADIENT OR APPROXIMATE AT OLD ITERATE
C A(N,N)      <--> ON ENTRY: APPROXIMATE HESSIAN AT OLD ITERATE
C                    IN UPPER TRIANGULAR PART (AND UDIAG)
C                  ON EXIT:  UPDATED APPROX HESSIAN AT NEW ITERATE
C                    IN LOWER TRIANGULAR PART AND DIAGONAL
C                  [LOWER TRIANGULAR PART OF SYMMETRIC MATRIX]
C UDIAG        --> ON ENTRY: DIAGONAL OF HESSIAN
C XPLS(N)      --> NEW ITERATE, X[K]
C GPLS(N)      --> GRADIENT OR APPROXIMATE AT NEW ITERATE
C EPSM         --> MACHINE EPSILON
C ITNCNT       --> ITERATION COUNT
C RNF          --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C IAGFLG       --> =1 IF ANALYTIC GRADIENT SUPPLIED, =0 OTHERWISE
C NOUPDT      <--> BOOLEAN: NO UPDATE YET
C                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C S(N)         --> WORKSPACE
C Y(N)         --> WORKSPACE
C T(N)         --> WORKSPACE
C
      DIMENSION X(N),G(N),XPLS(N),GPLS(N)
      DIMENSION A(NR,1)
      DIMENSION UDIAG(N)
      DIMENSION S(N),Y(N),T(N)
      LOGICAL NOUPDT,SKPUPD
C
C COPY HESSIAN IN UPPER TRIANGULAR PART AND UDIAG TO
C LOWER TRIANGULAR PART AND DIAGONAL
C
      DO 5 J=1,N
        A(J,J)=UDIAG(J)
        IF(J.EQ.N) GO TO 5
        JP1=J+1
        DO 4 I=JP1,N
          A(I,J)=A(J,I)
    4   CONTINUE
    5 CONTINUE
C
      IF(ITNCNT.EQ.1) NOUPDT=.TRUE.
      DO 10 I=1,N
        S(I)=XPLS(I)-X(I)
        Y(I)=GPLS(I)-G(I)
   10 CONTINUE
C+++
      DEN1=SDOT(N,S,1,Y,1)
C+++
      SNORM2=SNRM2(N,S,1)
C+++
      YNRM2=SNRM2(N,Y,1)
      IF(DEN1.LT.SQRT(EPSM)*SNORM2*YNRM2) GO TO 100
C     IF(DEN1.GE.SQRT(EPSM)*SNORM2*YNRM2)
C     THEN
        CALL MVMLTS(NR,N,A,S,T)
C+++
        DEN2=SDOT(N,S,1,T,1)
        IF(.NOT. NOUPDT) GO TO 50
C       IF(NOUPDT)
C       THEN
C
C         H <-- [(S+)Y/(S+)HS]H
C
          GAM=DEN1/DEN2
          DEN2=GAM*DEN2
          DO 30 J=1,N
            T(J)=GAM*T(J)
            DO 20 I=J,N
              A(I,J)=GAM*A(I,J)
   20       CONTINUE
   30     CONTINUE
          NOUPDT=.FALSE.
C       ENDIF
   50   SKPUPD=.TRUE.
C
C       CHECK UPDATE CONDITION ON ROW I
C
        DO 60 I=1,N
          TOL=RNF*MAX(ABS(G(I)),ABS(GPLS(I)))
          IF(IAGFLG.EQ.0) TOL=TOL/SQRT(RNF)
          IF(ABS(Y(I)-T(I)).LT.TOL) GO TO 60
C         IF(ABS(Y(I)-T(I)).GE.TOL)
C         THEN
            SKPUPD=.FALSE.
            GO TO 70
C         ENDIF
   60   CONTINUE
   70   IF(SKPUPD) GO TO 100
C       IF(.NOT.SKPUPD)
C       THEN
C
C         BFGS UPDATE
C
          DO 90 J=1,N
            DO 80 I=J,N
              A(I,J)=A(I,J)+Y(I)*Y(J)/DEN1-T(I)*T(J)/DEN2
   80       CONTINUE
   90     CONTINUE
C       ENDIF
C     ENDIF
  100 RETURN
      END
      SUBROUTINE SNDOFD(NR,N,XPLS,FCN,FPLS,A,SX,RNOISE,STEPSZ,ANBR)
C PURPOSE
C -------
C FIND SECOND ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A"
C TO THE SECOND DERIVATIVE (HESSIAN) OF THE FUNCTION DEFINED BY THE SUBP
C "FCN" EVALUATED AT THE NEW ITERATE "XPLS"
C
C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE
C 1) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION
C    IF NO ANALYTICAL USER FUNCTION HAS BEEN SUPPLIED FOR EITHER
C    THE GRADIENT OR THE HESSIAN AND IF THE OPTIMIZATION FUNCTION
C    "FCN" IS INEXPENSIVE TO EVALUATE.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C XPLS(N)      --> NEW ITERATE:   X[K]
C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C FPLS         --> FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C A(N,N)      <--  FINITE DIFFERENCE APPROXIMATION TO HESSIAN
C                  ONLY LOWER TRIANGULAR MATRIX AND DIAGONAL
C                  ARE RETURNED
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C RNOISE       --> RELATIVE NOISE IN FNAME [F(X)]
C STEPSZ(N)    --> WORKSPACE (STEPSIZE IN I-TH COMPONENT DIRECTION)
C ANBR(N)      --> WORKSPACE (NEIGHBOR IN I-TH DIRECTION)
C
C
      DIMENSION XPLS(N)
      DIMENSION SX(N)
      DIMENSION STEPSZ(N),ANBR(N)
      DIMENSION A(NR,1)
C
C FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION
C OF I-TH UNIT VECTOR.
C
      OV3 = 1.0/3.0
      DO 10 I=1,N
        STEPSZ(I)=RNOISE**OV3 * MAX(ABS(XPLS(I)),1./SX(I))
        XTMPI=XPLS(I)
        XPLS(I)=XTMPI+STEPSZ(I)
        CALL FCN(N,XPLS,ANBR(I))
        XPLS(I)=XTMPI
   10 CONTINUE
C
C CALCULATE COLUMN I OF A
C
      DO 30 I=1,N
        XTMPI=XPLS(I)
        XPLS(I)=XTMPI+2.0*STEPSZ(I)
        CALL FCN(N,XPLS,FHAT)
        A(I,I)=((FPLS-ANBR(I))+(FHAT-ANBR(I)))/(STEPSZ(I)*STEPSZ(I))
C
C CALCULATE SUB-DIAGONAL ELEMENTS OF COLUMN
        IF(I.EQ.N) GO TO 25
        XPLS(I)=XTMPI+STEPSZ(I)
        IP1=I+1
        DO 20 J=IP1,N
          XTMPJ=XPLS(J)
          XPLS(J)=XTMPJ+STEPSZ(J)
          CALL FCN(N,XPLS,FHAT)
          A(J,I)=((FPLS-ANBR(I))+(FHAT-ANBR(J)))/(STEPSZ(I)*STEPSZ(J))
          XPLS(J)=XTMPJ
   20   CONTINUE
   25   XPLS(I)=XTMPI
   30 CONTINUE
      RETURN
      END
      SUBROUTINE TREGUP(NR,N,X,F,G,A,FCN,SC,SX,NWTAKE,STEPMX,STEPTL,
     +     DLT,IRETCD,XPLSP,FPLSP,XPLS,FPLS,MXTAKE,IPR,METHOD,UDIAG)
C
C PURPOSE
C -------
C DECIDE WHETHER TO ACCEPT XPLS=X+SC AS THE NEXT ITERATE AND UPDATE THE
C TRUST REGION DLT.
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
C                  LOWER TRIANGULAR PART AND DIAGONAL.
C                  HESSIAN OR APPROX IN UPPER TRIANGULAR PART
C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C SC(N)        --> CURRENT STEP
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C NWTAKE       --> BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
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 XPLS ACCEPTED AS NEXT ITERATE;
C                       DLT TRUST REGION FOR NEXT ITERATION.
C                    =1 XPLS UNSATISFACTORY BUT ACCEPTED AS NEXT ITERATE
C                       BECAUSE XPLS-X .LT. SMALLEST ALLOWABLE
C                       STEP LENGTH.
C                    =2 F(XPLS) TOO LARGE.  CONTINUE CURRENT ITERATION
C                       WITH NEW REDUCED DLT.
C                    =3 F(XPLS) SUFFICIENTLY SMALL, BUT QUADRATIC MODEL
C                       PREDICTS F(XPLS) SUFFICIENTLY WELL TO CONTINUE
C                       CURRENT ITERATION WITH NEW DOUBLED DLT.
C XPLSP(N)    <--> WORKSPACE [VALUE NEEDS TO BE RETAINED BETWEEN
C                  SUCCESIVE CALLS OF K-TH GLOBAL STEP]
C FPLSP       <--> [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C XPLS(N)     <--  NEW ITERATE X[K]
C FPLS        <--  FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C MXTAKE      <--  BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C METHOD       --> ALGORITHM TO USE TO SOLVE MINIMIZATION PROBLEM
C                    =1 LINE SEARCH
C                    =2 DOUBLE DOGLEG
C                    =3 MORE-HEBDON
C UDIAG(N)     --> DIAGONAL OF HESSIAN IN A(.,.)
C
      DIMENSION X(N),XPLS(N),G(N)
      DIMENSION SX(N),SC(N),XPLSP(N)
      DIMENSION A(NR,1)
      LOGICAL NWTAKE,MXTAKE
      DIMENSION UDIAG(N)
C
      IPR=IPR
      MXTAKE=.FALSE.
      DO 100 I=1,N
        XPLS(I)=X(I)+SC(I)
  100 CONTINUE
      CALL FCN(N,XPLS,FPLS)
      DLTF=FPLS-F
C+++
      SLP=SDOT(N,G,1,SC,1)
C
C NEXT STATEMENT ADDED FOR CASE OF COMPILERS WHICH DO NOT OPTIMIZE
C EVALUATION OF NEXT "IF" STATEMENT (IN WHICH CASE FPLSP COULD BE
C UNDEFINED).
      IF(IRETCD.EQ.4) FPLSP=0.0
C$    WRITE(IPR,961) IRETCD,FPLS,FPLSP,DLTF,SLP
      IF(IRETCD.NE.3 .OR. (FPLS.LT.FPLSP .AND. DLTF.LE. 1.E-4*SLP))
     +                                                     GO TO 130
C     IF(IRETCD.EQ.3 .AND. (FPLS.GE.FPLSP .OR. DLTF.GT. 1.E-4*SLP))
C     THEN
C
C       RESET XPLS TO XPLSP AND TERMINATE GLOBAL STEP
C
        IRETCD=0
        DO 110 I=1,N
          XPLS(I)=XPLSP(I)
  110   CONTINUE
        FPLS=FPLSP
        DLT=.5*DLT
C$      WRITE(IPR,951)
        GO TO 230
C     ELSE
C
C       FPLS TOO LARGE
C
  130   IF(DLTF.LE. 1.E-4*SLP) GO TO 170
C       IF(DLTF.GT. 1.E-4*SLP)
C       THEN
C$        WRITE(IPR,952)
          RLN=0.
          DO 140 I=1,N
            RLN=MAX(RLN,ABS(SC(I))/MAX(ABS(XPLS(I)),1./SX(I)))
  140     CONTINUE
C$        WRITE(IPR,962) RLN
          IF(RLN.GE.STEPTL) GO TO 150
C         IF(RLN.LT.STEPTL)
C         THEN
C
C           CANNOT FIND SATISFACTORY XPLS SUFFICIENTLY DISTINCT FROM X
C
            IRETCD=1
C$          WRITE(IPR,954)
            GO TO 230
C         ELSE
C
C           REDUCE TRUST REGION AND CONTINUE GLOBAL STEP
C
  150       IRETCD=2
            DLTMP=-SLP*DLT/(2.*(DLTF-SLP))
C$          WRITE(IPR,963) DLTMP
            IF(DLTMP.GE. .1*DLT) GO TO 155
C           IF(DLTMP.LT. .1*DLT)
C           THEN
              DLT=.1*DLT
              GO TO 160
C           ELSE
  155         DLT=DLTMP
C           ENDIF
  160       CONTINUE
C$          WRITE(IPR,955)
            GO TO 230
C         ENDIF
C       ELSE
C
C         FPLS SUFFICIENTLY SMALL
C
  170     CONTINUE
C$        WRITE(IPR,958)
          DLTFP=0.
          IF (METHOD .EQ. 3) GO TO 180
C
C         IF (METHOD .EQ. 2)
C         THEN
C
          DO 177 I = 1, N
             TEMP = 0.0
             DO 173 J = I, N
                TEMP = TEMP + (A(J, I)*SC(J))
  173        CONTINUE
             DLTFP = DLTFP + TEMP*TEMP
  177     CONTINUE
          GO TO 190
C
C         ELSE
C
  180     DO 187 I = 1, N
             DLTFP = DLTFP + UDIAG(I)*SC(I)*SC(I)
             IF (I .EQ. N) GO TO 187
             TEMP = 0
             IP1 = I + 1
             DO 183 J = IP1, N
                TEMP = TEMP + A(I, J)*SC(I)*SC(J)
  183        CONTINUE
             DLTFP = DLTFP + 2.0*TEMP
  187     CONTINUE
C
C         END IF
C
  190     DLTFP = SLP + DLTFP/2.0
C$        WRITE(IPR,964) DLTFP,NWTAKE
          IF(IRETCD.EQ.2 .OR. (ABS(DLTFP-DLTF).GT. .1*ABS(DLTF))
     +         .OR. NWTAKE .OR. (DLT.GT. .99*STEPMX)) GO TO 210
C         IF(IRETCD.NE.2 .AND. (ABS(DLTFP-DLTF) .LE. .1*ABS(DLTF))
C    +         .AND. (.NOT.NWTAKE) .AND. (DLT.LE. .99*STEPMX))
C         THEN
C
C           DOUBLE TRUST REGION AND CONTINUE GLOBAL STEP
C
            IRETCD=3
            DO 200 I=1,N
              XPLSP(I)=XPLS(I)
  200       CONTINUE
            FPLSP=FPLS
            DLT=MIN(2.*DLT,STEPMX)
C$          WRITE(IPR,959)
            GO TO 230
C         ELSE
C
C           ACCEPT XPLS AS NEXT ITERATE.  CHOOSE NEW TRUST REGION.
C
  210       CONTINUE
C$          WRITE(IPR,960)
            IRETCD=0
            IF(DLT.GT. .99*STEPMX) MXTAKE=.TRUE.
            IF(DLTF.LT. .1*DLTFP) GO TO 220
C           IF(DLTF.GE. .1*DLTFP)
C           THEN
C
C             DECREASE TRUST REGION FOR NEXT ITERATION
C
              DLT=.5*DLT
              GO TO 230
C           ELSE
C
C             CHECK WHETHER TO INCREASE TRUST REGION FOR NEXT ITERATION
C
  220         IF(DLTF.LE. .75*DLTFP) DLT=MIN(2.*DLT,STEPMX)
C           ENDIF
C         ENDIF
C       ENDIF
C     ENDIF
  230 CONTINUE
C$    WRITE(IPR,953)
C$    WRITE(IPR,956) IRETCD,MXTAKE,DLT,FPLS
C$    WRITE(IPR,957)
C$    WRITE(IPR,965) (XPLS(I),I=1,N)
      RETURN
C
  951 FORMAT(55H TREGUP    RESET XPLS TO XPLSP. TERMINATION GLOBAL STEP)
  952 FORMAT(26H TREGUP    FPLS TOO LARGE.)
  953 FORMAT(38H0TREGUP    VALUES AFTER CALL TO TREGUP)
  954 FORMAT(54H TREGUP    CANNOT FIND SATISFACTORY XPLS DISTINCT FROM,
     +       27H X.  TERMINATE GLOBAL STEP.)
  955 FORMAT(53H TREGUP    REDUCE TRUST REGION. CONTINUE GLOBAL STEP.)
  956 FORMAT(21H TREGUP       IRETCD=,I3/
     +       21H TREGUP       MXTAKE=,L1/
     +       21H TREGUP       DLT   =,E20.13/
     +       21H TREGUP       FPLS  =,E20.13)
  957 FORMAT(32H TREGUP       NEW ITERATE (XPLS))
  958 FORMAT(35H TREGUP    FPLS SUFFICIENTLY SMALL.)
  959 FORMAT(54H TREGUP    DOUBLE TRUST REGION.  CONTINUE GLOBAL STEP.)
  960 FORMAT(50H TREGUP    ACCEPT XPLS AS NEW ITERATE.  CHOOSE NEW,
     +       38H TRUST REGION.  TERMINATE GLOBAL STEP.)
  961 FORMAT(18H TREGUP    IRETCD=,I5/
     +       18H TREGUP    FPLS  =,E20.13/
     +       18H TREGUP    FPLSP =,E20.13/
     +       18H TREGUP    DLTF  =,E20.13/
     +       18H TREGUP    SLP   =,E20.13)
  962 FORMAT(18H TREGUP    RLN   =,E20.13)
  963 FORMAT(18H TREGUP    DLTMP =,E20.13)
  964 FORMAT(18H TREGUP    DLTFP =,E20.13/
     +       18H TREGUP    NWTAKE=,L1)
  965 FORMAT(14H TREGUP       ,5(E20.13,3X))
      END
