* ************************************************************************* subroutine lnsrch( alpha, alphas, alphab, trial, dir, ds, x, + osol, rgra, rgnor, dbd, newp, fcalc, lfc, + fuval, fx, difval, oldfi, nefval, elvar, + elptr, xst, bound, ba, nba, su, nsu, nqa, + ftstem, bep, beptr, kdist, inform, ipdevc, + nflag, iflag ) * ************************************************************************* * Purpose : * --------- * Once the superbasic search direction DS is computed, this * routine computes a step of length ALPHA along this direction * such that the trial point is still feasible and the cost * function sufficiently reduced. * Since various bounds on the variables are present in the * calculation, the technique used in this linesearch procedure * is similar to that used by Bertsekas. * Parameters : * ------------ * alpha ( dble ) * input : meaningless. * output : the computed value of the steplength. * alphas ( dble ) * input : meaningless. * output : The maximum steplength such that the feasibility * on the superbasic variables remains feasible. * alphab ( dble ) * input : meaningless. * output : The maximum steplength such that the feasibility * on the basic variables remains feasible. * trial ( int ) * input : meaningless. * output : the number of linesearch iterations performed * before finding a new and better feasible point. * dir ( dble ) * input : the superbasic components of the direction. * output : the basic components of the projected direction * are computed. * ds ( dble ) * input : vector used as workspace. * output : meaningless. * x ( dble ) * input : the current iterate vector. * output : a new and better iterate is produced. * osol ( dble ) * input : array used as workspace. * output : it contains the previous iterate. * rgra ( dble ) * input : the reduced gradient vector. * output : unmodified. * rgnor ( dble ) * input : the euclidean norm of the reduced gradient vector. * output : unmodified. * dbd ( dble ) * input : array used as workspace. * output : vector containing for each superbasic variable * the absolute difference between its value and * its active bound value. * newp ( log ) * input : meaningless. * output : .true. iff a new and better feasible point * has been found. * .false. otherwise. * fcalc ( int ) * input : meaningless. * output : this vector contains the indices of the element * functions for which information is required. * lfc ( int ) * input : meaningless. * output : the number of element functions for which * information is required. * fuval ( dble ) * input : array used to store the function and derivative * values of the element functions. * output : the element functions whose indices are in FCALC(i) * for i=1 up to i=LFC are evaluated at new point X. * fx ( dble ) * input : the objective function value. * output : the objective function value evaluated at * new point X. * difval ( dble ) * input : meaningless. * output : the decrease in the objective function value. * oldfi ( dble ) * input : array used as workspace. * output : it contains the element function values * evaluated at the previous iterate. * nefval ( int ) * input : the number of element function evaluations. * output : this number is increased by one each time * an element function needs to be evaluated. * elvar ( int ) * input : array containing the indices of the varaibles * in the first element, followed by those in the * second element, etc. * output : unmodified. * elptr ( int ) * input : array whose kth value is the position of * the first variable of element k, in the * list ELVAR. * output : unmodified. * xst ( int ) * input : vector containing the status of the variables. * output : unmodified. * bound ( int ) * input : meaningless. * output : 0 iff no basic and superbasic variables hit a bound. * 1 iff at least one basic variable hits a bound. * 2 iff no basic variable hit a bound and a least one * superbasic variable hits a bound. * 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. * nqa ( int ) * input : the number of quasi-active 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. * inform ( int ) * input : must be set to zero the first time this routine * is called. * If it is not equal to zero, the routine is then * re-entered with some information that was needed * to proceed the calculation. * output : If it is not equal to zero, some element function * values are needed to proceed the calculation. * A return is then made to the main program. * Otherwise, no further information is needed. * ipdevc ( int ) * input : output device unit number for printing messages. * output : unmodified. * nflag ( log ) * input : meaningless. * output : .true. iff the difference between two successive * function values along the current linesearch * is insignificant compared to the noise on * these function values. * .false. otherwise. * iflag ( int ) * input : should be equal to zero when the routine is called. * output : = 18, if the directional derivative at the beginning * of the linesearch procedure is non-negative. * Otherwise, it remains unmodified. * Routines used : * --------------- * sqrt, abs, min, max. * dcopy, dsetcd, active, xpdbs, dcopcd, bding. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer ba(*), nba, su(*), nsu, nqa, ftstem(*), + bep(*), beptr(*), nefval, fcalc(*), lfc, + trial, elvar(*), elptr(*), oldfi(*), bound, + xst(*), kdist(*), ipdevc, iflag, inform logical newp, nflag double precision alpha, alphas, alphab, dir(*), ds(*), x(*), + osol(*), rgra(*), rgnor, dbd(*), fx, + difval, fuval(*) * Internal variables integer i, j, ik, k, kk, kd, kdis, inoise, ifc logical move, bfeas double precision sqrem, alpham, lambda, rd, rdini, ard, + bdact, f0, mindec, calfa, dsmax, + dsnor, slope * Save variables save i, j, ik, k, kk, kd, kdis, inoise, ifc save move, bfeas save sqrem, alpham, lambda, rd, rdini, ard, + bdact, f0, mindec, calfa, dsmax, + dsnor, slope double precision active external active * 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 * Parameter definition. integer ttrial parameter( ttrial = 10 ) * * Reverse communication test. * if( inform.ne.0 ) goto 100 * * Some initializations. * sqrem = sqrt(epsmch) alphas = huge calfa = half**ttrial rd = zero inoise = 0 bound = 0 newp = .false. move = .false. nflag = .false. * call dcopy( arcs, x, 1, osol, 1 ) call dsetcd( nba, ba, dir, zero ) * * We compute the maximum steplength ALPHAS, * such that the feasibility on the superbasic * variables remains feasible. * do 20 ik = 1 , nsu k = su(ik) * * The direction is set to minus the reduced * gradient value for the quasi active variables. * if( ik.gt.nsu-nqa ) dir(k) = -rgra(k) * * We test if superbasic component K of the * direction is too small to be considered. * if( abs(dir(k)).gt.sqrem ) then bdact = active(k,dir) dbd(k) = abs( x(k)-bdact ) alpha = dbd(k) / abs(dir(k)) * * We test if variable K is near its bound. * if( dbd(k).le.tol*(one+abs(bdact)) + .or. alpha.le.sqrem ) then dir(k) = zero else move = .true. rd = rd + rgra(k)*dir(k) rdini = rd alphas = min( alphas , alpha ) if( kd.eq.0 .or. abs(dir(k)).gt.abs(dir(kd)) ) kd = k endif * * We computes the basic components of * the direction. * kdis = kdist(ik) do 30 i = beptr(k) , beptr(k)+kdis-1 kk = bep(i) if( kk.lt.0 ) then dir(-kk) = dir(-kk) - dir(k) else dir(kk) = dir(kk) + dir(k) endif 30 continue endif 20 continue * * No move is possible. No * new point can be found. * if( .not.move ) return * * We compute the euclidean norm DSNOR * of the superbasic direction. * dsnor = zero dsmax = abs(dir(kd)) do 110 ik = 1 , nsu k = su(ik) if( k.ne.kd ) then dsnor = dsnor + (dir(k)/dsmax)**2 endif 110 continue dsnor = dsmax * sqrt(one+dsnor) * * We test if the superbasic direction is descent ? * slope = rd / (rgnor*dsnor) if( slope.ge.zero ) then * * The superbasic direction is not descent. * iflag = 18 return endif * * We compute the maximum steplength ALPHAB, * such that the feasibility on the basic * variables remains feasible. * alphab = huge do 40 ik = 1 , nba k = ba(ik) * * We test if basic component K of the * direction is too small to be considered. * if( abs(dir(k)).gt.sqrem ) then bdact = active(k,dir) alpha = abs( (x(k)-bdact)/dir(k) ) * * We test if variable K is near its bound. * if( abs(x(k)-bdact).le.tol*(one+abs(bdact)) + .or. alpha.le.epsmch ) then * * No move is possible. No * new point can be found. * alphab = zero return endif alphab = min( alpha , alphab ) else dir(k) = zero endif 40 continue * * Depending on the computed values * ALPHAS and ALPHAB, we define the * first value of ALPHA. * if( alphab.le.alphas ) then alpha = min( one , alphab ) bfeas = .true. else alpha = one if( alphas.ge.one ) then bfeas = .true. else bfeas = .false. endif endif alpham = min( alphab , alphas ) * * The linesearch iteration loop. * if ( ttrial .ge. 1 ) then trial = 1 50 continue * * We build the feasible basic-superbasic direction * corresponding to the value of steplength ALPHA. * call xpdbs( alpha, alpham, x, dir, ds, dbd, rgra, rd, + bfeas, ba, nba, su, nsu, ftstem, bep, beptr, + kdist ) * * We define the minimum function decrease required. * mindec = tenm2 * rd * alpha * * We select the element function that need * to be evaluated at new point X. * lfc = 0 do 60 i = 1 , elem do 70 j = elptr(i) , elptr(i+1)-1 k = elvar(j) if( x(k).ne.osol(k) ) then lfc = lfc + 1 fcalc(lfc) = i oldfi(i) = fuval(i) nefval = nefval + 1 go to 60 endif 70 continue 60 continue inform = 29 return * * We compute the new objective function value. * 100 continue inform = 0 f0 = fx fx = zero do 80 i = 1 , elem fx = fx + fuval(i) 80 continue * * We test if the decrease in the objective * function value is sufficient. * difval = fx - f0 if( difval.gt.mindec ) then * * No sufficient decrease. * fx = f0 do 90 ifc = 1 , lfc i = fcalc(ifc) fuval(i) = oldfi(i) 90 continue else * * A new feasible point is found. * newp = .true. endif * * We test if the decrease in the objective function value * is insignificant compared to the noise on the function * value. * if( abs(difval).le. elem*epsmch*(one+abs(f0)) ) then inoise = inoise + 1 if( inoise.ge.2 ) then write(ipdevc,1001) 1001 format(' Function value in the noise !!! ') nflag = .true. fx = f0 call dcopcd( nba, ba, osol, x ) call dcopcd( nsu, su, osol, x ) return endif endif * * We check if a variable hits a bound. * if( newp .or. ( inoise.eq.1 .and. + alpha.eq.alphab .and. alpha.le.calfa ) ) then call bding( bound, ba, nba, su, nsu, x, xst ) return endif * * No new point has been found, we reset * point X to its original values. * call dcopcd( nba, ba, osol, x ) call dcopcd( nsu, su, osol, x ) * * We update the steplength ALPHA. * if( trial.eq.3 .and. alpha.gt.alphas ) then * * We avoid a large number of halving steps. * alpha = alphas else if( alpha.gt.alphas ) then * * Bissection. * alpha = half * alpha else * * Safeguarded quadratic interpolation. * ard = alpha * rdini lambda = - half * ard / ( difval - ard ) alpha = min( max( tenm2 , lambda) , half ) * alpha endif endif if( alpha.le.alpham ) bfeas = .true. * trial = trial + 1 if ( trial .le. ttrial ) go to 50 endif * * The maximum number of linesearch iterations has been * performed without finding a better feasible point. * return end