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

      subroutine prod( x, fuval, gptr, lg, negval, hptr, fcalc, lfc,
     +                 inform,  curv,  ds,  gds,  dw1,  dw2,  elvar,
     +                 elptr, elst,  ba, nba,  su, nsu,  bep, beptr,
     +                 ftstem, kdist, bgx )

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

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

*     This routine computes the product of the reduced Hessian
*     matrix times the conjugate gradient direction ds and store
*     it in array gds.

*                                t
*                       gds  =  Z G Z ds

*     This matrix vector product can be obtained by calling a
*     special procedure that computes the product of a symmetric
*     partially separable matrix times a vector, or by computing
*     differences in the element gradient values.

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

*     x       ( dble )

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

*     fuval   ( dble )

*       input  : array used to store the function and derivative
*                values 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.

*     hptr    ( int )

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

*     fcalc   ( int )

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

*     lfc     ( int )

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

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

*     curv    ( dble )

*       input  : meaningless.
*       output : the curvature of the conjugate gradient
*                direction ds.

*     ds      ( dble )

*       input  : the superbasic components of the conjugate
*                gradient direction.
*       output : the basic corresponding components of the
*                conjugate gradient direction are computed.

*     gds     ( dble )

*       input  : meaningless
*       output : it contains the product of the reduced Hessian
*                matrix times the conjugate gradient direction ds.

*     dw1     ( dble )

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

*     dw2     ( dble )

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

*     elvar   ( int )

*       input  : array containing the indices of the varaibles
*                in the first element, followed by those in the
*                second element, etc.
*       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.

*     elst    ( int )

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

*     ba      ( int )

*       input  : vector containing the indices of the basic
*                variables.
*       output : unmodified.

*     nba     ( int )

*       input  : the number of basic variables.
*       output : unmodified.

*     su      ( int )

*       input  : vector containing the indices of the superbasic
*                variables.
*       output : unmodified.

*     nsu     ( int )

*       input  : the number of superbasic variables.
*       output : unmodified.

*     bep     ( int )

*       input  : array containing the flow augmenting paths of
*                the superbasic variables.
*       output : unmodified.

*     beptr   ( int )

*       input  : array whose kth value is the position of the
*                first element of the flow augmenting path of
*                superbasic arc k, in array BEP.
*       output : unmodified.

*     ftstem  ( int )

*       input  : vector containing for each superbasic variable,
*                the length from its origine node to the stem
*                node, and the length of its flow augmenting path.
*       output : unmodified.

*     kdist   ( int )

*       input  : vector containing for each superbasic variable
*                the length of its flow augmenting path.
*       output : unmodified.

*     bgx     ( dble )

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

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

*     mod.
*     dsetcd, dsetvl, dembo, psprdt, ddotcd.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           elvar(*),  elptr(*),  gptr(*),  lg,  negval,
     +                  hptr(*),  elst(*),  ba(*), nba,  su(*), nsu,
     +                  bep(*),   beptr(*),    ftstem(*),    inform,
     +                  fcalc(*), lfc,  kdist(*)
      double precision  x(*), fuval(*), curv, ds(*), gds(*), dw1(*),
     +                  dw2(*), bgx(*)

*     Internal variables

      integer           ik, k, kk,  kdis, i, j,  iel, ndmb
      double precision  ddotcd

*     Save variables

      save              ik, k, kk,  kdis, i, j,  iel, ndmb

*     Common specifications

      integer           arcs, nodes, elem
      common / prbdim / arcs, nodes, elem

      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
*
*     On input, vector ds contains only the superbasic
*     components of the conjugate gradient direction.
*     We now compute the corresponding basic components
*     of the conjugate gradient direction and store them
*     in vector ds (i.e.  ds = Z ds ).
*
      call dsetcd( nba, ba, ds, zero )
      do 10 ik = 1 , nsu
        k = su(ik)
        kdis = kdist(ik)
*
*       We compute the corresponding basic components
*       of the conjugate gradient direction.
*
        do 20 i = beptr(k) , beptr(k)+kdis-1
          kk = bep(i)
          if( kk.lt.0 ) then
            ds(-kk) = ds(-kk) - ds(k)
          else
            ds(kk)  = ds(kk)  + ds(k)
          endif
 20     continue
 10   continue
*
*     We search the indices of the element Hessians
*     that are involved in the matrix vector product.
*
      lfc  = 0
      ndmb = 0
      do 30 iel = 1 , elem
        do 40 j = elptr(iel) , elptr(iel+1)-1
          if( ds(elvar(j)).ne.zero ) then
            lfc = lfc + 1
            if( mod(elst(iel),8) .eq. 7 ) then
              ndmb       = ndmb + 1
              fcalc(lfc) = -iel
            else
              fcalc(lfc) = iel
            endif
            go to 30
          endif
 40     continue
 30   continue
*
*     We compute the product  G ds = G (Z ds),
*     and we store it in array  gds.
*
      call dsetvl( arcs, gds, 1, zero )
*
*     The matrix vector product is obtained by calling a special
*     procedure that computes the product of a symmetric partially
*     separable matrix times a vector.
*
      if ( ndmb .lt. lfc ) then
        call psprdt( fcalc, lfc, fuval,  hptr,  ds, gds,  elvar,
     +               elptr,  elst,  dw1,  dw2 )
      endif
*
*     The matrix vector product is obtained by computing
*     differences in the element gradient values.
*
      if ( ndmb .gt. 0 ) then
 100    continue
        call dembo( inform, elptr, elvar, elst, gptr, lg, negval,
     +              fcalc, lfc, fuval,  x, dw1(arcs+1),  ds, gds,
     +              dw1, dw2, bgx(arcs+1) )
        if( inform.ne.0 ) return
      endif
*
*                               t          t
*     We compute the product   Z G Z ds = Z gds
*     and we store it in array  gds.
*
      do 50 ik = 1 , nsu
        k = su(ik)
        kdis = kdist(ik)
        do 60 i = beptr(k) , beptr(k)+kdis-1
          kk = bep(i)
          if( kk.lt.0 ) then
            gds(k) = gds(k) - gds(-kk)
          else
            gds(k) = gds(k) + gds(kk)
          endif
 60     continue
 50   continue
*
*     We compute the curvate of the conjugate
*     gradient direction ds.
*
      curv = ddotcd( nsu, su, gds, ds )
*
      return
      end