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

      subroutine prcond( fuval, lg, hptr, diag, dpr, w1, w2, elptr,
     +                   elvar,  elst,   iss,  ip,  su,  nsu,  nqa,
     +                   ftstem, beptr,  bep,  kdist )

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

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

*     This routine computes the diagonal preconditioner that is
*     used in the preconditioned conjugate gradient method. This
*     preconditioner is stored in array dpr and is computed as
*     follows :

*                                   t
*                    dpr = | diag( Z diag(G) Z ) |
*

*     where diag(G) is the diagonal of matrix G.

*     When Dembo's method is used, the Hessian matrix is not
*     stored at all and the diagonal of this Hessian matrix
*     is estimated as follows :

*                                  (g(x) - g(x_))[i]
*                  (diag(G))[i] =  -----------------
*                                  (  x  -   x_ )[i]

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

*     fuval   ( dble )

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

*     diag    ( dble )

*       input  : meaningless.
*       output : it contains the diagonal elements of the
*                Hessian matrix, for the variables that are
*                included in the independent set IP.

*     dpr     ( dble )

*       input  : meaningless.
*       output : it contains the diagonal preconditioner
*                corresponding to independent set IP.

*     w1      ( dble )

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

*     w2      ( dble )

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

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

*     iss     ( int )

*       input  : vector containing for each variable, the indice
*                of the independent set in which it is included.
*                It is equal to zero when the variable does not
*                belongs to any independent set.
*       output : unmodified.

*     ip      ( int )

*       input  : the indice of the current independent set.
*       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.

*     nqa     ( int )

*       input  : the number of superbasic variables that are
*                quasi active.
*       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.

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

*     bep     ( int )

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

*     kdist   ( int )

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

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

*     mod, sqrt, abs.
*     dsetvl, elint, range, ltcmvp.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           lg, hptr(*),  elptr(*), elvar(*),  elst(*),
     +                  iss(*),  ip,  su(*), nsu,  nqa,  ftstem(*),
     +                  beptr(*),  bep(*),  kdist(*)
      double precision  fuval(*),  diag(*),  dpr(*),  w1(*),  w2(*)

*     Internal variables

      integer           i, j, jj, jp, ik, k, kk, iel, deb, intdim,
     +                  eldim, kdis, els
      logical           elint
      double precision  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
*
*     We initialize the diagonal values to zero.
*
      call dsetvl( arcs, diag, 1, zero )
*
*     We loop on the elements functions.
*
      do 10  iel = 1 , elem
*
*       We get some informations about
*       the element function IEL.
*
        els = mod(elst(iel),8)
        if( els.eq.7 ) then
          deb = arcs + lg
        else
          deb = hptr(iel)
        endif
        eldim = elptr(iel+1) - elptr(iel)
        i = 0
*
*       We loop on the variables that belong
*       to the element function IEL.
*
        do 20  j = elptr(iel), elptr(iel+1)-1
          i  = i + 1
          jj = elvar(j)
*
*         We test if the variable JJ is included
*         in the current independent set IP.
*
          if( iss(jj).eq.ip ) then
*
*           The variable JJ is included
*           in independent set IP.
*
            if( els.eq.7 ) then
*
*             Diagonal estimation of the Hessian matrix
*             when Dembo's method is applied.
*
              diag(jj) = diag(jj) + dpr(deb+j)
*
            else if( elint(elst(iel)) ) then
*
*             We estimate the JJ-th component of
*             the diagonal of the Hessian matrix.
*             The element Hessian matrix has an
*             internal representation.
*
              call range( iel, eldim, intdim, w1, w2, 0 )
              call dsetvl( eldim, w1, 1, zero )
              w1(i) = one
              call range( iel, eldim, intdim, w1, w2, 1 )
              w1(i) = zero
              call ltcmvp( intdim, fuval(deb), w2, w1,.true. )
              call range( iel, eldim, intdim, w1, w2, 2 )
              diag(jj) = diag(jj) + w2(i)
*
            else
*
*             We estimate the JJ-th component of
*             the diagonal of the Hessian matrix.
*             The element Hessian matrix has an
*             internal representation.
*
              if( eldim.eq.1 .or. i.eq.1 ) then
                diag(jj) = diag(jj) + fuval(deb)
              else
                jp = (i-1)*eldim - (( i**2 - 3*i + 2 ) / 2 )
                diag(jj) = diag(jj) + fuval(deb+jp)
              endif
            endif
          endif
 20     continue
 10   continue
*
*     We obtain the diagonal preconditioner
*     from the diagonal of the Hessian matrix.
*
      sqrem = sqrt(epsmch)
      do 30 ik = 1 , nsu-nqa
        k    = su(ik)
        kdis = kdist(ik)
        dpr(ik) = diag(k)
*
*       We go over the flow augmenting path
*       of the superbasic variable K.
*
        do 40 j = beptr(k) , beptr(k)+kdis-1
          kk      = abs(bep(j))
          dpr(ik) = dpr(ik) + diag(kk)
 40     continue
*
*       The components of the diagonal preconditioner
*       must be positive and not too small. Otherwise,
*       they are set to one.
*
        dpr(ik) = abs(dpr(ik))
        if( dpr(ik).le.sqrem ) dpr(ik) = one
 30   continue
*
      return
      end