*     *************************************************************************

      subroutine update( inform, x, xst, elst, elptr, elvar,
     +                   fuval,  fcalc,  lfc,  olgra,  gptr,
     +                   lg,  hptr,  nefval, negval, nehval,
     +                   dx, s, y, w4, recalc,   prec,  dpr,
     +                   nbfgs, nrk1, noupd )

*     *************************************************************************

*     Purpose :
*     ---------

*     This routine updates the element gradient vectors and the
*     element Hessians matrices that are stored in array FUVAL.

*     A considerable flexibility in the computation of the
*     derivatives is provided in LSNNO. If the element gradient
*     and/or the element Hessian are not supplied by the user,
*     the software contains an assortment of routines for
*     computing suitable approximations of these quantities.
*     The way these quantities will be computed is indicated
*     in array ELST that contains the status of the elemnt
*     functions.

*     Every time informations on some element function values or
*     on some element gradient values or on some element Hessian
*     values are needed, a return is made to the main program.
*     This routine is then re-entered with the information available
*     and the calculation proceeds.

*     Parameters :
*     ------------

*     inform  ( int )

*       input  : must be set to zero the first time this routine
*                is called.
*                If it is not equal to zero, the routine is then
*                re-entered with some information that was needed
*                to proceed the calculation.
*       output : If it is not equal to zero, some element function
*                values or some element gradient values or some
*                element Hessian values are needed to proceed the
*                calculation. A return is then made to the main
*                program.
*                Otherwise, no further information is needed.

*     x       ( dble )

*       input  : the current iterate vector.
*       output : unmodified.

*     xst     ( int )

*       input  : vector containing the status of the variables.
*       output : unmodified.

*     elst    ( int )

*       input  : vector containing the status of the element
*                functions.
*       output : the status of the element functions updated
*                by a quasi-Newtom formula may have been
*                changed.

*     elptr   ( int )

*       input  : array whose kth value is the position of
*                the first variable of element k, in the
*                list ELVAR.
*       output : unmodified.

*     elvar   ( int )

*       input  : array containing the indices of the varaibles
*                in the first element, followed by those in the
*                second element, etc.
*       output : unmodified.

*     fuval   ( dble )

*       input  : array used to store the function and derivative
*                values of the element functions.
*       output : the derivative values of the element functions
*                are updated.

*     fcalc   ( int )

*       input  : meaningless.
*       output : this vector contains the indices of the element
*                functions for which information is required.

*     lfc     ( int )

*       input  : meaningless.
*       output : the number of element functions for which
*                information is required.

*     olgra   ( dble )

*       input  : this array contains the element gradient
*                vectors evaluated at the previous iterate.
*       output : it contains the differences between the
*                element gradient vectors evaluated at point
*                X and the element gradient vectors evaluated
*                at the previous iterate.

*     gptr    ( int )

*       input  : array whose ith value is the position of the
*                first component of the ith element gradient
*                in FUVAL.
*       output : unmodified.

*     lg      ( int )

*       input  : the total lenght of the element gradient vectors.
*       output : unmodified.

*     hptr    ( int )

*       input  : array whose ith value is the position of the
*                first component of the ith element Hessian
*                in FUVAL.
*       output : unmodified.

*     nefval  ( int )

*       input  : the number of element function evaluations.
*       output : this number is increased by one each time
*                an element function needs to be evaluated.

*     negval  ( int )

*       input  : the number of element gradient evaluations.
*       output : this number is increased by one each time
*                an element gradient needs to be evaluated.

*     nehval  ( int )

*       input  : the number of element Hessian evaluations.
*       output : this number is increased by one each time
*                an element Hessian needs to be evaluated.

*     dx      ( dble )

*       input  : array used as workspace.
*       output : it contains the differences between the
*                current iterate X and the previous iterate.

*     s       ( dble )

*       input  : array used as workspace.
*       output : it contains for each element function the
*                differences between the current iterate X
*                and the previous iterate.

*     y       ( dble )

*       input  : array used as workspace.
*       output : it contains for each element function the
*                differences between the element gradient
*                vector evaluated at the current iterate X
*                and the element gradient vector evaluated
*                at the previous iterate.

*     w4      ( dble )

*       input  : array used as workspace.
*       output : meaningless.

*     recalc  ( int  )

*       input  : array containing the indices of the elements whose
*                derivatives need updating.
*       output : unmodified.

*     prec    ( log )

*       input  : .true.  if the user wants to use a preconditioned
*                        conjugate gradient scheme to solve the linear
*                        systems of equations which arise at each
*                        iteration of the minimization procedure.
*                .false. if the user does not want to use a
*                        preconditioner.
*       output : unmodified.

*     dpr     ( dble )

*       input  : array containing the information necessary to build
*                the preconditioner associated with Dembo's strategy.
*       output : the updated information.

*     nbfgs   ( int )

*       input  : meaningless.
*       output : the number of element Hessians that have been
*                updated by the BFGS formula.

*     nrk1    ( int )

*       input  : meaningless.
*       output : the number of element Hessians that have been
*                updated by the Rank-One formula.

*     noupd   ( int )

*       input  : meaningless.
*       output : the number of element Hessians for which
*                no update was needed.

*     Routines used :
*     ---------------

*     mod, sqrt, abs, max.
*     dcopy, ddot, analg, elint, fdingf, fdelgf,
*     fdinhg, fdelhg, fdinhf, fdelhf, range, gathr0,
*     hessid, selhqn, elcvx, scbfgs, scrk1.

*     Programming :
*     -------------

*     D. Tuyttens

*     ========================================================================

*     Routine parameters

      integer           inform,   xst(*),  elst(*),   elptr(*),
     +                  elvar(*), gptr(*), lg, hptr(*), nefval,
     +                  negval,  nehval,  fcalc(*), lfc, nbfgs,
     +                  nrk1,  noupd, recalc(*)
      logical           prec
      double precision  x(*),  fuval(*),  dx(*),   s(*),  y(*),
     +                  w4(*), olgra(*), dpr(*)

*     Internal variables

      integer           iel, els, dim, eldim, intdim, i, ij, k,
     +                  debx, debg, debh,  ifc, nrec
      logical           analg, elint, elhqn, elcvx,  newp,
     +                  upd, negcur, qnd
      double precision  sbs, yts, alpha, sqrem, ddot

*     Save variables

      save              iel, els, dim, eldim, intdim, i, ij, k,
     +                  debx, debg, debh,  ifc, nrec
      save              newp, upd, negcur, qnd
      save              sbs, yts, alpha, sqrem

*     Common specifications

      integer           arcs, nodes, elem
      double precision  epsmch, huge, tiny, tol
      common / prbdim / arcs, nodes, elem
      common / prbmch / epsmch, huge, tiny, tol

      double precision  zero, one, two, three, half, tenm1, tenm2, tenm4
      common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4
*
*     Reverse communication test.
*
      if( inform.ne.0 ) then
        inform = inform - 14
        go to( 200, 910, 920, 500, 500,
     +         300, 300, 700, 700, 700,
     +         700, 700, 700 ) inform
      endif
*
*     We obtain from the informations stored in array
*     ELST, the way the element gradients and element
*     Hessians will be estimated.
*
      nrec = lfc
      call dcopy( lg, fuval(gptr(1)), 1, olgra, 1 )
*
*
*     The gradient and Hessian are provided by the user.
*
      qnd = .false.
      lfc = 0
      do 100 ifc = 1, nrec
        iel = recalc( ifc )
        els = mod( elst(iel), 8 )
        if( analg(elst(iel)) .and. els.eq.0 ) then
          lfc = lfc + 1
          fcalc(lfc) = iel
        endif
        qnd = qnd .or. els .ge. 3
 100  continue
*
*     Element gradient and Hessian evaluations.
*
      if( lfc .gt. 0 ) then
        negval = negval + lfc
        nehval = nehval + lfc
        inform = 2
        return
      endif
 910  continue
*
*     Only the gradient is provided by the user.
*
      lfc = 0
      do 210 ifc = 1, nrec
        iel = recalc( ifc )
        els = mod( elst(iel), 8 )
        if( analg(elst(iel)) .and. els.ne.0 ) then
          lfc = lfc + 1
          fcalc(lfc) = iel
        endif
 210  continue
*
*     Element gradient evaluations.
*
      if( lfc .gt. 0 ) then
        negval = negval + lfc
        inform = 3
        return
      endif
 920  continue
*
*     Only the Hessian is  provided by the user.
*
      lfc = 0
      do 120 ifc = 1, nrec
        iel = recalc( ifc )
        els = mod( elst(iel), 8 )
        if( .not. analg(elst(iel)) .and. els.eq.0 ) then
          lfc = lfc + 1
          fcalc(lfc) = iel
        endif
 120  continue
*
*     Element Hessian evaluations.
*
      if( lfc .gt. 0 ) then
        nehval = nehval + lfc
        inform = 1
        return
      endif
 200  continue
      inform = 0
*
*     The gradient vector is estimated by
*     differences in the function values.
*
      lfc = 1
      if( nrec .ge. 1 ) then
        ifc = 1
 10     continue
        iel = recalc(ifc)
        if( .not. analg(elst(iel)) ) then
          fcalc(1) = iel
 300      continue
          if( elint(elst(iel)) ) then
*
*           The element gradient has an internal representation.
*
            call fdingf( iel, elptr, elvar, x, w4(arcs+1), s, y,
     +                   gptr,  fuval,  nefval,  inform )
          else
*
*           The element gradient has an elemental representation.
*
            call fdelgf( iel,  elptr, elvar,  xst, x,
     +                   gptr, fuval, nefval, inform )
          endif
          if( inform.ne.0 ) return
*
*         Element function evaluations.
*
        endif
        ifc = ifc + 1
        if ( ifc .le. nrec ) go to 10
      endif
*
*     The Hessian matrix is estimated by
*     differences in the gradient values.
*
      lfc = 1
      if( nrec .ge. 1 ) then
        ifc = 1
 40     continue
        iel = recalc(ifc)
        if( mod( elst(iel), 8 ) .eq. 1 ) then
          fcalc(1) = iel
 500      continue
          if( elint(elst(iel)) ) then
*
*           The element Hessian has an internal representation.
*
            call fdinhg( iel, elptr, elvar, x, w4(arcs+1), s, y, w4,
     +                   fuval,  gptr,  hptr,  negval,  inform )
          else
*
*           The element Hessian has an elemental representation.
*
            call fdelhg( iel, elptr, elvar, xst, x, gptr, hptr,
     +                   fuval,  negval,  w4(arcs+1),  inform )
          endif
          if( inform.ne.0 ) return
*
*         Element gradient evaluations.
*
        endif
        ifc = ifc + 1
        if ( ifc .le. nrec ) go to 40
      endif
*
*     The Hessian matrix is estimated by
*     differences in the function values.
*
      lfc = 1
      if( nrec .ge. 1 ) then
        ifc = 1
 50     continue
        iel = recalc(ifc)
        if( mod( elst(iel), 8 ) .eq. 2 ) then
          fcalc(1) = iel
 700      continue
          if( elint(elst(iel)) ) then
*
*           The element Hessian has an internal representation.
*
            call fdinhf( iel, elptr, elvar, x, w4(arcs+1), s, y,
     +                   w4,  hptr,  fuval,  nefval,  inform )
          else
*
*           The element Hessian has an elemental representation.
*
            call fdelhf( iel, elptr, elvar, xst, x, w4(arcs+1),
     +                   s,  hptr,  fuval,  nefval,  inform )
          endif
          if( inform.ne.0 ) return
*
*         Element function evaluations.
*
        endif
        ifc = ifc + 1
        if ( ifc .le. nrec ) go to 50
      endif
*
*     The element Hessian matrices are updated by a quasi-Newton formula
*     or Dembo's preconditioner needs computing
*
      if( qnd ) then
        sqrem = sqrt(epsmch)
        nbfgs = 0
        nrk1  = 0
        noupd = 0
*
*       We compute the differences between the
*       iterate X and the previous iterate.
*
        do 70 k = 1 , arcs
          dx(k) = x(k) - dx(k)
 70     continue
*
*       We compute the differences between the gradient
*       evaluated at point X and the gradient evaluated
*       at the previous iterate.
*
        ij = gptr(1)
        do 75 i = 1 , lg
          olgra(i) = fuval(ij) - olgra(i)
          ij = ij + 1
 75     continue
*
*       We loop on the element functions
*       that need to be updated.
*
        do 80 ifc = 1 , nrec
          iel   = recalc(ifc)
          els   = mod( elst(iel), 8 )
*
*         Dembo's strategy without preconditioning.
*         In this case, nothing needs updating.
*
          if( els.eq.7 .and. .not.prec ) goto 80
*
          eldim = elptr(iel+1) - elptr(iel)
          debx  = elptr(iel)
          debg  = gptr(iel) - elem
          debh  = hptr(iel)
*
*         We obtain the dimension of the
*         element Hessian matrix.
*
          if( elint(elst(iel)) ) then
            call range( iel, eldim, intdim, s, y , 0 )
            dim = intdim
          else
            dim = eldim
          endif
*
*         We compute the vectors S and Y.
*
          if( dim .lt. eldim ) then
            if( els.eq.7 .and. prec ) then
              call dcopy( dim, olgra(debg), 1, s, 1 )
              call range( iel, eldim, intdim, s, y, 2 )
              call gathr0( eldim, s, 1, dx, elvar(debx) )
            else
              call gathr0( eldim, y, 1, dx, elvar(debx) )
              call range( iel, eldim, intdim, y, s, 1 )
              call dcopy( dim, olgra(debg), 1, y, 1 )
            endif
          else
            call gathr0( dim, s, 1, dx, elvar(debx) )
            call dcopy( dim, olgra(debg), 1, y, 1 )
          endif
*
*         We test if this element Hessian matrix needs
*         really to be updated.
*
          newp = .false.
          do 90 i = 1 , dim
            newp = newp .or. abs(s(i)).gt.epsmch
 90       continue
*
          if( newp ) then
*
*           The element Hessian matrix needs
*           to be updated. We compute its
*           curvate estimate.
*
            yts = ddot( dim, s, 1, y, 1 )
*
*           The Hessian matrix is set to a scaling of
*           the identity matrix before its is first
*           updated by a quasi-Newton formula.
*
            if( els.eq.6 .and. .not.elhqn(elst(iel)) ) then
              sbs = ddot( dim, s, 1, s, 1 )
              if( sbs.gt.sqrem ) then
                alpha = max( sqrem , abs(yts / sbs) )
              else
                alpha = one
              endif
              call hessid( dim, fuval(debh), alpha )
              call selhqn( elst(iel) )
            endif
*
*           When Dembo's preconditioner is used,
*           we update the diagonal of the Hessian
*           matrix and we store it in array FUVAL.
*
            if( els.eq.7 .and. prec ) then
              debh = arcs + lg + elptr(iel)
              do 95 ij = 1 , eldim
                if( abs(s(ij)).gt.sqrem ) then
                  dpr(debh) = y(ij) / s(ij)
                else
                  dpr(debh) = one
                endif
                debh = debh + 1
 95           continue

            else if( dim.eq.1 ) then
*
*            The element function has
*            only one variable.
*
             fuval(debh) = y(1) / s(1)
*
            else if( elcvx(elst(iel)) ) then
*
*             We test if the element function is convex.
*
              if( yts.le.zero ) go to 110
*
*             The BFGS update is applied.
*
              call scbfgs( dim, fuval(debh), s, y, epsmch,
     +                     epsmch,  negcur,  w4 )
*
*             We test if a negative curvate has been observed.
*
              if( negcur ) then
*
*               A negative curvate has been detected.
*               We will used the Rank-One update.
*
                go to 110
              else
                nbfgs = nbfgs + 1
              endif
              else
*
*             The Rank-One update is applied.
*
 110          continue
              call scrk1( dim, fuval(debh), s, y, epsmch,
     +                    epsmch,  upd,  w4 )
*
              if( upd ) then
*
*               We update the status of the element
*               function to correspond to the status
*               of an element function that is non
*               convex.
*
                call snelcv( elst(iel) )
                nrk1 = nrk1 + 1
              else
*
*               No update has been done.
*
                noupd = noupd + 1
              endif
            endif
          endif
 80     continue
      endif
*
 900  continue
      inform = 0
*
      return
      end