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 * <g,d>
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, <s,s> and <s,y>

      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 * <g,d>);
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 <g,d>.

      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,'<g,d> = ',
     *             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

