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

      subroutine fdinhg( iel, elptr, elvar, x, w1, w2, w3, w4,
     +                   fuval,  gptr,  hptr,  negval, inform )

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

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

*     This routine estimates the IELth element Hessian matrix
*     by forward differences in the element gradient values.
*     This IELth element Hessian has an internal representation.

*     Every time information on the element gradient value is needed,
*     a return is made to the main program. This routine is then
*     re-entered with the information available and the calculation
*     proceeds.

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

*     iel     ( int )

*       input  : the indice of the element Hessian matrix
*                that will be estimated.
*       output : unmodified.

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

*     x       ( dble )

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

*     w1      ( dble )

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

*     w2      ( dble )

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

*     w3      ( dble )

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

*     w4      ( dble )

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

*     fuval   ( dble )

*       input  : array used to store the function and derivative
*                values of the element functions.
*       output : the components of indices HPTR(IEL) up to
*                HPTR(IEL+1)-1 will be set to the estimate
*                of the element Hessian matrix by forward
*                differences in the element gradient values.

*     gptr    ( int )

*       input  : array whose ith value is the position of the
*                first component of the ith element gradient
*                in FUVAL.
*       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.

*     negval  ( int )

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

*     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 re-entered
*                with an element gradient value that was needed
*                to proceed the calculation.
*       output : If it is not equal to zero, an element gradient
*                value is needed to proceed the calculation. A return
*                is then made to the main program.
*                Otherwise, no further information is needed.

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

*     sqrt.
*     range, dcopcd, dcopy, dsetvl, dxpycd.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           iel, elptr(*), elvar(*), gptr(*),
     +                  hptr(*), negval, inform
      double precision  x(*), w1(*), w2(*), w3(*), w4(*),
     +                  fuval(*)

*     Internal variables

      integer           i, j, ij, l, eldim, intdim, debg,
     +                  debh
      double precision  sqrem

*     Save variables

      save              i, j, ij, l, eldim, intdim, debg
      save              debh
      save              sqrem

*     Common specifications

      double precision  epsmch, huge, tiny, tol
      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 ) goto 100
*
*     Estimate of the IELth element Hessian matrix by
*     forward differences in the element gradient values.
*
      sqrem = sqrt(epsmch)
      debg  = gptr(iel)
      debh  = hptr(iel)
      eldim = elptr(iel+1) - elptr(iel)
      call range( iel, eldim, intdim, w1, w1, 0 )
      call dcopcd( eldim, elvar(elptr(iel)), x, w3 )
      call dcopy( intdim, fuval(debg), 1, w4, 1 )
      call dsetvl( intdim, w1, 1, zero )
*
      i = 1
 10   continue
      w1(i) = sqrem
      call range( iel, eldim, intdim, w1, w2, 3 )
      call dxpycd( eldim, elvar(elptr(iel)), x, w2 )
      negval = negval + 1
      inform = 5
      return
*
*     Element gradient evaluation.
*
 100  continue
      inform = 0
      call dcopcd( eldim, elvar(elptr(iel)), w3, x )
      l = debh
      do 20 j = 1 , i-1
        l  = l - intdim + i - j
        ij = i - j
        fuval(l) = half*(fuval(l) + (fuval(debg+ij-1)-w4(ij))/sqrem)
 20   continue
      l = debh
      do 30 j = i , intdim
        fuval(l) = (fuval(debg+j-1)-w4(j))/sqrem
        l = l + 1
 30   continue
      debh = debh + intdim - i + 1
      i = i + 1
      if( i .le. intdim ) go to 10
      call dcopy( intdim, w4, 1, fuval(debg), 1 )
*
      return
      end