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

      subroutine mbst( pr, thr, rthr, depth,  ba, nba,  nsu, su,  nb,
     +                 optmal,  ibding,   fr,  to,  x,   xst,   what,
     +                 ipdevc,  info,  iflag )

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

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

*     This routine builds the initial maximal basis spanning tree
*     and the four data structure vectors PR, THR, RTHR and DEPTH.
*     The tree is obtained by using a Greedy algorithm.

*     The indices of the basic variables are stored in BA and the
*     indices of the superbasic variables are stored in SU.

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

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

*     ba      ( int )

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

*     nba     ( int )

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

*     nsu     ( int )

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

*     su      ( int )

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

*     nb      ( int )

*       input  : meaningless.
*       output : the number of nonbasic variables.

*     optmal  ( int )

*       input  : meaningless.
*       output : = 2, if the superbasic set is empty,
*                = 0, otherwise.

*     ibding  ( int )

*       input  : meaningless.
*       output : = 0, if the basis is fully free,
*                = 1, otherwise.

*     fr      ( int )

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

*     to      ( int )

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

*     x       ( dble )

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

*     xst     ( int )

*       input  : vector containing the status of the variables.
*       output : this vector is updated.

*     what    ( int )

*       input  : if WHAT <  2, no output is generated,
*                        >= 2, some informations about
*                              the basis are printed.
*       output : unmodified.

*     ipdevc  ( int )

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

*     info    ( int )

*       input  : meaningless.
*       output : on exit under an error, it contains information
*                about the error.

*     iflag   ( int )

*       input  : should be equal to zero when the routine is called.
*       output : = 10, if the graph is not connected.
*                = 15, if the total number of variables partitioned in
*                      the fixed, superbasic and nonbasic sets differs
*                      from the total number of variables of the
*                      problem.
*                = 16, if the number of basic variables is incorrect.
*                Otherwise, it remains unmodified.

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

*     xlower, xupper, gxfix, sxnbl, sxnbu, gxsu,
*     sxbdg, isetvl.

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

*     D. Tuyttens

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

*     Routine parameters

      integer           pr(*), thr(*), rthr(*), depth(*), ba(*), nba,
     +                  su(*), nsu, nb, optmal, ibding, fr(*), to(*),
     +                  xst(*), what, ipdevc, info, iflag
      double precision  x(*)

*     Internal variables

      integer           ns, nsaux, nf, k, i, nt, frk, ttk,  mxt, mxp,
     +                  mxa, mxtf, mxpf, mxaf,  nd1,  ison,  end, ik,
     +                  mikf, mik
      logical           gxfix, gxsu
      double precision  xlower, xupper, lb, ub

*     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
*
*     Some initializations.
*
      nb  = 0
      nsu = 0
      nf  = 0
      end = arcs
*
*     The vector SU is organized in such a way that it
*     contains first the indices of the variables that
*     are strictly between their bounds followed by the
*     indices of the bounded and the fixed variables.
*
      do 10 k = 1 , arcs
*
*       Is variable K fixed ?
*
        if( gxfix(xst(k)) ) then
          nf = nf + 1
          su(end) = k
          end = end - 1
        else
*
*         Is variable K lower bounded ?
*
          lb = xlower(k)
          if( x(k)-lb .le. tol*(one+abs(lb)) ) then
            call sxnbl(xst(k))
          endif
*
*         Is variable K upper bounded ?
*
          ub = xupper(k)
          if( ub-x(k) .le. tol*(one+abs(ub)) ) then
            call sxnbu(xst(k))
          endif
*
*         We add indice K in vector SU.
*
          if( gxsu(xst(k)) ) then
            nsu     = nsu + 1
            su(nsu) = k
          else
            nb = nb + 1
            call sxbdg(xst(k))
            su(end) = k
            end = end - 1
          endif
        endif
*
 10   continue
*
*     We test if the number of fixed, non bounded and bounded
*     variables is equal to the total number of variables.
*
      if( nf+nsu+nb .ne. arcs ) then
        info  = nf+nsu+nb
        iflag = 15
        return
      endif
*
*     Initialisations of the vectors PR, THR, RTHR, DEPTH
*     and BA.
*
      nba = 0
      nd1 = nodes + 1
      call isetvl( nd1, pr, 1, 0 )
      pr(nd1)    = nd1
      depth(nd1) = 0
      thr(nd1)   = 1
      rthr(1)    = nd1
      thr(1)     = nd1
      rthr(nd1)  = 1
      depth(1)   = 1
      pr(1)      = nd1
      ba(1)      = 0
      optmal     = 0
      ibding     = 0
*
*     We build the maximal basis spanning tree by using
*     a Greedy algorithm.
*
      do 30 nt = 3 , nd1
*
*       Reset some variables to zero.
*
        mxt  = 0
        mxp  = 0
        mxa  = 0
        mxtf = 0
        mxpf = 0
        mxaf = 0
*
*       Loop on the arcs and try to find a new candidate
*       to enter the tree.
*
        do 40 ik = 1 , arcs
          k = su(ik)
*
*         Is arc K already in the tree ?
*
          if( k.lt.0 ) goto 40
          frk = fr(k)
          ttk = to(k)
*
          if( pr(frk).ne.0 ) then
            if( pr(ttk).eq.0 ) then
*
*             Arc K is a new candidate to enter the tree.
*
              if( gxfix(xst(k)) ) then
                mxtf = ttk
                mxpf = frk
                mxaf = k
                mikf = ik
              else
                mxt  = ttk
                mxp  = frk
                mxa  = k
                mik  = ik
                if( ik.le.nsu ) goto 100
              endif
            endif
*
          else
*
            if( pr(ttk).ne.0 ) then
*
*             Arc K is a new candidate to enter the tree.
*
              if( gxfix(xst(k)) ) then
                mxtf = frk
                mxpf = ttk
                mxaf = k
                mikf = ik
              else
                mxt  = frk
                mxp  = ttk
                mxa  = k
                mik  = ik
                if( ik.le.nsu ) goto 100
              endif
            endif
          endif
*
*       We try next arc.
*
 40     continue
*
*       We test if the graph is not connected.
*
        if( mxa.eq.0 .and. mxaf.eq.0 ) then
          iflag = 10
          return
        else
*
          if( mxa.eq.0 ) then
*
*           A fixed arc enters the tree.
*
            mxt = mxtf
            mxp = mxpf
            mxa = mxaf
            ik  = mikf
            nf  = nf - 1
          else
*
*           A bounded variable enters the tree.
*
            ik  = mik
            nb  = nb - 1
            ibding = 1
          endif
*
*         A new arc enters the tree. We update the
*         vectors PR, THR, RTHR, DEPTH and BA.
*
 100      continue
          pr(mxt)    = mxp
          depth(mxt) = depth(mxp) + 1
          ison       = thr(mxp)
          thr(mxp)   = mxt
          rthr(mxt)  = mxp
          thr(mxt)   = ison
          rthr(ison) = mxt
          call sxba(xst(mxa))
          call sxndis(xst(mxt))
          nba = nba + 1
          ba(mxt) = mxa
          su(ik)  = -su(ik)
        endif
*
 30   continue
*
*     We set in SU the indices of the remaining
*     superbasic variables.
*
      nsaux = nsu
      nsu   = 0
      ns    = 0
      do 50 i = 1 , arcs
        if( su(i).gt.0 ) then
          if( i.le.nsaux ) nsu = nsu + 1
          ns = ns + 1
          su(ns) = su(i)
        endif
 50   continue
*
*     We test if the number of fixed, non bounded, bounded and
*     basic variables is equal to the total number of variables.
*
      if( nf+nsu+nb+nba .ne. arcs ) then
        info  = nba
        iflag = 16
        return
      endif
*
*     We test if the maximal basis is fully free.
*
      if( what.ge.2 ) then
        if( ibding.eq.0 ) then
          write(ipdevc,1003) nb
 1003     format(' BASIS FULLY FREE  - ',I4,' NONBASIC VARIABLES.')
        else
          write(ipdevc,1004) nb
 1004     format(' BASIS NOT FULLY FREE  - ',I4,' NONBASIC VARIABLES.')
        endif
      endif
*
*     Is the superbasic set empty ?
*
      if( nsu.eq.0 ) then
        optmal = 2
        if( what.ge.2 ) write(ipdevc,1006)
 1006   format(/' SUPERBASIC SET IS EMPTY. ')
      endif
*
      return
      end