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                  <real-value> 
C     OPTIMALITY-TOLERANCE                   <real-value> 
C     MAX-OUTER-ITERATIONS                   <integer-value>
C     MAX-INNER-ITERATIONS                   <integer-value>
C     MAX-FUNCTION-EVALUATIONS               <integer-value>
C     OUTPUT-DETAIL                          <integer-value>
C     NCOMP-ARRAY                            <integer-value>
C     MLCOMP-ARRAY                           <integer-value>
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
