C ================================================================= C File: genlinma.f C ================================================================= C ================================================================= C Module: Main program C ================================================================= C Last update of any of the component of this module: C C April 30, 2008. 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 GENLIN solves problems of the form C ------------------------------------ C C min f(x) C C subject to C C Ai x = bi, 0 <= i <= mleq, C Ci x >= di, mleq+1 <= i <= ml, C l <= x <= u, C C with ml >= 1 linear constraints and n variable. C C GENLIN is a method to solve linearly constrained problems based on C GENCAN. C C GENLIN 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 GENLIN must cite: C C M. Andretta, E. G. Birgin and J. M. Martínez, "Partial Spectral C Projected Gradient Method with Active-Set Strategy for Linearly C Constrained Optimization", submitted, 2008. C C Every published work that uses ALGENCAN must cite: 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 and C C R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, C "Augmented Lagrangian methods under the Constant Positive Linear C Dependence constraint qualification", to appear in Mathematical C Programming, 2005. C C Every published work that uses GENCAN must cite: C C E. G. Birgin and J. M. Martínez, "Large-scale active-set C box-constrained optimization method with spectral projected C gradients", Computational Optimization and Applications 23, pp. C 101-125, 2002. C C (See other related works at the TANGO home page.) C ***************************************************************** C HOW TO USE GENLIN TO SOLVE YOUR OWN PROBLEM C -------------------------------------------- C C You will find below 7 subroutines and that YOU SHOULD 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 linear equality C constraints (mleq), the total number of linear constraints (ml) C (ml must be at least one), the coeficient matrix A and the right C hand vector b such that each row of A (lets say Ai) represent a C constraint Ai x >= b(i). The first mleq rows of A (e components C of b) must contain the coeficients of the equality linear C constraints. You must also set inform to 0 if the problem is C well defined. If an error occurs during the definition of the C problem, set inform to a non-zero value. C C 2) SUBROUTINE EVALF MUST also be MODIFIED to evaluate C the objective function. C C 3) SUBROUTINE EVALG, to compute the gradient vector of the C objective function is NOT MANDATORY but HIGHLY RECOMMENDED. It C must be provided by the user if he/she set gtype = 0. If the C user sets gtype = 1 then finite differences approximations will C be used. (See the description of gtype parameter in subroutine C param.) C C 4) OPTIONALY, you can modify SUBROUTINES EVALH and EVALSH to C evaluate the Hessian matrix of the objective function in the C dense or sparse form, respectively. EVALH is MANDATORY if the C user chooses to use BETRALIN (intype = 2) instead of GENLIN. C If GENLIN is chosen (intype = 1), then EVALSH is required when C hptype = 0 and optional otherwise. (See the description of hptype C parameter in subroutine param.) C C 5) OPTIONALY, another way of using second-order information when C GENLIN is chosen, is to provide, not individual subroutines to C compute the Hessian of the objective function, but to provide a C SUBROUTINE EVALHP to compute the product of an arbitrary vector C times the Hessian matrix. There is not a clear advantage in coding C EVALHP instead of coding EVALH. This alternative was incorporated C to be primarily used in the CUTEr interfaces. This subroutine is C optional and is required when hptype = 1. (See the description of C hptype parameter in subroutine param.) C C 6) OPTIONALLY, in the SUBROUTINE PARAM you may modify some C arguments like feasibility and optimality tolerances and maximum C number of iterations and functional evaluations. See the detailed C description of each argument in subroutine PARAM. If you prefer, C leave all these parameters with the suggested values. C C 7) Alternatively, if you are interested in modifying any of the C parameters described in (5) but you prefer not to modify subroutine C PARAM, you can use the genlin.dat SPECIFICATION FILE. The C advantage of using the specification file is that you can test C the method with different parameters without compiling it every C time. This file IS NOT MANDATORY. So, if you will not modify the C default parameters or if you fill more comfortable modifying C subroutine PARAM, you do not need to worry about genlin.dat C specification file. C C The specification file can have any number of lines, with at most C 80 characters per line. Blank lines are admitted and lines starting C with '*' or '#' are considered comment lines and will be ignored. C Each line must have just a "keyword" or a keyword followed by an C integer or a real value when required. As the interpreter is not C case-sensitive, you can write the keywords in lower case, upper C case or any combination of them. You can also insert blanks in any C place before or after the keywords. Moreover, just the first 10 C characters of the keyword are relevant and the rest will be ignored. C Available keywords are: C C ANALYTIC-GRADIENT C FINITE-DIFFERENCES-GRADIENT C HESSIANS-PROVIDED C LAGRHESS-PRODUCT-PROVIDED C INCREMENTAL-QUOTIENTS C BFGS-QN-APPROXIMATION C ADAPTIVE-HESSIAN C AUTO-BDSOLVER C GENCAN-BDSOLVER C BETRA-BDSOLVER C SPARSE-BETRA-BDSOLVER C UNPRECONDITIONED-CG C BFGS-QN-PRECONDITIONER C CHECK-DERIVATIVES C AVOID-LOSS-OF-FEASIBILITY C FEASIBILITY-TOLERANCE C OPTIMALITY-TOLERANCE C MAX-OUTER-ITERATIONS C MAX-INNER-ITERATIONS C MAX-FUNCTION-EVALUATIONS C OUTPUT-DETAIL C NCOMP-ARRAY C MLCOMP-ARRAY C C By default, GENLIN uses: C C ANALYTIC-GRADIENT C ADAPTIVE-HESSIAN C BFGS-QN-PRECONDITIONER C AUTO-BDSOLVER C BFGS-QN-PRECONDITIONER C FEASIBILITY-TOLERANCE 1.0d-08 C OPTIMALITY-TOLERANCE 1.0d-08 C MAX-INNER-ITERATIONS 200000 C MAX-FUNCTION-EVALUATIONS 1000000 C OUTPUT-DETAIL 1 C NCOMP-ARRAY 5 C MLCOMP-ARRAY 5 C C and derivatives are not tested. C C 8) Finally, SUBROUTINE ENDP can be modified by the user to write, C save or compute any type of extra information related to the C solution of the problem. Subroutine ENDP will be called just once C after algencan found a solution of the problem. You can modify C this subroutine for, for example, save the solution in a file C chosed by you and in a special format to be processed later by C another software. C C Besides the modifications that you must do to solve your own C problem using GENLIN, GENLIN 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 genlin. C ****************************************************************** C ****************************************************************** program genlinma implicit none C ****************************************************************** C NOW, WE SET THE MAXIMUM NUMBER OF VARIABLES (NMAX) AND THE MAXIMUM C NUMBER OF LINEAR CONSTRAINTS (MLMAX) THAT THE PROBLEM IS ALLOWED c TO HAVE. IF YOUR PROBLEM HAS MORE THAN THE DECLARED VALUES FOR C NMAX AND MLMAX THEN YOU MUST MODIFY THESE VALUES IN SUCH A WAY C THAT THE NUMBER OF VARIABLES (N) OF YOUR PROBLEM BE SMALLER THAN C OR EQUAL TO NMAX AND THE NUMBER OF LINEAR CONSTRAINTS (ML) OF YOUR C PROBLEM BE SMALLER THAN OR EQUAL TO MLMAX. C ****************************************************************** C PARAMETERS integer mlmax,nmax parameter ( mlmax = 2000 ) parameter ( nmax = 1550 ) integer ldh,lda,ldabar parameter ( ldh = nmax ) parameter ( lda = mlmax ) parameter ( ldabar = nmax ) integer lwd24 parameter ( lwd24 = 2500500 ) C LOCAL SCALARS logical checkder,lossfeas character * 6 precond integer gtype,hptype,inform,intype,iprint,maxtotfc,maxtotit,n, + ncomp,totcgcnt,totfcnt,totgcnt,tothcnt,totchcnt,totiter, + ml,mleq,mlcomp double precision epsfeas,epsopt,f,gpsupn,cnorm real time C LOCAL ARRAYS integer wi1(nmax),wi2(nmax),wi3(mlmax),wi4(2*nmax+mlmax), + wi5(nmax+mlmax) double precision l(nmax),u(nmax),wd1(nmax),wd2(nmax),wd3(nmax), + wd4(nmax),wd5(nmax),wd6(nmax),wd7(nmax),wd8(nmax), + wd9(nmax),wd10(nmax),wd11(nmax),wd12(nmax),wd13(nmax), + wd14(nmax),wd15(nmax),wd16(nmax),wd17(nmax),wd18(nmax), + wd19(nmax),wd20(nmax),wd21(nmax),wd22(mlmax),wd23(mlmax), + wd24(lwd24),x(nmax),A(lda,nmax),Abar(ldabar,mlmax), + H(ldh,nmax),b(mlmax) C EXTERNAL SUBROUTINES external solver C SET UP PROBLEM DATA call inip(n,x,l,u,ml,mleq,lda,A,b,inform) if ( inform .ne. 0 ) then write(*,*) 'An error occured during problem definition' stop end if C SET SOME SOLVER ARGUMENTS call param(gtype,hptype,intype,precond,checkder,lossfeas,epsfeas, +epsopt,maxtotit,maxtotfc,iprint,ncomp,mlcomp) C CALL OPTIMIZATION SOLVER call solver(n,x,l,u,ml,mleq,lda,A,b,gtype,hptype,intype,precond, +checkder,lossfeas,epsfeas,epsopt,maxtotit,maxtotfc,iprint,ncomp, +mlcomp,f,cnorm,gpsupn,totiter,totfcnt,totgcnt,tothcnt,totcgcnt, +totchcnt,time,inform,wi1,wi2,wi3,wi4,wi5,wd1,wd2,wd3,wd4,wd5,wd6, +wd7,wd8,wd9,wd10,wd11,wd12,wd13,wd14,wd15,wd16,wd17,wd18,wd19, +wd20,wd21,wd22,wd23,wd24,lwd24,ldabar,Abar,ldh,H) C WRITE ADDITIONAL OUTPUT INFORMATION CODED BY THE USER call endp(n,x,l,u,ml,mleq,lda,A,b) stop end C ****************************************************************** C ****************************************************************** subroutine param(gtype,hptype,intype,precond,checkder,lossfeas, +epsfeas,epsopt,maxtotit,maxtotfc,iprint,ncomp,mlcomp) C SCALAR ARGUMENTS character * 6 precond logical checkder,lossfeas integer gtype,hptype,iprint,intype,maxtotfc,maxtotit,ncomp,mlcomp double precision epsfeas,epsopt C Parameters of the subroutine: C ============================= C C On Return: C ========== C C GTYPE integer C ------------- C C Type of first derivatives calculation according to the following C convention: C C 0 means true first derivatives. In this case, subroutines evalg C and evaljac must be modified by the user to compute the C gradient of the objective function and the gradients of the C constraints, respectively. C C 1 means that a finite difference approximation will be used. In C this case, subroutines evalg and evaljac may have an empty body C but must be present. It is also recommended that those C empty-body subroutines set flag = - 1. Last but not least, the C option gtype = 1 is not cheap neither safe. C C HPTYPE integer C -------------- C C Type of Hessian-vector product according to the following C convention: C C The way in which the product of the Hessian matrix by a vector will C be done depends on the value of the parameter hptype in the following C way: C C 0 means that the subroutine to compute the Hessian of the C objective function (evalh) was provided by the user. So, the C product of the Hessian times a vector will be computed using the C Hessian provided by this subroutine. C C 1 means that a subroutine (evalhp) to compute the product of the C Hessian of the objective function times an arbitrary vector is C being provided by the user. C C 9 means that an incremental quotients approximation will be used. C This option requires the evaluation of an extra gradient at each C Conjugate Gradient iteration. If gtype = 0 then this gradient C evaluation will be done using the user supplied subroutine C evalg. On the other hand, if gtype = 1, the gradient calculation C will be done using just calls to the user provided subroutine C evalf. nind calls will be done, where nind is the dimension of C the current face of the active-set method. C C If you did not code subroutine evalg to compute the gradient of C the objective function then your options finished here. C C 4 means that the Hessian matrix will be approximated and then the C product of the Hessian approximation by the vector will be C computed exactly. The approximation is a BFGS approximation of C the Hessian (without memory) and using the spectral C approximation (instead of the identity) as initial C approximation. C C Numerical experiments suggested that this option is not C convenient for unconstrained or just linearly-constrained problems. C (Note that this option was developed to be used in the Augmented C Lagrangian framework.) C C This option does NOT require an extra gradient evaluation per C iteration and, in this sense, each CG iteration is C computationally cheaper than a CG iteration of the previous C choices. C C Quadratic subproblems are convex with this choice. C C In the linearly-constrained context, options hptype = 2,3,5 and 6 C are all identical to hptype = 9. C C INTYPE integer C -------------- C C This parameter can be used to select the algorithm that will be C used to solve the unconstrained subproblems within the faces. C C 1 means GENLIN will be used (that is, GENCAN will be used within the C face) C 2 means BETRALIN will be used (that is, BETRA will be used within C the face) C C PRECOND character * 6 C --------------------- C C Indicates the type of preconditioning that will be used for C Conjugates Gradients according to the following convention: C C 'NONE' means no preconditioner at all. C C 'QNCGNA' means Quasi-Newton Correction of the Gauss-Newton C approximation of the Hessian. The exact form is this C preconditioner is described in: C C E. G. Birgin and J. M. Martínez, "Structured minimal- C memory inexact quasi-Newton method and secant C preconditioners for Augmented Lagrangian Optimization", C submitted, 2005. C C CHECKDER logical C ---------------- C C If you are using finite differences aproximations for the C derivatives (gtype = 1) then checkder must be set to FALSE. C On the other hand, if you are using your own coded derivatives C you may would like to test them against a finite differences C approximation to verify if they are o.k. In this case, set C checkder = TRUE. C C LOSSFEAS logical C ------------------------ C C If a loss of feasibility up to sqrt(epsfeas) is allowed during the C computations, lossfeas must be set to TRUE (recommended). In this C case, if the final iterate does not satisfy the constraints with C the desired precision, we sugest that the user re-run the method C using the final point obtained as the initial point, using C lossfeas = FALSE. C C EPSFEAS double precision C ------------------------ C C Feasibility tolerance for the sup-norm of the linear constraints. C C EPSOPT double precision C ----------------------- C C Optimality tolerance for the sup-norm of the projected gradient of C the objective function. C C MAXTOTIT integer C ---------------- C C Maximum total number of iterations. C C MAXTOTFC integer C ---------------- C C Idem MAXTOTIT but for number of functional evaluations. C C IPRINT integer C -------------- C C Controls the ammount of information of the output according to the C following convention: C C 0 means no output at all. C C 1 means information at each outer iteration but without any C information of how the subproblems are being solved. C C 2 means the same as 1 plus information of each inner iteration. C C 3 means the same as 2 plus information of the line searches and C the calculation of the truncated Newton direction (using CG) of C each inner iteration. C C In all cases, an output file named solution.txt with the final C point, Lagrange mutipliers and penalty parameters will be C generated. Moreover, the same output of the screen will be saved C in a file named algencan.out. C C NCOMP integer C ------------- C C Every time a vector is printed, just its first ncomp component C will be displayed. C C MLCOMP integer C -------------- C C Every time the coeficients of a linear constraint is printed, just C its first mlcomp component will be displayed. C ****************************************************************** C FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET SOME GENLIN C ARGUMENTS RELATED TO DERIVATIVES, STOPPING CRITERIA AND OUTPUT: C ****************************************************************** gtype = 0 hptype = 0 intype = 1 precond = 'QNCGNA' checkder = .false. lossfeas = .true. epsfeas = 1.0d-08 epsopt = 1.0d-08 maxtotit = 200000 maxtotfc = 5 * maxtotit iprint = 1 ncomp = 5 mlcomp = 5 C ****************************************************************** C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE PARAM. C ****************************************************************** end