C ================================================================= C File: alspgma.f C ================================================================= C Last update: May 2nd, 2005. C Users are encouraged to download periodically updated versions of C this code at the TANGO home page: C C www.ime.usp.br/~egbirgin/tango/ C ***************************************************************** C ***************************************************************** C ALSPG solves problems of the form C --------------------------------- C C min f(x) C C subject to C C c_j(x) = 0, j in E, C c_j(x) <= 0, j in I, C l <= x <= u, C C where E is the set of indices of the equality constraints, I is C the set of indices of the inequality constraints, and there are C n variables and m constraints. C C ALSPG is an Augmented Lagrangian method that uses SPG to solve C the bound-constrained problems. C C ALSPG is part of the TANGO Project. C C Visit the TANGO home page in the web: C C www.ime.usp.br/~egbirgin/tango/ C ***************************************************************** C TANGO LICENSE: C -------------- C C The codes of TANGO are free for non-commercial use. They may be C distributed within a "Research Group" but not distributed to C third parties without the authorization of E. G. Birgin or J. M. C Martínez. The code must not be included in any commercial package C without previous agreement. C C Every published work that uses ALSPG or SPG must cite: C C E. G. Birgin, J. M. Martínez and M. Raydan, "Nonmonotone spectral C projected gradient methods on convex sets", SIAM Journal on C Optimization 10, pp. 1196-1211, 2000, C C and/or C C E. G. Birgin, J. M. Martínez and M. Raydan, "Algorithm 813: SPG - C software for convex-constrained optimization", ACM Transactions C on Mathematical Software 27, pp. 340-349, 2001. C C (See other related works at the TANGO home page.) C ***************************************************************** C HOW TO USE ALSPG TO SOLVE YOUR OWN PROBLEM C ------------------------------------------ C C You will find below 7 subroutines and that YOU MUST MODIFY to C solve your own problem. The modifications that you must do are: C C 1) In the SUBROUTINE INIP you must set the number of variables C of your problem (n), the initial point (x), the lower and upper C bounds on the variables (l and u), the number of constraints (m) C (set m equal to zero if your problems has no constraints), a C vector (equatn) that indicates, for each constraint, whether it C is and equality constraint or an inequality (least-or-equal) C constraint, the initial estimation of the Lagrange multipliers C (you can set it equal to zero if you do not have a better C estimative) and the initial penalty parameters (rho). C C 2) SUBROUTINES EVALF, EVALC, EVALG AND EVALJAC, MUST also be C MODIFIED to evaluate the objective function, the constraints, C the gradient vector of the objective function and the gradients C of the constraints, respectively. You will find the way in C which these subroutines must be modified below. C C 3) SUBROUTINE PROJ must be modified to compute the projection of C an arbitrary point onto the lower-level convex constraints. If the C lower-level are box constraints then you do not need to modify C the current version of SUBROUTINE PROJ. C C 4) OPTIONALLY, in the SUBROUTINE PARAM you may modify tolerances C and maximum of iterations related to the stopping criteria. There C are also two arguments related to the output in the screen and C a parameter related to the derivatives calculation. If you C prefer, leave them with the suggested values. C C Besides the modifications that you must do to solve your own C problem using ALSPG, ALSPG has many other arguments to C control a large variety of issues of the algorithmic behaviour. C To see how to control these other arguments, see their detailed C description at the begining of subroutine alspg (for arguments C related to Augmented Lagrangians) and at the begining of C subroutine spg (for arguments related to the bound-constraints C internal solver). C ****************************************************************** C ****************************************************************** subroutine inip(n,x,l,u,m,lambda,rho,equatn) C This subroutine must set some problem data. For achieving this C objective YOU MUST MODIFY it according to your problem. See below C where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C This subroutine has no input parameters. C C On Return C C n integer, C number of variables, C C x double precision x(n), C initial point, C C l double precision l(n), C lower bounds on x, C C u double precision u(n), C upper bounds on x, C C m integer, C number of constraints (excluding the bounds), C C lambda double precision lambda(m), C initial estimation of the Lagrange multipliers, C C rho double precision rho(m), C initial penalty parameters, C C equatn logical equatn(m), C for each constraint i, set equatn(i) = .true. if it is an C equality constraint of the form c_i(x) = 0, and set C equatn(i) = .false. if it is an inequality constraint of C the form c_i(x) <= 0. C SCALAR ARGUMENTS integer m,n C ARRAY ARGUMENTS logical equatn(m) double precision l(n),lambda(m),rho(m),u(n),x(n) C LOCAL SCALARS integer i C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR PROBLEM C DATA: C ****************************************************************** C Number of variables n = 3 C Initial point x(1) = 1.0d0 x(2) = 2.0d0 x(3) = 3.0d0 C Lower and upper bounds l(1) = -10.0d0 l(2) = -20.0d0 l(3) = -30.0d0 u(1) = 10.0d0 u(2) = 20.0d0 u(3) = 30.0d0 C Number of constraints (equalities plus inequalities) m = 2 C Lagrange multipliers approximation. Most users prefer to use the C null initial Lagrange multipliers estimates. However, if the C problem that you are solving is "slightly different" from a C previously solved problem of which you know the correct Lagrange C multipliers, we encourage you to set these multipliers as initial C estimates. Of course, in this case you are also encouraged to use C the solution of the previous problem as initial estimate of the C solution. Similarly, most users prefer to use rho = 10 as initial C penalty parameters. But in the case mentioned above (good C estimates of solution and Lagrange multipliers) larger values of C the penalty parameters (say, rho = 1000) may be more useful. More C warm-start procedures are being elaborated. do i = 1,m lambda(i) = 0.0d0 end do C Initial penalty parameters do i = 1,m rho(i) = 10.0d0 end do C For each constraint i, set equatn(i) = .true. if it is an equality C constraint of the form c_i(x) = 0, and set equatn(i) = .false. if C it is an inequality constraint of the form c_i(x) <= 0. equatn(1) = .true. equatn(2) = .false. C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE INIP. C ****************************************************************** end C ****************************************************************** C ****************************************************************** subroutine evalf(n,x,f,flag) C This subroutine must compute the objective function. For achieving C this objective YOU MUST MODIFY it according to your problem. See C below where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C f double precision, C objective function value at x, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the objective C function. (For example, trying to compute the square root C of a negative number, dividing by zero or a very small C number, etc.) If everything was o.k. you must set it C equal to zero. C SCALAR ARGUMENTS integer flag,n double precision f C ARRAY ARGUMENTS double precision x(n) C Objective function C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR OBJECTIVE C FUNCTION: C ****************************************************************** flag = 0 f = x(1) + x(2) + x(3) C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALF. C ****************************************************************** end C ****************************************************************** C ****************************************************************** subroutine evalg(n,x,g,flag) C This subroutine must compute the gradient vector of the objective C function. For achieving these objective YOU MUST MODIFY it in the C way specified below. However, if you decide to use numerical C derivatives (we dont encourage this option at all!) you dont need C to modify evalg. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C g double precision g(n), C gradient vector of the objective function evaluated at x, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of any component C of the gradient vector. (For example, trying to compute C the square root of a negative number, dividing by zero or C a very small number, etc.) If everything was o.k. you C must set it equal to zero. C SCALAR ARGUMENTS integer flag,n C ARRAY ARGUMENTS double precision g(n),x(n) C Gradient vector C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE GRADIENT C VECTOR OF YOUR OBJECTIVE FUNCTION: C ****************************************************************** flag = 0 g(1) = 1.0d0 g(2) = 1.0d0 g(3) = 1.0d0 C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALG. C ****************************************************************** end C ****************************************************************** C ****************************************************************** subroutine evalc(n,x,i,c,flag) C This subroutine must compute the i-th constraint of your problem. C For achieving this objective YOU MUST MOFIFY it according to your C problem. See below the places where your modifications must be C inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C i integer, C index of the constraint to be computed, C C On Return C C c double precision, C i-th constraint evaluated at x, C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C SCALAR ARGUMENTS integer i,flag,n double precision c C ARRAY ARGUMENTS double precision x(n) C Constraints C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR C CONSTRAINTS: C ****************************************************************** flag = 0 if ( i .eq. 1 ) then c = x(1) - x(2) return end if C If i = 2 ... c = - x(2) - x(3) + 3.0d0 C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALC. C ****************************************************************** end C ****************************************************************** C ****************************************************************** subroutine evaljac(n,x,i,indjac,valjac,nnzjac,flag) C This subroutine must compute the gradient of the constraint i. For C achieving these objective YOU MUST MODIFY it in the way specified C below. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C i integer, C index of the constraint whose gradient will be computed, C C On Return C C nnzjac integer, C number of perhaps-non-null elements of the computed C gradient, C C indjac integer indjac(nnzjac), C see below, C C valjac double precision valjac(nnzjac), C the non-null value of the partial derivative of the i-th C constraint with respect to the indjac(k)-th variable must C be saved at valjac(k). C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C SCALAR ARGUMENTS integer flag,i,m,n,nnzjac C ARRAY ARGUMENTS integer indjac(n) double precision x(n),valjac(n) C Sparse gradient vector of the i-th constraint C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE GRADIENTS C OF YOUR CONSTRAINTS: C ****************************************************************** flag = 0 if ( i .eq. 1 ) then nnzjac = 2 indjac(1) = 1 valjac(1) = 1.0d0 indjac(2) = 2 valjac(2) = - 1.0d0 return end if C If i = 2 ... nnzjac = 2 indjac(1) = 2 valjac(1) = - 1.0d0 indjac(2) = 3 valjac(2) = - 1.0d0 C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALJAC. C ****************************************************************** end C ****************************************************************** C ****************************************************************** subroutine param(gtype,epsfeas,epsopt,maxoutit,maxtotit,maxtotfc, +iprint,ncomp) C This subroutine must set some alspg arguments related to C stopping criteria and output. OPTIONALLY, you may modify the C suggested default values. C C Parameters of the subroutine: C C On Entry: C C This subroutine has no input parameters. C C On Return C C gtype integer, C type of derivatives calculation according to the C following convention: C C = 0, true derivatives. In this case, subrotuines evalg C and evaljac must be modified by the user to compute C the derivatives of the objective function and the C constraints. C = 1, finite differences approximation provided by C ALSPG. In this case, subroutines evalg and C evaljac may have an empty body but must be present. C It is also recommended that those empty-body C subroutines set flag = - 1. Last but not least, C the option gtype = 1 is not cheap neither safe. C C epsfeas double precision, C feasibility tolerance (for the sup-norm of the C constraints), C C epsopt double precision, C optimality tolerance (for the sup-norm of the projected C gradient of the Augmented Lagrangian), C C maxoutit integer, C maximum number of outer iterations, C C maxtotit integer, C maximum total number of inner iterations (total means C summing up the inner iterations of each outer iteration), C C maxtotfc integer, C maximum total number of Augmented Lagrangian functional C evaluations, C C iprint integer, C controls the ammount of information of the output C according to the following convention: C C = 0, no output at all, C = 1, information at each outer iteration but without any C information of how the subproblems are being solved, C = 2, the same as 1 plus information of each inner C iteration, C = 3, the same as 2 plus information of the line searches C and the calculation of truncated Newton direction C (using CG) of each inner iteration, C C (In all cases, an output file named solution.txt with C the final point, Lagrange mutipliers and penalty C parameters will be generated. Moreover, the same output C of the screen will be saved in a file named C alspg.out.) C C ncomp integer, C every time a vector is printed, just its first ncomp C component will be printed. C SCALAR ARGUMENTS integer gtype,iprint,maxoutit,maxtotfc,maxtotit,ncomp double precision epsfeas,epsopt C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET SOME ALSPG C ARGUMENTS RELATED TO DERIVATIVES, STOPPING CRITERIA AND OUTPUT: C ****************************************************************** gtype = 0 epsfeas = 1.0d-04 epsopt = 1.0d-04 maxoutit = 50 maxtotit = 1000000 maxtotfc = 5 * maxtotit iprint = 3 ncomp = 5 C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE PARAM. C ****************************************************************** end C ***************************************************************** C ***************************************************************** subroutine proj(n,x,flag) C This subroutine computes the projection of an arbitrary point onto C the feasible set. Since the feasible set can be described in many C ways, its information is in a common block. C C Parameters: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C point to be projected. C C On Return C C x double precision x(n), C projected point, C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the projection. (For example, C trying to compute the square root of a negative number, C dividing by zero or a very small number, etc.) If C everything was o.k. you must set it equal to zero. C PARAMETERS integer nmax parameter ( nmax = 100000 ) C COMMON ARRAYS double precision l(nmax),u(nmax) C SCALAR ARGUMENTS integer n,flag C ARRAY ARGUMENTS double precision x(n) C COMMON BLOCKS common /bounds/ l,u C LOCAL SCALARS integer i flag = 0 C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE PROJ TO PROJECT AN C ARBITRAY POINT ONTO THE LOWER-LEVEL SET: C ****************************************************************** do i = 1,n x(i) = max( l(i), min( x(i), u(i) ) ) end do C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE PARAM. C ****************************************************************** end C ****************************************************************** C EXCEPT IF YOU ARE AN EXPERIENT USER, YOUR MODIFICATIONS FINISH C HERE !!!! C ****************************************************************** C ****************************************************************** C ALL THE INFORMATION THAT ALSPG PRINTS IN THE SCREEN WILL ALSO C BE SAVED IN TEXT FILE CALLED alspg.out. THE SOLUTION, AS WELL C AS THE FINAL LAGRANGE MULTIPLIERS APPROXIMATION AND THE FINAL C PENALTY PARAMETERS WILL BE SAVED IN A TEXT FILE CALLED C solution.txt. C ****************************************************************** C ***************************************************************** C ***************************************************************** program alspgma C ****************************************************************** C NOW, WE SET THE MAXIMUM NUMBER OF VARIABLES (NMAX) AND THE MAXIMUM C NUMBER OF CONSTRAINTS (MMAX) THAT THE PROBLEM IS ALLOWED TO HAVE. C IF YOUR PROBLEM HAS MORE THAN THE DECLARED VALUES FOR NMAX AND C MMAX THEN YOU MUST MODIFY THESE VALUES IN SUCH A WAY THAT THE C NUMBER OF VARIABLES (N) OF YOUR PROBLEM BE SMALLER THAN OR EQUAL C TO NMAX AND THE NUMBER OF CONSTRAINTS (M) OF YOUR PROBLEM BE C SMALLER THAN OR EQUAL TO MMAX. C C ATTENTION: PARAMETERS NMAX AND MMAX ARE ALSO DECLARED IN TWO MORE C PLACES, SUBROUTINES EVALNAL AND EVALHD. IF YOU CHANGE THEIR VALUES C HERE THEN YOU MUST ALSO DO THE SAME MODIFICATION IN THOSE OTHER C PLACES TOO. THE VALUES OF NMAX AND MMAX MUST BE THE SAME IN ALL C PLACES! C ****************************************************************** C PARAMETERS integer mmax,nmax parameter ( mmax = 100000 ) parameter ( nmax = 100000 ) C COMMON ARRAYS double precision l(nmax),u(nmax) C LOCAL SCALARS integer gtype,i,inform,m,n,maxoutit,maxtotit,maxtotfc,outiter, + totfcnt,totgcnt,totiter,iprint,ncomp double precision f,nalpsupn,epsfeas,epsopt,snorm real time,dtime C LOCAL ARRAYS logical equatn(mmax) double precision lambda(mmax),rho(mmax),wd(6*nmax+4*mmax),x(nmax) real dum(2) C COMMON BLOCKS common /bounds/ l,u C EXTERNAL SUBROUTINES external dtime,alspg C SET UP PROBLEM DATA call inip(n,x,l,u,m,lambda,rho,equatn) C SET SOME ALSPG ARGUMENTS call param(gtype,epsfeas,epsopt,maxoutit,maxtotit,maxtotfc,iprint, +ncomp) C CALL ALSPG open(10, file='alspg.out') time = dtime(dum) call easyalspg(n,x,m,lambda,rho,equatn,epsfeas,epsopt,maxoutit, +maxtotfc,maxtotit,iprint,ncomp,gtype,f,snorm,nalpsupn,outiter, +totiter,totfcnt,totgcnt,inform,wd) time = dtime(dum) time = dum(1) close(10) C SAVE THE SOLUTION open(10, file='solution.txt') C Solution point write(10,9000) do i = 1,n write(10,9010) i,x(i) end do C Lagrange multipliers and penalty parameters write(10,9020) do i = 1,m write(10,9030) i,lambda(i),rho(i) end do close(10) C WRITE STATISTICS open(10,file='alspgma-tabline.out') write(10,9040) time,inform,n,m,outiter,totiter,totfcnt,totgcnt, + f,snorm,nalpsupn close(10) stop C NON-EXECUTABLE STATEMENTS 9000 format(/,'FINAL POINT:',//,2X,'INDEX',16X,'X(INDEX)') 9010 format(I7,1P,D24.16) 9020 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS AND ', + 'PENALTY PARAMETERS: ',//,2X,'INDEX',11X,'LAMBDA(INDEX)', + 15X,'RHO(INDEX)') 9030 format(I7,1P,D24.16,1X,1P,D24.16) 9040 format(F8.2,1X,I1,1X,I5,1X,I5,1X,I3,1X,I7,1X,I7,1X,I7,1X, +1P,D24.16,1X,1P,D7.1,1X,1P,D7.1) end C ****************************************************************** C STOP HERE ALL YOUR MODIFICATIONS !!!!!! C ****************************************************************** C ================================================================= C File: alspg.f C ================================================================= C Last update of EASYALSPG: March 16th, 2005. subroutine easyalspg(n,x,m,lambda,rho,equatn,epsfc,epsop,maxoutit, +maxtotfc,maxtotit,iprint,ncomp,gtype,f,snorm,nalpsupn,outiter, +totiter,totfcnt,totgcnt,inform,wd) C SCALAR ARGUMENTS integer gtype,inform,iprint,m,maxoutit,maxtotfc,maxtotit,n,ncomp, + outiter,totfcnt,totiter,totgcnt double precision epsfc,epsop,f,nalpsupn,snorm C ARRAY ARGUMENTS logical equatn(m) double precision lambda(m),rho(m),wd(6*n+4*m),x(n) C This subroutine aims to simplify the usage of alspg. C CONSTANTS FOR GENERAL USES C Steps: h = max( steabs, sterel * abs( x ) ) should be a number C such that h is small ( relatively to x ) and x + h is different C from x. So, h is something that can be used a a step for a finite C differences approximation of a partial derivative relative to x. C Epsilons: something smaller than max( epsabs, epsrel * abs( x ) ) C should be considered as ``zero'' when compared with x. It is used, C for example, to detect that a step taken during a line search is C too small. C Infinitys: infrel is a big number that may appear in the C calculations. infabs is a number that should never be reached in C the calculations and is used the represent ``infinite''. Detailed C explanations of how are they used are rather cumbersome. double precision steabs,sterel,epsabs,epsrel,infabs,infrel parameter ( steabs = 1.0d-10 ) parameter ( sterel = 1.0d-07 ) parameter ( epsabs = 1.0d-20 ) parameter ( epsrel = 1.0d-10 ) parameter ( infabs = 1.0d+99 ) parameter ( infrel = 1.0d+20 ) C LOCAL SCALARS character * 3 pptype integer maxitncp double precision epsopfin,epsopini,epsopcon,rhomult,rhofrac, + lammin,lammax lammin = - 1.0d+20 lammax = 1.0d+20 rhomult = 10.0d0 rhofrac = 0.5d0 pptype = 'ONE' epsopini = epsop epsopfin = epsop epsopcon = 1.0d-01 maxitncp = 9 C Call alspg call alspg(n,x,m,lambda,rho,equatn,lammin,lammax,rhomult,rhofrac, +pptype,epsfc,epsopini,epsopfin,epsopcon,maxitncp,maxoutit, +maxtotfc,maxtotit,iprint,ncomp,gtype,steabs,sterel,epsabs,epsrel, +infabs,infrel,f,wd(1),wd(m+1),snorm,wd(2*m+1),nalpsupn,outiter, +totiter,totfcnt,totgcnt,inform,wd(2*m+n+1),wd(3*m+n+1), +wd(4*m+n+1)) return end C ****************************************************************** C ****************************************************************** C Last update of ALSPG or any of its dependencies: C C May 2nd, 2005. C C See report of modifications at the end of this file. subroutine alspg(n,x,m,lambda,rho,equatn,lammin,lammax,rhomult, +rhofrac,pptype,epsfc,epsopini,epsopfin,epsopcon,maxitncp,maxoutit, +maxtotfc,maxtotit,iprint,ncomp,gtype,steabs,sterel,epsabs,epsrel, +infabs,infrel,f,c,sigma,snorm,nal,nalpsupn,outiter,totiter, +totfcnt,totgcnt,inform,lambar,sigpre,wd) C SCALAR ARGUMENTS character * 3 pptype integer gtype,inform,iprint,m,maxitncp,maxoutit,maxtotfc,maxtotit, + n,ncomp,outiter,totfcnt,totgcnt,totiter double precision epsabs,epsfc,epsopcon,epsopfin,epsopini,epsrel,f, + infabs,infrel,lammax,lammin,nalpsupn,rhofrac,rhomult, + snorm,steabs,sterel C ARRAY ARGUMENTS logical equatn(m) double precision c(m),lambar(m),lambda(m),nal(n),rho(m),sigma(m), + sigpre(m),wd(5*n),x(n) C Solves the general-constrained minimization problem C C Minimize f(x) C C subject to C C c_i(x) = 0, i \in E, C C c_i(x) <= 0, i \in I, C C l <= x <= u C C using the method of multipliers described in C C R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, C "On Augmented Lagrangian methods with general lower-level C constraints", Technical Report MCDO-050304, Department of Applied C Mathematics, UNICAMP, Brazil, 2005. C C Parameters of the subroutine: C C On entry: C C n integer C number of variables C C m integer C number of constraints (equations plus inequations C C x double precision x(n) C initial estimation of the solution C C lambda double precision lambda(m) C initial Lagrange multipliers estimation C C rho double precision rho(m) C initial penalty parameters C C equatn logical equatn(m) C equatn is logical array that, for each constraint, C indicates whether the constraint is an equality constraint C (.true.) or an inequality constraint (.false.) C C lammin double precision C see below C C lammax double precision C the Lagrange multipliers approximations used for the C definition of the subproblems must be inside the box C [lammin,lammax]. C C rhomult double precision C when increased, the penalty parameters will be multiplied C by rhomult C C rhofrac double precision C the penalty parameters will be, basically, increased when C the infeasibility is larger than a fraction of the C previous infeasibility. This fraction is rhofrac. C C epsfc double precision C required precision for feasibility and complementarity C C epsopini double precision C required precision for the first subproblem optimality C condition (sup-norm) C C epsopfin double precision C required precision for optimality C C epsopcon double precision C the precision used for each subproblem is allways C computed as epsopk = max( epsopfin, epsopcon * C epsop_k-1 ) C C maxitncp integer C maximum allowed number of outer iterations without C progress in feasibility C C maxoutit integer C maximum allowed number of outer iterations C C maxtotit integer C maximum allowed number of inner iterations C C maxtotfc integer C maximum allowed number of objective function evaluations C C lambar double precision lambar(m) C sigpre double precision sigpre(m) C working vectors C C On return: C C x double precision x(n) C final solution estimation C C lambda double precision lambda(m) C final Lagrange multipliers estimation C C rho double precision rho(m) C final penalty parameters C C f double precision C objective functional value at the solution C C c double precision c(m) C constraints at the solution C C s double precision c(m) C constraints values at the solution C C snorm double precision C constraints violation is defined as abs( c(i) for equality C constraints and abs( max( c(i), - lambar(i) / rho(i) ) ) C for inequalities, where lambar(i) is the projection of the C Lagrange multiplier approximation lambda(i) onto the C interval [lammin,lammax]. snorm represents the maximum of C these violations. C C nal double precision nal(n) C gradient vector of the Augmented Lagrangian at the C solution C C nalpsupn double precision C sup-norm of the continuous projected gradient of the C Augmented Lagrangian at the solution C C outiter integer C number of outer iterations C C totiter integer C total number of inner iterations C C totfcnt integer C total number of objective functional evaluations C C totgcnt integer C total number of gradient evaluations C C inform integer C C This output parameter tells what happened in this C subroutine, according to the following conventions: C C 0 = convergence with feasibility, optimality and C complementarity; C C 1 = it was achieved the maximum allowed number of C outer iterations (maxoutit); C C 2 = it was achieved the maximum allowed total number of C inner iterations (maxtotit); C C 3 = it was achieved the maximum allowed total number of C functional evaluations (maxtotfc); C C 4 = the algorithm stopped by ``lack of feasibility C progress'', that means that the current point is C infeasible (the constraints violation is larger than C the tolerance epsfc) and the constraint violation has C not even simple decrease during maxitncp consecutive C iterations; in this case, the problem may be C infeasible; C C -90 = error in evalfc subroutine; C C -91 = error in evalgjc subroutine; C C -92 = error in evalhd subroutine. C LOCAL SCALARS integer fcnt,flag,gcnt,i,icrit,iter,mprint,nprint,maxit,maxfc, + rhoincr double precision al,snormprev,epsopk,nalpi C SET UP SOME CONSTANTS C just for printing nprint = min0(n,ncomp) mprint = min0(m,ncomp) C INTIALIZATION icrit = 0 outiter = 0 totiter = 0 totfcnt = 0 totgcnt = 0 if ( m .ne. 0 ) then epsopk = epsopini else epsopk = epsopfin end if C COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS call evalfc(n,x,f,m,c,inform) totfcnt = totfcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 1 ) then write(*, 2000) inform write(10,2000) inform end if go to 500 end if C COMPUTE LAGRANGE MULTIPLIERS do i = 1,m lambar(i) = max( lammin, min( lammax, lambda(i) ) ) end do C COMPUTE COMPLEMENTARITY AND FEASIBILITY VIOLATIONS snormprev = infabs snorm = 0.0d0 do i = 1,m if ( equatn(i) ) then sigma(i) = c(i) else sigma(i) = max( c(i), - lambar(i) / rho(i) ) end if snorm = max( snorm, abs( sigma(i) ) ) end do C COMPUTE CONTINUOUS PROJECTED GRADIENT NORM if ( gtype .eq. 0 ) then call evalnal(n,x,m,lambda,rho,equatn,nal,inform) else ! if ( gtype .eq. 1 ) then call evalnaldiff(n,x,m,lambda,rho,equatn,nal,sterel,steabs, + inform) end if totgcnt = totgcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 1 ) then write(*, 2000) inform write(10,2000) inform end if go to 500 end if do i = 1,n wd(i) = x(i) - nal(i) end do call proj(n,wd,flag) if ( flag .ne. 0 ) then inform = - 92 if ( iprint .ge. 1 ) then write(*, 2000) inform write(10,2000) inform end if go to 500 end if nalpsupn = 0.0d0 do i = 1,n nalpi = wd(i) - x(i) nalpsupn = max( nalpsupn, abs( nalpi ) ) end do C PRINT INITIAL INFORMATION if ( iprint .ge. 1 ) then write(*, 900) n,m write(10, 900) n,m end if C MAIN LOOP 100 continue C PRINT ITERATION INFORMATION if ( iprint .ge. 1 ) then write(*, 970) outiter write(*, 980) nprint,(x(i),i=1,nprint) write(*, 990) mprint,(lambda(i),i=1,mprint) write(*, 1000) mprint,(rho(i),i=1,mprint) write(*, 1020) f,snorm,nalpsupn,totiter,totfcnt,totgcnt write(10, 970) outiter write(10, 980) nprint,(x(i),i=1,nprint) write(10, 990) mprint,(lambda(i),i=1,mprint) write(10,1000) mprint,(rho(i),i=1,mprint) write(10,1020) f,snorm,nalpsupn,totiter,totfcnt,totgcnt end if C SAVING INTERMEDIATE DATA FOR CRASH REPORT open(20, file='alspg-tabline.out') write(20,3000) n,m,outiter,totiter,totfcnt,totgcnt,f,snorm, + nalpsupn close(20) C TEST STOPPING CRITERIA if ( snorm .le. epsfc .and. nalpsupn .le. epsopfin ) then ! THE POINT IS FEASIBLE AND OPTIMAL inform = 0 if ( iprint .ge. 1 ) then write(*, 1060) inform,epsfc,epsfc,epsopfin write(10,1060) inform,epsfc,epsfc,epsopfin end if go to 500 end if if ( outiter .ge. maxoutit ) then ! MAXIMUM NUMBER OF OUTER ITERATIONS EXCEEDED, STOP inform = 1 if ( iprint .ge. 1 ) then write(*, 1070) inform,maxoutit write(10,1070) inform,maxoutit end if go to 500 end if if ( totiter .gt. maxtotit ) then ! MAXIMUM NUMBER OF INNER ITERATIONS EXCEEDED, STOP inform = 2 if ( iprint .ge. 1 ) then write(*, 1080) inform,maxtotit write(10,1080) inform,maxtotit end if go to 500 end if if ( totfcnt .gt. maxtotfc ) then ! MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED, STOP inform = 3 if ( iprint .ge. 1 ) then write(*, 1090) inform,maxtotfc write(10,1090) inform,maxtotfc end if go to 500 end if C TEST IF THE INFEASIBILITY DECREASED SUFFICIENTLY if ( snorm .gt. epsfc .and. snorm .ge. snormprev ) then icrit = icrit + 1 if ( icrit .ge. maxitncp ) then ! LACK OF FEASIBILITY PROGRESS, STOP inform = 4 if ( iprint .ge. 1 ) then write(*, 1100) inform,maxitncp write(10,1100) inform,maxitncp end if go to 500 end if else icrit = 0 end if C DO AN OUTER ITERATION outiter = outiter + 1 C OPTIMALITY REQUERIMENT FOR THE SUBPROBLEM if ( m .ne. 0 ) then if ( outiter .eq. 1 ) then epsopk = epsopini else epsopk = max( epsopfin, epsopk * epsopcon ) end if else epsopk = epsopfin end if maxit = 5000 maxfc = 10 * maxit C CALL THE INNER-SOLVER call easyspg(n,x,m,lambar,rho,equatn,epsopk,maxit,maxfc,iprint, +ncomp,gtype,sterel,steabs,epsrel,epsabs,infrel,infabs,al,nal, +nalpsupn,iter,fcnt,gcnt,inform,wd) totiter = totiter + iter totfcnt = totfcnt + fcnt totgcnt = totgcnt + gcnt if ( inform .lt. 0 ) then if ( iprint .ge. 1 ) then write(*, 2000) inform write(10,2000) inform end if go to 500 end if C COMPUTE OBJECTIVE FUNCTION AND CONSTRAINTS call evalfc(n,x,f,m,c,inform) totfcnt = totfcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 1 ) then write(*, 2000) inform write(10,2000) inform end if go to 500 end if C UPDATE LAGRANGEAN MULTIPLIERS APPROXIMATION do i = 1,m call evaldpdy(c(i),rho(i),lambar(i),equatn(i),lambda(i)) lambar(i) = max( lammin, min( lammax, lambda(i) ) ) end do C COMPUTE COMPLEMENTARITY AND FEASIBILITY VIOLATIONS snormprev = snorm snorm = 0.0d0 do i = 1,m sigpre(i) = sigma(i) if ( equatn(i) ) then sigma(i) = c(i) else sigma(i) = max( c(i), - lambar(i) / rho(i) ) end if snorm = max( snorm, abs( sigma(i) ) ) end do C UPDATE PENALTY PARAMETERS if ( pptype .eq. 'ONE' ) then if ( snorm .gt. rhofrac * snormprev ) then do i = 1,m rho(i) = rhomult * rho(i) end do if ( iprint .ge. 1 ) then write(*, 1030) rhofrac * snormprev,rhomult write(10,1030) rhofrac * snormprev,rhomult end if else if ( iprint .ge. 1 ) then write(*, 1050) write(10,1050) end if end if else ! if ( pptype .eq. 'DIF' ) then rhoincr = 0 do i = 1,m if ( abs(sigma(i)) .gt. rhofrac * abs(sigpre(i)) ) then rho(i) = rhomult * rho(i) rhoincr = rhoincr + 1 end if end do if ( iprint .ge. 1 ) then if ( rhoincr .ne. 0 ) then write(*, 1040) rhoincr,m,rhomult write(10,1040) rhoincr,m,rhomult else write(*, 1050) write(10,1050) end if end if end if C ITERATE go to 100 C STOP 500 continue return C NON-EXECUTABLE STATEMENTS 900 format( /,' Entry to ALSPG.',/, + /,' Number of variables: ',I7, + ' Number of constraints: ',I7) 970 format( /,' ALSPG outer iteration: ',I7) 980 format( /,' Current point (first ',I7, ' components): ', + /,6(1X,1PD11.4)) 990 format( /,' Updated Lagrange multipliers (first ',I7, + ' components): ', + /,6(1X,1PD11.4)) 1000 format( /,' Updated penalty parameters (first ',I7, + ' components): ', + /,6(1X,1PD11.4)) 1020 format( /,' Objective function value',38X,' = ',1PD11.4, + /,' Sup-norm of the constraints',35X,' = ',1PD11.4, + /,' Sup-norm of the projected gradient of the Augmented', + ' Lagrangian = ',1PD11.4,/, + /,' Up-to-now number of inner iterations',26X,' = ',4X,I7, + /,' Up-to-now number of Augmented Lagrangian function', + ' evaluations',1X,' = ',4X,I7, + /,' Up-to-now number of Augmented Lagrangian gradient', + ' evaluations',1X,' = ',4X,I7) 1030 format( /,' The desired infeasibility improvement was not', + ' achieved (',1PD11.4,').', + /,' The penalty parameter will be increased multiplying', + ' by rhofrac = ',1PD11.4,'.') 1040 format( /,' The desired infeasibility improvement was not', + ' achieved in ',I7,' over ',I7, + /,' constraints.', + /,' Penalty parameters will be increased multiplying', + ' by rhofrac = ',1PD11.4,'.') 1050 format( /,' Desired feasibility improvement was achieved,', + /,' so, penalty parameter(s) will not be modified.') 1060 format( /,' Flag of ALSPG = ',I2, + ' (Convergence with Sup-norm of the Augmented', + ' Lagrangian', + /,' projected gradient smaller than ',1PD11.4,',', + ' constraints Sup-norm smaller than', + /,' ',1PD11.4,', and largest complementarity violation', + ' smaller than ',1PD11.4,')') 1070 format( /,' Flag of ALSPG = ',I2, + ' (It was exceeded the maximum allowed number of outer', + /,' iterations (maxoutit =',I7,').)') 1080 format( /,' Flag of ALSPG = ',I2, + ' (It was exceeded the maximum allowed number of inner', + /,' iterations (maxinnit =',I7,').)') 1090 format( /,' Flag of ALSPG = ',I2, + ' (It was exceeded the maximum allowed number of', + /,' functional evaluations (maxfc =',I7,').)') 1100 format( /,' Flag of ALSPG = ',I2, + ' (Constraint violations have not decreased', + ' substantially over',I2,' outer iterations. ', + /,' Problem possibly infeasible.)') 2000 format( /,' Flag of ALSPG = ',I3,' Fatal Error',/, + /,' The following codes means: ',/, + /,' -90 : error in evalf subroutine', + /,' -91 : error in evalc subroutine', + /,' -92 : error in evalg subroutine', + /,' -93 : error in evaljac subroutine',/) 3000 format(I5,1X,I5,1X,I3,1X,I7,1X,I7,1X,I7,1X,1P,D24.16,1X,1PD7.1,1X, + 1PD7.1) end C ***************************************************************** C ***************************************************************** subroutine evalfc(n,x,f,m,c,inform) C This subroutine computes the objective function and the C constraints. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C m integer, C number of constraints, C C On Return C C f double precision, C objective function value at x, C C c double precision c(m), C constraints at x, C C inform integer C = 0 means that the evaluation was successfuly done, C < 0 means that some error occured during the evaluation. C SCALAR ARGUMENTS double precision f integer inform,m,n C ARRAY ARGUMENTS double precision c(m),x(n) C LOCAL SCALARS integer flag,i inform = 0 C COMPUTE OBJECTIVE FUNTION call evalf(n,x,f,flag) if ( flag .ne. 0 ) then inform = - 90 return end if C COMPUTE CONSTRAINTS do i = 1,m C COMPUTE THE i-TH CONSTRAINT call evalc(n,x,i,c(i),flag) if ( flag .ne. 0 ) then inform = - 91 return end if end do end C ***************************************************************** C ***************************************************************** subroutine evalal(n,x,m,lambda,rho,equatn,al,inform) C This subroutine computes the Augmented Lagrangian function. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C point at which the Augmented Lagrangian function will be C evaluated. C C m integer, C number of constraints (equalities plus inequalities), C C lambda double precision lambdae(m), C current estimation of the Lagrange multipliers, C C rho double precision rho(m) C penalty parameters, C C equatn logical equatn(m) C equatn is logical array that, for each constraint, C indicates whether the constraint is an equality C constraint (.true.) or an inequality constraint C (.false.) C C On Return C C al double precision al, C value of the Augmented Lagrangian function, C C inform integer C = 0 means that the evaluation was successfuly done, C < 0 means that some error occured during the evaluation. C SCALAR ARGUMENTS double precision al integer inform,m,n C ARRAY ARGUMENTS logical equatn(m) double precision lambda(m),rho(m),x(n) C LOCAL SCALARS integer flag,i double precision c,f,p inform = 0 C COMPUTE OBJECTIVE FUNTION call evalf(n,x,f,flag) if ( flag .ne. 0 ) then inform = - 90 return end if C COMPUTES AL(x) = f(x) + \sum_j=1^m P(c(j),rho(j),lambda(j)) al = f do i = 1,m C COMPUTE i-TH CONSTRAINT call evalc(n,x,i,c,flag) if ( flag .ne. 0 ) then inform = - 91 return end if C ADD P(c(i),rho(i),lambda(i)) call evalp(c,rho(i),lambda(i),equatn(i),p) al = al + p end do end C ***************************************************************** C ***************************************************************** subroutine evalp(y,rho,lambda,equatn,p) C SCALAR ARGUMENTS logical equatn double precision lambda,p,rho,y if ( equatn ) then p = y * ( lambda + 0.5d0 * rho * y ) else if ( lambda + rho * y .ge. 0.0d0 ) then p = y * ( lambda + 0.5d0 * rho * y ) else p = - 0.5d0 * lambda ** 2 / rho end if end if return end C ***************************************************************** C ***************************************************************** subroutine evalnal(n,x,m,lambda,rho,equatn,nal,inform) C This subroutine computes the gradient of the Augmented Lagrangian C function. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C point at which the gradient will be evaluated. C C m integer, C number of constraints, C C lambda double precision lambda(m), C current estimation of the Lagrange multipliers, C C rho double precision rho(m) C penalty parameters, C C equatn logical equatn(m) C equatn is logical array that, for each constraint, C indicates whether the constraint is an equality C constraint (.true.) or an inequality constraint C (.false.) C C On Return C C nal double precision nal(n), C gradient of the Augmented Lagrangian function, C C inform integer C = 0 means that the evaluation was successfuly done, C < 0 means that some error occured during the evaluation. C ****************************************************************** C HERE WE SET THE MAXIMUM NUMBER OF VARIABLES (NMAX) AND THE MAXIMUM C NUMBER OF CONSTRAINTS (MMAX) THAT THE PROBLEM IS ALLOWED TO HAVE. C IF YOUR PROBLEM HAS MORE THAN THE DECLARED VALUES FOR NMAX AND C MMAX THEN YOU MUST MODIFY THESE VALUES IN SUCH A WAY THAT THE C NUMBER OF VARIABLES (N) OF YOUR PROBLEM BE SMALLER THAN OR EQUAL C TO NMAX AND THE NUMBER OF CONSTRAINTS (M) OF YOUR PROBLEM BE C SMALLER THAN OR EQUAL TO MMAX. C C ATTENTION: PARAMETERS NMAX AND MMAX ARE ALSO DECLARED IN TWO MORE C PLACES, THE MAIN PROGRAM AND SUBROUTINE EVALHD. IF YOU CHANGE C THEIR VALUES HERE THEN YOU MUST ALSO DO THE SAME MODIFICATION IN C THOSE OTHER PLACES TOO. THE VALUES OF NMAX AND MMAX MUST BE THE C SAME IN ALL PLACES! C ****************************************************************** C PARAMETERS integer mmax,nmax parameter ( mmax = 100000 ) parameter ( nmax = 100000 ) C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF PARAMETERS NMAX AND MMAX. C ****************************************************************** C COMMON ARRAYS double precision dpdc(mmax) C SCALAR ARGUMENTS integer inform,m,n C ARRAY ARGUMENTS logical equatn(m) double precision lambda(m),nal(n),rho(m),x(n) C LOCAL SCALARS integer flag,i,j,jcnnz double precision c C LOCAL ARRAYS integer jcvar(nmax) double precision jcval(nmax) C COMMON BLOCKS common /graddata/ dpdc inform = 0 C COMPUTE THE GRADIENT OF THE OBJECTIVE FUNCTION call evalg(n,x,nal,flag) if ( flag .ne. 0 ) then inform = - 92 return end if C COMPUTE \nabla L(x) = \nabla f(x) + \sum_{j=1}^m dPdcj * dcjdx do j = 1,m C COMPUTE THE j-TH CONSTRAINT call evalc(n,x,j,c,flag) if ( flag .ne. 0 ) then inform = - 91 return end if C COMPUTE dP/dc call evaldpdy(c,rho(j),lambda(j),equatn(j),dpdc(j)) if ( dpdc(j) .ne. 0.0d0 ) then C COMPUTE THE GRADIENT OF THE j-TH CONSTRAINT call evaljac(n,x,j,jcvar,jcval,jcnnz,flag) if ( flag .ne. 0 ) then inform = - 93 return end if C ADD dPdc * dcdx do i = 1,jcnnz nal(jcvar(i)) = nal(jcvar(i)) + dpdc(j) * jcval(i) end do end if end do end C ***************************************************************** C ***************************************************************** subroutine evaldpdy(y,rho,lambda,equatn,dpdy) C SCALAR ARGUMENTS logical equatn double precision y,rho,lambda,dpdy if ( equatn ) then dpdy = lambda + rho * y else dpdy = max( 0.0d0, lambda + rho * y ) end if return end C ================================================================= C File: spg.f C ================================================================= C Last update of EASYSPG: March 16th, 2005. subroutine easyspg(n,x,m,lambda,rho,equatn,epsgpsn,maxit,maxfc, +iprint,ncomp,gtype,sterel,steabs,epsrel,epsabs,infrel,infabs, +f,g,gpsupn,iter,fcnt,gcnt,inform,wd) C SCALAR ARGUMENTS integer fcnt,gcnt,gtype,m,maxfc,maxit,n,ncomp,inform,iprint,iter double precision epsabs,epsgpsn,epsrel,f,gpsupn,infabs,infrel, + steabs,sterel C ARRAY ARGUMENTS logical equatn(m) double precision g(n),lambda(m),rho(m),wd(5*n),x(n) double precision gamma,sigma1,sigma2 parameter ( gamma = 1.0d-04 ) parameter ( sigma1 = 0.1d0 ) parameter ( sigma2 = 0.9d0 ) integer mininterp parameter ( mininterp = 4 ) double precision nint,lspgma,lspgmi parameter ( nint = 2.0d0 ) parameter ( lspgma = 1.0d+10 ) parameter ( lspgmi = 1.0d-10 ) integer pmax parameter ( pmax = 100 ) integer tmax parameter ( tmax = 10000 ) integer maxitnfp,maxitngp,p double precision epsnfp,fmin,epsgpen,gpeucn2 double precision lastfv(pmax),lastgpns(tmax) maxitngp = tmax maxitnfp = maxit epsnfp = 0.0d0 fmin = - infrel p = 10 epsgpen = epsgpsn call spg(n,x,m,lambda,rho,equatn,p,epsgpen,epsgpsn,maxitnfp, +epsnfp,maxitngp,fmin,maxit,maxfc,nint,mininterp,gtype,iprint, +ncomp,f,g,gpeucn2,gpsupn,iter,fcnt,gcnt,inform,wd(1),wd(n+1), +wd(2*n+1),wd(3*n+1),wd(4*n+1),lastfv,lastgpns,lspgma,lspgmi, +gamma,sigma1,sigma2,sterel,steabs,epsrel,epsabs,infrel,infabs) end C ***************************************************************** C ***************************************************************** subroutine spg(n,x,m,lambda,rho,equatn,p,epsgpen,epsgpsn,maxitnfp, +epsnfp,maxitngp,fmin,maxit,maxfc,nint,mininterp,gtype,iprint, +ncomp,f,g,gpeucn2,gpsupn,iter,fcnt,gcnt,inform,s,y,xbest,d,wd, +lastfv,lastgpns,lspgma,lspgmi,gamma,sigma1,sigma2,sterel,steabs, +epsrel,epsabs,infrel,infabs) C SCALAR ARGUMENTS integer fcnt,gcnt,gtype,inform,iprint,iter,m,maxfc,maxit,maxitnfp, + maxitngp,mininterp,n,ncomp,p double precision epsabs,epsgpen,epsgpsn,epsnfp,epsrel,f,fmin, + gamma,gpeucn2,gpsupn,infabs,infrel,lspgma,lspgmi,nint, + sigma1,sigma2,steabs,sterel C ARRAY ARGUMENTS logical equatn(m) double precision d(n),g(n),lambda(m),lastfv(0:p-1), + lastgpns(0:maxitngp-1),rho(m),s(n),wd(n),x(n),xbest(n), + y(n) C Subroutine SPG implements the Spectral Projected Gradient C Method (Version 2: "feasible continuous projected path") to C find the local minimizers of a given function with convex C constraints, described in C C E. G. Birgin, J. M. Martinez and M. Raydan, "Nonmonotone spectral C projected gradient methods for convex sets", SIAM Journal on C Optimization 10, pp. 1196-1211, 2000. C C The user must supply the external subroutines evalf, evalg C and proj to evaluate the objective function and its gradient C and to project an arbitrary point onto the feasible region. C C This version 17 JAN 2000 by E.G.Birgin, J.M.Martinez and M.Raydan. C Reformatted 03 OCT 2000 by Tim Hopkins. C Final revision 04 MAY 2001 by E.G.Birgin, J.M.Martinez and M.Raydan. C C On Entry C C n integer C number of variables C C x double precision x(n) C initial estimation of the solution C C m integer C lambda double precision lambda(m) C rho double precision rho(m) C These three parameters are not used nor modified by C SPG and they are passed as arguments to the user- C defined subroutines evalal and evalnal to compute the C objective function and its gradient, respectively. C Clearly, in an Augmented Lagrangian context, if SPG is C being used to solve the convex-constrained subproblems, m C would be the number of constraints, lambda the Lagrange C multipliers approximation and rho the penalty parameters C C p integer C number of previous functional values to be considered C in the nonmonotone line search C C RECOMMENDED: p = 10 C C CONSTRAINTS: p >= 0 C C epsgpen double precision C epsgpen means EPSilon for the Projected Gradient Euclidian C Norm. It is a small positive number for declaring C convergence when the Euclidian norm of the continuous C projected gradient is less than or equal to epsgpen C C RECOMMENDED: epsgpen = 1.0d-05 C C CONSTRAINTS: epsgpen >= 0.0 C C epsgpsn double precision C epsgpsn means EPSilon for the Projected Gradient Sup Norm. C It is a small positive number for declaring convergence C when the sup norm of the continuous projected gradient is C less than or equal to epsgpsn C C RECOMMENDED: epsgpsn = 1.0d-05 C C CONSTRAINTS: epsgpsn >= 0.0 C C maxitnfp integer C maxitnfp means MAXimum of ITerations with No Function C Progress. See below for more details. C C epsnfp double precision C epsnfp means EPSilon for No Function Progress. It is a C small positive number for declaring ''lack of progress in C the objective function value'' if f(x_k) - f(x_{k+1}) <= C epsnfp * max{ f(x_j) - f(x_{j+1}, j < k } during maxitnfp C consecutive iterations. This stopping criterion may be C inhibited setting maxitnfp equal to maxit. C C RECOMMENDED: maxitnfp = 5 and epsnfp = 1.0d-02 C C CONSTRAINTS: maxitnfp >= 1 and epsnfp >= 0.0 C C maxitngp integer C maxitngp means MAXimum of ITerations with No Gradient C Progress. If the order of the Euclidian norm of the C continuous projected gradient did not change during C maxitngp consecutive iterations then the execution stops. C C RECOMMENDED: maxitngp = 10 C C CONSTRAINTS: maxitngp >= 1 C C fmin double precision C function value for the stopping criteria f <= fmin C C There is a stopping criterion that stops SPG if a C point with a functional value smaller than fmin is found. C The idea behind this stopping criterion is to stop the C method if the objective function is not bounded from C bellow. C C RECOMMENDED: fmin = - infabs C C CONSTRAINTS: there are no constraints for this argument C C maxit integer C maximum number of allowed iterations C C RECOMMENDED: maxit = 1000 C C CONSTRAINTS: maxit >= 0 C C maxfc integer C maximum allowed number of functional evaluations C C RECOMMENDED: maxfc = 5 * maxit C C CONSTRAINTS: maxfc >= 1 C C nint double precision C Constant for the interpolation. See the description of C sigma1 and sigma2 above. Sometimes, in a line search, we C take the new trial step as the previous one divided by C nint C C RECOMMENDED: nint = 2.0 C C CONSTRAINTS: nint > 1.0. C C mininterp integer C Constant for testing if, after having made at least C mininterp interpolations, the steplength is too small. In C that case, failure of the line search is declared (may be C the direction is not a descent direction due to an error C in the gradient calculations). Use mininterp greater C than or equal to maxfc for inhibit this stopping C criterion C C RECOMMENDED: mininterp = 4 C C CONSTRAINTS: mininterp >= 1 C C gtype integer C gtype indicates in which way the gradient of the C objective function will be computed. If the user have C been implemented the user-supplied evalnal subroutine to C compute the gradient of the objective function then C gtype argument must be set to 0 (ZERO) and the user- C supplied evalnal subroutine will be called by SPG any C time the gradient would be required. C C subroutine evalnal(n,x,m,lambda,rho,equatn,g,flag) C C C On Entry: C C C n integer, C C number of variables, C C C C x double precision x(n), C C current point, C C C C m integer, C C number of constraints (equalities plus C C inequalities), C C C C lambda double precision lambda(m), C C current estimation of the Lagrange C C multipliers, C C C C rho double precision rho(m) C C penalty parameters, C C C C NOTE: arguments m, lambda and rho are useful when C C SPG is being used for solving the box- C C constrained subproblems of an Augmented Lagrangian C C framework. When SPG is being used stand-alone C C for solving a bound-constrained problem, these C C arguments are dummy arguments. C C C C On Return C C C C g double precision g(n), C C gradient of the objective function at x, C C C C flag integer C C 0 means ''no errors'', C C 1 means ''some error occurs in the gradient C C evaluation''. C C C C SCALAR ARGUMENTS C integer flag,m,n C C C ARRAY ARGUMENTS C double precision g(n),lambda(m),rho(m),x(n) C C C ''Here it should be the body of evalnal subroutine C C that saves in g the gradient vector of the C C objective at x. Moreover, it sets flag equal to 0 C C if the calculation was successfully done and sets C C flag equal to any other value different from 0 if C C the gradient vector is not well defined at the C C current point x. If SPG gtype argument was C C setted to 1, i.e., the finite difference C C approximation provided by SPG will be used, then C C this subroutine must even be present for C C compilation purposes but it will never be called.'' C C end C C If, on the other hand, the user is not able to provide C evalnal subroutine, gtype argument must be set to 1 C (ONE). In this case, every time SPG needs to compute C the gradient of the objective function, an internal C subroutine that approximates it by finite-differences C will be used (be aware that it maybe very time C consuming). Moreover, note that the evalnal subroutine C must still be present (with an empty body). C C RECOMMENDED: gtype = 0 (provided you have the evalg C subroutine) C C CONSTRAINTS: allowed values are just 0 or 1. C C iprint integer C Commands printing. Nothing is printed if iprint is C smaller than 2. If iprint is greater than or equal to C 2, SPG iterations information is printed. If iprint C is greater than or equal to 3, line searches and C Conjugate Gradients information is printed. C C RECOMMENDED: iprint = 2 C C CONSTRAINTS: allowed values are just 2 or 3. C C ncomp integer C This constant is just for printing. In a detailed C printing option, ncomp component of some vectors will be C printed C C RECOMMENDED: ncomp = 5 C C CONSTRAINTS: ncomp >= 0 C C d double precision d(n) C lastfv double precision lastfv(p) C lastgpns double precision lastgpns(maxitngp) C s double precision s(n) C y double precision y(n) C xbest double precision xbest(n) C wd double precision wd(n) C working vectors C C lspgmi double precision C See below C C lspgma double precision C The spectral steplength, called lamspg, is projected onto C the box [lspgmi,lspgma] C C RECOMMENDED: lspgmi = 1.0d-10 and lspgma = 1.0d+10 C C CONSTRAINTS: lspgma >= lspgmi > 0.0 C C gamma double precision C Constant for the Armijo criterion C f(x + alpha d) <= f(x) + gamma * alpha * C C RECOMMENDED: gamma = 1.0d-04 C C CONSTRAINTS: 0.0 < gamma < 0.5. C C sigma1 double precision C See below C C sigma2 double precision C Constant for the safeguarded interpolation. If alpha_new C is not inside the interval [sigma1, sigma * alpha] then C we take alpha_new = alpha / nint C C RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9 C C CONSTRAINTS: 0 < sigma1 < sigma2 < 1. C C sterel double precision C See below C C steabs double precision C This constants mean a ''relative small number'' and ''an C absolute small number'' for the increments in finite C difference approximations of derivatives C C RECOMMENDED: epsrel = 1.0d-07 and epsabs = 1.0d-10 C C CONSTRAINTS: sterel >= steabs > 0 C C epsrel double precision C See below C C epsabs double precision C See below C C infrel double precision C See below C C infabs double precision C This four constants mean a ''relative small number'', C ''an absolute small number'', ''a relative large number'' C and ''an absolute large number''. Basically, a quantity A C is considered negligible with respect to another quantity C B if |A| < max ( epsrel * |B|, epsabs ) C C RECOMMENDED: epsrel = 1.0d-10, epsabs = 1.0d-20, C infrel = 1.0d+20, infabs = 1.0d+99 C C CONSTRAINTS: epsrel >= epsabs >= 0.0 C infabs >= infrel >= 0.0 C C On Return C C x double precision x(n) C Final estimation to the solution C C f double precision C Function value at the final estimation C C g double precision g(n) C Gradient at the final estimation C C gpeucn2 double precision C Squared Euclidian norm of the continuous projected C gradient at the final estimation C C gpsupn double precision C the same as before but with sup-norm C C iter integer C number of iterations C C fcnt integer C number of function evaluations C C gcnt integer C number of gradient evaluations C C inform integer C This output parameter tells what happened in this C subroutine, according to the following conventions: C C 0 = convergence with small Euclidian norm of the C continuous projected gradient (smaller than epsgpen); C C 1 = convergence with small sup-norm of the continuous C projected gradient (smaller than epsgpsn); C C 2 = the algorithm stopped by ''lack of progress'', that C means that f(xk) - f(x_{k+1}) <= epsnfp * C max{ f(x_j) - f(x_{j+1}, j < k } during maxitnfp C consecutive iterations. If desired, set maxitnfp C equal to maxit to inhibit this stopping criterion. C C 3 = the algorithm stopped because the order of the C Euclidian norm of the continuous projected gradient C did not change during maxitngp consecutive C iterations. Probably, we are asking for an C exaggerated small norm of continuous projected C gradient for declaring convergence. If desired, set C maxitngp equal to maxit to inhibit this stopping C criterion. C C 4 = the algorithm stopped because the functional value c is very small (smaller than fmin). If desired, set C fmin equal to minus infabs to inhibit this stopping C criterion. C C 6 = too small step in a line search. After having made at C least mininterp interpolations, the steplength C becames small. ''small steplength'' means that we are C at point x with direction d and step alpha, and C C alpha * ||d||_infty < max( epsabs, epsrel * C ||x||_infty ). C C In that case failure of the line search is declared C (may be the direction is not a descent direction due C to an error in the gradient calculations). If C desired, set mininterp equal to maxfc to inhibit this C stopping criterion. C C 7 = it was achieved the maximum allowed number of C iterations (maxit); C C 8 = it was achieved the maximum allowed number of C function evaluations (maxfc); C C < 0 = error in evalal, evalnal or proj subroutines. C LOCAL SCALARS integer i,itnfp,flag,nprint double precision bestprog,currprog,epsgpen2,fbest,fprev,gpi, + gpnmax,lamspg,sts,sty,xnorm C EXTERNAL SUBROUTINES external spgls,evalal,evalnal,proj C INTRINSIC FUNCTIONS intrinsic abs,max,min,mod C ================================================================== C Initialization C ================================================================== C Set some initial values: C counters, iter = 0 fcnt = 0 gcnt = 0 C just for printing, nprint = min0( n, ncomp ) C for testing convergence, epsgpen2 = epsgpen ** 2 C for the nonmonotone line-search, do i = 0,p - 1 lastfv(i) = - infabs end do C for testing progress in f, and fprev = infabs bestprog = 0.0d0 itnfp = 0 C for testing progress in the projected gradient norm. do i = 0,maxitngp - 1 lastgpns(i) = infabs end do C Print problem information if( iprint .ge. 2 ) then write(*, 977) n write(*, 980) nprint,(x(i),i=1,nprint) write(10,977) n write(10,980) nprint,(x(i),i=1,nprint) end if C Project initial guess. If the initial guess is infeasible, C projection puts it into the feasible convex set. call proj(n,x,flag) if ( flag .ne. 0 ) then inform = - 92 if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if C Compute x Euclidian norm xnorm = 0.0d0 do i = 1,n xnorm = xnorm + x(i) ** 2 end do xnorm = sqrt( xnorm ) C Initialize best solution do i = 1,n xbest(i) = x(i) end do C Compute objective function and gradient at the initial point call evalal(n,x,m,lambda,rho,equatn,f,inform) fcnt = fcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if if ( gtype .eq. 0 ) then call evalnal(n,x,m,lambda,rho,equatn,g,inform) else ! if ( gtype .eq. 1 ) then call evalnaldiff(n,x,m,lambda,rho,equatn,g,sterel,steabs, + inform) end if gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if C Store functional value for the nonmonotone line search lastfv(0) = f C Initialize best functional value fbest = f C Compute continuous-project-gradient Euclidian and Sup norms do i = 1,n wd(i) = x(i) - g(i) end do call proj(n,wd,flag) if ( flag .ne. 0 ) then inform = - 92 if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if gpsupn = 0.0d0 gpeucn2 = 0.0d0 do i = 1,n gpi = wd(i) - x(i) gpsupn = max( gpsupn, abs( gpi ) ) gpeucn2 = gpeucn2 + gpi ** 2 end do C Print initial information if( iprint .ge. 2 ) then write(*, 981) iter write(*, 985) nprint,(x(i),i=1,nprint) write(*, 986) nprint,(g(i),i=1,nprint) write(*, 1002) f,sqrt(gpeucn2),gpsupn,fcnt,gcnt write(10,981) iter write(10,985) nprint,(x(i),i=1,nprint) write(10,986) nprint,(g(i),i=1,nprint) write(10,1002) f,sqrt(gpeucn2),gpsupn,fcnt,gcnt end if C ================================================================== C Main loop C ================================================================== 100 continue C ================================================================== C Test stopping criteria C ================================================================== C Test whether the continuous-projected-gradient Euclidian norm C is small enough to declare convergence if ( gpeucn2 .le. epsgpen2 ) then inform = 0 if ( iprint .ge. 2 ) then write(*, 990) inform,epsgpen write(10,990) inform,epsgpen end if go to 500 end if C Test whether the continuous-projected-gradient Sup norm C is small enough to declare convergence if ( gpsupn .le. epsgpsn ) then inform = 1 if ( iprint .ge. 2 ) then write(*, 991) inform,epsgpsn write(10,991) inform,epsgpsn end if go to 500 end if C Test whether we performed many iterations without good progress C of the functional value currprog = fprev - f bestprog = max( currprog, bestprog ) if ( currprog .le. epsnfp * bestprog ) then itnfp = itnfp + 1 if ( itnfp .ge. maxitnfp ) then inform = 2 if ( iprint .ge. 2 ) then write(*, 992) inform,epsnfp,maxitnfp write(10,992) inform,epsnfp,maxitnfp end if go to 500 endif else itnfp = 0 endif C Test whether we have performed many iterations without good C reduction of the euclidian-norm of the projected gradient gpnmax = 0.0d0 do i = 0,maxitngp - 1 gpnmax = max( gpnmax, lastgpns(i) ) end do lastgpns(mod( iter, maxitngp )) = gpeucn2 if ( gpeucn2 .ge. gpnmax ) then inform = 3 if ( iprint .ge. 2 ) then write(*, 993) inform,maxitngp write(10,993) inform,maxitngp end if go to 500 endif C Test whether the functional value is very small if ( f .le. fmin ) then inform = 4 if ( iprint .ge. 2 ) then write(*, 994) inform,fmin write(10,994) inform,fmin end if go to 500 end if C Test whether the number of iterations is exhausted if ( iter .ge. maxit ) then inform = 7 if ( iprint .ge. 2 ) then write(*, 997) inform,maxit write(10,997) inform,maxit end if go to 500 end if C Test whether the number of functional evaluations is exhausted if ( fcnt .ge. maxfc ) then inform = 8 if ( iprint .ge. 2 ) then write(*, 998) inform,maxfc write(10,998) inform,maxfc end if go to 500 end if C ================================================================== C The stopping criteria were not satisfied, a new iteration will be C made C ================================================================== iter = iter + 1 C ================================================================== C Save current values, f, x and g C ================================================================== fprev = f do i = 1,n s(i) = x(i) y(i) = g(i) end do C ================================================================== C Compute new iterate C ================================================================== C Compute spectral steplength if ( iter .eq. 1 .or. sty .le. 0.0d0 ) then lamspg = max( 1.0d0, xnorm ) / sqrt( gpeucn2 ) else lamspg = sts / sty end if lamspg = min( lspgma, max( lspgmi, lamspg ) ) C Perform a line search with safeguarded quadratic interpolation C along the direction of the spectral continuous projected C gradient call spgls(n,x,m,lambda,rho,equatn,f,g,lamspg,p,lastfv,nint, +mininterp,fmin,maxfc,iprint,fcnt,inform,wd,d,gamma,sigma1,sigma2, +epsrel,epsabs) if ( inform .lt. 0 ) then if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if C Save the functional value for the nonmonotone line search lastfv(mod( iter, p )) = f C Compare the functional value against the best one and, if smaller C update the best functional value and the corresponding best point if ( f .lt. fbest ) then fbest = f do i = 1,n xbest(i) = x(i) end do end if C Compute the gradient at the new iterate if ( gtype .eq. 0 ) then call evalnal(n,x,m,lambda,rho,equatn,g,inform) else ! if ( gtype .eq. 1 ) then call evalnaldiff(n,x,m,lambda,rho,equatn,g,sterel,steabs, + inform) end if gcnt = gcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if C ================================================================== C Prepare for the next iteration C ================================================================== C Compute x Euclidian norm xnorm = 0.0d0 do i = 1,n xnorm = xnorm + x(i) ** 2 end do xnorm = sqrt( xnorm ) C Compute s = x_{k+1} - x_k, y = g_{k+1} - g_k, and sts = 0.0d0 sty = 0.0d0 do i = 1,n s(i) = x(i) - s(i) y(i) = g(i) - y(i) sts = sts + s(i) ** 2 sty = sty + s(i) * y(i) end do C Compute continuous-project-gradient Euclidian and Sup norms do i = 1,n wd(i) = x(i) - g(i) end do call proj(n,wd,flag) if ( flag .ne. 0 ) then inform = - 92 if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if gpsupn = 0.0d0 gpeucn2 = 0.0d0 do i = 1,n gpi = wd(i) - x(i) gpsupn = max( gpsupn, abs( gpi ) ) gpeucn2 = gpeucn2 + gpi ** 2 end do C Print information of this iteration if( iprint .ge. 2 ) then write(*, 983) iter write(*, 985) nprint,(x(i),i=1,nprint) write(*, 986) nprint,(g(i),i=1,nprint) write(*, 1002) f,sqrt(gpeucn2),gpsupn,fcnt,gcnt write(10,983) iter write(10,985) nprint,(x(i),i=1,nprint) write(10,986) nprint,(g(i),i=1,nprint) write(10,1002) f,sqrt(gpeucn2),gpsupn,fcnt,gcnt end if C ================================================================== C Test some stopping criteria that may occur inside the line C searches C ================================================================== if ( inform .eq. 6 ) then if ( iprint .ge. 2 ) then write(*, 996) inform,mininterp,epsrel,epsabs write(10,996) inform,mininterp,epsrel,epsabs end if go to 500 end if C ================================================================== C Iterate C ================================================================== go to 100 C ================================================================== C End of main loop C ================================================================== C ================================================================== C Report output status and return C ================================================================== 500 continue C Set f and x with the best solution and its functional value C f = fbest C do i = 1,n C x(i) = xbest(i) C end do C Print final information if( iprint .ge. 2 ) then write(*, 982) iter write(*, 985) nprint,(x(i),i=1,nprint) write(*, 986) nprint,(g(i),i=1,nprint) write(*, 1002) f,sqrt(gpeucn2),gpsupn,fcnt,gcnt write(10,982) iter write(10,985) nprint,(x(i),i=1,nprint) write(10,986) nprint,(g(i),i=1,nprint) write(10,1002) f,sqrt(gpeucn2),gpsupn,fcnt,gcnt end if return C Non-executable statements 977 format(/1X, 'Entry to SPG. Number of variables: ',I7) 980 format(/1X,'Initial point (first ',I6, ' components): ', */,6(1X,1PD11.4)) 981 format(/1X,'SPG iteration: ',I6, ' (Initial point)') 982 format(/1X,'SPG iteration: ',I6, ' (Final point)') 983 format(/,1X,'SPG iteration: ',I6) 985 format(1X,'Current point (first ',I6, ' components): ', */,6(1X,1PD11.4)) 986 format(1X,'Current gradient (first ',I6, ' components): ', */,6(1X,1PD11.4)) 990 format(/1X,'Flag of SPG = ',I3, *' (convergence with Euclidian-norm of the projected gradient', */,1X,'smaller than ',1PD11.4,')') 991 format(/1X,'Flag of SPG = ',I3, *' (convergence with sup-norm of the projected gradient', */,1X,'smaller than ',1PD11.4,')') 992 format(/1X,'Flag of SPG= ',I3, *' (The algorithm stopped by lack of enough progress. This means', */,1X,'that f(x_k) - f(x_{k+1}) .le. ',1PD11.4, *' * max [ f(x_j)-f(x_{j+1}, j < k ]',/,1X,'during ',I7, *' consecutive iterations') 993 format(/1X,'Flag of SPG = ',I3, *' (The algorithm stopped because the order of the', */,1X,'Euclidian-norm of the continuous projected gradient did', *' not change during ',/,1X,I7,' consecutive iterations.', *' Probably, an exaggerated small norm of the',/,1X,'continuous', *' projected gradient is required for declaring convergence') 994 format(/1X,'Flag of SPG = ',I3, *' (The algorithm stopped because the functional value is', */,1X,'smaller than ',1PD11.4) 996 format(/1X,'Flag of SPG = ',I3, *' (Too small step in a line search. After having made at ', */,1X,'least ',I7,' interpolations, the steplength becames small.', *' Small means that',/,1X,'we were at point x with direction d', *' and took a step alpha such that',/,1X,'alpha * |d_i| .lt.', *' max [',1PD11.4,' * |x_i|,',1PD11.4,' ] for all i') 997 format(/1X,'Flag of SPG = ',I3, *' (It was exceeded the maximum allowed number of iterations', */,1X,'(maxit=',I7,')') 998 format(/1X,'Flag of SPG = ',I3, *' (It was exceeded the maximum allowed number of functional', */,1X,'evaluations (maxfc=',I7,')') 1002 format(1X,'Functional value: ', 1PD11.4, */,1X,'Euclidian-norm of the continuous projected gradient: ', *1PD11.4, */,1X,'Sup-norm of the continuous projected gradient: ',1PD11.4, */,1X,'Functional evaluations: ',I7, */,1X,'Gradient evaluations: ',I7) 1000 format(/1X,'Flag of SPG = ',I3,' Fatal Error') end C ****************************************************************** C ****************************************************************** subroutine spgls(n,x,m,lambda,rho,equatn,f,g,lamspg,p,lastfv,nint, +mininterp,fmin,maxfc,iprint,fcnt,inform,xtrial,d,gamma,sigma1, +sigma2,epsrel,epsabs) C SCALAR ARGUMENTS integer fcnt,m,maxfc,mininterp,n,inform,iprint,p double precision epsabs,epsrel,f,fmin,gamma,lamspg,nint,sigma1, + sigma2 C ARRAY ARGUMENTS logical equatn(m) double precision d(n),g(n),lambda(m),lastfv(0:p-1),rho(m),x(n), + xtrial(n) C Safeguarded quadratic interpolation, used in the Spectral C Projected Gradient directions. C C On Entry C C n integer C the order of the x C C x double precision x(n) C current point C C m integer C lambda double precision lambda(m) C rho double precision rho(m) C These three parameters are not used nor modified by C SPG and they are passed as arguments to the user- C defined subroutines evalal and evalnal to compute the C objective function and its gradient, respectively. C Clearly, in an Augmented Lagrangian context, if SPG is C being used to solve the bound-constrainted subproblems, m C would be the number of constraints, lambda the Lagrange C multipliers approximation and rho the penalty parameters C C f double precision C function value at the current point C C g double precision g(n) C gradient vector at the current point C C lamspg double precision C spectral steplength C C nint double precision C constant for the interpolation. See the description of C sigma1 and sigma2 above. Sometimes we take as a new C trial step the previous one divided by nint C C RECOMMENDED: nint = 2.0 C C mininterp integer C constant for testing if, after having made at least C mininterp interpolations, the steplength is so small. In C that case failure of the line search is declared (may be C the direction is not a descent direction due to an error C in the gradient calculations) C C RECOMMENDED: mininterp = 4 C C fmin double precision C functional value for the stopping criterion f <= fmin C C maxfc integer C maximum number of functional evaluations C C iprint integer C Commands printing. Nothing is printed if iprint is C smaller than 2. If iprint is greater than or equal to C 2, SPG iterations information is printed. If iprint C is greater than or equal to 3, line searches information C is printed. C C RECOMMENDED: iprint = 2 C C CONSTRAINTS: allowed values are just 2 or 3. C C xtrial double precision xtrial(n) C d double precision d(n) C working vectors C C gamma double precision C constant for the Armijo criterion C f(x + alpha d) <= f(x) + gamma * alpha * <\nabla f(x),d> C C RECOMMENDED: gamma = 10^{-4} C C sigma1 double precision C sigma2 double precision C constant for the safeguarded interpolation C if alpha_new \notin [sigma1, sigma*alpha] then we take C alpha_new = alpha / nint C C RECOMMENDED: sigma1 = 0.1 and sigma2 = 0.9 C C sterel double precision C steabs double precision C this constants mean a ``relative small number'' and ``an C absolute small number'' for the increments in finite C difference approximations of derivatives C C RECOMMENDED: epsrel = 10^{-7}, epsabs = 10^{-10} C C epsrel double precision C epsabs double precision C infrel double precision C infabs double precision C this constants mean a ``relative small number'', ``an C absolute small number'', and ``infinite or a very big C number''. Basically, a quantity A is considered C negligible with respect to another quantity B if |A| < C max ( epsrel * |B|, epsabs ) C C RECOMMENDED: epsrel = 10^{-10}, epsabs = 10^{-20}, C infrel = 10^{+20}, infabs = 10^{+99} C C On Return C C x double precision C final estimation of the solution C C f double precision C functional value at the final estimation C C fcnt integer C number of functional evaluations used in the line search C C inform integer C This output parameter tells what happened in this C subroutine, according to the following conventions: C C 0 = convergence with an Armijo-like criterion C (f(xnew) <= f(x) + gamma * alpha * ); C C 4 = the algorithm stopped because the functional value C is smaller than fmin; C C 6 = too small step in the line search. After having made C at least mininterp interpolations, the steplength C becames small. ''small steplength'' means that we are C at point x with direction d and step alpha, and, for C all i, C C | alpha * d(i) | <= max ( epsrel * |x(i)|, epsabs ). C C In that case failure of the line search is declared C (maybe the direction is not a descent direction due C to an error in the gradient calculations). Use C mininterp > maxfc to inhibit this criterion; C C 8 = it was achieved the maximum allowed number of C function evaluations (maxfc); C C < 0 = error in evalf subroutine. C LOCAL SCALARS logical samep integer i,interp,flag double precision alpha,atmp,fmax,ftrial,gtd C EXTERNAL SUBROUTINES external evalal,proj C INTRINSIC FUNCTIONS intrinsic abs,max C Initialization interp = 0 C Compute fmax fmax = lastfv(0) do i = 1,p - 1 fmax = max( fmax, lastfv(i) ) end do C Compute first trial point, spectral projected gradient direction, C and directional derivative . alpha = 1.0d0 do i = 1,n xtrial(i) = x(i) - lamspg * g(i) end do call proj(n,xtrial,flag) if ( flag .ne. 0 ) then inform = - 92 if ( iprint .ge. 2 ) then write(*, 1000) inform write(10,1000) inform end if return end if gtd = 0.0d0 do i = 1,n d(i) = xtrial(i) - x(i) gtd = gtd + g(i) * d(i) end do C Print presentation information if ( iprint .ge. 3 ) then write(*, 980) lamspg,gtd write(10,980) lamspg,gtd end if call evalal(n,xtrial,m,lambda,rho,equatn,ftrial,inform) fcnt = fcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if C Print information of the first trial if ( iprint .ge. 3 ) then write(*, 999) alpha,ftrial,fcnt write(10,999) alpha,ftrial,fcnt end if C Main loop 100 continue C Test Armijo stopping criterion if ( ftrial .le. fmax + gamma * alpha * gtd ) then f = ftrial do i = 1,n x(i) = xtrial(i) end do inform = 0 if ( iprint .ge. 3 ) then write(*, 990) inform write(10,990) inform end if go to 500 end if C Test whether f is very small if ( ftrial .le. fmin ) then f = ftrial do i = 1,n x(i) = xtrial(i) end do inform = 4 if ( iprint .ge. 3 ) then write(*, 994) inform write(10,994) inform end if go to 500 end if C Test whether the number of functional evaluations is exhausted if ( fcnt .ge. maxfc ) then if ( ftrial .lt. f ) then f = ftrial do i = 1,n x(i) = xtrial(i) end do end if inform = 8 if ( iprint .ge. 3 ) then write(*, 998) inform write(10,998) inform end if go to 500 end if C Compute new step (safeguarded quadratic interpolation) interp = interp + 1 if ( alpha .lt. sigma1 ) then alpha = alpha / nint else atmp = ( - gtd * alpha ** 2 ) / + ( 2.0d0 * ( ftrial - f - alpha * gtd ) ) if ( atmp .lt. sigma1 .or. atmp .gt. sigma2 * alpha ) then alpha = alpha / nint else alpha = atmp end if end if C Compute new trial point do i = 1,n xtrial(i) = x(i) + alpha * d(i) end do call evalal(n,xtrial,m,lambda,rho,equatn,ftrial,inform) fcnt = fcnt + 1 if ( inform .lt. 0 ) then if ( iprint .ge. 3 ) then write(*, 1000) inform write(10,1000) inform end if return end if C Print information of the current trial if ( iprint .ge. 3 ) then write(*, 999) alpha,ftrial,fcnt write(10,999) alpha,ftrial,fcnt end if C Test whether at least mininterp interpolations were made and two C consecutive iterates are close enough samep = .true. do i = 1,n if ( abs( alpha * d(i) ) .gt. + max( epsrel * abs( x(i) ), epsabs ) ) then samep = .false. end if end do if ( interp .ge. mininterp .and. samep ) then if ( ftrial .lt. f ) then f = ftrial do i = 1,n x(i) = xtrial(i) end do end if inform = 6 if ( iprint .ge. 3 ) then write(*, 996) inform write(10,996) inform end if go to 500 end if C Iterate go to 100 C Return 500 continue return C Non-executable statements 980 format(/,6x,'SPG (spectral steplength ',1PD11.4,1x,' = ', * 1PD11.4,')',/, * /,6x,'SPG Line search') 999 format(6x,'Alpha= ',1PD11.4,' F= ',1PD11.4,' FE= ',I5) 990 format(6x,'Flag of SPG Line search = ',I3, * ' (Convergence with an Armijo-like criterion)') 994 format(6x,'Flag of SPG Line search = ',I3, * ' (Small functional value, smaller than ',/, * 6X,'parameter fmin)') 996 format(6x,'Flag of SPG Line search = ',I3, * ' (Too small step in the interpolation)') 998 format(6x,'Flag of SPG Line search = ',I3, * ' (Too many functional evaluations)') 1000 format(6x,'Flag of SPG Line search = ',I3,' Fatal Error') end C ****************************************************************** C ****************************************************************** subroutine evalnaldiff(n,x,m,lambda,rho,equatn,g,sterel,steabs, +inform) C SCALAR ARGUMENTS integer n,m,inform double precision sterel,steabs C ARRAY ARGUMENTS logical equatn(m) double precision x(n),lambda(m),rho(m),g(n) C Approximates the gradient vector g(x) of the objective function by C central finite differences. This subroutine, which works in the C full space, is prepared to replace the subroutine evalnal (to C evaluate the gradient vector) in the case of the lastest have not C being provided by the user. C C Parameters of the subroutine: C C On entry: C C n integer C number of variables C C x double precision x(n) C current point C C m integer C lambda double precision lambda(m) C rho double precision rho(m) C equatn logical equatn(m) C These four parameters are not used nor modified by C GENCAN and they are passed as arguments to the user- C defined subroutines evalal and evalnal to compute the C objective function and its gradient, respectively. C Clearly, in an Augmented Lagrangian context, if GENCAN is C being used to solve the bound-constrainted subproblems, m C would be the number of constraints, lambda the Lagrange C multipliers approximation and rho the penalty parameters. C equatn is logical array that, for each constraint, C indicates whether the constraint is an equality constraint C (.true.) or an inequality constraint (.false.) C C sterel double precision C See below C C steabs double precision C This constants mean a ''relative small number'' and ''an C absolute small number'' for the increments in finite C difference approximations of derivatives C C RECOMMENDED: epsrel = 1.0d-07 and epsabs = 1.0d-10 C C CONSTRAINTS: sterel >= steabs > 0 C C On Return C C g double precision g(n) C approximation of the gradient vector at x C C inform integer C 0 = no errors, C < 0 = there was an error in the gradient calculation. C LOCAL SCALARS integer j double precision tmp,step,fplus,fminus inform = 0 do j = 1,n tmp = x(j) step = max( steabs, sterel * abs( tmp ) ) x(j) = tmp + step call evalal(n,x,m,lambda,rho,equatn,fplus,inform) if ( inform .lt. 0 ) then return end if x(j) = tmp - step call evalal(n,x,m,lambda,rho,equatn,fminus,inform) if ( inform .lt. 0 ) then return end if g(j) = ( fplus - fminus ) / ( 2.0d0 * step ) x(j) = tmp end do return end C ****************************************************************** C ****************************************************************** double precision function drand(ix) C This is the random number generator of Schrage: C C L. Schrage, A more portable Fortran random number generator, ACM C Transactions on Mathematical Software 5 (1979), 132-138. double precision ix double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/ xhi= ix/b16 xhi= xhi - dmod(xhi,1.d0) xalo= (ix-xhi*b16)*a leftlo= xalo/b16 leftlo= leftlo - dmod(leftlo,1.d0) fhi= xhi*a + leftlo k= fhi/b15 k= k - dmod(k,1.d0) ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k if (ix.lt.0) ix= ix + p drand= ix*4.656612875d-10 return end