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