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

      subroutine bcycle( x, fr, tt, pr, depth, lib, bep, ncycle,
     +                   rgmx, alpha, pivot, outarc, isa, inode,
     +                   jnode,  nc,  iflag )

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

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

*     This routine checks whether or not the flow can be pushed
*     around the cycle generated by the nonbasic arc ISA.
*     It computes the maximal steplength ALPMAX for feasibility.
*     It also stores the flow augmenting path of the nonbasic
*     arc ISA in array BEP.

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

*     x       ( dble )

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

*     fr      ( int )

*       input  : vector containing the origine nodes of the arcs.
*       output : unmodified.

*     tt      ( int )

*       input  : vector containing the end nodes of the arcs.
*       output : unmodified.

*     pr      ( int )

*       input  : the predecessor vector.
*       output : unmodified.

*     depth   ( int )

*       input  : the depth vector.
*       output : unmodified.

*     lib     ( int )

*       input  : vector containing the indices of the basic
*                variables. Its ith component corresponds to
*                the indice of the basic arc having the
*                following endnodes i and PR(i).
*       output : unmodified.

*     bep     ( int )

*       input  : meaningless.
*       output : this array contains the flow augmenting
*                path of the nonbasic arc ISA.

*     ncycle  ( int )

*       input  : meaningless.
*       output : the length of the flow augmenting path
*                of the nonbasic arc ISA.

*     rgmx    ( dble )

*       input  : the reduced cost of the nonbasic arc ISA.
*       output : unmodified.

*     alpha   ( dble )

*       input  : meaningless.
*       output : the maximum steplength allowed for keeping
*                feasibility. When ALPHA = 0.0, no change
*                of flow is possible.

*     pivot   ( log )

*       input  : meaningless.
*       output : .true.  iff a pivoting is needed,
*                .false. otherwise.

*     outarc  ( int )

*       input  : meaningless.
*       output : the indice of the basic arc candidate for
*                pivoting.

*     isa     ( int )

*       input  : the indice of the nonbasic arc.
*       output : unmodified.

*     inode   ( int )

*       input  : meaningless.
*       output : one of the endnodes of the nonbasic arc ISA.

*     jnode   ( int )

*       input  : meaningless.
*       output : the other endnode of the nonbasic arc ISA.
*                The path connecting nodes JNODE and NC is
*                the pivot' path.

*     nc      ( int )

*       input  : meaningless.
*       output : one of the endnodes of the basic arc candidate
*                for pivoting. The other endnode is PR(NC).

*     iflag   ( int )

*       input  : should be equal to zero when the routine is
*                called.
*       output : = 9, iff the linear problem is unbounded,
*                Otherwise, it remains unmodified.

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

*     max, min, abs.
*     xlower, xupper, mxflow.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           fr(*), tt(*), lib(*), pr(*), depth(*),
     +                  bep(*), ncycle,  outarc,  isa,  inode,
     +                  jnode,  nc,  iflag
      double precision  x(*),  alpha,  rgmx,  xlower,  xupper
      logical           pivot

*     Internal variables

      integer           iarc,  kifr,  kit,  isig, isig1,
     +                  ifr,  it,  ndif,  ndit,  i,  imx,  imn,
     +                  kimx, kimn, ndmx, ndmn, nxtimx, nxtimn
      double precision  upflow,  mxflow,  bigfl,  tolrf

*     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
*
*     Maximum steplength ALPHA for keeping
*     feasibility on nonbasic arc ISA.
*
      if( rgmx.lt.zero ) then
        alpha = xupper(isa) - x(isa)
      else
        alpha = x(isa) - xlower(isa)
      endif
*
*     Some initializations.
*
      tolrf  = epsmch
      bigfl  = one + max( alpha, tolrf )
      outarc = isa
      pivot  = .false.
      kifr   = -1
      kit    = 1
      ifr    = fr(isa)
      it     = tt(isa)
*
*     We start by tracing the cycle from the node
*     with the biggest depth.
*
      ndif = depth (ifr)
      ndit = depth (it)
      ndmx = max( ndif , ndit )
      ndmn = min( ndif , ndit )
*
      if( ndif.ge.ndit ) then
        imx  = ifr
        kimx = kifr
        imn  = it
        kimn = kit
        isig = -1
      else
        imx  = it
        kimx = kit
        imn  = ifr
        kimn = kifr
        isig = 1
      endif
*
      ncycle = 0
*
*     We climb the deepest branch until finding a
*     node IMX having the same depth as node IMN.
*
      do 10 i = 1 , ndmx-ndmn
*
*       We obtain the indice IARC of the basic arc.
*       Its orientation in the cycle can be deduced
*       from the sign of NXTIMX.
*
        iarc   = lib(imx)
        nxtimx = pr(imx) * kimx
        if( imx.ne.fr(iarc) ) nxtimx = -nxtimx
*
*       We compute the maxixmum steplength for keeping
*       feasibility on basic arc IARC.
*
        upflow = mxflow( x(iarc), xupper(iarc), xlower(iarc),
     +                   iarc,  nxtimx,  rgmx,  bigfl )
*
        if( upflow.le.tolrf ) then
*
*         The basic arc IARC is blocking.
*
          alpha  = zero
          pivot  = .true.
          outarc = iarc
          isig1  = -isig
          nc     = imx
          go to 200
        endif
*
        if( alpha.ge.upflow ) then
*
*         The basic arc IARC has the smallest maximum
*         steplength for keeping feasibility.
*
          alpha  = upflow
          pivot  = .true.
          outarc = iarc
          isig1  = -isig
          nc     = imx
        endif
*
*       We store the basic arc IARC and its orientation
*       in array BEP.
*
        imx    = abs(nxtimx)
        ncycle = ncycle + 1
        if( nxtimx.gt.0 ) then
          bep(ncycle) = iarc
        else
          bep(ncycle) =-iarc
        endif
*
 10   continue
*
*     We climb the two branches of the tree at the same time
*     until reaching the stem node.
*
 100  if( imx .ne. imn ) then
*
*       First branch of the tree.
*
*       We obtain the indice IARC of the basic arc.
*       Its orientation in the cycle can be deduced
*       from the sign of NXTIMX.
*
        iarc   = lib(imx)
        nxtimx = pr(imx) * kimx
        if( imx.ne.fr(iarc) ) nxtimx = -nxtimx
*
*       We compute the maxixmum steplength for keeping
*       feasibility on basic arc IARC.
*
        upflow = mxflow( x(iarc), xupper(iarc), xlower(iarc),
     +                   iarc,  nxtimx,  rgmx,  bigfl )
*
        if( upflow.le.tolrf) then
*
*         The basic arc IARC is blocking.
*
          alpha  = zero
          pivot  = .true.
          outarc = iarc
          isig1  = -isig
          nc     = imx
          go to 200
        endif
*
        if( alpha.ge.upflow ) then
*
*         The basic arc IARC has the smallest maximum
*         steplength for keeping feasibility.
*
          alpha  = upflow
          pivot  = .true.
          outarc = iarc
          isig1  = -isig
          nc     = imx
        endif
*
*       We store the basic arc IARC and its orientation
*       in array BEP.
*
        imx    = abs(nxtimx)
        ncycle = ncycle + 1
        if( nxtimx.gt.0 ) then
          bep(ncycle) = iarc
        else
          bep(ncycle) =-iarc
        endif
*
*       Second branch of the tree.
*
*       We obtain the indice IARC of the basic arc.
*       Its orientation in the cycle can be deduced
*       from the sign of NXTIMX.
*
        iarc   = lib(imn)
        nxtimn = pr(imn) * kimn
        if( imn.ne.fr(iarc) ) nxtimn = -nxtimn
*
*       We compute the maxixmum steplength for keeping
*       feasibility on basic arc IARC.
*
        upflow = mxflow( x(iarc), xupper(iarc), xlower(iarc),
     +                   iarc,  nxtimn,  rgmx,  bigfl )
*
        if( upflow.le.tolrf ) then
*
*         The basic arc IARC is blocking.
*
          alpha  = zero
          pivot  = .true.
          outarc = iarc
          isig1  = isig
          nc     = imn
          go to 200
        endif
*
        if( alpha.ge.upflow ) then
*
*         The basic arc IARC has the smallest maximum
*         steplength for keeping feasibility.
*
          alpha  = upflow
          pivot  = .true.
          outarc = iarc
          isig1  = isig
          nc     = imn
        endif
*
*       We store the basic arc IARC and its orientation
*       in array BEP.
*
        imn    = abs(nxtimn)
        ncycle = ncycle + 1
        if( nxtimn.gt.0 ) then
          bep(ncycle) = iarc
        else
          bep(ncycle) =-iarc
        endif
        go to 100
*
      endif
*
*     We test if the linear problem is unbounded.
*
 200  if( alpha.gt.huge ) then
*
*       The linear problem in unbounded.
*
        iflag = 9
        return
      else
*
*       The path connecting nodes JNODE to NC
*       is the pivot' path.
*
        inode = it
        jnode = ifr
        if( isig1.eq.-1 ) then
          inode = ifr
          jnode = it
        endif
      endif
*
      return
      end