      SUBROUTINE DMLC01 (DMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X,ACC,
     *                   IPRINT, NFMAX, IW, LIW, W, LW)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c     File: DMLC.for    Minimization with linear constraints.
C>> 2001-10-05 DMLC Krogh  Minor format changes -- 1PE => 1P,E
C>> 2001-06-18 DMLC Krogh  Changed ". LT." to " .LT."
C>> 1996-07-08 DMLC Krogh  Multiple changes for C conversion.
c>> 1996-03-30 DMLC Krogh  Removed MAX from type stmt., added external.
c>> 1995-11-15 DMLC Krogh  Moved formats up for C conversion.
C>> 1994-11-11 DMLC Krogh  Declared all vars.
c>> 1994-11-02 DMLC Krogh  Changes to use M77CON
c>> 1992-04-27 DMLC CLL
c>> 1992-04-17 CLL
c>> 1992-03-27 CLL
c>> 1992-02-05 CLL Declared all floating point variables.
c>> 1992-01-16 CLL
c>> 1991-06-10 CLL & FTK Editing comments.
c>> 1991-04-30 CLL & FTK
c>> 1991-04-19 CLL
c>> 1991-04-15 FTK & CLL  Made FIRST an argument of DMLC20
c>> 1990-09-12 C. L. Lawson and F. T. Krogh, JPL
c>> 1990-07-25 C. L. Lawson, JPL
c>> 1990-07-12 C. L. Lawson, JPL
c>> 1990-04-03 C. L. Lawson, JPL
c--D replaces "?": ?MLC,?MLC01,?MLC02,?MLCFG,?MLC20,?MLC21,?MLC04,?MLC05
c--&      ?MLC13,?MLC15,?MLC03,?MLC12,?MLC06,?MLC11,?MLC09,?MLC16,?MLC19
c--&      ?MLC07,?MLC18,?MLC08,?MLC14,?MLC10,?MLC17,?ERV1
c     Also calls     ERMOR, ERMSG, IERM1, IERV1
C     ------------------------------------------------------------------
c     This algorithm and the original code is due to
c     M. J. D. Powell, Cambridge Univ., 1989.
c     1990-04-03 This version adapted for usage at JPL by
c     C. L. Lawson.  Instrumented with "c--" lines to automate type
c     conversions.
c     1990-07-11 CLL Added new argument NFMAX to DMLC01.  Added new
c     subroutine DMLC20 to compute gradient by finite differences.
c     Added argument HAVEG in the user-provided subr DMLCFG.
c     Added args XL and XU to subr _MLC09 in order to be able to pass
c     them on to DMLC20.
c     Added reference to [D/R]1MACH in _MLC15.
c     1990-07-12 CLL  Transposed indices in the array A(,).
c     The constraints are now Ax .le. b rather than (A**t)x .le. b.
c     1990-07-19 CLL  Reorganized the arg list of DMLC01.
c        Prev args INFO, NACT, IACT() are embedded in new arg IW().
c        New item, FVAL, and previous args W() and PAR() are combined
c        into new arg W().  The new item FVAL is the final value of the
c        objective function.
c     1990-07-24 CLL  Setting INFO = 0 initially.  It keeps this value
c        till it is changed for some reason.  Previously it was set to
c        4 initially.
c        Changed specification of IPRINT.  It now contains 5 print
c        parameters.  This gives more flexibility, such as requesting
c        printing only at termination.
c        Moved the main block of intermediate printing into a new subr,
c        DMLC21.  This allows cleaner logic for testing when to print.
c        DMLC21 is called both from DMLC02 and from _MLC05.  In revising
c        the print control we removed variables NFVALK and ITERP, and
c        added variables NFVPR and TOLPRT.
c        Printing of error messages is now done using JPL MATH77 error
c        printing subroutines.  Requires linking of
c        DERV1, ERFIN, ERMOR, ERMSG, IERM1, and IERV1.  The calls are in
c        DMLC01, DMLC02, and _MLC15.  ERFIN is not called directly.
c        Setting INFO = 9 in DMLC02 when terminate due to solution being
c        determined just by the constraints.  Also compute and return
c        the func value in this case.
c     1990-09-12 CLL and FTK Added arg. W20 for _MLC20.
c     1991-04-19 CLL Change to have NFVALS count every call to _MLCFG.
c        Move test of IPFREQ to be sure last iteration gets printed if
c        it is a multiple of IPFREQ.  Changed usage of NFVPR.
c 1992 March/April CLL Additions to the code in MINFUN (C05) and
c    LSRCH (C09).  The code will now return INFO = 2 in some cases
c    in which it previously returned INFO = 3.  It will also finish in
c    fewer iterations in some cases.  Added more to error msg on quit-
c    ting with INFO = 3.  Added HAVEG into arg list of _MLC05.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c--   Begin Mask code changes.
c     Names of subroutines have been changed as follows:
c     Original names          New DP names         New SP names
c     GETMIN                  DMLC01               SMLC01
c     MINFLC                  DMLC02               SMLC02
c     EQCONS                  DMLC03               SMLC03
c     GETFES                  DMLC04               SMLC04
c     MINFUN                  DMLC05               SMLC05
c     CONRES                  DMLC06               SMLC06
c     GETD                    DMLC07               SMLC07
c     SDEGEN                  DMLC08               SMLC08
c     LSRCH                   DMLC09               SMLC09
c     NEWCON                  DMLC10               SMLC10
c     SATACT                  DMLC11               SMLC11
c     ADDCON                  DMLC12               SMLC12
c     ADJTOL                  DMLC13               SMLC13
c     DELCON                  DMLC14               SMLC14
c     INITZU                  DMLC15               SMLC15
c     KTVEC                   DMLC16               SMLC16
c     SDIRN                   DMLC17               SMLC17
c     STEPBD                  DMLC18               SMLC18
c     ZBFGS                   DMLC19               SMLC19
c     FGCALC                  DMLCFG               SMLCFG
c     New code by FTK -->     DMLC20               SMLC20
c     New code by CLL -->     DMLC21               SMLC21
c--   End Mask code changes.

C
C  This is main entry point to a package of subroutines that calculate
C  the least value of a differentiable function of several variables
C  subject to linear constraints on the values of the variables.
c     ------------------------------------------------------------------
c                  Subroutine Arguments
c
C  DMLCFG [in]  Name of subroutine having an interface of the form
C                SUBROUTINE DMLCFG (N, X, F, G, HAVEG)
c                INTEGER N
c                [DOUBLE PRECISION or REAL] X(N), F, G(N)
C                LOGICAL HAVEG
c     provided by the user to compute values of the objective function,
c     F, and optionally its gradient vector, G(), at the argument value
c     given by X().  The user may choose either to write code to compute
c     G() or not.  If DMLCFG computes G() it must also set HAVEG = .true
c     If DMLCFG does not compute G() it must set HAVEG = .false.
c     In this latter case the package will make additional calls to
c     DMLCFG to compute a finite difference approximation to the
c     gradient vector.
c     The vectors X() and G() will be of length N.  The subroutine must
c     not change the values of N or X().
c     DMLCFG may use COMMON block(s) if neccessary to communicate
c     data between the user's main program and DMLCFG.
C  N is the number of variables and must be set by the user.
C  M is the number of linear constraints (excluding simple bounds) and
C     must be set by the user.
C  MEQ is the number of constraints that are equalities and must be set
C     by the user.
C  A(.,.) is a 2-dimensional array whose rows are the gradients of
C     the M constraint functions.  Its entries must be set by the user
C     and its dimensions must be at least M by N.
C  LDA is the actual first dimension of the array A that is supplied by
C     the user, so its value may not be less than M.
C  B(.) is a vector of constraint right hand sides that must also be set
C     by the user.  Specifically the constraints on the variables X(I)
C     I=1(1)N are
C          A(K,1)*X(1)+...+A(K,N)*X(N) .EQ. B(K)  K=1,...,MEQ
C          A(K,1)*X(1)+...+A(K,N)*X(N) .LE. B(K)  K=MEQ+1,...,M  .
C     Note that the data that define the equality constraints come
C     before the data of the inequalities.
C  XL(.) and XU(.) are vectors whose components must be set to lower and
C     upper bounds on the variables.  Choose very large negative and
C     positive entries if a component should be unconstrained, or set
C     XL(I)=XU(I) to freeze the I-th variable.  Specifically these
C     simple bounds are
C          XL(I) .LE. X(I) and X(I) .LE. XU(I)  I=1,...,N  .
C  X(.) is the vector of variables of the optimization calculation.  Its
C     initial elements must be set by the user to an estimate of the
C     required solution.  The subroutines can usually cope with poor
C     estimates, and there is no need for X(.) to be feasible initially.
C     These variables are adjusted automatically and the values that
C     give the least feasible calculated value of the objective function
C     are available in X(.) on the return from DMLC01.
C  ACC is a tolerance on the first order conditions at the calculated
C     solution of the optimization problem.  These first order
C     conditions state that, if X(.) is a solution, then there is a set
C     of active constraints with indices IACT(K) K=1(1)NACT, say, such
C     that X(.) is on the boundaries of these constraints, and the
C     gradient of the objective function can be expressed in the form
C          GRAD(F)=PAR(1)*GRAD(C(IACT(1)))+...
C                        ...+PAR(NACT)*GRAD(C(IACT(NACT)))  .
C     Here PAR(K) K=1(1)NACT are Lagrange multipliers that are
C     nonpositive
C     for inequality constraints, and GRAD(C(IACT(K))) is the gradient
C     of the IACT(K)-th constraint function, so it is A(IACT(K),.) if
C     IACT(K) .LE. M, and it is minus or plus the J-th coordinate vector
c     if the constraint is the lower or upper bound on X(J)
c     respectively.  The normal return from the calculation occurs when
c     X(.) is feasible and the sum of squares of components of the
c     vector RESKT(.) is at most ACC**2, where RESKT(.) is the
c     N-component vector of residuals of the first order condition that
c     is displayed above.  Sometimes the package cannot satisfy this
c     condition, because noise in the function values can prevent a
c     change to the variables, no line search being allowed to increase
c     the objective function.
C  IPRINT  [in, integer]  Contains print control parameters, packed as
c     follows:
c     IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST
c     Printing only begins after a first feasible point is found.
c     The items to be printed are selected by IPMORE.  The other
c     parameters determine when to print.
c     IPFREQ is zero or positive.  If positive, printing will be done on
c        iterations 0, IPFREQ, 2*IPFREQ, etc.
c     IPTOL = 0 or 1.  1 means to print the new TOL value
c        and the standard items each time TOL is changed.
c     IPFRST = 0 or 1.  1 means to print on the first iteration, i.e. on
c        iteration No. 0.
c     IPMORE = 0 or 1.  0 means the items to be printed are ITERC,
c        NFVALS, F, X(1:N), and G(1:N).  1 means to also print
c        IACT(1:NAXT), PAR(1:NACT) and RESKT(1:N).
c     IPLAST = 0 or 1.  1 means to print the final results, and the
c        reason for termination.
c  NFMAX  [in, integer]  If positive this sets an upper limit on the
c     number of function evaluations.  When gradients are estimated by
c     differences, the count includes the function evaluations for this
c     purpose also.  If zero there is no upper limit.
c  IW(.)  [work, out, integer]  Integer work space.  Also used to return
c     status information.  Its length must be at least LIW.
c     On return:
c        IW(1) contains INFO.  Indicates reason for termination.
c              See below for explanation of values.
c        IW(2) contains ITERC.  No. of iterations.
c        IW(3) contains NFVALS.  No. of function evaluations, not
c              counting extra func evals done to estimate the gradient
c              when the gradient is not computed explicitly.  In this
c              latter case the actual No. of func evals will be
c              (K+1)*NFVALS, where K is the number of solution
c              components whose lower and upper bounds are not equal.
c        IW(4) contains NACT.  The number of active constraints at the
c              solution point.  Will be in the interval [0, N].
c        IW(5:4+NACT) contains IACT(1:NACT).  IACT() is used as work
c              space of length M + 2*N.  On return the first NACT
c              locations contain the indices of the constraints that
c              are active at the final point.  Indices [1:M] refer to
c              rows of the system Ax .le. b.  Indices [M+1:M+N]
c              refer to component lower bounds.  Indices [M+N+1:M+2*N]
c              refer to component upper bounds.
C
c  LIW  [in, integer]  Dimension of the IW() array.  Require
c       LIW .ge. 4 + M + 2*N.
c
C  W(.)  [work, out, float]  Working space array of float variables.
c     Also used to return results.
c     Its length must be at least LW.
c     On return:
c     W(1) contains FVAL, the final value of the objective function.
c     W(2:1+N) contains the final gradient vector, GRAD(1:N).
c     W(2+N:1+2*N) contains the Kuhn-Tucker residual vector, RESKT(1:N).
c     W(2+2*N:1+2*N+NACT) contains the Lagrange multipliers, PAR(1:NACT)
c        where NACT has a value in [0,N].  GRAD(), RESKT(), and PAR()
c        are mentioned above in the description of ACC.
c  LW  [in, integer]  Dimension of the W() array.  Require
c      LW .ge. 3 + M + 16*N + N**2.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c                      Return values of INFO
c
c     INFO indicates the reason for termination.
c
C          INFO=1   X(.) is feasible and the condition that depends on
C          ACC is satisfied.
C          INFO=2   X(.) is feasible and rounding errors are preventing
C          further progress.
C          INFO=3   X(.) is feasible but the objective function fails to
C     decrease although a decrease is predicted by the current gradient
C     vector.  If this return occurs and RESKT(.) has large components
C     then the user's calculation of the gradient of the objective
C     function may be incorrect.  One should also question the coding of
C     the gradient when the final rate of convergence is slow.
C          INFO=4   In this case the calculation cannot begin because IA
C     is less than N or because the lower bound on a variable is greater
C     than the upper bound.
C          INFO=5   This value indicates that the equality constraints
C     are inconsistent.   These constraints include any components of
C     X(.) that are frozen by setting XL(I)=XU(I).
C          INFO=6   In this case there is an error return because the
C     equality constraints and the bounds on the variables are found to
C     be inconsistent.
C          INFO=7   This value indicates that there is no vector of
C     variables that satisfies all of the constraints.  Specifically,
C     when this return or an INFO=6 return occurs, the current active
C     constraints (whose indices are IACT(K) K=1(1)NACT) prevent any
C     change in X(.) that reduces the sum of constraint violations,
C     where only bounds are included in this sum if INFO=6.
C          INFO=8   In this case the limit on the number of calls of
C     subroutine DMLCFG (see below) has been reached, and there would
C     have been further calculation otherwise.
c          INFO = 9  The solution is determined just by the constraints.
c     In this case PAR() and RESKT() are not defined on return.
c     FVAL will be defined on return.  G() will only be defined on
c     return if DMLCFG returns HAVEG = .true.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  Note 1.   The variables N, M, MEQ, LDA, ACC, IPRINT, NFMAX and the
C     elements of the arrays A(.,.), B(.), XL(.) and XU(.) are not
c     altered by the optimization procedure.  Their values, and the
c     initial components of X(.) must be set on entry to DMLC01.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C  Note 2.   A paper on the method of calculation and a report on the
C     main features of the computer code are available from the author
C     M.J.D.Powell (D.A.M.T.P., University of Cambridge, Silver Street,
C     Cambridge CB3 9EW, England).
c     ------------------------------------------------------------------
      external DMLCFG
      integer LDA, LIW, LW
      integer I, IBRES, ID
      integer IFVAL, IG, IGM, IGS, IIACT, IINFO, IITERC, INACT, INFO
      integer INFVAL, IPAR, IPRINT, IRESKT, IU, IW(LIW), IW20, IXBIG
      integer IXS, IZ, IZTG, KTNORM, M, MEQ, N, NFMAX
      double precision A(LDA,*),ACC, B(*),TEMP, XL(*),XU(*),X(*),W(LW)
      parameter(IFVAL = 1, KTNORM = 2, IG = 3)
      parameter(IINFO = 1, IITERC = 2, INFVAL = 3, INACT = 4, IIACT = 5)
C     ------------------------------------------------------------------
c                                   Check sizes of LIW and LW
      INFO = 0
      if(LIW .lt. 4 + M + 2*N) then
         INFO = 4
         call IERM1('DMLC01', INFO, 0,'Dimension LIW is too small.',
     *        'LIW',LIW, ',')
         call IERV1('Need',4 + M + 2*N,'.')
      endif
c
      if(LW .lt. 3 + M + N*(16 + N)) then
         INFO = 4
         call IERM1('DMLC01', INFO, 0,'Dimension LW is too small.',
     *        'LW',LW, ',')
         call IERV1('Need',3 + M + N*(16 + N),'.')
      endif
      if(INFO .ne. 0) then
         IW(1) = INFO
         return
      endif
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C                         Partition the array W()
c                         parameter(IFVAL = 1, KTNORM = 2, IG = 3)
      IRESKT=IG+N
      IPAR = IRESKT + N
      IZ=IPAR + N
      IU=IZ+N*N
      IXBIG=IU+N
      IBRES=IXBIG+N
      ID=IBRES+M+N+N
      IZTG=ID+N
      IGM=IZTG+N
      IXS=IGM+N
      IGS=IXS+N
C                  Need 4*N+1 locations beginning at IW20.
      IW20=IGS+N
C                         Partition the array IW()
c     parameter(IINFO = 1, IITERC = 2, INFVAL = 3, INACT = 4, IIACT = 5)
C
c     Zero the RESKT vector so it will be safe to compute its norm on
c     return from DMLC02 even in cases when DMLC02 does not compute
c     RESKT.
c
      do 10 I = IRESKT, IRESKT + N -1
         W(I) = 0.0D0
   10 continue
c
C     Call the optimization package.
C
      CALL DMLC02 (DMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X,ACC,
     *  IW(IIACT),IW(INACT),W(IPAR),IPRINT, NFMAX,
     *  IW(IINFO), IW(IITERC), IW(INFVAL),
     *  W(IG),W(IZ),W(IU),W(IXBIG),W(IRESKT),W(KTNORM),W(IBRES),W(ID),
     *  W(IZTG),W(IGM),W(IXS),W(IGS), W(IFVAL), W(IW20))
c
      TEMP = 0.0D0
      do 20 I = IRESKT, IRESKT+N-1
         TEMP = TEMP + W(I)**2
   20 continue
      W(KTNORM) = sqrt(TEMP)
      RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC02 (DMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X,
     *  ACC,IACT,NACT,PAR, IPRINT, NFMAX, INFO, ITERC, NFVALS,
     *  G,Z,U,XBIG,RESKT, ENRMKT, BRES,D,ZTG,GM,XS,GS, FVAL, W20)
c
c     Main subroutine to control the Powell constrained minimization
c     package.
c     1990-07-19 CLL  Added NFMAX, ITERC, NFVALS, and FVAL to arg list.
c        Also added FVAL to arg list of DMLC05 that is called
c        from this subr.
c        NFMAX provides limit on no. of func evaluations, but zero means
c        no limit.  This was previously overloaded onto INFO.
c        ITERC & NFVALS give output of counts of iterations and
c        function evaluations.
c        FVAL provides output of the final objective function value.
c     ------------------------------------------------------------------
c                                 Arguments
c
c  DMLCFG, N,M,MEQ,A(),LDA,B(),XL(),XU(),X(),
c  ACC,IACT(),NACT,PAR(),
c  IPRINT [in, integer]  Contains print parameters packed as follows:
c           IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST
c  NFMAX [in, integer]  Limit on No. of function evaluations.  But zero
c         means no limit.
c  INFO [out, integer]  Termination status flag.
c  ITERC [out, integer]  Count of No. of iterations.  Initialized to 0
c        in this subr.  Incremented in DMLC05.
c  NFVALS [out, integer]  Count of number of function evaluations.
c  G(),Z(),U(),XBIG(),
c  RESKT()
c  ENRMKT  Euclidean norm of RESKT().
c  BRES(),D(),
c  ZTG(),GM(),XS(),GS(),
c  FVAL [out, float]  Final (best) value of F.
c  W20  [out, work, float] Working storage for DMLC20.
c     ------------------------------------------------------------------
c                    Descriptions of some of the internal variables.
c
c  GMODE [integer]  Initialized to 0 here.  Sent to DMLC05 for eventual
c        use by DMLC20.
c  HAVEG [logical]  Arg in call to DMLCFG, not used in this subroutine.
c  NFVPR [integer]  Used to avoid printing results from the same
c        function evaluation more than once.  Initialized to -1 in this
c        subr.  Modified in DMLC05.  Used in DMLC02 and DNLC05.
c     ------------------------------------------------------------------
      external DMLCFG

      integer GMODE, I, IACT(*), INFO
      integer IPFREQ, IPFRST, IPLAST, IPMORE, IPRINT
      integer IPTOL, ITERC, K, LDA, M, MEQ, MEQL, MP, MSAT
      integer MTOT, N, NACT, NFMAX, NFVALS, NFVPR

      double precision A(LDA,*),ACC, B(*),BRES(*),D(*),ENRMKT, FVAL
      double precision G(*), GM(*), GS(*), PAR(*), RELACC, RESKT(*)
      double precision SSQKT, TOL, U(*)
      double precision W20(0:*), XBIG(*), XL(*), XU(*), X(*), XS(*)
      double precision Z(*), ZTG(*), ZZNORM
      logical HAVEG
c     ------------------------------------------------------------------
C     Initialize ZZNORM, ITERC, NFVPR, NFVALS, INFO
C
      ZZNORM=-1.0D0
      W20(0) = -1.0d0
      W20(1) = -1.0d0
      GMODE = 0
      ITERC=0
      NFVPR = -1
      NFVALS=0
      INFO = 1
c
c                         Decompose IPRINT into 5 integers:
c           IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST
C
      IPFREQ = IPRINT/100
      I = IPRINT - IPFREQ * 100
      IPTOL = I/8
      I = I - IPTOL * 8
      IPFRST = I/4
      I = I - IPFRST * 4
      IPMORE = I/2
      IPLAST = I - IPMORE * 2
c
C                      Check the bounds on N, M and MEQ.
C
      if(N .lt. 1 .or. M .lt. 0 .or. MEQ .lt. 0 .or. M .lt. MEQ) then
         INFO = 4
         call IERM1('DMLC02', INFO, 0,'Bad value for M, N, or MEQ.',
     *        'M',M, ',')
         call IERV1('N',N,',')
         call IERV1('MEQ',MEQ,'.')
         GOTO 40
      END IF
C
C                      Initialize RELACC, Z, U and TOL.
C
      CALL DMLC15 (N,M,XL,XU,X,IACT,MEQL,INFO,Z,U,XBIG,RELACC)
      IF (INFO .EQ. 4) go to 40
      TOL=max(0.01D0,10.0D0*RELACC)
      if(IPFREQ .ne. 0 .or. IPFRST .ne. 0 .or.
     *   IPLAST .ne. 0 .or. IPTOL .ne. 0) then
         print'(/6x,''Beginning subroutine DMLC02, called from DMLC01:''
     *      /6x,''N ='',i4,'',  M ='',i4,'',  MEQ  ='',i4,'',  ACC ='',
     *      g11.3/6x,''RELACC ='',g11.3,'',  TOL = '',g11.3)',
     *      N, M, MEQ, ACC,  RELACC, TOL
      endif
C
C     Add any equality constraints to the active set.
C
      IF (MEQ .GT. 0) THEN
          CALL DMLC03 (N,M,MEQ,A,LDA,B,XU,IACT,MEQL,INFO,Z,U,RELACC,XS,
     1      GS)
          IF (INFO .EQ. 5) THEN
            call ERMSG('DMLC02',INFO,0,
     *           'Equality constraints are inconsistent.','.')
            GOTO 40
          END IF
      END IF
      NACT=MEQL
      MSAT=MEQL
C
C     Add the bounds to the list of constraints.
C
      MTOT=NACT
      DO 10 I=1,N
      IF (XL(I) .LT. XU(I)) THEN
          MTOT=MTOT+2
          IACT(MTOT-1)=M+I
          IACT(MTOT)=M+N+I
      END IF
   10 CONTINUE
C
C     Try to satisfy the bound constraints.
C
      CALL DMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z,U,XBIG,
     1  RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,RESKT,XS,GS,W20)
      IF (MSAT .LT. MTOT) THEN
         INFO=6
         call ERMSG('DMLC02',INFO,0,
     *           'Equalities and bounds are inconsistent.','.')
         GOTO 40
      END IF
C
C     Add the ordinary inequalities to the list of constraints.
C
      IF (M .GT. MEQ) THEN
          MP=MEQ+1
          DO 20 K=MP,M
          MTOT=MTOT+1
   20     IACT(MTOT)=K
      END IF
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   30 continue
c                                    Begin main loop
C
C                          Correct any constraint violations.
C
      CALL DMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z,U,XBIG,
     1  RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,RESKT,XS,GS,W20)
      IF (MSAT .LT. MTOT) THEN
         INFO=7
         call ERMSG('DMLC02',INFO,0,
     *           'Constraints are inconsistent.','.')
         GOTO 40
      ELSE IF (MEQL .EQ. N) THEN
c
c        Here the soln is determined just by the constraints.
c        PAR() and RESKT() are not defined on return in this case.
c        FVAL has not yet been evaluated.  We evaluate FVAL here in
c        order to return its value.  We leave G() not computed if HAVEG
c        is false, however it would be easy to compute G by calling
c        DMLC20 if that seemed preferable.
c
C%%        (*dmlcfg)( n, x, fval, g, &haveg );
         call DMLCFG(N, X, FVAL, G, HAVEG)
         NFVALS = NFVALS+1
         ENRMKT = 0.0d0
         INFO = 9

         if(IPFREQ .ne. 0 .or. IPFRST .ne. 0 .or.
     *      IPLAST .ne. 0 .or. IPTOL .ne. 0) then
            print '(/'' [1] DMLC02..''/
     *    '' The solution is determined by the equality constraints.'')'
            call DMLC21(IPMORE, .false., N, ITERC, NFVALS,
     *              FVAL, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT)
         endif
         GOTO 40
      END IF
C
C     Minimize the objective function in the case when constraints are
C       treated as degenerate if their residuals are less than TOL.
C
c          DMLC05 is MINFUN
      CALL DMLC05 (DMLCFG, N,M,A,LDA,B,XL,XU,X,ACC,IACT,NACT,
     *  PAR, IPFREQ, IPTOL, IPFRST, IPMORE, INFO, G, Z,
     *  U,XBIG,RELACC,ZZNORM,TOL,MEQL,MTOT,ITERC, NFVPR, NFVALS,
     *  NFMAX,RESKT, SSQKT, BRES,D,ZTG,GM,XS,GS, FVAL, GMODE,
     *  HAVEG, W20)
C
C     Reduce TOL if necessary.
C
      IF (TOL .GT. RELACC .AND. NACT .GT. 0) THEN
          if (NFMAX .le. 0 .or. NFVALS .lt. NFMAX) then
C                  DMLC13 is ADJTOL
              CALL DMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC,
     *                     TOL, MEQL)
              GOTO 30
          else
              INFO=8
          endif
      END IF
c
c           Here when INFO = 1, 2, 3, or 8.
c           We treat 1 and 2 as normal conditions and 3 and 8 as errors.
c
      if (INFO .EQ. 1) then
         if (IPLAST .ne. 0) print*,
     *       '[1]  DMLC01 quitting because requested accuracy',
     *       ' has been attained.'

      elseif (INFO .EQ. 2) then
         if (IPLAST .ne. 0) print*,
     *       '[2]  DMLC01 quitting due to limitation of',
     *       ' computational precision.'

      elseif (INFO .EQ. 3) then
         call ERMSG('DMLC02', INFO, 0,
     *   'F not decreasing, although gradient indicates it should.',',')
         if(HAVEG) then
            call ERMOR('Could be due to limitation of precision',',')
            call ERMOR('or to incorrect code for the gradient.','.')
         else
            call ERMOR('Could be due to limitation of precision.','.')
         endif

      elseif (INFO .EQ. 8) then
         call IERM1('DMLC02', INFO, 0,
     *   'Reached specified max number of function evaluations.',
     *   'Count',NFVALS, ',')
         call IERV1('NFMAX',NFMAX,'.')
      endif
c
c                           Test to print results before quitting.
c
      ENRMKT = sqrt(SSQKT)
      if(IPLAST .ne. 0 .and. NFVALS .ne. NFVPR) then
c        print*,'DMLC02.. Debug.. Call C21 on leaving C02.'
         call DMLC21(IPMORE, .true., N, ITERC, NFVALS,
     *           FVAL, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT)
      endif
c
   40 return
      end
c     ==================================================================
      SUBROUTINE DMLC03 (N,M,MEQ,A,LDA,B,XU,IACT,MEQL,INFO,Z,U,RELACC,
     1  AM,CGRAD)
C
C     Try to add the next equality constraint to the active set.
c     ------------------------------------------------------------------
      integer N, M, MEQ, LDA, IACT(*), MEQL, INFO
      integer I, IZ, J, JM, K, KEQ, NP
      double precision A(LDA,*),B(*),XU(*)
      double precision RELACC, RHS, SUM, SUMABS, VMULT
      double precision Z(*),U(*),AM(*),CGRAD(*)
c     ------------------------------------------------------------------
      DO 50 KEQ=1,MEQ
      IF (MEQL .LT. N) THEN
          NP=MEQL+1
          IACT(NP)=KEQ
          CALL DMLC12 (N,M,A,LDA,IACT,MEQL,Z,U,RELACC,NP,AM,CGRAD)
          IF (MEQL .EQ. NP) GOTO 50
      END IF
C
C     If linear dependence occurs then find the multipliers of the
C       dependence relation and apply them to the right hand sides.
C
      SUM=B(KEQ)
      SUMABS=abs(B(KEQ))
      IF (MEQL .GT. 0) THEN
          DO 10 I=1,N
   10     AM(I)=A(KEQ,I)
          K=MEQL
   20     VMULT=0.0D0
          IZ=K
          DO 30 I=1,N
          VMULT=VMULT+Z(IZ)*AM(I)
   30     IZ=IZ+N
          VMULT=VMULT*U(K)
          J=IACT(K)
          IF (J .LE. M) THEN
              DO 40 I=1,N
   40         AM(I)=AM(I)-VMULT*A(J,I)
              RHS=B(J)
          ELSE
              JM=J-M-N
              AM(JM)=AM(JM)-VMULT
              RHS=XU(JM)
          END IF
          SUM=SUM-RHS*VMULT
          SUMABS=SUMABS+abs(RHS*VMULT)
          K=K-1
          IF (K .GE. 1) GOTO 20
      END IF
C
C     Error return if the constraints are inconsistent.
C
      IF (abs(SUM) .GT. RELACC*SUMABS) THEN
          INFO=5
          GOTO 60
      END IF
   50 CONTINUE
   60 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z,
     1  U,XBIG,RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,GMNEW,PARNEW,
     2  CGRAD,W20)
C
C     Make the correction to X for the active constraints.
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, INFO, MEQL, MSAT, MTOT
      integer I, ITEST, MSATK, INDXBD
      double precision A(LDA,*),B(*),XL(*),XU(*),X(*),PAR(*),G(*)
      double precision U(*),XBIG(*),BRES(*),D(*)
      double precision GM(*),GMNEW(*),PARNEW(*), CGRAD(*)
      double precision W20(0:*), RELACC, STEPCB, SUMRES, SUMRSK, TOL
      double precision Z(*), ZTG(*)
c     ------------------------------------------------------------------
      INFO=0
   10 CALL DMLC11 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,INFO,Z,U,XBIG,RELACC,
     1  TOL,MEQL)
      IF (INFO .GT. 0) MSAT=NACT
      IF (MSAT .EQ. MTOT) GOTO 60
C
C     Try to correct the infeasibility.
C
   20 MSATK=MSAT
      SUMRSK=0.0D0
   30 CALL DMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG,BRES,
     1  D,ZTG,RELACC,TOL,STEPCB,SUMRES,MEQL,MSAT,MTOT,INDXBD,GM,GMNEW,
     2  PARNEW,CGRAD,W20)
C
C     Include the new constraint in the active set.
C
      IF (STEPCB .GT. 0.0D0) THEN
          DO 40 I=1,N
          X(I)=X(I)+STEPCB*D(I)
   40     XBIG(I)=max(XBIG(I),abs(X(I)))
          CALL DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,GMNEW,
     *                 CGRAD)
      END IF
C
C     Test whether to continue the search for feasibility.
C
      IF (MSAT .LT. MTOT) THEN
          IF (STEPCB .EQ. 0.0D0) GOTO 50
          IF (MSATK .LT. MSAT) GOTO 20
          IF (SUMRSK .EQ. 0.0D0 .OR. SUMRES .LT. SUMRSK) THEN
              SUMRSK=SUMRES
              ITEST=0
          END IF
          ITEST=ITEST+1
          IF (ITEST .LE. 2) GOTO 30
C
C     Reduce TOL if it may be too large to allow feasibility.
C
   50     IF (TOL .GT. RELACC) THEN
              CALL DMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC,
     1          TOL,MEQL)
              GOTO 10
          END IF
      END IF
   60 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC05 (DMLCFG, N,M,A,LDA,B,XL,XU,X,ACC,IACT,NACT,
     *  PAR, IPFREQ, IPTOL, IPFRST, IPMORE,
     *  INFO,G,Z,U,XBIG,RELACC,ZZNORM,TOL,MEQL,MTOT,ITERC,NFVPR,NFVALS,
     *  NFMAX,RESKT, SSQKT, BRES,D,ZTG,GM,XS,GS, FVAL, GMODE,
     *  HAVEG, W20)
C
C     MINFUN: Solve the minimization problem using the current value of
c             TOL.
c     1990-07-19 CLL Added FVAL to argument list to return final value
c     of F.
c     1990-07-23 CLL Removed IPRINT from arg list and added IPFREQ,
c        IPTOL, IPFRST, IPMORE, and NFVPR.  NFVPR keeps a record of
c        last function evaluation which printing has been done.
c     1990-07-23 CLL Added logical variable TOLPRT.
c     ------------------------------------------------------------------
c                     Arguments
c
c  DMLCFG, N,M,A(),LDA,B(),XL(),XU(),X(),ACC,IACT(),NACT,
c  PAR()
c  IPFREQ [in, integer]  Zero or positive.  If positive, printing of
c        intermediate results using DMLC21 will be done on iterations
c        0, IPFREQ, 2*IPFREQ, etc.
c  IPTOL [in, integer]  = 0 or 1.  1 means to print the new TOL value
c        and items printed by DMLC21 each time this subr is entered.
c  IPFRST [in, integer]  = 0 or 1.  1 means to print using DMLC21 on
c        the first iteration, i.e. iteration No. 0.
c  IPMORE [in, integer]  = 0 or 1.  Will be passed to DMLC21 to cause
c        printing to include more items when = 1 than when = 0.
c  INFO [inout, integer]  Status flag.
c  G(),Z(),U(),XBIG(),RELACC,ZZNORM,TOL,MEQL,MTOT,
c  NFVALS [inout, integer]  Count of calls for function evaluation.
c        Only counts calls to DMLCFG when the algorithm directly needs a
c        function value, not calls to DMLCFG that occur in DMLC20 for
c        estimation of the gradient.
c  NFMAX [in, integer]  Limit on No. of function evaluations, but zero
c        means no limit.
c  RESKT()
c  SSQKT    Sum of Squares of elts of RESKT().
c  BRES(),D(),ZTG(),GM(),XS(),GS(),
c  FVAL [out, float]  On return set to current (best) value of F.
c  GMODE [inout, integer]  For use by DMLC20.
c  HAVEG [out, logical]  Primarily an internal variable in this subr,
c     but returned so calling subr can use it to select proper error msg
c     when INFO = 3.
c     Set true by DMLCFG if DMLCFG provides the gradient
c     vector and false otherwise.  If false, this subr calls DMLC20 to
c     compute a finite difference approximation for the gradient.
c  W20  [out, work, float] Working storage for DMLC20.
c     ------------------------------------------------------------------
c                     Description of some of the internal variables.
c
c  DOPRNT [logical]  When IPTOL = 1 this subr prints TOL on entry and
c     sets DOPRNT = true.  Also if IPFRST = 1 and ITERC = 0 this subr
c     sets DOPRNT = true.
c     Having DOPRNT = true causes subsequent printing of intermediate
c     results for the current iteration, after PAR() and RESKT() have
c     been computed.
c  HAVEG [logical]  Set true by DMLCFG if DMLCFG provides the gradient
c     vector and false otherwise.  If false, this subr calls DMLC20 to
c     compute a finite difference approximation for the gradient.
c     ------------------------------------------------------------------
      external D1MACH, DMLCFG
      integer N, M, LDA, IACT(*), NACT, IPFREQ, IPTOL, IPFRST, IPMORE,
     *  INFO, MEQL, MTOT, ITERC, NFVPR, NFVALS, NFMAX, GMODE
      integer I, ITERK, K, MSAT, INDXBD
      double precision D1MACH
      double precision A(LDA,N), ACC, B(*), BRES(*), DDOTG, DIFF
      double precision D(*), EPS, F, FPREV, FVAL, G(N), GM(*), GS(*)
      double precision PAR(*), RELACC, RELAXF, RESKT(*)
      double precision SSQKT, STEP, STEPCB, SUM, TOL
      double precision U(*), W20(0:*)
      double precision X(*), XBIG(*), XL(N), XS(*), XU(N)
      double precision Z(*), ZTG(*), ZZNORM
      logical DOPRNT, HAVEG
      save EPS, F
      data EPS / 0.0d0 /
c     ------------------------------------------------------------------
c                  Initialize the minimization calculation.
c
*     print'(/''DMLC05.. Debug.. On entry, INFO = '',i6,'',  ITERC ='',
*        i6)', INFO, ITERC
      if(EPS .eq. 0.0d0) then
         EPS = D1MACH(4)
      endif
      MSAT=MTOT
      ITERK=ITERC
      IF (NFVALS .EQ. 0 .OR. INFO .EQ. 1) THEN
c++ CODE for ~.C. is active
          CALL DMLCFG (N,X,F,G, HAVEG)
c++ CODE for .C. is inactive
C%%        (*dmlcfg)( n, x, &f, g, haveg );
C          HAVEG = HAVEG
c++ END
*         print'('' DMLC05 798.. X()='',2g23.15/(19x,2g23.15))',
*    *      (X(I),I=1,N)
*         print'('' .. F ='',g23.15)', F
*         print'('' ............ G()='',2g23.15/(19x,2g23.15))',
*    *      (G(I),I=1,N)

          NFVALS=NFVALS+1
          if(.not. HAVEG) CALL DMLC20(DMLCFG,N,X,F,G,XL,XU,
     *                         GMODE, NFVALS, W20)
      END IF
      FPREV=abs(F+F+1.0D0)
c                         Setting DIFF here is not needed for the
c                         algorithm, but is convenient when doing
c                         debug printing.  CLL 3/26/92
      DIFF = 1.0d0
      IF (IPTOL .ne. 0) THEN
          print'(/'' New value of TOL = '',g13.5)', TOL
          DOPRNT = .true.
      elseif(IPFRST .ne. 0 .and. ITERC .eq. 0) then
          DOPRNT = .true.
      else
          DOPRNT = .false.
      END IF
C
C     Calculate the next search direction.
C          DMLC06 is CONRES
   10 CALL DMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG,BRES,D,
     1  ZTG,RELACC,TOL,STEPCB,DDOTG,MEQL,MSAT,MTOT,INDXBD,GM,RESKT,XS,
     2  GS,W20)
C
C     Calculate the Kuhn Tucker residual vector.
C          DMLC16 is KTVEC
      CALL DMLC16 (N,M,A,LDA,IACT,NACT,PAR,G,RESKT,Z,U,BRES,RELAXF,MEQL,
     1  SSQKT,XS,GS)
c
c                     Test for printing results every IPFREQ iterations.
c
      if(IPFREQ .ne. 0) then
         if( mod(ITERC, IPFREQ) .eq. 0) DOPRNT = .true.
      endif
c
      if(DOPRNT) then
         call DMLC21(IPMORE, .true., N, ITERC, NFVALS,
     *           F, TOL, X, G, NACT, IACT, PAR, RESKT, sqrt(SSQKT))
         NFVPR = NFVALS
         DOPRNT = .false.
      endif
C
C     Test for convergence.
C
*     print'(/'' DMLC05 849..''/'' DMLC05 849..'',3g13.5/
*    *   ''             F, FPREV, FPREV-F='',3g13.5/
*    *   ''                          DIFF='',g13.5/
*    *   ''             TOL, RELACC, NACT='',2g13.5,i5)',
*    *   ACC**2, SSQKT, DDOTG, F, FPREV, FPREV-F, DIFF,
*    *   TOL, RELACC, NACT
      IF (SSQKT .LE. ACC*ACC) THEN
          INFO=1
          GOTO 70
      END IF
      IF (DDOTG .GE. 0.0D0) THEN
          INFO=2
          GOTO 70
      END IF
C
C     Test for termination due to no decrease in F.
C
      IF (F .GE. FPREV) THEN
          IF (TOL .EQ. RELACC .OR. NACT .EQ. 0) THEN
              IF (DIFF .GT. 0.0D0) GOTO 20
          END IF
*         print'(/'' DMLC05 869.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F
          if(abs(DDOTG) .le. EPS * abs(F)) then
             INFO = 2
          else
             INFO=3
          endif
          GOTO 70
      END IF
   20 continue
      DIFF=FPREV-F
      FPREV=F
C
C     Test that more calls of DMLCFG are allowed.
C
      IF (NFMAX .gt. 0 .and. NFVALS .ge. NFMAX) THEN
          INFO=8
          GOTO 70
      END IF
C
C                            Test whether to reduce TOL.
C
      IF (TOL .GT. RELACC .AND. ITERC .GT. ITERK .AND.
     1  0.1D0*RELAXF .GE. max(DIFF,-0.5D0*DDOTG)) then
*        print'('' DMLC05.. Debug.. return to reduce TOL.''/
*    *   '' ........ TOL='',g13.5,'',   RELACC='',g13.5,'',  INFO='',
*    *   i5)', TOL, RELACC, INFO
         GOTO 70
      endif
c
C     Calculate the step along the search direction.
C
      ITERC=ITERC+1
*     print'(/'' DMLC05 895.. ITERC = '',i6)', ITERC
C          DMLC09 is LSRCH
      CALL DMLC09 (DMLCFG, N,X,G,D,XS,GS,RELACC,STEPCB,DDOTG,F,
     *             STEP,NFVALS, NFMAX,BRES, XL, XU, GMODE, W20)
      IF (STEP .EQ. 0.0D0) THEN
*         print'(/'' DMLC05 918.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F
          INFO=3
          if(abs(DDOTG) .le. EPS * abs(F)) then
             INFO = 2
          else
             SUM=0.0D0
             DO 50 I=1,N
   50        SUM=SUM+abs(D(I)*GS(I))
*            print'(/'' DMLC05 890..      DDOTG, SUM='',2g13.5/
*    *         '' RELACC*SUM, DDOTG+RELACC*SUM='',2g13.5)',
*    *         DDOTG, SUM, RELACC*SUM, DDOTG+RELACC*SUM
             IF (DDOTG+RELACC*SUM .GE. 0.0D0) INFO=2
          endif
          GOTO 70
      END IF
C
C     Revise XBIG.
C
      DO 60 I=1,N
   60 XBIG(I)=max(XBIG(I),abs(X(I)))
C
C     Revise the second derivative approximation.
C
      CALL DMLC19 (N,X,NACT,G,Z,ZTG,XS,GS,ZZNORM)
C
C     Add a constraint to the active set if it restricts the step.
C
      IF (STEP .EQ. STEPCB) THEN
          K=IACT(INDXBD)
          IF (K .GT. M) THEN
              K=K-M
              IF (K .LE. N) THEN
                  X(K)=XL(K)
              ELSE
                  X(K-N)=XU(K-N)
              END IF
          END IF
          CALL DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,XS,GS)
      END IF
      GOTO 10
C
   70 continue
      FVAL = F
      return
      end
c     ==================================================================
      SUBROUTINE DMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG,
     1  BRES,D,ZTG,RELACC,TOL,STEPCB,SUMRES,MEQL,MSAT,MTOT,INDXBD,
     2  GM,GMNEW,PARNEW,CGRAD,W20)
C
c     CONRES
C     Calculate and partition the residuals of the inactive constraints,
C       and set the gradient vector when seeking feasibility.
C     W20 has length 4*N+1.  Only locations 1 to N are used here.
c     ------------------------------------------------------------------
      integer N,M,LDA,IACT(*),NACT,MEQL,MSAT,MTOT,INDXBD
      integer I, IDIFF, J, JM, K, KL, MDEG, MSATK
      double precision A(LDA,*),B(*), DDOTG, PAR(*)
      double precision G(*), U(*),XBIG(*),BRES(*)
      double precision D(*), GM(*),GMNEW(*),PARNEW(*), CGRAD(*)
      double precision RELACC, RES, RESABS, STEPCB, SUM, SUMRES
      double precision TEMP, TOL
      double precision W20(0:*), XL(*),XU(*),X(*)
      double precision Z(*), ZTG(*)
c     ------------------------------------------------------------------
      IDIFF=MTOT-MSAT
C
      IF (IDIFF .GT. 0) THEN
          DO 10 I=1,N
   10     G(I)=0.0D0
          SUMRES=0.0D0
      END IF
      MSATK=MSAT
      MDEG=NACT
      MSAT=NACT
      KL=MEQL+1
      DO 50 K=KL,MTOT
      J=IACT(K)
C
C     Calculate the residual of the current constraint.
C
      IF (J .LE. M) THEN
          RES=B(J)
          RESABS=abs(B(J))
          DO 20 I=1,N
          RES=RES-X(I)*A(J,I)
   20     RESABS=RESABS+abs(XBIG(I)*A(J,I))
      ELSE
          JM=J-M
          IF (JM .LE. N) THEN
              RES=X(JM)-XL(JM)
              RESABS=abs(XBIG(JM))+abs(XL(JM))
          ELSE
              JM=JM-N
              RES=XU(JM)-X(JM)
              RESABS=abs(XBIG(JM))+abs(XU(JM))
          END IF
      END IF
      BRES(J)=RES
C
C     Set TEMP to the relative residual.
C
      TEMP=0.0D0
      IF (RESABS .NE. 0.0D0) TEMP=RES/RESABS
      IF (K .GT. MSATK .AND. TEMP .LT. 0.0D0) THEN
          IF (TEMP+RELACC .GE. 0.0D0) THEN
              IF (J .LE. M) THEN
                  SUM=abs(B(J))
                  DO 30 I=1,N
   30             SUM=SUM+abs(X(I)*A(J,I))
              ELSE
                  JM=J-M
                  IF (JM .LE. N) THEN
                      SUM=abs(X(JM))+abs(XL(JM))
                  ELSE
                      SUM=abs(X(JM-N))+abs(XU(JM-N))
                  END IF
              END IF
              IF (abs(RES) .LE. SUM*RELACC) TEMP=0.0D0
          END IF
      END IF
C
C     Place the residual in the appropriate position.
C
      IF (K .LE. NACT) GOTO 50
      IF (K .LE. MSATK .OR. TEMP .GE. 0.0D0) THEN
          MSAT=MSAT+1
          IF (MSAT .LT. K) THEN
              IACT(K)=IACT(MSAT)
          END IF
          IF (TEMP .GT. TOL) THEN
              IACT(MSAT)=J
          ELSE
              MDEG=MDEG+1
              IACT(MSAT)=IACT(MDEG)
              IACT(MDEG)=J
          END IF
C
C     Update the gradient and SUMRES if the constraint is violated when
C       seeking feasibility.
C
      ELSE
          IF (J .LE. M) THEN
              DO 40 I=1,N
   40         G(I)=G(I)+A(J,I)
          ELSE
              J=J-M
              IF (J .LE. N) THEN
                  G(J)=G(J)-1.0D0
              ELSE
                  G(J-N)=G(J-N)+1.0D0
              END IF
          END IF
          SUMRES=SUMRES+abs(RES)
      END IF
   50 CONTINUE
C
C     Seek the next search direction unless DMLC06 was called from
C     DMLC04 [GETFES] and feasibility has been achieved.
C
      STEPCB=0.0D0
      IF (IDIFF .GT. 0 .AND. MSAT .EQ. MTOT) GOTO 60
c          GETD
      CALL DMLC07 (N,M,A,LDA,IACT,NACT,PAR,G,Z,U,D,ZTG,RELACC,DDOTG,
     *             MEQL, MDEG,GM,GMNEW,PARNEW,CGRAD,W20)
C
C     Calculate the (bound on the) step-length due to the constraints.
C
      IF (DDOTG .LT. 0.0D0) THEN
          CALL DMLC18 (N,M,A,LDA,IACT,BRES,D,STEPCB,DDOTG,MDEG,MSAT,
     1      MTOT,INDXBD)
      END IF
      IF (IDIFF .EQ. 0) SUMRES=DDOTG
   60 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC07 (N,M,A,LDA,IACT,NACT,PAR,G,Z,U,D,ZTG,RELACC,
     1  DDOTG,MEQL,MDEG,GM,GMNEW,PARNEW,CGRAD,W20)
C
c     GETD
C     Initialize GM and cycle backwards through the active set.
C     W20 has length 4*N+1.  Only locations 0 to N are used here.
c     If W20(1) is nonnegative it means [D/S]MLC20 has been called and
c     has stored into W20(1:N) estimates of the errors in the gradient
c     values stored in G(1:N).
c     If W20(1) is negative the user is computing the gradients and
c     W20(1:N) is not being otherwise used.
c     This subr stores a value into W20(0).  This will only be used if
c     [D/S]MLC20 is called.
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, MEQL, MDEG
      integer I, IZ, J, JM, K
      double precision A(LDA,*), CGRAD(*), D(*), DDOTG, DDOTGM
      double precision G(*), GM(*), GMNEW(*)
      double precision PAR(*), PARNEW(*), RELACC, SIZE, TEMP, U(*)
      double precision W20(0:*), Z(*), ZTG(*)
c     ------------------------------------------------------------------
   10 DO 20 I=1,N
   20 GM(I)=G(I)
      K=NACT
   30 IF (K .GT. 0) THEN
C
C     Set TEMP to the next multiplier, but reduce the active set if
C       TEMP has an unacceptable sign.
C
          TEMP=0.0D0
          IZ=K
          DO 40 I=1,N
          TEMP=TEMP+Z(IZ)*GM(I)
   40     IZ=IZ+N
          TEMP=TEMP*U(K)
          IF (K .GT. MEQL .AND. TEMP .GT. 0.0D0) THEN
              CALL DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,K)
              GOTO 10
          END IF
C
C     Update GM using the multiplier that has just been calculated.
C
          J=IACT(K)
          IF (J .LE. M) THEN
              DO 50 I=1,N
   50         GM(I)=GM(I)-TEMP*A(J,I)
          ELSE
              JM=J-M
              IF (JM .LE. N) THEN
                  GM(JM)=GM(JM)+TEMP
              ELSE
                  GM(JM-N)=GM(JM-N)-TEMP
              END IF
          END IF
          PAR(K)=TEMP
          K=K-1
          GOTO 30
      END IF
C
C     Calculate the search direction and DDOTG.
C
      DDOTG=0.0D0
      IF (NACT .LT. N) THEN
c             SDEGEN
         CALL DMLC08 (N,M,A,LDA,IACT,NACT,PAR,Z,U,D,ZTG,GM,RELACC,
     1      DDOTGM,MEQL,MDEG,GMNEW,PARNEW,CGRAD)
         IF (DDOTGM .LT. 0.0D0) THEN
            SIZE = 0.0D0
            DO 60 I=1,N
               TEMP = D(I) * G(I)
               DDOTG=DDOTG + TEMP
               if (W20(1) .ge. 0.0D0) then
                  SIZE = SIZE + abs(D(I)*W20(I))
               else
                  SIZE = SIZE + abs(TEMP)
               endif
   60       continue

            if(W20(1) .lt. 0.0d0) SIZE = RELACC * SIZE
            if (DDOTG .ge. -SIZE) then
               DDOTG = 0.0D0
               W20(0) = 0.0D0
            else
               W20(0) = SIZE/abs(DDOTG)
            endif
         ELSE
            W20(0) = 0.0D0
         END IF
      END IF
      RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC08 (N,M,A,LDA,IACT,NACT,PAR,Z,U,D,ZTG,GM,RELACC,
     1  DDOTGM,MEQL,MDEG,GMNEW,PARNEW,CGRAD)
C
c     SDEGEN
C     Calculate the search direction and branch if it is not downhill.
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, MEQL, MDEG
      integer I, IDROP, ITEST, IZ, J, JM, K, KU, MP, NP
      double precision A(LDA,*), CGRAD(*), D(*), DDOTGM, DTEST
      double precision GM(*), GMNEW(*)
      double precision PAR(*), PARNEW(*), U(*)
      double precision RATIO, RELACC, SUM, TEMP, THETA
      double precision Z(*), ZTG(*)
c     ------------------------------------------------------------------
      MP=MEQL+1
      DTEST=0.0D0
C
   10 CALL DMLC17 (N,NACT,Z,D,ZTG,GM,RELACC,DDOTGM)
      IF (DDOTGM .EQ. 0.0D0) GOTO 120
C
C     Branch if there is no need to consider any degenerate constraints.
C     The test gives termination if two consecutive additions to the
C       active set fail to increase the predicted new value of F.
C
      IF (NACT .EQ. MDEG) GOTO 120
      NP=NACT+1
      SUM=0.0D0
      DO 20 J=NP,N
   20 SUM=SUM+ZTG(J)**2
      IF (DTEST .GT. 0.0D0 .AND. SUM .GE. DTEST) THEN
          IF (ITEST .EQ. 1) GOTO 120
          ITEST=1
      ELSE
          DTEST=SUM
          ITEST=0
      END IF
C
C     Add a constraint to the active set if there are any significant
C       violations of degenerate constraints.
C
      K=NACT
      CALL DMLC10 (N,M,A,LDA,IACT,NACT,Z,U,D,RELACC,MDEG,GMNEW,PARNEW,
     1  CGRAD)
      IF (NACT .EQ. K) GOTO 120
      PAR(NACT)=0.0D0
C
C     Calculate the new reduced gradient and Lagrange parameters.
C
   30 DO 40 I=1,N
   40 GMNEW(I)=GM(I)
      K=NACT
   50 TEMP=0.0D0
      IZ=K
      DO 60 I=1,N
      TEMP=TEMP+Z(IZ)*GMNEW(I)
   60 IZ=IZ+N
      TEMP=TEMP*U(K)
      PARNEW(K)=PAR(K)+TEMP
      IF (K .EQ. NACT) PARNEW(K)=min(PARNEW(K),0.0D0)
      J=IACT(K)
      IF (J .LE. M) THEN
          DO 70 I=1,N
   70     GMNEW(I)=GMNEW(I)-TEMP*A(J,I)
      ELSE
          JM=J-M
          IF (JM .LE. N) THEN
              GMNEW(JM)=GMNEW(JM)+TEMP
          ELSE
              GMNEW(JM-N)=GMNEW(JM-N)-TEMP
          END IF
      END IF
      K=K-1
      IF (K .GT. MEQL) GOTO 50
C
C     Set RATIO for linear interpolation between PAR and PARNEW.
C
      RATIO=0.0D0
      IF (MP .LT. NACT) THEN
          KU=NACT-1
          DO 80 K=MP,KU
          IF (PARNEW(K) .GT. 0.0D0) THEN
              RATIO=PARNEW(K)/(PARNEW(K)-PAR(K))
              IDROP=K
          END IF
   80     CONTINUE
      END IF
C
C     Apply the linear interpolation.
C
      THETA=1.0D0-RATIO
      DO 90 K=MP,NACT
   90 PAR(K)=min(THETA*PARNEW(K)+RATIO*PAR(K),0.0D0)
      DO 100 I=1,N
  100 GM(I)=THETA*GMNEW(I)+RATIO*GM(I)
C
C     Drop a constraint if RATIO is positive.
C
      IF (RATIO .GT. 0.0D0) THEN
          CALL DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP)
          DO 110 K=IDROP,NACT
  110     PAR(K)=PAR(K+1)
          GOTO 30
      END IF
C
C     Return if there is no freedom for a new search direction.
C
      IF (NACT .LT. N) GOTO 10
      DDOTGM=0.0D0
  120 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC09 (DMLCFG, N,X,G,D,XS,GS,RELACC,STEPCB, DDOTG,
     *                   F,STEP, NFVALS,NFMAX,GOPT, XL, XU, GMODE, W20)
c     LSRCH  Line search
c     4/18/91 CLL changed handling of ICOUNT and NFVALS.  Added GMODE.
c     ------------------------------------------------------------------
C                                  Subroutine Arguments
c
c     Input: DMLCFG, N, D(), RELACC, STEPCB, DDOTG, NFMAX, XL(), XU()
c     Inout: X(), G(), F, NFVALS, GMODE, W20()
c     Out:   GS(), STEP
c     Work space: XS(), GOPT()
c
c     DMLCFG, N,X(),G(),D()
c     XS()  [out]  Copy of the initial X().
c     GS()  [out]  Copy of the initial G().
c     RELACC,STEPCB, DDOTG,
c     F,STEP, NFVALS,NFMAX,GOPT(), XL(), XU(),
c     GMODE [integer,inout]  For use by DMLC20.
c     W20() [float,inout]  For use by DMLC20.
c     ------------------------------------------------------------------
      external DMLCFG
      integer N, NFVALS, NFMAX, GMODE
      integer I, ICOUNT
      double precision D(*)
      double precision DDOTG, DDOTGB, DGHGH, DGKNOT, DGLOW, DGMID, DGOPT
      double precision F, FBASE, FHGH, FLOW, FOPT
      double precision G(*), GOPT(*), GS(*)
      double precision RATIO, RELACC, RELINT, SBASE, SIZE, STEP, STEPCB
      double precision STPHGH, STPLOW, STPMIN, STPOPT, TEMP
      double precision W20(0:*), X(*), XL(*), XS(*), XU(*)
      logical HAVEG, KEEP, NOBACK
c     ------------------------------------------------------------------
C
C     Initialization.
c
*     print'(/'' DMLC09 1290.. STEPCB='',g13.5/
*    *   ''               .. D()='',3g13.5/(22x,3g13.5))',
*    *   STEPCB, (D(I),I=1,N)
      RELINT=0.9D0
      ICOUNT=0
      RATIO=-1.0D0
      DO 10 I=1,N
      XS(I)=X(I)
      GS(I)=G(I)
      GOPT(I)=G(I)
      IF (D(I) .NE. 0.0D0) THEN
          TEMP=abs(X(I)/D(I))
          IF (RATIO .LT. 0.0D0 .OR. TEMP .LT. RATIO) RATIO=TEMP
      END IF
   10 CONTINUE
      STEP=min(1.0D0,STEPCB)
      STPMIN=max(RELACC*RATIO,1.0D-12*STEP)
      STEP=max(STPMIN,STEP)
      SBASE=0.0D0
      FBASE=F
      DDOTGB=DDOTG
      STPLOW=0.0D0
      FLOW=F
      DGLOW=DDOTG
      STPHGH=0.0D0
      STPOPT=0.0D0
      FOPT=F
      DGOPT=abs(DDOTG)
      NOBACK = .false.
*     print'('' DMLC09 1318.. RATIO,STPMIN='',2g13.5/
*    *   ''                     DDOTG ='',g13.5)',
*    *   RATIO, STPMIN, DDOTG
C
C     Calculate another function and gradient value.
C
   20 continue
      if(NOBACK .and. STEP .lt. STPOPT) then
*        print*,' DMLC09 1355.. Ending DMLC09 on NOBACK.'
         go to 70
      endif
      NOBACK = .false.
      DO 30 I=1,N
   30    X(I)=XS(I)+STEP*D(I)
C%%        (*dmlcfg)( n, x, f, g, &haveg );
         call DMLCFG(N, X, F, G, HAVEG)
      NFVALS = NFVALS+1
      if(.not. HAVEG) CALL DMLC20(DMLCFG,N,X,F,G,XL,XU,GMODE,NFVALS,W20)
      ICOUNT=ICOUNT+1
*         print'(/'' DMLC09 1321..  ICOUNT='',i3/
*    *      ''   .. STEP,STPOPT ='',2g23.15)', ICOUNT, STEP, STPOPT
*         print'(''            .. X()='',2g23.15/(19x,2g23.15))',
*    *      (X(I),I=1,N)
*         print'(''             .. F ='',g23.15)',F
*         print'(''            .. G()='',2g23.15/(19x,2g23.15))',
*    *      (G(I),I=1,N)
      DGMID=0.0D0
      SIZE = 0.0d0
      DO 40 I=1,N
         TEMP = D(I)*G(I)
         DGMID =  DGMID + TEMP
         SIZE = SIZE + abs(TEMP)
   40 continue
      SIZE = SIZE * RELACC
*     print'(/'' DMLC09 1335..       FOPT - F ='',g13.5/
*    *   ''            ..          DGMID ='',g13.5/
*    *   ''            ..           SIZE ='',g13.5/
*    *   ''        .. DGOPT - abs(DGMID) ='',g13.5)',
*    *    FOPT - F, DGMID, SIZE, DGOPT - abs(DGMID)
      IF (F .LE. FOPT) THEN
         if(ICOUNT .eq. 1 .and. DGMID .le. 0.0d0) then
*           print '('' DMLC09 1381.. Setting NOBACK on new test.'')'
            NOBACK = .true.
            KEEP = .true.
         else
            KEEP = F .LT. FOPT .OR. abs(DGMID) .LT. DGOPT
         endif
         if(KEEP) then
              STPOPT=STEP
*             print'(/'' DMLC09 1384.. New STPOPT='',g13.5)', STPOPT
              FOPT=F
              DO 50 I=1,N
   50         GOPT(I)=G(I)
              DGOPT=abs(DGMID)
          END IF
      END IF
      IF (NFMAX .gt. 0 .and. NFVALS .ge. NFMAX) GOTO 70
C
C      Modify the bounds on the steplength or convergence.
C
      IF (F .GE. FBASE+0.1D0*(STEP-SBASE)*DDOTGB) THEN
          IF (STPHGH .GT. 0.0D0 .OR. F .GT. FBASE .OR. DGMID .GT.
     1      0.5D0*DDOTG) THEN
              STPHGH=STEP
              FHGH=F
              DGHGH=DGMID
              GOTO 60
          END IF
          SBASE=STEP
          FBASE=F
          DDOTGB=DGMID
      END IF
      IF (DGMID .GE. 0.7D0*DDOTGB) GOTO 70
      STPLOW=STEP
      FLOW=F
      DGLOW=DGMID
   60 IF (STPHGH .GT. 0.0D0 .AND. STPLOW .GE. RELINT*STPHGH) GOTO 70
C
C     Calculate the next step length or end the iterations.
C
*     print'(/'' DMLC09 1366.. ICOUNT, STEP='',i23,g23.15/
*    *   ''             STPLOW, STPHGH='',2g23.15)',
*    *   ICOUNT, STEP, STPLOW, STPHGH
      IF (STPHGH .EQ. 0.0D0) THEN
          IF (STEP .EQ. STEPCB) GOTO 70
          TEMP=10.0D0
          IF (DGMID .GT. 0.9D0*DDOTG) TEMP=DDOTG/(DDOTG-DGMID)
          STEP=min(TEMP*STEP,STEPCB)
*         print'('' DMLC09 1374.. To 20 with STEP='',g23.15)', STEP
          GOTO 20
      ELSE IF (ICOUNT .EQ. 1 .OR. STPLOW .GT. 0.0D0) THEN
          DGKNOT=2.0D0*(FHGH-FLOW)/(STPHGH-STPLOW)-0.5D0*(DGLOW+DGHGH)
          IF (DGKNOT .GE. 0.0D0) THEN
              RATIO=max(0.1D0,0.5D0*DGLOW/(DGLOW-DGKNOT))
          ELSE
              RATIO=(0.5D0*DGHGH-DGKNOT)/(DGHGH-DGKNOT)
          END IF
          STEP=STPLOW+RATIO*(STPHGH-STPLOW)
*         print'('' DMLC09 1384.. To 20 with STEP='',g23.15)', STEP
          GOTO 20
      ELSE
          STEP=0.1D0*STEP
          if (STEP .GE. STPMIN) then
*         print'('' DMLC09 1389.. To 20 with STEP='',g23.15)', STEP
             GOTO 20
          endif
      END IF
C
C     Return from subroutine.
C
   70 continue
      IF (STEP .NE. STPOPT) THEN
          STEP=STPOPT
          F=FOPT
          DO 80 I=1,N
             X(I)=XS(I)+STEP*D(I)
   80     G(I)=GOPT(I)
      END IF
      RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC10 (N,M,A,LDA,IACT,NACT,Z,U,D,RELACC,MDEG,ZZDIAG,
     1  GMNEW,CGRAD)
c     NEWCON
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, MDEG
      integer I, IADD, IZ, J, JM, JMV, K, KHIGH, NP
      double precision A(LDA,*), CGRAD(*), CVIOL, CVMAX, D(*), GMNEW(*)
      double precision RELACC, SAVABS, SAVSUM, SUM, SUMABS, SUMD, TEMP
      double precision U(*)
      double precision Z(*), ZZDIAG(*)
c     ------------------------------------------------------------------
C
C     Initialization.
C
      NP=NACT+1
      KHIGH=MDEG
      IZ=0
      DO 20 I=1,N
      ZZDIAG(I)=0.0D0
      DO 10 J=NP,N
   10 ZZDIAG(I)=ZZDIAG(I)+Z(IZ+J)**2
   20 IZ=IZ+N
C
C     Calculate the scalar products of D with its constraints.
C
   30 CVMAX=0.0D0
      DO 50 K=NP,KHIGH
      J=IACT(K)
      IF (J .LE. M) THEN
          SUM=0.0D0
          SUMABS=0.0D0
          SUMD=0.0D0
          DO 40 I=1,N
          TEMP=D(I)*A(J,I)
          SUM=SUM+TEMP
          SUMABS=SUMABS+abs(TEMP)
   40     SUMD=SUMD+ZZDIAG(I)*A(J,I)**2
      ELSE
          JM=J-M
          IF (JM .LE. N) THEN
              SUM=-D(JM)
          ELSE
              JM=JM-N
              SUM=D(JM)
          END IF
          SUMABS=abs(SUM)
          SUMD=ZZDIAG(JM)
      END IF
C
C     Pick out the most violated constraint, or return if the
C       violation is negligible.
C
      IF (SUM .GT. RELACC*SUMABS) THEN
          CVIOL=SUM*SUM/SUMD
          IF (CVIOL .GT. CVMAX) THEN
              CVMAX=CVIOL
              IADD=K
              SAVSUM=SUM
              SAVABS=SUMABS
          END IF
      END IF
   50 CONTINUE
      IF (CVMAX .LE. 0.0D0) GOTO 140
      IF (NACT .EQ. 0) GOTO 120
C
C     Set GMNEW to the gradient of the most violated constraint.
C
      J=IACT(IADD)
      IF (J .LE. M) THEN
          JMV=0
          DO 60 I=1,N
   60     GMNEW(I)=A(J,I)
      ELSE
          JMV=J-M
          DO 70 I=1,N
   70     GMNEW(I)=0.0D0
          IF (JMV .LE. N) THEN
              GMNEW(JMV)=-1.0D0
          ELSE
              JMV=JMV-N
              GMNEW(JMV)=1.0D0
          END IF
      END IF
C
C     Modify GMNEW for the next active constraint.
C
      K=NACT
   80 TEMP=0.0D0
      IZ=K
      DO 90 I=1,N
      TEMP=TEMP+Z(IZ)*GMNEW(I)
   90 IZ=IZ+N
      TEMP=TEMP*U(K)
      J=IACT(K)
      IF (J .LE. M) THEN
          DO 100 I=1,N
  100     GMNEW(I)=GMNEW(I)-TEMP*A(J,I)
      ELSE
          JM=J-M
          IF (JM .LE. N) THEN
              GMNEW(JM)=GMNEW(JM)+TEMP
          ELSE
              GMNEW(JM-N)=GMNEW(JM-N)-TEMP
          END IF
      END IF
C
C     Revise the values of SAVSUM and SAVABS.
C
      SUM=0.0D0
      SUMABS=0.0D0
      DO 110 I=1,N
      TEMP=D(I)*GMNEW(I)
      SUM=SUM+TEMP
  110 SUMABS=SUMABS+abs(TEMP)
      SAVSUM=min(SAVSUM,SUM)
      SAVABS=max(SAVABS,SUMABS)
      K=K-1
      IF (K .GE. 1) GOTO 80
C
C     Add the new constraint to the active set if the constraint
C       violation is still significant.
C
      IF (JMV .GT. 0) D(JMV)=0.0D0
      IF (SAVSUM .LE. RELACC*SAVABS) GOTO 130
  120 K=NACT
      CALL DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IADD,GMNEW,CGRAD)
      IF (NACT .GT. K) GOTO 140
C
C     Seek another constraint violation.
C
      IADD=NP
  130 IF (NP .LT. KHIGH) THEN
          K=IACT(KHIGH)
          IACT(KHIGH)=IACT(IADD)
          IACT(IADD)=K
          KHIGH=KHIGH-1
          GOTO 30
      END IF
  140 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC11 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,INFO,Z,U,XBIG,
     1  RELACC,TOL,MEQL)
c     SATACT
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, INFO
      integer I, IDROP, IZ, J, JX, K, MEQL
      double precision A(LDA,*),B(*)
      double precision RELACC, RES, RESABS, RESBIG, SAVEX, SCALE
      double precision TEMP, TEMPA, TOL, U(*)
      double precision X(*), XBIG(*), XL(*),XU(*), Z(*)
c     ------------------------------------------------------------------
      IF (NACT .EQ. 0) GOTO 50
      DO 30 K=1,NACT
C
C     Calculate the next constraint residual.
C
      J=IACT(K)
      IF (J .LE. M) THEN
          RES=B(J)
          RESABS=abs(B(J))
          RESBIG=RESABS
          DO 10 I=1,N
          TEMPA=A(J,I)
          TEMP=TEMPA*X(I)
          RES=RES-TEMP
          RESABS=RESABS+abs(TEMP)
   10     RESBIG=RESBIG+abs(TEMPA)*XBIG(I)
      ELSE
          JX=J-M
          IF (JX .LE. N) THEN
              RES=X(JX)-XL(JX)
              RESABS=abs(X(JX))+abs(XL(JX))
              RESBIG=XBIG(JX)+abs(XL(JX))
              SAVEX=XL(JX)
          ELSE
              JX=JX-N
              RES=XU(JX)-X(JX)
              RESABS=abs(X(JX))+abs(XU(JX))
              RESBIG=XBIG(JX)+abs(XU(JX))
              SAVEX=XU(JX)
          END IF
      END IF
C
C     Shift X if necessary.
C
      IF (RES .NE. 0.0D0) THEN
          TEMP=RES/RESABS
          IF (K .LE. MEQL) TEMP=-abs(TEMP)
          IF (TOL .EQ. RELACC .OR. TEMP+RELACC .LT. 0.0D0) THEN
              INFO=1
              SCALE=RES*U(K)
              IZ=K
              DO 20 I=1,N
              X(I)=X(I)+SCALE*Z(IZ)
              IZ=IZ+N
   20         XBIG(I)=max(XBIG(I),abs(X(I)))
              IF (J .GT. M) X(JX)=SAVEX
C
C     Else flag a constraint deletion if necessary.
C
          ELSE IF (RES/RESBIG .GT. TOL) THEN
              IACT(K)=-IACT(K)
          END IF
      END IF
   30 CONTINUE
C
C     Delete any flagged constraints and then return.
C
      IDROP=NACT
  40  IF (IACT(IDROP) .LT. 0) THEN
          IACT(IDROP)=-IACT(IDROP)
          CALL DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP)
      END IF
      IDROP=IDROP-1
      IF (IDROP .GT. MEQL) GOTO 40
   50 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,ZTC,
     1  CGRAD)
c     ADDCON
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, INDXBD
      integer I, ICON, INEWBD, IPIV, IZ, IZNBD, J, JP, NP
      double precision A(LDA,*), CGRAD(*)
      double precision RELACC, SUM, SUMABS, TEMP, TEMPA, TEMPB
      double precision U(*), WCOS, WPIV, WSIN
      double precision Z(*), ZTC(*)
c     ------------------------------------------------------------------
      NP=NACT+1
      ICON=IACT(INDXBD)
      IACT(INDXBD)=IACT(NP)
      IACT(NP)=ICON
C
C     Form ZTC when the new constraint is a bound.
C
      IF (ICON .GT. M) THEN
          INEWBD=ICON-M
          IF (INEWBD .LE. N) THEN
              TEMP=-1.0D0
          ELSE
              INEWBD=INEWBD-N
              TEMP=1.0D0
          END IF
          IZNBD=INEWBD*N-N
          DO 10 J=1,N
   10     ZTC(J)=TEMP*Z(IZNBD+J)
C
C     Else form ZTC for an ordinary constraint.
C
      ELSE
          DO 20 I=1,N
   20     CGRAD(I)=A(ICON,I)
          DO 30 J=1,N
          ZTC(J)=0.0D0
          IZ=J
          DO 30 I=1,N
          ZTC(J)=ZTC(J)+Z(IZ)*CGRAD(I)
   30     IZ=IZ+N
      END IF
C
C     Find any Givens rotations to apply to the last columns of Z.
C
      J=N
   40 JP=J
      J=J-1
      IF (J .GT. NACT) THEN
          IF (ZTC(JP) .EQ. 0.0D0) GOTO 40
          IF (abs(ZTC(JP)) .LE. RELACC*abs(ZTC(J))) THEN
              TEMP=abs(ZTC(J))
          ELSE IF (abs(ZTC(J)) .LE. RELACC*abs(ZTC(JP))) THEN
              TEMP=abs(ZTC(JP))
          ELSE
              TEMP=abs(ZTC(JP))*sqrt(1.0D0+(ZTC(J)/ZTC(JP))**2)
          END IF
          WCOS=ZTC(J)/TEMP
          WSIN=ZTC(JP)/TEMP
          ZTC(J)=TEMP
C
C     Apply the rotation when the new constraint is a bound.
C
          IZ=J
          IF (ICON .GT. M) THEN
              DO 50 I=1,N
              TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ)
              Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1)
              Z(IZ+1)=TEMP
   50         IZ=IZ+N
              Z(IZNBD+JP)=0.0D0
C
C     Else apply the rotation for an ordinary constraint.
C
          ELSE
              WPIV=0.0D0
              DO 60 I=1,N
              TEMPA=WCOS*Z(IZ+1)
              TEMPB=WSIN*Z(IZ)
              TEMP=abs(CGRAD(I))*(abs(TEMPA)+abs(TEMPB))
              IF (TEMP .GT. WPIV) THEN
                  WPIV=TEMP
                  IPIV=I
              END IF
              Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1)
              Z(IZ+1)=TEMPA-TEMPB
   60         IZ=IZ+N
C
C     Ensure orthogonality of Z(.,JP) to CGRAD.
C
              SUM=0.0D0
              IZ=JP
              DO 70 I=1,N
              SUM=SUM+Z(IZ)*CGRAD(I)
   70         IZ=IZ+N
              IF (SUM .NE. 0.0D0) THEN
                  IZ=IPIV*N-N+JP
                  Z(IZ)=Z(IZ)-SUM/CGRAD(IPIV)
              END IF
          END IF
          GO TO 40
      END IF
C
C     Test for linear independence in the proposed new active set.
C
      IF (ZTC(NP) .EQ. 0.0D0) GOTO 90
      IF (ICON .LE. M) THEN
          SUM=0.0D0
          SUMABS=0.0D0
          IZ=NP
          DO 80 I=1,N
          TEMP=Z(IZ)*CGRAD(I)
          SUM=SUM+TEMP
          SUMABS=SUMABS+abs(TEMP)
   80     IZ=IZ+N
          IF (abs(SUM) .LE. RELACC*SUMABS) GOTO 90
      END IF
C
C     Set the new diagonal element of U() and return.
C
      U(NP)=1.0D0/ZTC(NP)
      NACT=NP
   90 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC,TOL,
     1  MEQL)
C
C     ADJTOL:  Change the tolerance TOL to a smaller value, if possible.
c     May also recompute XBIG().
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, MEQL
      integer I, J, JM, K, KL
      double precision A(LDA,*), B(*)
      double precision RELACC, RES, RESABS, TOL, VIOL
      double precision X(*), XBIG(*), XL(*), XU(*)
c     ------------------------------------------------------------------
C     Set VIOL to the greatest relative constraint residual of the first
C       NACT constraints.
      VIOL=0.0D0
      IF (NACT .GT. MEQL) THEN
          KL=MEQL+1
          DO 20 K=KL,NACT
          J=IACT(K)
          IF (J .LE. M) THEN
              RES=B(J)
              RESABS=abs(B(J))
              DO 10 I=1,N
              RES=RES-A(J,I)*X(I)
   10         RESABS=RESABS+abs(A(J,I)*XBIG(I))
          ELSE
              JM=J-M
              IF (JM .LE. N) THEN
                  RES=X(JM)-XL(JM)
                  RESABS=XBIG(JM)+abs(XL(JM))
              ELSE
                  JM=JM-N
                  RES=XU(JM)-X(JM)
                  RESABS=XBIG(JM)+abs(XU(JM))
              END IF
          END IF
          IF (RES .GT. 0.0D0) VIOL=max(VIOL,RES/RESABS)
   20     CONTINUE
      END IF
C
C     Adjust TOL.
C
      TOL=0.1D0*min(TOL,VIOL)
      IF (TOL .LE. RELACC+RELACC) THEN
          TOL=RELACC
          DO 30 I=1,N
   30     XBIG(I)=abs(X(I))
      END IF
      RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP)
C
C     DELCON: Cycle through the constraint exchanges that are needed.
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, IDROP
      integer I, IBD, ICON, IPIV, ISAVE, IZ, IZBD, J, JP, NM
      double precision A(LDA,*), DENOM
      double precision RELACC, RJJP, SUM, TEMP, TEMPA, TEMPB
      double precision U(*), UJP, WCOS,  WPIV, WSIN
      double precision Z(*)
c     ------------------------------------------------------------------
      NM=NACT-1
      IF (IDROP .EQ. NACT) GOTO 60
      ISAVE=IACT(IDROP)
C
C     Cycle through the constraint exchanges that are needed.
C
      DO 50 J=IDROP,NM
      JP=J+1
      ICON=IACT(JP)
      IACT(J)=ICON
C
C     Calculate the (J,JP) element of R.
C
      IF (ICON .LE. M) THEN
          RJJP=0.0D0
          IZ=J
          DO 10 I=1,N
          RJJP=RJJP+Z(IZ)*A(ICON,I)
   10     IZ=IZ+N
      ELSE
          IBD=ICON-M
          IF (IBD .LE. N) THEN
              IZBD=IBD*N-N
              RJJP=-Z(IZBD+J)
          ELSE
              IBD=IBD-N
              IZBD=IBD*N-N
              RJJP=Z(IZBD+J)
          END IF
      END IF
C
C     Calculate the parameters of the next rotation.
C
      UJP=U(JP)
      TEMP=RJJP*UJP
      DENOM=abs(TEMP)
      IF (DENOM*RELACC .LT. 1.0D0) DENOM=sqrt(1.0D0+DENOM*DENOM)
      WCOS=TEMP/DENOM
      WSIN=1.0D0/DENOM
C
C     Rotate Z when a bound constraint is promoted.
C
      IZ=J
      IF (ICON .GT. M) THEN
          DO 20 I=1,N
          TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ)
          Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1)
          Z(IZ+1)=TEMP
   20     IZ=IZ+N
          Z(IZBD+JP)=0.0D0
C
C     Rotate Z when an ordinary constraint is promoted.
C
      ELSE
          WPIV=0.0D0
          DO 30 I=1,N
          TEMPA=WCOS*Z(IZ+1)
          TEMPB=WSIN*Z(IZ)
          TEMP=abs(A(ICON,I))*(abs(TEMPA)+abs(TEMPB))
          IF (TEMP .GT. WPIV) THEN
              WPIV=TEMP
              IPIV=I
          END IF
          Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1)
          Z(IZ+1)=TEMPA-TEMPB
   30     IZ=IZ+N
C
C     Ensure orthogonality to promoted constraint.
C
          SUM=0.0D0
          IZ=JP
          DO 40 I=1,N
          SUM=SUM+Z(IZ)*A(ICON,I)
   40     IZ=IZ+N
          IF (SUM .NE. 0.0D0) THEN
              IZ=IPIV*N-N+JP
              Z(IZ)=Z(IZ)-SUM/A(ICON,IPIV)
          END IF
      END IF
C
C     Set the new diagonal elements of U.
C
      U(JP)=-DENOM*U(J)
      U(J)=UJP/DENOM
   50 CONTINUE
C
C     Return.
C
      IACT(NACT)=ISAVE
   60 NACT=NM
      RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC15 (N,M,XL,XU,X,IACT,MEQL,INFO,Z,U,XBIG,RELACC)
c     INITZU:  Initialize RELACC, Z(), U().  Adjust X().  Set XBIG().
c     1990-07-20 CLL Changed error handling and setting of INFO.
c     Expect INFO to be set to 1 before entering this subr.  This subr
c     sets INFO = 4 if it finds XL(j) > XU(j) for some j.
c     ------------------------------------------------------------------
      external D1MACH
      integer N, M, IACT(*), MEQL, INFO
      integer I, IZ, J, JACT, NN
      double precision D1MACH
      double precision RELACC
      double precision U(*)
      double precision XL(*),XU(*),X(*), XBIG(*), Z(*)
c     ------------------------------------------------------------------
C
C                                                Set RELACC.
c     July 1990 CLL Commented out following 6 lines and replaced with
c     reference to D1MACH(4) which is the smallest x such that the
c     stored value of 1+x is greater than 1.  Efforts to determine
c     such an x with portable code such as the following 6 lines have
c     generally eventually failed on some new computer.  Also using
c     D1MACH allows adjustments to be made for known deficiencies
c     in particular computers, for example for the Cray X/MP & Y/MP.
c
*     ZTPAR=100.0D0
*     RELACC=1.0D0
*  10 RELACC=0.5D0*RELACC
*     TEMPA=ZTPAR+0.5D0*RELACC
*     TEMPB=ZTPAR+RELACC
*     IF (ZTPAR .LT. TEMPA .AND. TEMPA .LT. TEMPB) GOTO 10
      RELACC = 100.0D0 * D1MACH(4)
C
C     Seek bound inconsistencies and bound equality constraints.
C
      MEQL=0
      DO 20 J=1,N
         IF (XL(J) .GT. XU(J)) then
            INFO = 4
            call IERM1('DMLC15',INFO,0,'Bad bounds: XL(j) > XU(j)',
     *                 'j',J,',')
            call DERV1('XL(j)',XL(J),',')
            call DERV1('XU(j)',XU(J),'.')
            return
         endif
      IF (XL(J) .EQ. XU(J)) MEQL=MEQL+1
   20 CONTINUE
C
C     Initialize U, Z and XBIG.
C
      JACT=0
      NN=N*N
      DO 30 I=1,NN
   30 Z(I)=0.0D0
      IZ=0
      DO 40 I=1,N
      IF (XL(I) .EQ. XU(I)) THEN
          X(I)=XU(I)
          JACT=JACT+1
          U(JACT)=1.0D0
          IACT(JACT)=I+M+N
          J=JACT
      ELSE
          J=I+MEQL-JACT
      END IF
      Z(IZ+J)=1.0D0
      IZ=IZ+N
   40 XBIG(I)=abs(X(I))
      RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC16 (N,M,A,LDA,IACT,NACT,PAR,G,RESKT,Z,U,BRES,
     *  RELAXF, MEQL,SSQKT,PARW,RESKTW)
C
C     Calculate the Lagrange parameters and the residual vector.
c
c     Input:  N, M, A(), LDA, IACT, NACT, G(), Z(), U(), BRES(), MEQL
c     Output: PAR(), RESKT, RELAXF, SSQKT
c     Work space: PARW(), RESKTW()
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), NACT, MEQL
      integer I, ICASE, IZ, J, JM, K, KK, KL
      double precision A(LDA,*), BRES(*), G(*)
      double precision PAR(*), PARW(*)
      double precision RELAXF, RESKT(*), RESKTW(*)
      double precision SSQKT, SSQKTW, TEMP, U(*), Z(*)
c     ------------------------------------------------------------------
      DO 10 I=1,N
   10 RESKT(I)=G(I)
*     print'(/'' DMLC16 1993.. G()='',2g23.15/(19x,2g23.15))',
*    *   (G(I),I=1,N)
*     print'('' .. U()='', 3g23.15/(8x,3g23.15))', (U(I),I=1,3)
*     print'('' .. Z()='', 2g23.15/(8x,2g23.15))', (Z(I),I=1,16)
*     print'('' .. A(,)='',2g23.15/(9x,2g23.15))',((A(I,J),J=1,4),I=1,3)
      IF (NACT .GT. 0) THEN
          ICASE=0
   20     DO 50 KK=1,NACT
          K=NACT+1-KK
          J=IACT(K)
          TEMP=0.0D0
          IZ=K
          DO 30 I=1,N
          TEMP=TEMP+Z(IZ)*RESKT(I)
*       print'('' DMLC16 1991.. Z(),RESKT(),TEMP='',3g13.5)',
*    *      Z(IZ),RESKT(I),TEMP
   30     IZ=IZ+N
          TEMP=TEMP*U(K)
*       print'('' DMLC16 1995.. U(K),TEMP='', 2g13.5)', U(K), TEMP
          IF (ICASE .EQ. 0) PAR(K)=0.0D0
          IF (K .LE. MEQL .OR. PAR(K)+TEMP .LT. 0.0D0) THEN
              PAR(K)=PAR(K)+TEMP
*           print'('' DMLC16 1999.. TEMP,PAR(K)='',2g13.5)',
*    *         TEMP,PAR(K)
          ELSE
              TEMP=-PAR(K)
              PAR(K)=0.0D0
          END IF
          IF (TEMP .NE. 0.0D0) THEN
              IF (J .LE. M) THEN
                  DO 40 I=1,N
                     RESKT(I)=RESKT(I)-TEMP*A(J,I)
*          print'('' DMLC16 2009.. TEMP,A='',2g13.5/
*    *       ''                  RESKT(I) ='',g13.5)',
*    *       TEMP, A(J,I), RESKT(I)
   40             continue
              ELSE
                  JM=J-M
                  IF (JM .LE. N) THEN
                      RESKT(JM)=RESKT(JM)+TEMP
*          print'('' DMLC16 2018.. TEMP,RESKT(JM)='',2g13.5)',
*    *       TEMP,RESKT(JM)
                  ELSE
                      RESKT(JM-N)=RESKT(JM-N)-TEMP
*          print'('' DMLC16 2022.. TEMP,RESKT(JM-N)='',2g13.5)',
*    *       TEMP,RESKT(JM-N)
                  END IF
              END IF
          END IF
   50     CONTINUE
C
C     Calculate the sum of squares of the KT residual vector.
C
          SSQKT=0.0D0
          IF (NACT .EQ. N) GOTO 130
*         print'(/'' DMLC16 2014.. RESKT()='',4g12.4/(23x,4g12.4))',
*    *       (RESKT(I),I=1,N)
          DO 60 I=1,N
   60     SSQKT=SSQKT+RESKT(I)**2
*        print'('' DMLC16 2018.. SSQKT='',g13.5)', SSQKT
C
C     Apply iterative refinement to the residual vector.
C
          IF (ICASE .EQ. 0) THEN
              ICASE=1
              DO 70 K=1,NACT
   70         PARW(K)=PAR(K)
              DO 80 I=1,N
   80         RESKTW(I)=RESKT(I)
              SSQKTW=SSQKT
              GOTO 20
          END IF
C
C     Undo the iterative refinement if it does not reduce SSQKT.
C
          IF (SSQKTW .LT. SSQKT) THEN
              DO 90 K=1,NACT
   90         PAR(K)=PARW(K)
              DO 100 I=1,N
  100         RESKT(I)=RESKTW(I)
              SSQKT=SSQKTW
          END IF
C
C     Calculate SSQKT when there are no active constraints.
C
      ELSE
          SSQKT=0.0D0
          DO 110 I=1,N
  110     SSQKT=SSQKT+G(I)**2
      END IF
C
C     Predict the reduction in F if one corrects any positive residuals
C       of active inequality constraints.
C
      RELAXF=0.0D0
      IF (MEQL .LT. NACT) THEN
          KL=MEQL+1
          DO 120 K=KL,NACT
          J=IACT(K)
          IF (BRES(J) .GT. 0.0D0) THEN
              RELAXF=RELAXF-PAR(K)*BRES(J)
          END IF
  120     CONTINUE
      END IF
  130 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC17 (N,NACT,Z,D,ZTG,GM,RELACC,DDOTGM)
c     SDIRN:  Compute a search direction.
c     ------------------------------------------------------------------
      integer N, NACT
      integer I, IZ, J, NP
      double precision D(*), DDOTGM, GM(*)
      double precision RELACC, SUM, SUMABS, TEMP
      double precision Z(*), ZTG(*)
c     ------------------------------------------------------------------
      DDOTGM=0.0D0
      IF (NACT .GE. N) GOTO 60
C
C     Premultiply GM by the transpose of Z.
C
      NP=NACT+1
      DO 20 J=NP,N
      SUM=0.0D0
      SUMABS=0.0D0
      IZ=J
      DO 10 I=1,N
      TEMP=Z(IZ)*GM(I)
      SUM=SUM+TEMP
      SUMABS=SUMABS+abs(TEMP)
   10 IZ=IZ+N
      IF (abs(SUM) .LE. RELACC*SUMABS) SUM=0.0D0
   20 ZTG(J)=SUM
C
C     Form D by premultiplying ZTG by -Z.
C
      IZ=0
      DO 40 I=1,N
      SUM=0.0D0
      SUMABS=0.0D0
      DO 30 J=NP,N
      TEMP=Z(IZ+J)*ZTG(J)
      SUM=SUM-TEMP
   30 SUMABS=SUMABS+abs(TEMP)
      IF (abs(SUM) .LE. RELACC*SUMABS) SUM=0.0D0
      D(I)=SUM
   40 IZ=IZ+N
C
C     Test that the search direction is downhill.
C
      SUMABS=0.0D0
      DO 50 I=1,N
      TEMP=D(I)*GM(I)
      DDOTGM=DDOTGM+TEMP
   50 SUMABS=SUMABS+abs(TEMP)
      IF (DDOTGM+RELACC*SUMABS .GE. 0.0D0) DDOTGM=0.0D0
   60 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC18 (N,M,A,LDA,IACT,BRES,D,STEPCB,DDOTG,MDEG,MSAT,
     1  MTOT,INDXBD)
c     STEPBD:
C     Set steps to constraint boundaries and find the least positive
C     one.
c
c>> 1990-04-02 CLL Changes to avoid jumping into scope of Block If.
c     ------------------------------------------------------------------
      integer N, M, LDA, IACT(*), MDEG, MSAT, MTOT, INDXBD
      integer I, IFLAG, J, JM, K, KL
      double precision A(LDA,*),BRES(*),D(*), DDOTG, SP, STEPCB, TEMP
c     ------------------------------------------------------------------
      IFLAG=0
      STEPCB=0.0D0
      INDXBD=0
      K=MDEG
   10 K=K+1
      IF (K .gt. MTOT) go to 30
c
c                    Use remote code to compute scalar product.
      go to 80
   20 continue
C
C     Set BRES(J) to indicate the status of the j-th constraint.
C
          IF (SP*BRES(J) .LE. 0.0D0) THEN
              BRES(J)=0.0D0
          ELSE
              BRES(J)=BRES(J)/SP
              IF (STEPCB .EQ. 0.0D0 .OR. BRES(J) .LT. STEPCB) THEN
                  STEPCB=BRES(J)
                  INDXBD=K
              END IF
          END IF
          GO TO 10
   30 continue
C
C     Try to pass through the boundary of a violated constraint.
C
   40 continue
      IF (INDXBD .le. MSAT) go to 70
          IFLAG=1
          K=INDXBD
c                    Use remote code to compute scalar product.
          go to 80
   50     continue
          MSAT=MSAT+1
          IACT(INDXBD)=IACT(MSAT)
          IACT(MSAT)=J
          BRES(J)=0.0D0
          INDXBD=MSAT
          DDOTG=DDOTG-SP
          IF (DDOTG .LT. 0.0D0 .AND. MSAT .LT. MTOT) THEN
C
C     Seek the next constraint boundary along the search direction.
C
              TEMP=0.0D0
              KL=MDEG+1
              DO 60 K=KL,MTOT
              J=IACT(K)
              IF (BRES(J) .GT. 0.0D0) THEN
                  IF (TEMP .EQ. 0.0D0 .OR. BRES(J) .LT. TEMP) THEN
                      TEMP=BRES(J)
                      INDXBD=K
                  END IF
              END IF
   60         CONTINUE
              IF (TEMP .GT. 0.0D0) THEN
                  STEPCB=TEMP
                  GOTO 40
              END IF
          END IF
   70 continue
      RETURN
c     ------------------------------------------------------------------
c         This is a remote block of code to compute the
C         scalar product of D with the current constraint normal.
C
   80     J=IACT(K)
          IF (J .LE. M) THEN
              SP=0.0D0
              DO 90 I=1,N
   90         SP=SP+D(I)*A(J,I)
          ELSE
              JM=J-M
              IF (JM .LE. N) THEN
                  SP=-D(JM)
              ELSE
                  SP=D(JM-N)
              END IF
          END IF
C                         Return from remote block is selected by IFLAG.
          IF (IFLAG .EQ. 0) then
             go to 20
          ELSE
             go to 50
          ENDIF
      END
c     ==================================================================
      SUBROUTINE DMLC19 (N,X,NACT,G,Z,ZTG,XS,GS,ZZNORM)
C
C     Test if there is sufficient convexity for the update.
c     ------------------------------------------------------------------
      integer N,NACT
      integer I, IZ, K, KP, KM, NP
      double precision DD, DG, G(*), GS(*), SUM, TEMP, WCOS, WSIN
      double precision X(*), XS(*), Z(*),ZTG(*), ZZNORM
c     ------------------------------------------------------------------
      DD=0.0D0
      DG=0.0D0
      TEMP=0.0D0
      DO 10 I=1,N
      XS(I)=X(I)-XS(I)
      DD=DD+XS(I)**2
      TEMP=TEMP+GS(I)*XS(I)
      GS(I)=G(I)-GS(I)
   10 DG=DG+GS(I)*XS(I)
      IF (DG .LT. 0.1D0*abs(TEMP)) GOTO 90
C
C     Transform the Z matrix.
C
      K=N
   20 KP=K
      K=K-1
      IF (K .GT. NACT) THEN
          IF (ZTG(KP) .EQ. 0.0D0) GOTO 20
          TEMP=abs(ZTG(KP))*sqrt(1.0D0+(ZTG(K)/ZTG(KP))**2)
          WCOS=ZTG(K)/TEMP
          WSIN=ZTG(KP)/TEMP
          ZTG(K)=TEMP
          IZ=K
          DO 30 I=1,N
          TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ)
          Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1)
          Z(IZ+1)=TEMP
   30     IZ=IZ+N
          GOTO 20
      END IF
C
C     Update the value of ZZNORM.
C
      IF (ZZNORM .LT. 0.0D0) THEN
          ZZNORM=DD/DG
      ELSE
          TEMP=sqrt(ZZNORM*DD/DG)
          ZZNORM=min(ZZNORM,TEMP)
          ZZNORM=max(ZZNORM,0.1D0*TEMP)
      END IF
C
C     Complete the updating of Z.
C
      NP=NACT+1
      TEMP=sqrt(DG)
      IZ=NP
      DO 40 I=1,N
      Z(IZ)=XS(I)/TEMP
   40 IZ=IZ+N
      IF (NP .LT. N) THEN
          KM=NP+1
          DO 80 K=KM,N
          TEMP=0.0D0
          IZ=K
          DO 50 I=1,N
          TEMP=TEMP+GS(I)*Z(IZ)
   50     IZ=IZ+N
          TEMP=TEMP/DG
          SUM=0.0D0
          IZ=K
          DO 60 I=1,N
          Z(IZ)=Z(IZ)-TEMP*XS(I)
          SUM=SUM+Z(IZ)**2
   60     IZ=IZ+N
          IF (SUM .LT. ZZNORM) THEN
              TEMP=sqrt(ZZNORM/SUM)
              IZ=K
              DO 70 I=1,N
              Z(IZ)=TEMP*Z(IZ)
   70         IZ=IZ+N
          END IF
   80     CONTINUE
      END IF
   90 RETURN
      END
c     ==================================================================
      SUBROUTINE DMLC20 (DMLCFG, N, X, F, G, XL, XU, GMODE, NFVALS, W)
c>> 1991-04-15 FTK & CLL  Made FIRST an argument of DMLC20
c>> 1990-09-11 Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA.
c
c This subr computes a finite difference approximation to a gradient
c vector.  Values of F are not computed outside the bounds given by XL()
c and XU().  For any I for which XL(I) = XU(I), G(I) will be set to zero
c and computation to estimate G(I) will be skipped.
c On entry F must have already been evaluated at X().
c The algorithm used is based on the following observations.
c 1. Since the gradient vector is being used in an iteration, there is
c    no point to computing it more accurately than is useful in the
c    iteration.  An extra N function values to get better error
c    estimates is not likely to pay for itself.
c 2. If the delta x used for the difference is in the direction of the
c    next move, then discretization error in approximating the gradient
c    is likely to help rather than hurt.
c 3. The sign of the change in X(I) for the next move is likely to be
c    the same as it was on the last move.
c 4. Because of 3 and 4, the increment used is larger than would be
c    suggested from a consideration of round off and discretization
c    errors.
c 5. When the distance between values on successive iterations gets very
c    small it is useful to have second derivative information in order
c    to balance the effects of round off and discretization errors.
c
c *************************** Variables ********************************
c
c D      The difference between the current X(I) and the value on the
c        previous iteration.
c DXN    L2 norm of the move just taken.
c DMLCFG [in, subroutine name]  Name of user-provided subroutine for
c    computing function values.  When called from this subr, DMLCFG must
c    compute F, set HAVEG = .false.,  and not store into G().
c EPSR   Save variable = relative machine precision.
c F      [in, float]  On entry F will contain the function value
c    evaluated at the given X().
c FAC    Save variable = sqrt(relative machine precision)
c FACX   Save variable, 4. * FAC.
c G()    [out, float]  On return G() will contain an approximation of
c    the N-dimensional gradient of f with respect to x at X() computed
c    using one-sided finite differences.
c GMODE [inout, integer] Has value 0, or 1.
c       0 means initial entry.  Compute 1-sided differences, assuming no
c         saved info, and set GMODE to 1.
c       1 Compute differences, using saved info.
c H      Increment in X(I) used for computing gradient.
c HB     Lower bound desired for value of abs(H).
c KOUNT  Counts interations.  If multiple of 10, two sided differences
c    are computed.
c LGL    W(LGL+I) contains the value of G(I) from previous iteration.
c LNOPP  Value of NOPP for start of next iteration.
c LPP    W(LPP+I) contains estimate of the second derivative (d/dv)G(I),
c    where v is the variable consisting of the linear combinations of
c    X's that gave the last move.
c LXL    W(LXL+I) contains X() from the previous iteration.
c MODE1  .true. if have 1-sided difference, = .false. if have 2-sided.
c N      [in, integer]  Number of components in X(), G(), XL(), & XU().
c NFVALS [inout, integer]  Counter of number of calls to DMLCFG.
c NOPP   Logical variable set .true. if second derivative information is
c    not available.
c TP     Temporary storage.
c X()    [in, work, float]  On entry contains the current N-dimensional
c    parameter vector.  This subr will alter components in X() but will
c    reset X() to its original contents on return.
c XI     Current base value for X(I) when computing F for a gradient.
c XL()   [in, float]  Lower bounds for components of X().
c XNEW   Current value for X(I) in computing F for a gradient.  Also
c    used for temporary storage of 1-sided gradient.
c XU()   [in, float]  Upper bounds for components of X().
c W()    [inout, float] Work array of dimension at least 4N+1.
c    W(0) and W(1) are initialized to -1.0 in [D/S]MLC02.  If this
c    subr is never called, i.e., if the user provides computed
c    gradients, then W(1) will remain at this initial value.
c    If this subr is called, it will (generallly) set W(1) positive.
c    [D/S]MLC07 will (generally) set W(0) positive for subsequent
c    testing when and if this subr is called.
c    W(1:N) contains an estimate of the error in G(I) for use in
c    [D/S]MLC07.
c    W(N+1 : 4*N) is used as temporary work space in this subr.  See the
c    usage of the indices LGL, LPP, LXL.
c     ------------------------------------------------------------------
c
c ********************** Variable Declarations *************************
c
      external D1MACH, DMLCFG
      integer GMODE, LGL, LPP, LXL, N
      integer I, KDEBUG, KOUNT, NFVALS
      double precision D1MACH
      double precision D, DXN, EPSR
      double precision F, FAC, FACX, FNEW, G(N), H, HB, TP
      double precision X(N), XI, XL(N), XNEW, XU(N), W(0:*)
      logical MODE1, NOPP, LNOPP
c%%     LOGICAL32 haveg;
      logical HAVEG
      save EPSR, FAC, FACX, KOUNT, LNOPP
      save KDEBUG
      data KDEBUG / 0 /
      data FAC / 0.0D0 /, KOUNT / 0 /
 1000    format ('0Gradient #=', I3, '  F=', 1p,e24.17, '  DXN=',
     1      1p,e15.8 / '  I', 11x, 'X', 14X, 'DX', 10x, 'G', 10x, 'ERR',
     2      9X, 'H', 9X, 'dG/dV')
c
c ********************** Start of Executable Code **********************
c
      if(FAC .eq. 0.0D0) then
c                             Get machine constants
         EPSR = D1MACH(4)
         FAC = sqrt(EPSR)
         FACX = 4.D0*FAC
      endif
c                             Set up base locations in W()
      LXL = N
      LGL = LXL + LXL
      LPP = LGL + LXL
      DXN = 0.D0
      if (GMODE .eq. 0) then
         GMODE = 1
         KOUNT = 1
         NOPP = .true.
      else
         KOUNT = KOUNT + 1
         if (mod(KOUNT,10) .eq. 1) then
            if (W(0) .gt. 1.D-5) KOUNT = 0
         else
            if (W(0) .gt. 0.01D0) KOUNT = 0
         end if
         do 10 I = 1, LXL
            DXN = DXN + (X(I) - W(LXL+I))**2
   10    continue
         DXN = sqrt(DXN)
         NOPP = LNOPP
      end if
      LNOPP = .false.
      if (KDEBUG .ne. 0) then
         print 1000, KDEBUG, F, DXN
         KDEBUG = KDEBUG + 1
      end if
      do 40 I = 1, LXL
         XI = X(I)
         D = XI - W(LXL+I)
         if (NOPP) then
            HB = FAC * (1.D0 + abs(XI))
         else
            HB = FAC * (sqrt(abs(F/W(LPP+I))) + FACX*(abs(X(I))+FACX))
         end if
   20    if ((W(LPP+I) .lt. 0.D0) .or. (mod(KOUNT, 10) .eq. 0)) then
            if ((XI + HB .le. XU(I)) .and. (XI - HB .ge. XL(I))) then
               MODE1 = .false.
               H = sign(min(100.D0*HB, min(XU(I)-XI, XI-XL(I))), D)
               go to 30
            end if
         end if
         MODE1 = .true.
         H = sign(max(HB, abs(1.D-4*D)), D)
   30    XNEW = XI + H
c                                 Adjust H to keep XNEW within bounds.
         if (XNEW .gt. XU(I)) then
            XNEW = XU(I)
            if (XNEW - XI .lt. HB) then
               if (XI - XL(I) .gt. XNEW - XI)  XNEW = max(XL(I), XI-HB)
            end if
            H = XNEW - XI
         else if (XNEW .lt. XL(I)) then
            XNEW = XL(I)
            if (XI - XNEW .lt. HB) then
               if (XU(I) - XI .gt. XI - XNEW)  XNEW = min(XU(I), XI+HB)
            end if
            H = XNEW - XI
         end if
         if(H .eq. 0.D0) then
            G(I) = 0.D0
            W(I) = 0.D0
         else
            X(I) = XNEW
C%%           (*dmlcfg)( lxl, x, &fnew, g, &haveg );
            CALL DMLCFG(LXL, X, FNEW, G, HAVEG)
            NFVALS = NFVALS+1
            G(I) = (FNEW - F) / H
            W(I) = EPSR * (abs(F) + abs(FNEW)) / abs(H)
            if (.not. MODE1) then
               XNEW = XI - H
               X(I) = XNEW
C%%              (*dmlcfg)( lxl, x, &fnew, g, &haveg );
               CALL DMLCFG(LXL, X, FNEW, G, HAVEG)
               NFVALS = NFVALS+1
               XNEW = (FNEW - F) / (XNEW - XI)
               W(I) = W(I) + 5.D-4 * abs(XNEW - G(I))
               G(I) = .5D0 * (G(I) + XNEW )
            end if
            X(I) = XI
         end if
         if (DXN .eq. 0.D0) then
            LNOPP = .true.
            if (MODE1) W(I) = 0.D0
         else
            W(LPP+I) = max(0.D0, abs((G(I) - W(LGL+I)) - 10.D0 *
     1         abs(W(I))) / DXN) + FAC * abs(F)
            TP = 10.D0 * abs(H) * W(LPP+I)
            if (MODE1) then
               if (TP .gt. abs(G(I))) then
                  if (mod(KOUNT, 10) .ne. 0) then
                     KOUNT = 0
                     go to 20
                  end if
               end if
               W(I) = W(I) + TP
            else
               if (abs(XNEW-G(I)) .gt. min(10.D0*TP, 1.D-3 * abs(G(I))))
     1             W(LPP+I) = -W(LPP+I)
            end if
         end if
         if (KDEBUG .ne. 0) then
            print '(I3, 1P,E22.14, E10.2, E12.4, E12.4, E11.3, E10.2)',
     *            I,X(I),X(I)-W(LXL+I),G(I),W(I),H,W(LPP+I)
         end if
         W(LXL+I) = XI
         W(LGL+I) = G(I)
   40 continue
      return
      end
c     ==================================================================
      SUBROUTINE DMLC21(IPMORE, ALL, N, ITERC, NFVALS,
     *           F, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT)
c>> 1990-07-11 CLL
c     This subr prints results (first feasible, intermediate, and
c     final).
c     ------------------------------------------------------------------
c                                 Arguments
c
c  IPMORE [in, integer]  = 0 or 1.  0 means print ITERC, NFVALS, F, X(),
c        and G().  1 means also print IACT(), PAR(), and RESKT().
c        However see ALL for an exception.
c  ALL [in, logical]  False means G(), PAR(), and RESKT() might not
c        contain valid data and therefore will not be printed.
c  N, ITERC, NFVALS,
c  F, TOL, X(), G(), NACT, IACT(), PAR(), RESKT(),
c  ENRMKT [in]  Euclidian norm of RESKT().
c     ------------------------------------------------------------------
c%%   long int i, k;
      integer I
      integer IPMORE, N, ITERC, NFVALS, NACT, IACT(*)
      double precision ENRMKT, F, X(N), G(N), PAR(*), RESKT(N), TOL
      logical ALL
c     ------------------------------------------------------------------
      print'(/1x,i6,'' iterations. '',i6,'' function evals.    F = '',
     * g13.5/1x,''  TOL ='',g13.5,''   Norm of RESKT ='',g13.5)',
     * ITERC, NFVALS, F, TOL, ENRMKT
C%%     printf( "\n     X =" );
C%%     for (i = 0; i < n; i+=5){
C%%         for (k = i; k < (i < n - 4 ? i+5 : n); k++)
C%%            printf( "%13.5g", x[k] );
C%%         if (i < n - 4) printf ( "\n        " );}
      print'(''     X ='',5g13.5/(8x,5g13.5))', (X(I),I=1,N)
      if (IPMORE .gt. 0) then
         if(ALL) then
C%%         printf( "\n     G =" );
C%%         for (i = 0; i < n; i+=5){
C%%             for (k = i; k < (i < n - 4 ? i+5 : n); k++)
C%%                printf( "%13.5g", g[k] );
C%%             if (i < n - 4) printf ( "\n        " );}
            print'(a,5g13.5/(8x,5g13.5))','     G =',(G(I),I=1,N)
            if (NACT .eq. N) then
               print*,'    Kuhn-Tucker residual vector is zero.'
            else
C%%            printf( "\n RESKT =" );
C%%            for (i = 0; i < n; i+=5){
C%%                for (k = i; k < (i < n - 4 ? i+5 : n); k++)
C%%                   printf( "%13.5g", reskt[k] );
C%%                if (i < n - 4) printf ( "\n        " );}
             print'('' RESKT ='',5g13.5/(8x,5g13.5))', (RESKT(I),I=1,N)
            end if
         end if
         if (NACT .eq. 0) then
            print*,'    No active constraints.'
         else
C%%         printf( "\n  IACT =" );
C%%         for (i = 0; i < nact; i+=12){
C%%             for (k = i; k < (i < nact - 11 ? i+12 : nact); k++)
C%%                printf( "%5ld", iact[k] );
C%%             if (i < nact - 11) printf ( "\n        " );}
            print'(a,12i5/(8x,12i5))','  IACT =',(IACT(I),I=1,NACT)
            if (ALL) then
C%%           printf( "\n   PAR =" );
C%%           for (i = 0; i < nact; i+=5){
C%%               for (k = i; k < (i < nact - 4 ? i+5 : nact); k++)
C%%                  printf( "%13.5g", par[k] );
C%%               if (i < nact - 4) printf ( "\n        " );}
C%%           printf( "\n" );
            print'(''   PAR ='',5g13.5/(8x,5g13.5))', (PAR(I),I=1,NACT)
            end if
         end if
      end if
      return
      end






