* ************************************************************************* subroutine optim( elvar, elptr, elst, dbstr, xst, fr, to, maxit, + itmaj, itmin, minit, itcg, itcgmx, gptr,lg,hptr, + pr, thr, rthr, depth, naibs, lib, nb, iss, ip, + ftstem, nro, kandis, npiv, nacti, tdeac, ba, su, + kdist, beptr, bep, optmal, inform, nefval, + negval, nehval, fcalc, lfc, prec, nipst, x, fx, + fuval, epsf, dir, olgra, rgra, rgmax, rgnor, + wpre, dw1, dw2, dw3, dw4, iw3, ipdevc, what, + freq, info, bound, nflag, iflag ) * ************************************************************************* * Purpose : * --------- * This routine minimizes a partially separable objective function * in the independent set IP. The method used is iterative. At each * iteration, a quadratic model of objective function is constructed. * The step is determined by using a truncated conjugate gradient * algorithm, followed by a linesearch. * The strategy for treating bound constraints is based on the usual * projection device. This technique is similar to that used by * Bertsekas. * The method combines specialized data structure, MINOS-like * partitionning of the variables and the concept of independent * superbasic sets for handling the network constraints. * Parameters : * ------------ * elvar ( int ) * input : array containing the indices of the variables 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 function k, in the list ELVAR. * output : unmodified. * elst ( int ) * input : the status of the element functions. * output : this status may be modified in routine UPDATE. * dbstr ( log ) * input : .true. if at least one element function is * treated with Dembo's strategy * .false. otherwise. * output : unmodified. * xst ( int ) * input : the status of the variables. * output : these status may be modified in routine STATUS. * 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. * maxit ( int ) * input : the maximum number of minor iterations. * output : unmodified. * itmaj ( int ) * input : the major iteration number. * output : unmodified. * itmin ( int ) * input : the minor iteration number. * output : this number is increased by one when * a new minor iteration is started. * minit ( int ) * input : the minor iteration number in the * current major iteration. * output : this number is increased by one when * a new minor iteration is started. * itcg ( int ) * input : the total number of conjugate gradient iterations * executed before routine OPTIM is called. * output : this number is updated by the truncated conjugate * gradient procedure PCG or CG. * itcgmx ( int ) * input : the maximal number of conjugate gradient * iterations needed if all the linear system * were solved exactly before routine OPTIM * is called. * We consider here that a NxN linear system * requires N conjugate gradient iterations to * be solved exactly. * This quantity is used to estimate the average * size of the linear systems solved during the * minimization. * output : this number is updated by the truncated conjugate * gradient procedure PCG or CG. * gptr ( int ) * input : array whose ith value is the position of the the first * component of the ith element gradient in FUVAL, with * respect to its internal variables. * output : unmodified. * lg ( int ) * input : the total length of the element gradient vectors. * output : unmodified. * hptr ( int ) * input : array whose ith value is the position of the first * component of the ith element Hessian in FUVAL, with * respect to its internal variables. * output : unmodified. * pr ( int ) * input : the predecessor data structure vector. * output : this vector is updated after each pivot. * thr ( int ) * input : the traversal data structure vector. * output : this vector is updated after each pivot. * rthr ( int ) * input : the reverse data structure vector. * output : this vector is updated after each pivot. * depth ( int ) * input : the depth data structure vector. * output : this vector is updated after each pivot. * naibs ( int ) * input : vector whose IPth component contains the number * of superbasic variables in independent set IP. * output : the component decreases by one each time a * superbasic variable is detected as bounded. * lib ( int ) * input : this vector contains the indices of the * basic variables composing the maximal * basis spanning tree. The ith component * corresponds to the indice of the basic * arc whose endnodes are i and PR(i). * output : this vector is updated each time a * pivoting is performed. * nb ( int ) * input : the number of nonbasic variables. * output : this number increases by one each time a * superbasic variable is detected as bounded. * iss ( int ) * input : vector whose kth component contains the indice of * independent superbasic set in which the kth variable * is included. If the kth component is zero, then * variable k is not included in one independent set. * output : the components corresponding to the indices of the * superbasic variables that are bounded are set to zero. * If the superbasic set becomes empty, the components * corresponding to the indices of the basic variables are * also set to zero. * If a basic variable is bounded and no non bounded * superbasic variable is found to pivot with it, then * the component corresponding to the indice of the basic * variable is set to zero. * ip ( int ) * input : the indice of the current independent set * that is optimized. * output : unmodified. * ftstem ( int ) * input : vector containing for each superbasic variable * the length of its origine node to the stem node * and the length of its flow augmenting path. * output : this vector is updated after each pivoting. * nro ( int ) * input : the sum of the lengths of the flow augmenting * paths of all the superbasic variables. * output : this value is updated after each pivoting. * kandis ( int ) * input : the length of the longest flow augmenting path. * output : this length is updated after each pivoting. * npiv ( int ) * input : the number of minor iterations involving a pivoting. * output : this number is increased by one if at least one pivoting * is performed during one minor iteration. * nacti ( int ) * input : the number of variables that have been activated. * output : this number increases by one each time a * superbasic variable is detected as bounded. * tdeac ( log ) * input : .true. iff a deactivation has been done before * this routine is called. * .false. otherwise. * output : unmodified. * ba ( int ) * input : vector containing the indices of the basic variables. * output : this vector is updated when the status of the * variables are updated in routine STATUS. * su ( int ) * input : vector containing the indices of the superbasic * variables. * output : this vector is updated when the status of the * variables are updated in routine STATUS. * kdist ( int ) * input : vector containing for each superbasic variable * the length of its flow augmenting path. * output : only the components of the superbasic variables * that are not bounded are updated and maintained. * beptr ( int ) * output : array whose kth value is the position of the * first element of the the flow augmenting path * of the superbasic variable K, in array BEP. * output : when a pivoting has been performed, this array * is updated. * bep ( int ) * input : it contains the flow augmenting paths of all * the superbasic variables in the current subspace. * output : when a pivoting has been performed, this vector * is updated. * optmal ( int ) * input : it is always equal to 2. * output : = 1, iff a quasi optimal solution is found * for independent set IP. * = 2, iff an optimal solution is found for * independent set IP. * 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 or some element gradient values or some * element Hessian values are needed to proceed the * calculation. A return is then made to the main * program. * Otherwise, no further information is needed. * nefval ( int ) * input : the number of element function evaluations. * output : this number is increased by one each time an * element function is evaluated. * negval ( int ) * input : the number of element gradient evaluations. * output : this number is increased by one each time an * element gradient is evaluated. * nehval ( int ) * input : the number of element Hessian evaluations. * output : this number is increased by one each time an * element Hessian is evaluated. * fcalc ( int ) * input : meaningless. * output : when INFORM > 0, this vector contains the indices of * the element functions for which information is required. * lfc ( int ) * input : meaningless. * output : when INFORM > 0, the number of element functions for * which information is required. * prec ( log ) * input : .true. if the user wants to use a preconditioned * conjugate gradient scheme to solve the linear * systems of equations which arise at each * iteration of the minimization procedure. * .false. if the user does not want to use a * preconditioner. * output : unmodified. * nipst ( log ) * input : .true. if the user whishes to use the concept * of independent superbasic sets in order * to exploit the parrallelism in both the * constraints and the cost function structure. * .false. if no decomposition is achieved. * output : unmodified. * x ( dble ) * input : the current iterate vector. * output : this vector is updated during each minor iteration. * fx ( dble ) * input : the objective function value evaluated at iterate X. * output : this value is updated during the linesearch procedure. * fuval ( dble ) * input : array used to store the function and derivative * values of the element functions. * output : the element function values are updated in the * linesearch routine LNSRCH. The derivatives values * of the element functions are updated in routine * UPDATE. * epsf ( dble ) * input : measure of the accuracy required to stop the minimization * procedure. * output : unmodified. * dir ( dble ) * input : array used to store the search direction computed * in routine PCG or CG. * output : the search direction of the last minor iteration. * olgra ( dble ) * input : this array is used to contain the element gradient * vectors evaluated at the previous iterate. * output : it contains the differences between the * element gradient vectors evaluated at point * X and the element gradient vectors evaluated * at the previous iterate. * rgra ( dble ) * input : the reduced gradient vector. * output : this vector is updated after each minor iteration. * rgmax ( dble ) * input : the infinity norm of the reduced gradient vector. * output : this quantity is updated after each minor iteration. * rgnor ( dble ) * input : the euclidean norm of the reduced gradient vector. * output : this quantity is updated after each minor iteration. * wpre ( dble ) * input : array used to store the diagonal preconditioner. * output : the diagonal preconditioner of the last minor * iteration. * dw1 ( dble ) * input : array used as workspace. * output : meaningless. * dw2 ( dble ) * input : array used as workspace. * output : meaningless. * dw3 ( dble ) * input : array used as workspace. * output : meaningless. * dw4 ( dble ) * input : array used as workspace. * output : meaningless. * iw3 ( int ) * input : array used as workspace. * output : meaningless. * ipdevc ( int ) * input : output device unit number for printing messages. * output : unmodified. * what ( int ) * input : amount of output generated. See User's Guide * for the possible values and the corresponding * amount of output. * output : unmodified. * freq ( int ) * input : frequency of output. * <= 0 : no iteration output is generated, * except warning messages. * > 0 : output every FREQ iteration. * output : unmodified. * info ( int ) * input : meaningless. * output : on exit under un error condition, it contains * information about the error. * 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. * 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 : it must be equal to zero when this routine is called. * output : exit condition of the routine. * If this flag is positive, an error has been * detected by the routine. * Routines used : * --------------- * max, min, mod. * gxba, snxbdg, normrg, robs, qatst, * intprt, prcond, dcopy, pcg, cg, * second, lnsrch, update, status, * hrgra, rgrad. * Programming : * ------------- * D. Tuyttens * ======================================================================== * Routine parameters integer elvar(*), elptr(*), elst(*), xst(*), fr(*), + to(*), itmaj, itmin, minit, itcg, itcgmx, + gptr(*), lg, hptr(*), pr(*), thr(*), rthr(*), + depth(*), naibs(*), ba(*), nba, su(*), nsu, + nb, iss(*), ip, ftstem(*), nro, kdist(*), + lib(*), beptr(*), bep(*), optmal, inform, + nefval, negval, nehval, ipdevc, what, freq, + info, bound, iflag, kandis, npiv, nacti,maxit, + fcalc(*), iw3(*), lfc logical prec, nipst, nflag, tdeac, dbstr double precision x(*), fx, fuval(*), epsf, dir(*), olgra(*), + rgra(*), rgmax, rgnor, wpre(*), dw1(*),dw2(*), + dw3(*), dw4(*) * Internal variables integer k, nroot, nqa, trial, it, cgdim, nbfgs, nrk1, + noupd, ij logical tprt, pivot, newp, gxba double precision ergmax, alpha, alphas, alphab, difval real time1, time2 * Save variables save k, nroot, nqa, trial, it, cgdim, nbfgs, nrk1 save noupd, ij save tprt, pivot, newp save ergmax, alpha, alphas, alphab, difval save time1, time2 * Common specifications integer arcs, nodes, elem double precision epsmch, huge, tiny, tol real xocpu, totcpu, cgcpu, prdcpu, lnscpu, updcpu common / prbdim / arcs, nodes, elem common / prbmch / epsmch, huge, tiny, tol common / prbcpu / xocpu, totcpu, cgcpu, prdcpu, lnscpu, updcpu * * Reverse communication test. * if( inform.ne.0 ) then if( inform.ge.15 .and. inform.le.27 ) goto 200 if( inform.eq.28 ) goto 300 if( inform.eq.29 ) goto 100 endif * * We store in BA the indices of the basic variables * in independent set IP. We store in SU the indices * of the superbasic variables in independent set IP. * nro = 0 minit = 0 nba = 0 nsu = 0 * do 10 k = 1 , arcs if( iss(k).eq.ip ) then if( gxba(xst(k)) ) then nba = nba + 1 ba(nba) = k else nsu = nsu + 1 su(nsu) = k endif call snxbdg(xst(k)) endif 10 continue * * We compute the reduced gradient norm. * call normrg( rgmax, rgnor, rgra, su, nsu ) * * We define the subspace tolerance. * ergmax = max( epsf , min( 1.0d-1 , rgmax ) * rgmax ) * * The optimality test is checked. * if( rgmax.gt.epsf ) then * * The quasi-optimality test is checked. * if( rgmax.gt.ergmax ) then * * We store in array BEP the flow augmenting paths for * all the superbasic variables whose indices are in SU. * The root of the subtree defined by the basic variables * whose indices are in BA is also produced. * call robs( nroot, fr, to, lib, ba, nba, su, nsu, + bep, beptr, ftstem, pr, depth, xst, kdist ) * * The minor iterations loop. * if ( maxit .gt. itmin ) then it = 1 50 continue itmin = itmin + 1 minit = minit + 1 * * We determine the quasi-active variables. * call qatst( nqa, su, nsu, x, rgra, dir, kdist ) * * We print out the intermediate results. * call intprt( .true., what, ipdevc, freq, tprt, tdeac, + itmin, minit, itmaj, itcg, cgdim, ip, nba, + nsu, nefval, negval, nehval, fx, x, rgnor, + rgmax, ergmax, nqa, nbfgs, nrk1, noupd, + pivot, npiv, alpha, alphab, alphas, trial, + difval ) * * We compute the diagonal preconditioner. * call second(time1) if( prec ) then call prcond( fuval, lg, hptr, dw1, wpre, dw2, dw3, + elptr, elvar, elst, iss, ip, su, nsu, + nqa, ftstem, beptr, bep, kdist ) endif * * We compute a descent direction using * a truncated conjugate gradient method. * if( dbstr ) then * * A special treatment for Dembo's method. * call dcopy( arcs, x, 1, dw4(arcs+1), 1 ) call dcopy( lg, fuval(gptr(1)), 1, wpre(arcs+1), 1 ) endif * 300 continue if( prec ) then * * The preconditioned conjugate gradient method is called. * call pcg( itcg, itcgmx, wpre, rgra, dir, dw1, dw2, dw3, + dw4, olgra, x, fuval, gptr, lg, negval, hptr, + fcalc, lfc, nqa, nsu, su, nba, ba, elst, elvar, + elptr, bep, beptr, ftstem, kdist, inform, + ipdevc, what, tprt ) else * * The conjugate gradient method is called. * call cg( itcg, itcgmx, rgra, dir, dw1, dw2, dw3, dw4, + olgra, x, fuval, gptr, lg, negval, hptr, fcalc, + lfc, nqa, nsu, su, nba, ba, elst, elvar, elptr, + bep, beptr, wpre, ftstem, kdist, + inform, ipdevc, what, tprt ) endif if( inform.ne.0 ) return * if( dbstr ) then * * A special treatment for Dembo's method. * call dcopy( arcs, dw4(arcs+1), 1, x, 1 ) call dcopy( lg, wpre(arcs+1), 1, fuval(gptr(1)), 1 ) endif * call second(time2) cgcpu = cgcpu + (time2-time1) * * Once the descent direction DIR is computed, a feasible and * better iterate is produced using a linesearch procedure. * call second(time1) 100 continue call lnsrch( alpha, alphas, alphab, trial, dir, dw4, x, + dw1, rgra, rgnor, dw2, newp, iw3, lfc, + fuval, fx, difval, dw3, nefval, elvar, elptr, + xst, bound, ba, nba, su, nsu, nqa, ftstem, + bep, beptr, kdist, inform, ipdevc, nflag, + iflag ) if( iflag.ne.0 ) return * if( inform.ne.0 ) then do 20 ij = 1 , lfc fcalc(ij) = iw3(ij) 20 continue return endif call second(time2) lnscpu = lnscpu + (time2-time1) * * We test if a better iterate is produced. * if( newp ) then * * A better iterate is produced. * call second(time1) 200 continue * * We update the element gradient vectors and the * element Hessian matrices. * call update( inform, x, xst, elst, elptr, elvar, + fuval, fcalc, lfc, olgra, gptr, + lg, hptr, nefval, negval, nehval, + dw1, dw2, dw3, dw4, iw3, prec, + wpre, nbfgs, nrk1, noupd ) if( inform.ne.0 ) then inform = inform + 14 return endif call second(time2) updcpu = updcpu + (time2-time1) else * * No better iterate is produced. * if( nflag .and. what.gt.0 ) then write(ipdevc,1030) 1030 format(/,' Function values in the noise !!! ') endif optmal = 1 return * * The minor iteration is stopped. * endif * * We update the status of the variables. * call status( bound, pr, thr, rthr, depth, 0, npiv, pivot, + lib, ba, nba, su, nsu, x, xst, nroot, nro, + kandis, ftstem, fr, to, bep, beptr, kdist, + iss, naibs, ip, nb, nacti ) * * We update the reduced gradient vector. * if( 2*(nba+nsu) .le. nro ) then call hrgra( fuval, gptr, dw1, rgra, dw2, dw3, pr, thr, + depth, fr, to, elvar, elptr, lib, su, nsu, + elst, xst, nroot ) else call rgrad( rgra, dw1, dw2, fuval, gptr, elvar, elptr, + elst, su, nsu, bep, beptr, kdist ) endif * * We update the norm of the reduced gradient vector. * call normrg( rgmax, rgnor, rgra, su, nsu ) * * We print out the intermediate results. * call intprt( .false., what, ipdevc, freq, tprt, tdeac, + itmin, minit, itmaj, itcg, cgdim, ip, nba, + nsu, nefval, negval, nehval, fx, x, rgnor, + rgmax, ergmax, nqa, nbfgs, nrk1, noupd, + pivot, npiv, alpha, alphab, alphas, trial, + difval ) * * The optimality test is checked. * if( rgmax.le.epsf ) return * * The quasi-optimality test is checked. * if( rgmax.le.ergmax ) then optmal = 1 return endif * * A new minor iteration is started. * it = it + 1 if ( it .le. maxit - itmin ) go to 50 endif * * The maximum number of iterations is reached. * info = maxit iflag = 13 return endif * * The quasi-optimality test is verified. * optmal = 1 endif * return end