* ************************************************************************* subroutine lsnno( na, nn, nel, elvar, elptr, elst, fr, to, + x, xst, fx, + fuval, lfuval, gptr, hptr, inform, fcalc, + lfc, ifc, it, maxit, epsf, prec, nipst, + iwk, liwk, wk, lwk, ipdevc, what, freq, + info, iflag ) * ************************************************************************* * * Purpose : * --------- * Framework of the Large Scale Nonlinear Network Optimization algorithm. * This subroutine solves problems like : * Minimize f(x) * subject to A x = b * l < x < u * where f(x) is a nonlinear partially separable function and A is a * node-arc incidence matrix associated with a given network. * Parameters : * ------------ * na ( int ) * input : number of arcs of the network. It is also the * number of variables of the considered problem. * output : unmodified. * nn ( int ) * input : number of nodes of the network. * output : unmodified. * nel ( int ) * input : number of element functions. * output : unmodified. * 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 : array of length NEL. * It contains the informations about the technique * to use for updating the derivatives of the * element functions. It is given by the user. * See the User's Guide for the possible values. * output : It contains the status of the element functions * at the current point. * fr ( int ) * input : array containing the origine nodes of the arcs. * output : unmodified. * to ( int ) * input : array containing the destination nodes of the arcs. * output : unmodified. * range ( ext ) * : supplied subroutine that performs various operations * related to the true range of the element function * second derivative matrices. * xlower ( ext ) * : supplied double precision function routine that returns * the lower bounds on the bounded variables. * xupper ( ext ) * : supplied double precision function routine that returns * the upper bounds on the bounded variables. * rhs ( ext ) * : supplied double precision function routine that returns * the suuply/demand on the nodes. * x ( dble ) * input : array containing the values of the variables at the * starting point. The feasibility of this point is * checked and a new feasible starting point is computed * if necessary. * output : array containing the values of the variables at the * best point found (usually the solution). * xst ( int ) * input : XST(k) = -1 iff the variable k is unconstrained. * XST(k) = 0 iff the variable k is fixed. * XST(k) = +1 iff the variable k is bounded. * output : it contains the status of the variables at the * current point. * fx ( dble ) * input : meaningless. * output : value of the objective function at the current point X. * fuval ( dble ) * input : array which is used to store function and derivative * information for the nonlinear element functions evaluated * at the current point X. * output : contains the same informations. * lfuval ( int ) * input : length of array FUVAL. * output : unmodified. * 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. * 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. * inform ( int ) * input : Initially, it must be set to zero. * output : it is set in LSNNO and passed back to the user * to prompt further action (reverse communication). * Possible s values of INFORM are detailled in the * User's Guide. * fcalc ( int ) * input : meaningless. * output : integer array of length NEL. * The LFC first components are set in LSNNO to the * indices of the element functions for which the * user must provide the value and/or derivatives at X. * lfc ( int ) * input : meaningless. * output : it is set by LSNNO to the number of element * functions for which information is needed at X. * ifc ( int ) * input : meaningless. * output : ifc(1) = number of element function evaluations, * ifc(2) = number of element gradient evaluations, * ifc(3) = number of element Hessian evaluations. * it ( int ) * input : meaningless. * output : it(1) = Major iterations needed to solve the problem * it(2) = Minor iterations needed to solve the problem * it(3) = Inner CG iterations needed to solve the problem * maxit ( int ) * input : maximum number of minor iterations that the algorithm * is allowed to perform. * output : unmodified. * epsf ( dble ) * input : measure of the accuracy required to stop the minimization * procedure. * output : unmodified. * 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. * iwk ( int ) * input : integer array used as workspace. * output : meaningless. * liwk ( int ) * input : length of integer array IWK. * output : unmodified. * wk ( dble ) * input : double precision array used as workspace. * output : meaningless. * lwk ( int ) * input : length of double precision array WK. * output : unmodified. * 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. * iflag ( int ) * input : meaningless. * output : exit condition of the routine. * If this flag is positive, an error has been * detected by the routine. * Routines used : * --------------- * mod. * wkorg, dcmach, inptst, second, simplx, * fgh0, mbst, dkstem, inset0, inschk, hrgra, * normrg, optim, deact, exflag, endprt. * Programming : * ------------- * D. Tuyttens * ========================================================================= * Routines parameters integer na, nn, nel, elvar(*), elptr(*), elst(*), + fr(*), to(*), xst(*), lfuval, gptr(*), + hptr(*), inform, lfc, fcalc(*), ifc(3), + it(3), maxit, liwk, iwk(*), lwk, + ipdevc, what, freq, info, iflag logical prec, nipst double precision x(*), fx, epsf, fuval(*), + wk(*) * Internal variables integer nefval, negval, nehval, optmal, bound, + nba, nsu, nb, ncan, opt, s, i, k, + nro, kandis, ndacti, ndacts, nip, ip, + itmaj, itmin, minit, itcg, itcgmx, + npiv, nacti, manimx, ok, opti, opnstp, + opdea, noisyf, lfuvmx, lg, info2, niptot logical fnoise, nflag, tdeac, dbstr real time1, time2 double precision rgmax, rgnor, finit * Pointers of the successive vectors stored in "IWK" and "WK". integer pr, thr, rthr, depth, naibs, ba, su, + iss, ftstem, iw1, iw2, beptr, bep, dir, + olgra, rgra, wpre, dw1, dw2, dw3, dw4, + iw3 * Save variables save nefval, negval, nehval, optmal, bound, + nba, nsu, nb, ncan, opt, s, i, k, + nro, kandis, ndacti, ndacts, nip, ip, + itmaj, itmin, minit, itcg, itcgmx, + npiv, nacti, manimx, ok, opti, opnstp, + opdea, noisyf, lfuvmx, lg, info2, niptot save fnoise, nflag, tdeac, dbstr save time1, time2 save rgmax, rgnor, finit save pr, thr, rthr, depth, naibs, ba, su, + iss, ftstem, iw1, iw2, beptr, bep, dir, + olgra, rgra, wpre, dw1, dw2, dw3, dw4, + iw3 * Internal variables brief description * nefval : (int) the number of element function evaluations. * negval : (int) the number of element gradient evaluations. * nehval : (int) the number of element Hessian evaluations. * manimx : (int) the maximum number of major iterations. * itmaj : (int) the number of major iterations. * itmin : (int) the number of minor iterations. * minit : (int) the current minor iteration number. * itcg : (int) the number of conjugate gradient iterations. * itcgmx : (int) the maximum number of conjugate iterations * needed if all the linear systems were solved * exactly. We suupose here that an NxN linear * system needs exactly N conjugate gradient * iterations to be solved exactly. * npiv : (int) the number of minor iteration involving a pivoting. * nba : (int) the number of basic variables. * nsu : (int) the number of superbasic variables. * nb : (int) the number of nonbasic variables. * ncan : (int) the number of nonbasic variables candidate to * be de-activated. * nip : (int) the number of independent sets created. * niptot : (int) the number of non empty independent sets. * ip : (int) indice of the current independent set. * b : (int) the number of basic variables in independent set ip. * s : (int) the number of superbasic variables in independent set ip. * optmal : (int) = 0, no optimization is already been done. * = 1, the quasi optimality test is verified. * = 2, the optimality test is verified. * bound : (int) = 0, no basic and superbasic variables hit a bound. * = 1, at least one basic variable hits a bound. * = 2, no basic variable hit a bound and a least one * superbasic variable hits a bound. * nro : ( int ) the sum of the lengths of the flow augmenting * paths of all the superbasic variables. * kandis : ( int ) the length of the longest flow augmenting path. * nacti : (int) the number of variables that have been activated. * ndacti : (int) the number of variables that have been deactivated * during the last deactivation. * ndacts : (int) the total number of variables that have been deactivated. * opt : (int) the state of the minimization procedure. * The following cases are possible. * ok : (int) every thing is OK. We continue the minimization. * opti : (int) the optimality test is verified. An optimal * solution is found. * opnstp : (int) a new iteration needs to be started. * opdea : (int) a deactivation is needed. * lfuvmx : (int) the minimum space required in array FUVAL. * lg : (int) the total length of the element gradient vectors. * info2 : (int) variable used to save the value of INFORM * before a return is made to the main program. * i : (int) a loop indice. * k : (int) a variable indice. * noisyf : (int) the number of consecutive major iterations * detecting that the objevctive function is noisy. * fnoise : (log) .true if the objective function is noisy in * all the independent sets. * nflag : (log) .true. if the objective function is noisy in * the current independent set. This value is * produced by routine OPTIM. * tdeac : (log) .true. means that some nonbasic variables have * been deactivated. * time1 : (real) the CPU time used from the begin of the execution. * time2 : (real) the CPU time used from the begin of the execution. * These two variables are used to determine the CPU * time used in the different parts of the method. * rgmax : (dble) the reduced gradient infinity norm. * rgnor : (dble) the reduced gradient euclidean norm. * finit : (dble) objective function value evaluated at * the starting point. * Common specifications. integer arcs, nodes, elem double precision epsmch, huge, tiny, tol double precision zero, one, two, three, half, tenm1, tenm2, tenm4 real xocpu, totcpu, cgcpu, prdcpu, lnscpu, updcpu integer mfc, imfc common / prbdim / arcs, nodes, elem common / prbmch / epsmch, huge, tiny, tol common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4 common / prbcpu / xocpu, totcpu, cgcpu, prdcpu, lnscpu, updcpu common / mfcom / mfc, imfc * * Reverse communication test. * if( inform.gt.0 ) then inform = info2 if( inform.le.14 ) then go to 100 else go to 500 endif endif * * Some initializations. * arcs = na nodes = nn elem = nel * ndacts = 0 optmal = 0 npiv = 0 nacti = 0 itmaj = 0 itmin = 0 minit = 0 itcg = 0 itcgmx = 0 nefval = 0 negval = 0 nehval = 0 noisyf = 0 cgcpu = 0. prdcpu = 0. lnscpu = 0. updcpu = 0. * * We organize the integer workspace IWK and the double * precision workspace DWK and we check if the reserved * length LIWK for the integer workspace and the reserved * length LWK for the double precision workspace are high * enough. * call wkorg( liwk, pr, thr, rthr, depth, naibs, ba, su, + iss, ftstem, iw1, iw2, iw3, beptr, bep, elptr, + lwk, dir, olgra, rgra, prec, wpre, dw1, + dw2, dw3, dw4, iflag, info ) if( iflag.ne.0 ) go to 999 if( mfc.eq.0 ) mfc = 1 * * We define the machine dependent constants. * call dcmach( epsmch, huge, tiny ) tol = epsmch**(two/three) * * We print out the problem specifications. * if( what.ge.1 ) then write(ipdevc,2000) (arcs-mfc+1)/mfc,nodes/mfc,elem, + bep+4*arcs-1,dw4+2*arcs-1,epsmch endif * * We build the gradient and the Hessian pointers lists. * We also initialize the status of the element functions * and the status of the variables. * call inptst( elptr, gptr, hptr, elst, dbstr, xst, wk, + info, iflag ) if( iflag.ne.0 ) go to 999 * * We determine the minimum space needed to store * the element function values, the element gradient * vectors and the element Hessian matrices. * lg = gptr(elem+1) - gptr(1) lfuvmx = hptr(elem+1) - 1 * * We test if the total available space * LFUVAL is large enough. * if( lfuvmx .gt. lfuval ) then info = hptr(elem+1)-1 iflag = 3 go to 999 endif * * We check the feasibility of the initial solution. * If it is not feasible, we compute a feasible one. * if( what.ge.2 ) write(ipdevc,2001) call second(time1) * call dcopy( arcs, x, 1, wk(dw1), 1 ) k = bep + arcs + arcs do 10 i = 1, arcs iwk(k+i-1) = fr(i) iwk(iss+i-1) = to(i) 10 continue call simplx( iwk(k), iwk(iss), wk(dw1), xst, wk(1), + iwk(pr), iwk(thr), iwk(rthr), iwk(depth), iwk(iw1), + iwk(iw2), iwk(bep), iflag, info, ipdevc, what ) if( iflag.ne.0 ) go to 999 call dcopy( arcs, wk(dw1), 1, x, 1 ) * call second(time2) xocpu = time2 - time1 * * Possible meanings for variable OPT. * ok = 1 opti = 2 opnstp = 3 opdea = 4 * * The major iteration loop. * if( what.ge.2 ) write(ipdevc,2002) call second(time1) * if( maxit .gt. 0 ) then manimx = 1 70 continue opt = ok * * First major iteration special treatment. * if( manimx.eq.1 ) then * * We evaluate the objective function, gradient, * and Hessian at the starting point X. * 100 continue call fgh0( inform, x, xst, elst, elptr, elvar, + fx, fuval, fcalc, lfc, gptr, lg, hptr, + nefval, negval, nehval, wk(dw1), wk(dw2), + wk(dw3), wk(dw4), prec, wk(wpre) ) * * We test if some informations about the element * function, gradient or Hessian are needed to * proceed the calculation. A return to the main * program is the made. * if( inform.ne.0 ) then info2 = inform if( info2.gt.5 ) then inform = 1 else if( info2.gt.2 ) then inform = 2 else inform = inform + 2 endif return endif finit = fx * * We obtain the initial maximal basis spanning tree * and the corresponding Basic-Superbasic-Nonbasic * decomposition. * call mbst( iwk(pr), iwk(thr), iwk(rthr), iwk(depth), + iwk(ba), nba, nsu, iwk(su), nb, optmal, + bound, fr, to, x, xst, what, ipdevc, info, + iflag ) if( iflag.ne.0 ) go to 999 * * We test if the superbasic set empty (OPTMAL=2). * if( nb.eq.0 .and. optmal.eq.2 ) then * * The superbasic set is empty and the * nonbasic set is empty. An optimal * solution is found. * opt = opti else if( optmal.eq.2 ) then * * The superbasic set is empty but the * nonbasic set is not empty. A deactivation * is needed. * opt = opdea else * * The superbasic set is not empty. * opt = ok endif endif * if( opt.eq.ok ) then * * For each superbasic arc K, we evaluate the * lengths from its endnodes to the stem node. * nro = 0 kandis = 0 do 20 i = su , su+nsu-1 k = iwk(i) call dkstem( k, fr, to, iwk(pr), iwk(depth), + iwk(ftstem), nro, kandis ) 20 continue * * We partition the variables into independent sets. * call inset0( nipst, nip, iwk(naibs), iwk(iss), optmal, + elvar, elptr, iwk(iw1), iwk(iw2), xst, + iwk(ftstem), iwk(pr), iwk(bep), fr, to, + iwk(ba), nba, info, iflag ) if( iflag.ne.0 ) go to 999 * * We check the consistency of the independent sets. * call inschk( nip, niptot, iwk(naibs), nsu, info, iflag ) if( iflag.ne.0 ) go to 999 * * We compute the reduced gradient vector. * call hrgra( fuval, gptr, wk(dw2), wk(rgra), wk(dw1), + wk(dw3), iwk(pr), iwk(thr), iwk(depth), + fr, to, elvar, elptr, iwk(ba), iwk(su), + nsu, elst, xst, 1 ) * * When the nonbasic set is not empty, * a decativation is done before starting * the first major iteration. * if( nb.ne.0 ) opt = opdea endif endif * * End of first major iteration special treatment. * if( opt.eq.ok ) then * * We compute the reduced gradient norm. * call normrg( rgmax, rgnor, wk(rgra), iwk(su), nsu ) * * The optimality test is checked. * if( rgmax.le.epsf ) then * * The optimality test is verified. * if( nb.eq.0 ) then * * The nonbasic set is empty, the * current point is an optimal solution. * optmal = 2 opt = opti else * * The nonbasic set is not empty. * if( optmal.eq.2 ) then * * The optimality test was verified before * the last deactivation and is still verified * after the deactivation. The current point * is an optimal solution. * opt = opti else * * A deactivation is needed. * if( optmal.eq.0 ) then optmal = 2 opt = opdea else optmal = 2 ndacti = 0 rgmax = zero opt = opdea endif endif endif endif * * Optimization of the different independent sets. * if( opt.eq.ok ) then * * Some initializations before * starting the optimization. * optmal = 2 fnoise = .true. * * We test if a new major iteration is started. * if( ndacti.ne.0 .or. itmaj.eq.0 ) then tdeac = .true. itmaj = itmaj + 1 if( what.ge.2 ) write(ipdevc,2003) itmaj,ndacti,niptot endif * * We loop on the independent sets. * ip = 1 30 continue s = iwk(naibs+ip-1) if( s.ne.0 ) then 500 continue call optim(elvar, elptr, elst, dbstr, xst, fr, to,maxit, + itmaj, itmin, minit, itcg, itcgmx, gptr, + lg, hptr, iwk(pr), iwk(thr), iwk(rthr), + iwk(depth), iwk(naibs), iwk(ba), nb, + iwk(iss), ip, iwk(ftstem), nro, kandis, + npiv, nacti, tdeac, iwk(iw1), iwk(iw1+nodes), + iwk(iw2), iwk(beptr), iwk(bep), optmal, + inform, nefval, negval, nehval, fcalc, lfc, + prec, nipst, x, fx, fuval, epsf, wk(dir), + wk(olgra), wk(rgra), rgmax, rgnor, wk(wpre), + wk(dw1), wk(dw2), wk(dw3), wk(dw4), iwk(iw3), + ipdevc, what,freq, info, bound, nflag,iflag ) if( iflag.ne.0 ) go to 999 * * We test if some informations about the element * function, gradient or Hessian are needed to * proceed the calculation. A return to the main * program is the made. * if( inform.ne.0 ) then info2 = inform if( info2.gt.19 .and. info2.ne.28 ) then inform = 1 else if( info2.gt.16 .or. info2.eq.28 ) then inform = 2 else inform = inform - 12 endif return endif * * Whe test if the objective function is noisy. * fnoise = fnoise .and. nflag endif * * Next independent set. * ip = ip + 1 if ( ip .le. nip ) go to 30 * * The optimization in the current major iteration * is terminated. We test the noise of the objective * function and we determine if a deactiavtion is needed. * if( fnoise ) then * * The objective function is noisy. * noisyf = noisyf + 1 if( noisyf.ge.2 ) then * * Two consecutive major iterations detect * that the objective function is noisy. * The algorithm is stopped. * iflag = 12 go to 999 else * * A deactivation is needed. * opt = opdea endif else * * The objective function is not too noisy. * noisyf = 0 endif * * An optimality test is checked before * a deactivation is done. * if( opt.eq.ok ) then if( optmal.eq.2 .and. nb.eq.0 ) then opt = opti else opt = opdea endif endif endif endif * * We test if a deactivation is needed. * if( opt.eq.opdea ) then * * A deactivation procedure is started. * call deact( fr, to, iwk(pr), iwk(depth), iwk(ftstem), + nro, kandis, iwk(bep+arcs), elvar, elptr, + iwk(su), nsu, iwk(ba), nba, iwk(bep), ncan, + optmal, xst, elst, ndacti, ndacts, nb, opt, + fnoise, wk(rgra), gptr, wk(dw1), x, fuval, + wk(dw2), nipst, nip, iwk(naibs), iwk(iss), + itmaj, iwk(iw1), iwk(iw2), info, iflag ) if( iflag.ne.0 ) go to 999 * * When the function is noisy and no variable * can be deactivated, the algorithm is stopped. * if( fnoise .and. ndacti.eq.0 ) then iflag = 12 go to 999 endif * * We test the consistency of the independent sets. * call inschk( nip, niptot, iwk(naibs), nsu, info, iflag ) if( iflag.ne.0 ) go to 999 endif * * We test if an optimal solution is found. * if( opt.eq.opti ) go to 999 * * A new major iteration is started. * manimx = manimx + 1 if( manimx .le. maxit ) go to 70 endif * * The maximum number iterations allowed is reached. * info = maxit iflag = 13 * * We prints out the message corresponding to the error * code IFLAG. * 999 continue if( iflag.ne.0 ) call exflag( ipdevc, iflag, info ) * * We obtain the CPU time used by the optimization. * call second(time2) totcpu = time2-time1 if( what.ge.2 ) write(ipdevc,2004) * * We verify the feasibility of the end solution. * call dcopy( arcs, x, 1, wk(dw1), 1 ) k = bep + arcs + arcs do 11 i = 1, arcs iwk(k+i-1) = fr(i) iwk(iss+i-1) = to(i) 11 continue call simplx( iwk(k), iwk(iss), wk(dw1), xst, wk(1), + iwk(pr), iwk(thr), iwk(rthr), iwk(depth), iwk(iw1), + iwk(iw2), iwk(bep), iflag, info, ipdevc, what ) call dcopy( arcs, wk(dw1), 1, x, 1 ) * * We prints out the results of the execution. * it(1) = itmaj it(2) = itmin it(3) = itcg ifc(1) = nefval ifc(2) = negval ifc(3) = nehval * call endprt( ipdevc, what, x, xst, iwk(su), nsu, iwk(ba), nba, nb, + itmaj, itmin, itcg, itcgmx, nefval, negval, nehval, + npiv, ndacts, nacti, wk(rgra), epsf, finit, fx ) * * Output formats. * 2000 format(/,' ********** PROBLEM SPECIFICATIONS **********',//, + ' NUMBER OF ARCS ',i4,/, + ' NUMBER OF NODES ',i4,/, + ' NUMBER OF ELEMENTS ',i4,//, + ' INTEGER SPACE USED ',i6,/, + ' DOUBLE SPACE USED ',i6,//, + ' MACHINE PRECISION ',d12.5,/) 2001 format(/,' ********** FLOW FEASIBILITY CHECKING **********') 2002 format(//,1x,24('>'),2x,'LSNNO - OPTIMIZATION PHASE ',24('<'),/) 2003 format(//' >> MAJOR ITERATION ',i3,' - DEACTIVATION OF ',i3, + ' VARIABLE(S) - ',i3,' IN-SET(S) <<') 2004 format(/,1x,24('>'),2x,'LSNNO - END OF OPTIMIZATION ',24('<'),//, * ' * OPTIMAL SOLUTION FOUND *') * return end