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

      subroutine fdelhf( iel, elptr, elvar, xst, x, w1, w2,
     +                   hptr,  fuval,  nefval,  inform )

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

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

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

*     Every time information on the element function 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.

*     xst     ( int )

*       input  : vector containing the status of the variables.
*       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.

*     hptr    ( int )

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

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

*     nefval  ( int )

*       input  : the number of element function evaluations.
*       output : this number is increased by one each time
*                one element function 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 function value that was needed
*                to proceed the calculation.
*       output : If it is not equal to zero, an element function
*                value is needed to proceed the calculation. A return
*                is then made to the main program.
*                Otherwise, no further information is needed.

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

*     max, abs, sign.
*     gxfix, dsetvl.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           iel, elptr(*), elvar(*), xst(*),
     +                  hptr(*), nefval, inform
      double precision  x(*), w1(*), w2(*), fuval(*)

*     Internal variables

      integer           i, j, k, ik, kk, eldim, debh,
     +                  isout, ieout, isin, iein
      logical           gxfix
      double precision  thrdem, tempf, tempx, tempxx

*     Save variables

      save              i, j, k, ik, kk, eldim, debh,
     +                  isout, ieout, isin, iein
      save              thrdem, tempf, tempx, tempxx

*     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.eq.8  ) goto 100
      if( inform.eq.9  ) goto 200
      if( inform.eq.10 ) goto 300
*
*     Estimate of the IELth element Hessian matrix by
*     forward differences in the element function values.
*
      thrdem = epsmch**(one/three)
      eldim  = elptr(iel+1) - elptr(iel)
      tempf  = fuval(iel)
      ik = 0
      isout = elptr(iel)
      ieout = elptr(iel+1)-1
      if ( ieout .ge. isout ) then
        i = isout
 10     continue
        k  = elvar(i)
        ik = ik + 1
        if( gxfix(xst(k)) ) then
          w2(ik) = tempf
        else
          tempx  = x(k)
          w1(ik) = thrdem * sign( max(abs(tempx),one) , tempx )
          x(k)   = tempx + w1(ik)
          nefval = nefval + 1
          inform = 8
          return
*
*         Element function evaluation.
*
 100      continue
          x(k) = tempx
          inform = 0
          w2(ik) = fuval(iel)
        endif
        i = i + 1
        if ( i .le. ieout ) go to 10
      endif
*
      debh = hptr(iel)
      ik = 0
      if ( ieout .ge. isout ) then
        i = isout
 20     continue
        k  = elvar(i)
        ik = ik + 1
        if( gxfix(xst(k)) ) then
          call dsetvl( eldim-ik+1, fuval(debh), 1, zero )
        else
          tempx = x(k)
          x(k)  = tempx + w1(ik)
          isin = ik+1
          iein = eldim
          if ( iein .ge. isin ) then
            j = isin
 30         continue
            kk = elvar(elptr(iel)+j-1)
            if( gxfix(xst(kk)) ) then
              fuval(debh+j-ik) = zero
            else
              tempxx = x(kk)
              x(kk)  = tempxx + w1(j)
              nefval = nefval + 1
              inform = 9
              return
*
*             Element function evaluation.
*
 200          continue
              x(kk)  = tempxx
              inform = 0
              fuval(debh+j-ik) = ((tempf-w2(j))+(fuval(iel)-w2(ik))) /
     +                           (w1(ik)*w1(j))
            endif
            j = j + 1
            if ( j .le. iein ) go to 30
          endif
          x(k)   = tempx + 2.0d0*w1(ik)
          nefval = nefval + 1
          inform = 10
          return
*
*         Element function evaluation.
*
 300      continue
          x(k)   = tempx
          inform = 0
          fuval(debh) = ((tempf-w2(ik))+(fuval(iel)-w2(ik))) /
     +                  (w1(ik)*w1(ik))
        endif
        debh = debh + eldim - ik + 1
        i = i + 1
        if ( i .le. ieout ) go to 20
      endif
      fuval(iel) = tempf
*
      return
      end