* ************************************************************************* subroutine xpdbs( alpha, alpham, x, dir, ds, dbd, rgra, rd, + bfeas, ba, nba, su, nsu, ftstem, bep, + beptr, kdist ) * ************************************************************************* * Purpose : * --------- * Given the steplength ALPHA, this routine computes a basic- * superbasic feasible direction. The superbasic direction is * projected on the feasible set of constraints. * Parameters : * ------------ * alpha ( dble ) * input : the current value of the steplength. * output : this value is reduced when the basic * part of the direction is not feasible * with respect to the bounds. * alpham ( dble ) * input : the minimum value between APLHAS and ALPHAB. * (see LNSRCH subroutine). * output : unmodified. * x ( dble ) * input : the current iterate vector. * output : a new iterate is computed. * dir ( dble ) * input : the superbasic components of the direction. * output : the basic components of the projected direction * are computed. * ds ( dble ) * input : meaningless. * output : it contains the projection of the supoerbasic * direction on the constraints. * dbd ( dble ) * input : vector containing for each superbasic variable * the absolute difference between its value and * its active bound value. * output : unmodified. * rgra ( dble ) * input : the reduced gradient vector. * output : unmodified. * rd ( dble ) * input : the slope of the direction. * output : this value is updated. * bfeas ( log ) * input : .true. iff ALPHA <= ALPHAB. * In this case, the basic part of * the direction is always feasible. * .false. otherwise. * output : it is set to .true. when ALPHA becomes * <= ALPHAB. * Otherwise, it remains unmodified. * ba ( int ) * input : vector containing the indices of the basic * variables. * output : unmodified. * nba ( int ) * input : the number of basic variables. * output : unmodified. * su ( int ) * input : vector containing the indices of the superbasic * variables. * output : unmodified. * nsu ( int ) * input : the number of superbasic variables. * output : unmodified. * ftstem ( int ) * input : vector containing for each superbasic variable, * the length from its origine node to the stem * node, and the length of its flow augmenting path. * output : unmodified. * bep ( int ) * input : array containing the flow augmenting paths of * the superbasic variables. * output : unmodified. * beptr ( int ) * input : array whose kth value is the position of the * first element of the flow augmenting path of * superbasic arc k, in array BEP. * output : unmodified. * kdist ( int ) * input : vector containing for each superbasic variable * the length of its flow augmenting path. * output : unmodified. * Routines used : * --------------- * min, max, sqrt. * dsetcd, active. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer ba(*), nba, su(*), nsu, ftstem(*), bep(*), + beptr(*), kdist(*) logical bfeas double precision alpha, alpham, x(*), dir(*), ds(*), + dbd(*), rgra(*), rd * Internal variables integer i, ik, k, kk, kdis, itfeas logical nofeas double precision alphax, aa, bdact, active * Common specifications double precision epsmch, huge, tiny, tol 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. * itfeas = 0 1 continue itfeas = itfeas + 1 nofeas = .false. * * We compute a feasible basic-superbasic direction. * call dsetcd( nba, ba, dir, zero ) do 10 ik = 1 , nsu k = su(ik) * * We project the superbasic direction * on the feasible set of constraints. * if( bfeas ) then ds(k) = alpha * dir(k) else if( dir(k).gt.zero ) then ds(k) = min( alpha*dir(k) , dbd(k) ) else ds(k) = max( alpha*dir(k) , -dbd(k) ) endif endif * * We obtain the basic direction. * kdis = kdist(ik) do 20 i = beptr(k) , beptr(k)+kdis-1 kk = bep(i) if( kk.lt.0 ) then dir(-kk) = dir(-kk) - ds(k) else dir(kk) = dir(kk) + ds(k) endif 20 continue 10 continue * * BFEAS = .true. means that ALPHA <= ALPHAB. * In this case, the basic part of the direction * is always feasible. * if( .not.bfeas ) then * * We test if the basic part of the direction is * still feasible with respect to the bounds. * alphax = huge do 30 ik = 1 , nba k = ba(ik) if( abs(dir(k)).le.sqrt(epsmch) ) then dir(k) = zero else bdact = active(k,dir) aa = abs( (x(k)-bdact) / dir(k) ) if( aa.lt.one ) then nofeas = .true. alphax = min( alphax , aa ) endif endif 30 continue * * We test if the basic feasibility is verified. * if( nofeas ) then * * We reduce steplength ALPHA until obtaining * basic feasibility. * alpha = alphax * alpha if( itfeas.ge.3 ) alpha = alpha - tenm4 if( alpha .le. alpham ) then alpha = alpham bfeas = .true. endif * * We compute a new basic-superbasic direction. * go to 1 endif rd = zero endif * * We update the basic components * of current iterate X. * do 40 ik = 1 , nba k = ba(ik) x(k) = x(k) + dir(k) 40 continue * * We update the superbasic components * of current iterate X and the slope. * do 50 ik = 1 , nsu k = su(ik) x(k) = x(k) + ds(k) if( .not.bfeas ) then rd = rd + rgra(k)*ds(k) endif 50 continue * return end