C     =================================================================
C     Module: Linear-constraints solver GENLIN
C     =================================================================

C     Last update of any of the component of this module: 
C
C     September 19, 2006.

C     ******************************************************************
C     ******************************************************************

      subroutine easygenlin(n,x,l,u,mp,lambdap,rhop,equatnp,linearp,ml,
     +meq,alin,acol,aval,annz,b,epsgpsn,fmin,iprint,ncomp,macheps,
     +bignum,f,g,gpsupn,ldh,H,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer annz,ldh,ml,meq,mp,n,ncomp,inform,iprint
      double precision bignum,epsgpsn,f,fmin,gpsupn,macheps

C     ARRAY ARGUMENTS
      logical equatnp(mp),linearp(mp)
      integer alin(annz),acol(annz)
      double precision x(n),l(n),u(n),lambdap(mp),rhop(mp),b(ml),g(n),
     +        H(ldh,n),aval(annz)

C     This subroutine aims to simplify the use of GENLIN. For this 
C     purpose it gives values to most of the GENLIN arguments and 
C     leaves to the user those arguments which he/she may would like to 
C     set by him/herself.
C
C     The arguments of EASYGENLIN are the input and output arguments of 
C     GENLIN that are supposed to be useful for a common user. The input 
C     arguments are mostly related to basic problem information, like 
C     dimension and bounds, and the initial point. There are also input 
C     arguments related to simple stopping criteria (like norm of the 
C     projected gradient, and maximum number of iterations and 
C     functional evaluations). There are also two input arguments 
C     related to control the amount of information written into the 
C     screen. The output arguments are related to information of the 
C     solution and some few performance measurements. Basically, on 
C     return, EASYGENLIN gives to the user the solution, the objective 
C     functional value and its gradient at the solution, Euclidian and 
C     sup-norm of the projected gradient at the solution, the number of 
C     iterations, functional and gradient evaluations, and Conjugate 
C     Gradient iterations used to reach the solution, and, finally, a 
C     flag that indicates the stopping criterion that was satisfied.
C
C     All the other arguments of GENLIN are setted with its default 
C     values by EASYGENLIN. EASYGENLIN divides the arguments of GENLIN 
C     in two sets. Those that are related to the behaviour of GENLIN are 
C     declared as Fortran parameters (constants). The other arguments of 
C     GENLIN, most of them related to alternative stopping criteria, and 
C     that may depend of, for example, maxit, are declared as local 
C     variables of EASYGENLIN.
C
C     GENLIN arguments that are defined as Fortran parameters in this 
C     subroutine are GENLIN arguments that should not be modified by a 
C     common user. They are arguments that modify the behaviour of 
C     GENLIN and whos values were selected because they are classical 
C     values in some cases or because some numerical experiments seemed 
C     to indicate that they are the best choices.
C
C     GENLIN arguments that are declared as local variables in this 
C     subroutine are GENLIN arguments that may be modified if, with 
C     their suggested values, GENLIN does not give the desired result. 
C     Most of them are related to Conjugate Gradients or to disabled 
C     stopping criteria that may be useful in bad-scaled problems or 
C     problems with not trustable derivatives.
C
C     Finally, this subroutine declares as local variables some 
C     arguments of GENLIN which in fact are output arguments. Most of 
C     them are related to quantities that can be used for statistics 
C     related to the GENLIN performance, like number Spectral Projected 
C     Gradient iterations, Truncated Newton iterations, Conjugate 
C     Gradient iterations, etc. As we assume that this values are not 
C     useful for the common user, this subroutine throw all of them 
C     away.

C     PARAMETERS
      include "dim.f"

      integer ldwd17
      parameter ( ldwd17 = 1000000 )

C     COMMON SCALARS
      integer m

C     COMMON ARRAYS
      logical equatn(mmax),linear(mmax)
      double precision lambda(mmax),rho(mmax)

C     COMMON BLOCK
      common /aldata/ lambda,rho,equatn,linear,m

C     This parameters should be arguments of easygenlin, as a common
C     user probably would like to modify them. However, for this
C     particular application, they will be fixed here.

      integer maxit,maxfc,intype,gtype,hptype
      parameter ( maxit     =      1000 )
      parameter ( maxfc     = 5 * maxit )
      parameter ( intype    =         0 )
      parameter ( gtype     =         0 )
      parameter ( hptype    =         1 )

C     These parameters might be related to the inner BETRA algorithm.
C     So, Marina will describe them here.

      logical scale
      integer msmaxit,ktest,promaxit
      double precision epsgzsn,scmaxcnd,scmaxal,rhoms,mseps,phieps,
     +        sigms,trgamma,delini,trdelmin,trdelmax,projeps,ceps      
      parameter ( scale     =   .false. )
      parameter ( msmaxit   =       100 )
      parameter ( ktest     =         1 )
      parameter ( promaxit  =      5000 )
      parameter ( epsgzsn   =   1.0d-05 )
      parameter ( scmaxcnd  =   0.2d+00 )
      parameter ( scmaxal   =   1.0d+05 )
      parameter ( rhoms     =   9.0d-01 )
      parameter ( mseps     =   1.0d-05 )
      parameter ( phieps    =   1.0d-05 )
      parameter ( sigms     =   2.0d-01 )
      parameter ( trgamma   =   1.0d-01 )
      parameter ( delini    =   1.0d+02 )
      parameter ( trdelmin  =   1.0d-08 )
      parameter ( trdelmax  =   1.0d+15 )
      parameter ( projeps   =   1.0d-08 )
      parameter ( ceps      =   1.0d-05 )

C     CONSTANTS FOR CLASSICAL LINE-SEARCH CONDITIONS

C     beta is the constant for the ''beta condition''. We use this 
C     condition to test whether is promising to extrapolate or not.

C     gamma is the constant for the sufficient decrease ''Armijo 
C     condition''.

C     theta is the constant for the ''angle condition''.

C     sigma1 and sigma2 are the constants for the safeguarding quadratic 
C     interpolations. We use them in a rather unusual way. Instead of 
C     discarding a new step anew if it does not belong to the interval 
C     [ sigma1 * aprev, sigma2 * aprev ], we discard it if it does not 
C     belong to the interval [ sigma1, sigma2 * aprev ]. In such a case 
C     we take something similar to ''anew = aprev / 2''.

      double precision beta,gamma,theta,sigma1,sigma2
      parameter ( beta   =   0.5d0 )
      parameter ( gamma  = 1.0d-04 )
      parameter ( theta  = 1.0d-06 )
      parameter ( sigma1 =   0.1d0 )
      parameter ( sigma2 =   0.9d0 )

C     CONSTANTS FOR SPECIFIC PROCEDURES (NOT SO CLASSICAL)

C     In line searches, when interpolating, the step may become so 
C     small that we should declare a line search failure indicating that 
C     direction may not be a descent direction. This decision is never 
C     take before doing at least mininterp interpolations.

C     In line searches, the beta condition (see above) may recommend to
C     extrapolate. We never do more than maxextrap extrapolations.

C     In the line searches, when we need to interpolate and the result 
C     of the quadratic interpolation is rejected, the new step is 
C     computed as anew = aprev / etaint. When the beta condition 
C     recommends to extrapolate, we compute anew = aprev * etaext.

C     When computing the Newton direction by Conjugate Gradients we 
C     never go further an artificial ''trust region''. This ''trust 
C     radius'' is never smaller than delmin.

C     In active set strategies, constants eta is used to decide whether 
C     the current face should be abandoned or not. In particular, the 
C     current face is abandoned when the norm of the internal to face 
C     component of the continuous projected gradient is smaller than 
C     ( 1 - eta ) times the norm of the continuous projected gradient. 
C     In this way, values of eta near 1 makes the method to work hard 
C     inside the faces and values of eta near 0 makes the method to 
C     abandon the faces very quickly.

C     We always use as a first step in a line search procedure along a
C     first order direction the spectral steplength. This steplength 
C     must belong to the interval [lspgmi,lspgma].

      integer maxextrap,mininterp
      parameter ( maxextrap =     0 )
      parameter ( mininterp =     4 )

      double precision etaint,etaext,delmin,eta,lspgma,lspgmi
      parameter ( etaint  =   2.0d0 )
      parameter ( etaext  =   2.0d0 )
      parameter ( delmin  = 1.0d+04 )
      parameter ( eta     =   0.1d0 )
      parameter ( lspgma  = 1.0d+10 )
      parameter ( lspgmi  = 1.0d-10 )

C     GENLIN INPUT ARGUMENTS THAT WILL BE SETTED BELOW
      integer cgscre,maxitnfp,maxitngp,maxitnqmp,trtype
      double precision cgepsf,cgepsi,cggpnf,cgmia,cgmib,delta0,epsgpen,
     +        epsnfp,epsnqmp

C     GENLIN OUTPUT ARGUMENTS THAT WILL BE DISCARDED
      integer nit,nef,neg,neh,ncg,nch,spgnit,spgnef,tnnit,tnnef,
     +        tnintcnt,tnexgcnt,tnexbcnt,tnstpcnt,tnintfe,tnexgfe,
     +        tnexbfe,spginit,spginef,trnit,trnef,trexgcnt,trexbcnt,
     +        trexgfe,trexbfe
      double precision gpeucn2

C     LOCAL ARRAYS (WORKING ARRAYS FOR GENLIN)
      integer wi1(nmax),wi2(nmax),wi3(mlmax),wi4(nmax+mlmax),
     +        wi5(mlmax+2*nmax)
      double precision A(mlmax,nmax),Abar(nmax,mlmax),wd1(nmax),
     +        wd2(nmax),wd3(nmax),wd4(nmax),wd5(nmax),wd6(nmax),
     +        wd7(nmax),wd8(mlmax),wd9(nmax),wd10(mlmax),wd11(nmax),
     +        wd12(nmax),wd13(nmax),wd14(nmax+mlmax),wd15(nmax),
     +        wd16(nmax),wd17(ldwd17),wd18(nmax),wd19(nmax),wd20(nmax),
     +        wd21(mlmax)

C     LOCAL SCALARS
      integer i,j

C     ARGUMENTS RELATED TO STOPPING CRITERIA

C     Besides the stopping criterion related to the sup-norm of the 
C     continuous projected gradient, there is another stopping criterion 
C     related to its Euclidian norm. So, GENLIN stops the process if it 
C     finds a point at which the Euclidian norm of the continuous 
C     projected gradient is smaller than epsgpen.

      epsgpen   =    0.0d0

C     Sometimes, is the problem is bad scaled, to request a small 
C     gradient norm at the solution may be inadequate. For this reason, 
C     if this norm is not decreasing during maxitngp (MAXimum of 
C     ITerations with No Gradient Progress) consecutive iterations then 
C     we stop the method with a warning. 
     
      maxitngp  =    maxit

C     maxitnfp means MAXimum of allowed number of iterations with No 
C     Progress in the objective functional value. ''Progress'' from one 
C     iteration to the next one refers to ( fnew - fprev ). Since the 
C     begining of the algorithm we save the ''best progress'' and 
C     consider that there was no progress in an iteration if the 
C     progress of this iterations was smaller than epsnfp times the best 
C     progress. Finally, the algorithm stops if there was no progress 
C     during maxitnfp consecutive iterations.

      maxitnfp  =    maxit
      epsnfp    =    0.0d0

C     ARGUMENTS RELATED TO CONJUGATE GRADIENTS

C     When computing the Truncated Newton direction by Conjugate 
C     Gradients there is something similar to a ''trust-region radius''. 
C     This trust radius is updated from iteration to iteration depending 
C     on the agreement of the objective function and its quadratic 
C     model. But an initial value for the trust radius is required. If 
C     the user has a good guess for this initial value then it should be 
C     passed to GENLIN using the delta0 arguments. On the other hand, if 
C     delta0 is set to -1, a default value depending on the norm of the 
C     current point will be used.

      delta0    =  - 1.0d0

C     The ''trust-region'' can be like a ball (using Euclidian norm) or 
C     like a box (using sup-norm). This choice can be made using trtype 
C     (TRust region TYPE) argument. trtype equal to 0 means Euclidian 
C     norm and trtype equal to 1 means sup-norm.

      trtype    =        0

C     When the method is far from the solution, it may be not useful to 
C     do a very large effort in computing the Truncated Newton direction 
C     precisely. To avoid it, a fixed maximum number of iterations for 
C     Conjugate Gradients can be given to GENLIN. If the user would like 
C     to choose this maximum number of iterations for Conjugate 
C     Gradient then it should use cgmia and cgmib arguments. On the other 
C     hand, if he/she prefers to leave this task to GENLIN then he/she 
C     should set cgmia = -1.0d0 and cgmib = -1.0d0.
 
      cgmia     =  - 1.0d0
      cgmib     =  - 1.0d0

C     If the task of deciding the accuracy for computing the Truncated 
C     Newton direction is leaved to GENLIN then a default strategy based 
C     on increasing accuracies will be used. The proximity to the 
C     solution is estimated observing the norm of the projected gradient 
C     at the current point and locating it between that norm at the 
C     initial point and the expected value of that norm at the solution. 
C     Then the accuracy for the Truncated Newton direction of the 
C     current iteration will be computed taking a precision located in 
C     the same relative position with respect to two given values for 
C     the accuracies for the first and the last Truncated Newton 
C     direction calculations. These two accuracies (cgepsi and cgepsf, 
C     respectively) must be given by the user. Moreover, the expected 
C     value of the projected gradient norm at the solution (cggpnf) must
C     also be given by the user who must indicate setting argument 
C     cgscre to 1 or 2 if that norm is the Euclidian or the sup-norm.
      
      cggpnf    =  max( 1.0d-04, max( epsgpen, epsgpsn ) ) 
      cgscre    =        2
      cgepsi    =  1.0d-01
      cgepsf    =  1.0d-05

C     The next two arguments are used for an alternative stopping 
C     criterion for Conjugate Gradients. Conjugate Gradients method is 
C     stopped if the quadratic model makes no progress during maxitnqmp 
C     (MAXimum of ITerations with No Quadratic Model Progress) 
C     consecutive iterations. In this context, ''no progress'' means 
C     that the progress is smaller than epsnqmp (EPSilon to measure the 
C     No Quadratic Model Progress) times the best progress obtained 
C     during the previous iterations.

      epsnqmp   =  1.0d-04
      maxitnqmp =        5

C     CONVERT A TO DENSE FORMAT

      do j = 1,n
          do i = 1,ml
              A(i,j) = 0.0d0
          end do
      end do

      do i = 1,annz
          A(alin(i),acol(i)) = A(alin(i),acol(i)) + aval(i)
      end do

      do i = 1,ml
          wd21(i) = b(i)
      end do

C     COPY AUGMENTED LAGRANGIAN ARGUMENTS INTO A COMMON BLOCK

      m = mp
      do i = 1,m
          lambda(i) = lambdap(i)
          rho(i)    = rhop(i)
          equatn(i) = equatnp(i)
          linear(i) = linearp(i)
      end do

C     FINALLY, CALL GENLIN

      call genlin(n,x,l,u,ml,meq,mlmax,A,wd21,epsgpen,epsgpsn,epsgzsn,
     +maxitnfp,epsnfp,maxitngp,fmin,maxit,maxfc,delta0,cgmia,cgmib,
     +cgscre,cggpnf,cgepsi,cgepsf,epsnqmp,maxitnqmp,etaint,etaext,
     +mininterp,maxextrap,intype,gtype,hptype,trtype,scale,scmaxcnd,
     +scmaxal,iprint,ncomp,wd1,wd2,wd3,wd4,wi1,wi2,wd5,wd6,wd7,wd8,wd9,
     +wi3,wi4,nmax,Abar,ldh,H,wd10,wi5,wd11,wd12,wd13,wd14,wd15,wd16,
     +wd17,ldwd17,wd18,wd19,wd20,ktest,eta,delmin,lspgmi,lspgma,theta,
     +gamma,beta,sigma1,sigma2,rhoms,mseps,msmaxit,phieps,sigms,trgamma,
     +delini,trdelmin,trdelmax,promaxit,projeps,ceps,macheps,bignum,f,g,
     +gpeucn2,gpsupn,nit,nef,neg,neh,ncg,nch,spgnit,spgnef,tnnit,tnnef,
     +tnintcnt,tnexgcnt,tnexbcnt,tnstpcnt,tnintfe,tnexgfe,tnexbfe,
     +spginit,spginef,trnit,trnef,trexgcnt,trexbcnt,trexgfe,trexbfe,
     +inform)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalfb(n,x,f,flag)

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      integer m

C     COMMON ARRAYS
      logical equatn(mmax),linear(mmax)
      double precision lambda(mmax),rho(mmax)

C     COMMON BLOCK
      common /aldata/ lambda,rho,equatn,linear,m

      call evalun(n,x,m,lambda,rho,equatn,linear,f,flag)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalgb(n,x,g,flag)

C     SCALAR ARGUMENTS
      integer flag,n

C     ARRAY ARGUMENTS
      double precision g(n),x(n)

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      integer m

C     COMMON ARRAYS
      logical equatn(mmax),linear(mmax)
      double precision lambda(mmax),rho(mmax)

C     COMMON BLOCK
      common /aldata/ lambda,rho,equatn,linear,m

      call evalnun(n,x,m,lambda,rho,equatn,linear,g,flag)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalhb(n,x,ldh,H,flag)

C     SCALAR ARGUMENTS
      integer flag,ldh,n

C     ARRAY ARGUMENTS
      double precision H(ldh,n),x(n)

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      integer m

C     COMMON ARRAYS
      logical equatn(mmax),linear(mmax)
      double precision lambda(mmax),rho(mmax)

C     COMMON BLOCK
      common /aldata/ lambda,rho,equatn,linear,m

C     LOCAL SCALARS
      integer i,j

      call evalhun(n,x,m,lambda,rho,equatn,linear,ldh,H,flag)

      do j = 1,n
          do i = 1,j - 1
              H(i,j) = H(j,i)
          end do
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalhdb(nind,ind,n,x,d,hd,flag)

C     SCALAR ARGUMENTS
      integer flag,n,nind
     
C     ARRAY ARGUMENTS
      integer ind(nind)
      double precision d(n),hd(n),x(n)

      flag = - 1
      
      end