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

      subroutine dembo( inform,  elptr, elvar, elst, gptr, lg, negval,
     +                  fcalc, lfc, fuval, x, xold, ds, gds, dw1, dw2,
     +                  bgx  )

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

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

*     This routine computes the product of a partially separable
*     matrix times the vector ds, by evaluating differences in
*     the element gradient values. The result is stored in array
*     gds.

*                          g(x + sigma ds) - g(x)
*                  gds  =  ----------------------
*                                   sigma

*     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 gradient
*                values are needed to proceed the calculation. A 
*                return is then made to the main program.
*                Otherwise, no further information is needed.

*     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.

*     elst    ( int )

*       input  : vector containing the status of the element
*                functions.
*       output : unmodified.

*     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.

*     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.

*     fcalc   ( int )

*       input  : this vector contains the indices of the element
*                Hessian matrices involved in the matrix vector
*                product.
*       output : unmodified.

*     lfc     ( int )

*       input  : the number of element Hessian matrices involved
*                in the matrix vector product.
*       output : unmodified.

*     fuval   ( dble )

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

*     x       ( dble )

*       input  : the current iterate vector.
*       output : it contains the vector values at which the element
*                gradients need to be evaluated.

*     xold    ( dble )

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

*     ds      ( dble )

*       input  : the vector we want to multiply times the
*                Hessian matrix.
*       output : unmodified.

*     gds     ( dble )

*       input  : meaningless
*       output : it contains the product of the Hessian matrix
*                times the vector ds.

*     dw1     ( dble )

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

*     dw2     ( dble )

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

*     bgx     ( dble )

*       input  : array with the gradient at x.
*       output : unmodified (only used for Dembo's strategy).

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

*     sqrt.
*     elint, range, scattr.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           inform, elptr(*), elvar(*), elst(*),  gptr(*), 
     +                  lg,  negval,  fcalc(*),  lfc
      double precision  fuval(*),   x(*),   xold(*),  ds(*),   gds(*),
     +                  dw1(*), dw2(*), bgx(*)

*     Internal variables

      integer           i, j, ij, iel, eldim, intdim, debg, debg2, ifc,
     +                  lx
      logical           elint
      double precision  sigma

*     Save variables

      save              i, j, ij, iel, eldim, intdim, debg, debg2, ifc,
     +                  lx
      save              sigma

*     Common specifications

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

*
*     Reverse communication test.
*
      if( inform.ne.0 ) goto 100
*
*     We compute the vector   x + sigma ds.
*
      lx = 0
      sigma = sqrt(epsmch)
      do 10 ifc = 1 , lfc
        iel = fcalc(ifc)
        if( iel .lt. 0 ) then
          lx        = lx + 1
          iel       = -iel
          fcalc(lx) = iel
          do 20 j = elptr(iel) , elptr(iel+1)-1
            ij    = elvar(j)
            x(ij) = xold(ij) + sigma*ds(ij)
 20       continue
        endif
 10   continue
*
*     We evaluate the gradient at point  x + sigma ds.
*
      lfc    = lx
      negval = negval + lfc
      inform = 28
      return
*
*     Element gradient evaluations.
*
 100  continue
      inform = 0
*
*     We loop on the element Hessians involved
*     in the matrix vector product.
*
      do 30 ifc = 1 , lfc
        iel   = fcalc(ifc)
        eldim = elptr(iel+1) - elptr(iel)
        debg  = gptr(iel)-1
        debg2 = debg - elem
        if( eldim.eq.1 ) then
*
*         The one dimensional case.
*
          ij      = elvar(elptr(iel))
          gds(ij) = gds(ij) + (fuval(debg+1)-bgx(debg2+1)) / sigma
        else
          if( elint(elst(iel)) ) then
*
*           The element gradient has an internal representation.
*
            call range( iel, eldim, intdim, dw1, dw2, 0 )
            do 40 i = 1 , intdim
              dw1(i) = (fuval(debg+i)-bgx(debg2+i)) / sigma
 40         continue
            call range( iel, eldim, intdim, dw1, dw2, 2 )
          else
*
*           The element gradient has an elemental representation.
*
            do 50 i = 1 , eldim
              dw2(i) = (fuval(debg+i)-bgx(debg2+i)) / sigma
 50         continue
          endif
*
*         The result of the computation is accumulated in array gds.
*
          call scattr( eldim, dw2, 1, gds, elvar(elptr(iel)) )
        endif
 30   continue
*
      return
      end