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