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

      subroutine simplx( fr, tt, x, xst, w,  pr, thr, rthr,
     +                   depth, lin, lib, bep, iflag, info,
     +                   ipdevc, what )

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

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

*     This routine tests if the flow vector X is feasible.

*     If it is not feasible, an "all artificial Phase 1" is used
*     to find a feasible flow vector X. This procedure introduces
*     artificial arcs. Their orientation and flows are chosen to
*     enforce the feasibility of the extended flow vector and the
*     nonnegativity of their associated variables. The following
*     linear network problem is considered :

*                           a         a
*                 Min  f(x,x ) = sum x
*                            a
*                 s.t. Ax + x  = b
*
*                      l <= x <= u
*                            a
*                      0 <= x

*             a
*     If f(x,x ) = 0 at the solution of this problem, then the
*                    flow vector X is a feasible flow vector for
*                    the original problem.
*                > 0 at the solution of this problem, then the
*                    original problem is infeasible.

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

*     fr      ( int )

*       input  : vector containing the origine nodes of the
*                arcs.
*       output : it contains also the origine nodes of the
*                artificial arcs.

*     tt      ( int )

*       input  : vector containing the end nodes of the arcs.
*       output : it contains also the end nodes of the artificial
*                arcs.

*     x       ( dble )

*       input  : the flow vector. The routine tests if it is
*                feasible.
*       output : if it was not feasible on input, the routine
*                computes a feasible flow vector X (if the
*                problem is not infeasible).

*     xst     ( int )

*       input  : vector containing the status of the variables.
*       output : unmodified.

*     w       ( dble )

*       input  : array used as workspace.
*       output : it contains the dual variables.

*     pr      ( int )

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

*     thr     ( int )

*       input  : meaningless.
*       output : the traversal vector.

*     rthr    ( int )

*       input  : meaningless.
*       output : the reverse vector.

*     depth   ( int )

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

*     lin     ( int )

*       input  : meaningless.
*       output : vector containing the indices of the
*                nonbasic variables.

*     lib     ( int )

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

*     bep     ( int )

*       input  : array used as workspace.
*       output : this array is used to store the flow
*                augmenting paths of the nonbasic arcs.

*     iflag   ( int )

*       input  : should be equal to zero when the routine
*                is called.
*       output : = 4, if the bounds are inconsistent,
*                = 5, if the value of a fixed variable
*                     is not between the bounds.
*                = 7, if the maximum number of iterations
*                     has been performed before obtaining
*                     a feasible solution,
*                = 8, the problem is infeasible,
*                = 9, if the linear problem is unbounded,
*                = 14,if the value of the linear objective
*                     function does not decrease.
*                Otherwise, it remains unmodified.

*     info    ( int )

*       input  : meaningless.
*       output : when IFLAG > 0, it contains some
*                information about the error detected.

*     ipdevc  ( int )

*       input  : output device unit number for printing messages.
*       output : unmodified.

*     what    ( int )

*       input  : >= 2, the routine prints out the message saying
*                      if the flow is feasible or not.
*                <  2, no output is produced.
*       output : unmodified.

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

*     abs, max, min.
*     gxfix, xlower, xupper, testfs, basis, canli, bcycle,
*     updtfl, pivot.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           fr(*), tt(*),  lin(*),  lib(*),   xst(*),
     +                  pr(*), thr(*), rthr(*), depth(*), bep(*),
     +                  ipdevc, iflag, info, what
      double precision  x(*), w(*)

*     Internal variables

      integer           it, mxiter,  ndegn,  npts,  npoint, itmp,
     +                  outarc,  imove,   inode, jnode,  nc,  ns,
     +                  ncycle, i, k, nd1
      double precision  alpha,  rgmx,   obj,
     +                  xupper, xlower, rhs, sumbi, tolfsb
      logical           pivout, opt, feasbl
      logical           gxfix

*     Common specifications

      integer           arcs, nodes, elem
      double precision  epsmch, huge, tiny, tol
      integer           mfc, imfc
      common / prbdim / arcs, nodes, elem
      common / prbmch / epsmch, huge, tiny, tol
      common / mfcom  / mfc, imfc

      double precision  zero, one, two, three, half, tenm1, tenm2, tenm4
      common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4
*
*     Some initializations.
*
      if( mfc.eq.0 ) mfc = 1
      nd1    = nodes + 1
      iflag  = 0
      ndegn  = 0
      npts   = 0
      npoint = 0
      tolfsb = nd1 * epsmch * mfc * mfc
      alpha  = huge
      mxiter = 10000
      sumbi  = zero
*
*     We compute the error on the right hand side vector.
*
      do 20 i = 1 , nodes
        sumbi  = sumbi + rhs(i)
 20   continue
      sumbi  = abs(sumbi)
*
*     We check if the bounds are consistent, and we
*     set the flow vector X in the feasible domein.
*
      do 21 k = 1 , arcs
*
*       We test if the bounds are consistent.
*
        if( xlower(k).gt.xupper(k) ) then
          info  = k
          iflag = 4
          return
        else if( gxfix(xst(k)) ) then
*
*         We check if the fixed variable K is between the bounds.
*
          if( x(k).lt.xlower(k) .or. x(k).gt.xupper(k) ) then
            info  = k
            iflag = 5
            return
          endif
        else
*
*         We set the flow value of arc K in the feasible domein.
*
          x(k) = max( xlower(k) , min( x(k) , xupper(k)) )
        endif
 21   continue
*
*     We check if the flow vector X is feasible.
*
      call testfs( fr, tt, x, feasbl, tolfsb,
     +             obj, sumbi, ipdevc, what )
*
*     If the flow vector X is feasible, then nothing is to do.
*
      if( feasbl ) return
*
*     We build the first maximal basis spanning tree and we compute
*     the dual variables vector W.
*
      call basis( fr, pr, thr, rthr, depth, lib, lin, ns, w )

*      write(ipdevc,2221) sumbi,tolfsb,obj
* 2221 format(' sumbi = ',d12.5,' tolfsb = ',d12.5,' obj = ',d12.5)
*
*     The iteration loop.
*
      do 10 it = 1 , mxiter
*
*       We price out the nonbasics variables and
*       we find the nonbasic arc IMOVE candidate
*       to enter the basis.
*
        call canli ( npoint,  imove,  rgmx, opt,
     +               fr, tt, lin, ns, x, xst, w )
*
*       OPT = .true. means that the flow vector X
*                    is optimal for the linear problem.
*
        if( opt ) then
          go to 500
        endif
*
*       The routine BCYCLE tests if the nonbasic arc IMOVE
*       is blocked. It computes the maximal steplength ALPHA,
*       and it stores the flow augmenting path of the nonbasic
*       arc IMOVE in array BEP.
*
        itmp = lin(imove)
*
        call bcycle( x, fr, tt, pr, depth,  lib,  bep, ncycle,
     +               rgmx, alpha, pivout, outarc, itmp, inode,
     +               jnode,  nc,  iflag )
*
        if( iflag .ne. 0 ) return
*
*       ALPHA =  0.0      means that the nonbasic arc ITMP is blocked.
*                         We pivot it into the basis (degenerate pivot).
*       PIVOUT = .false.  means that the nonbasic arc ITMP hits a bound
*                         before a basic arc does, so no pivoting is
*                         required .
*
        if( alpha.ne.zero ) then
*
*         We update the flow vector X and the objective
*         function value OBJ.
*
          call updtfl( x, obj, alpha, rgmx, itmp,
     +                 bep, ncycle, iflag )

          if( iflag.ne.0 ) return
*
        else
*
*         Degenerate pivot.
*
          ndegn = ndegn + 1
        endif
*
*       We test if a pivot is required.
*
        if( pivout ) then
*
*         We pivot the basic arc OUTARC with the
*         nonbasic arc ITMP. The vector LIB is
*         also updated.
*
          call pivot( itmp, pr, thr, rthr, depth,
     +                lib, inode, jnode,  nc, w )

          npts = npts + 1
*
*         We update the vector LIN containing the
*         indices of the nonbasic variables. The
*         old basic variable OUTARC is added to the
*         nonbasic variables if it is a non artificial
*         variable.
*
          lin(imove) = outarc
          if( outarc.gt.arcs ) then
            lin(imove) = lin(ns)
            ns         = ns - 1
          endif
        endif
*
*       We test if the objective function value
*       is small enough to consider that the
*       optimal solution of the linear problem
*       is found and that the flow vector X is
*       feasible for the original problem.
*
*        write(ipdevc,2222) it,obj,rgmx,alpha,itmp,outarc
* 2222   format(' It = ',i4,' obj = ',d12.5,2x,d12.5,2x,d12.5,2x,i4,i4)

        if( obj.le.(sumbi + tolfsb) ) then
          go to 500
        endif
*
*     Next iteration.
*
 10   continue
*
*     The maximum number of iterations has been performed
*     before reaching an optimal solution.
*
      iflag = 7
      return
*
*     We verify the feasibility of the computed flow vector X.
*
500   continue

      call testfs( fr, tt, x, feasbl, tolfsb,
     +             obj, sumbi, ipdevc, what )
*
*     We test if the problem is infeasible.
*
      if( .not.feasbl ) then
        iflag = 8
      endif
*
      return
      end