      program alspgma

C     ALSPG Test driver. This master file uses ALSPG to solve a variant 
C     of the location problem described in the Test Problems section of
C
C     E. G. Birgin, J. M. Martinez and M. Raydan, "Algorithm 813: SPG - 
C     software for convex-constrained optimization", ACM Transactions on 
C     Mathematical Software 27, pp. 340-349, 2001.
C 
C     It generates the problems, calls the optimizer to solve them and 
C     writes the reports (tables).
C
C     This variant of the location problem will be part of the Numerical
C     Experiments of a working paper by R.Andreani, E.G.Birgin, 
C     J.M.Martinez and L.Schuverdt.
C
C     This version 07 DEC 2004 by R.Andreani, E.G.Birgin, J.M.Martinez 
C     and L.Schuverdt.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax    =                630000 )
      integer npmax
      parameter ( npmax    =               1000000 )
      integer nvsmax
      parameter ( nvsmax   =            13 * npmax )
      integer nmax
      parameter ( nmax     = 2 * ( npmax + ncmax ) )
      integer mmax
      parameter ( mmax     =         npmax + ncmax )

C     COMMON SCALARS
      integer nc,np,totnvs
      double precision xdisp,xscal,ydisp,yscal

C     COMMON ARRAYS
      logical equatn(mmax)
      integer nvs(npmax)
      double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax),
     +        pol(nvsmax),radii(ncmax),vert(2*nvsmax)

C     LOCAL SCALARS
      integer i,inform,n,nvmapp,nvmipp,nx,ny,m,pnum,maxoutit,maxtotit,
     +        maxtotfc,outiter,totfcnt,totgcnt,totiter,iprint,ncomp
      double precision procit,propol,rmax,rmin,secmar,step,xstep,ystep,
     +        f,nalpsupn,epsfc,epsop,cnorm
      real time

C     LOCAL ARRAYS
      double precision c(mmax),lambda(mmax),nal(nmax),rho(mmax),x(nmax),
     +        wd(5*nmax)
      real dum(2)

C     EXTERNAL SUBROUTINES
      external dtime,genpro,alspg

C     COMMON BLOCKS
      common /circl/ ccent,radii,nc
      common /poly1/ nvs,np,totnvs
      common /poly2/ vert
      common /poly3/ edges
      common /poly4/ pcent
      common /poly5/ pol
      common /ellip/ xdisp,xscal,ydisp,yscal
      common /pdata/ equatn

C     DEFINE THE PROBLEM

C     GRID STEP (THIS IS JUST FOR SCALING THE PROBLEM)

      step   = 1.0d0

C     SECURITY MARGIN (BETWEEN 0 AND 1) TO AVOID CITIES OVERLAPPING

      secmar = 0.1d0

C     DIFFERENT STEPS MAY BE USED IN THE GRID

      xstep  =  step
      ystep  =  step

C     MINIMUM AND MAXIMUM RATIOS FOR THE CIRCLE- AND POLYGON-CITIES 
C     GENERATION

      rmin   =   0.5d0            * ( step / 2.0d0 )
      rmax   = ( 1.0d0 - secmar ) * ( step / 2.0d0 )

C     READ PROBLEM DATA

C     THIS PROBLEM DATA DEFINES DIFFERENT INSTANCIES OF THE PROBLEM AND 
C     IS RELATED TO THE NUMBER OF VARIABLES AND CONSTRAINTS WHICH 
C     APPROXIMATELY ARE:
C
C     - THE NUMBER OF POLYGON-CITIES IS APPROXIMATELY:
C       NPC = NX * NY * PROCIT * PROPOL
C
C     - THE NUMBER OF CIRCLE-CITIES IS APPROXIMATELY:
C       NCC = MX 8 NY * PROCIT * ( 1 - PROPOL )
C
C     - THE TOTAL NUMBER IF CITIES IS APPROXIMATELY:
C       NC = NPC + NCC
C
C     - NUMBER OF VARIABLES: 
C       2 * ( NC + 1 )
C
C     - NUMBER OF LINEAR CONSTRAINTS:
C       NPC * (NVMAPP - NVMIPP ) / 2 + 4
C
C     - NUMBER OF NONLINEAR CONSTRAINTS:
C       NCC + NC + 1 
C
C     WHERE:
C
C     - PNUM IS A NUMBER THAT IDENTIFIES THE PROBLEM INSTANCE,
C     - NX IS THE NUMBER OF POINTS IN THE GRID ABSCISSA,
C     - NY IS THE NUMBER OF POINTS IN THE GRID ORDINATE,
C     - PROCIT IS THE PROBABILITY OF HAVING A CITY AT A GRID POINT,
C     - PROPOL IS THE PROBABILITY OF A CITY TO BE A POLYGON (THE CITIES
C       WHICH ARE NOT POLYGONS ARE CIRCLES),
C     - NVMIPP IS THE MINIMUM NUNBER OF VERTICES OF A POLYGON,
C     - NVMAPP IS THE MAXIMUM NUMBER OF VERTICES OF A POLYGON.

      write(*,fmt=9000)

      read (*,fmt=*) pnum,nx,ny,procit,propol,nvmipp,nvmapp

      if ( procit .lt. 0.0d0 .or. procit .gt. 1.0d0 ) then
          write (*,fmt=*) 'PROCIT MUST BE BETWEEN ZERO AND ONE'
          stop
      end if

      if (propol .lt. 0.0d0 .or. propol .gt. 1.0d0 ) then
          write (*,fmt=*) 'PROPOL MUST BE BETWEEN ZERO AND ONE'
          stop
      end if

      if ( nvmipp .lt. 3 ) then
          write (*,fmt=*) 'NVMIPP MUST BE GREATER THAN OR EQUAL TO 3'
          stop
      end if

      if ( nvmapp .lt. nvmipp ) then
          write (*,fmt=*) 'NVMAPP MUST BE GREATER THAN OR EQUAL TO'
          write (*,fmt=*) 'NVMIPP. INCREASE NVMAPP OR REDUCE NVMIPP.'
          stop
      end if

C     GENERATE THE PROBLEM

      call genpro(nx,ny,xstep,ystep,procit,propol,nvmipp,nvmapp,rmin,
     +rmax,inform)

      if ( inform .ne. 0 ) then
          write (*,fmt=*) 'THERE WAS AN ERROR IN THE PROBLEM GENERATION'
          write (*,fmt=*) '(PROBABLY DUE TO MEMORY SPACE AVAILABILITY).'
          stop
      end if

      write (*,fmt=*) 'NUMBER OF CIRCLES: ',nc
      write (*,fmt=*) 'NUMBER OF POLYGONS: ',np
      write (*,fmt=*) 'NUMBER OF VERTICES (AND EDGES): ',totnvs

C     call drawsol(pnum,n,x,nx,ny,xstep,ystep,rmax) 
C     stop

C     SET NUMBER OF VARIABLES AND CONSTRAINTS (THAT WILL BE PENALIZED)

      n = 2 * ( np + nc )
      m = np + nc

      do i = 1,m
          equatn(i) = .false.
      end do

C     DEFINE INITIAL POINT

C     do i = 1,n
C         x(i) = 0.0d0
C     end do

      do i = 1,np
          x(2*i-1) = pcent(2*i-1)
          x(2*i  ) = pcent(2*i  )
      end do

      do i = 1,nc
          x(2*(np+i)-1) = ccent(2*i-1)
          x(2*(np+i)  ) = ccent(2*i  )
      end do

C     SET UP THE ARGUMENTS OF THE OPTIMIZATION ALGORITHM

      epsfc    =   1.0d-04
      epsop    =   1.0d-04

      maxoutit =        50
      maxtotfc =   5000000
      maxtotit =   1000000

      iprint   =         1
      ncomp    =         5

      do i = 1,m
          lambda(i) =  0.0d0
      end do

      do i = 1,m
          rho(i)    = 10.0d0
      end do

C     CALL THE OPTIMIZER

      call dtime(dum)

      call easyalspg(n,x,m,lambda,rho,epsfc,epsop,maxoutit,maxtotit,
     +maxtotfc,iprint,ncomp,f,c,cnorm,nal,nalpsupn,outiter,totiter,
     +totfcnt,totgcnt,inform,wd)
      
      call dtime(dum)

      time = dum(1) + dum(2)

C     WRITE STATISTICS

      write (*,fmt=9010) f,cnorm,nalpsupn,inform
      write (*,fmt=9020) outiter,totiter,totfcnt,totgcnt
      write (*,fmt=9030) time

      open(30, file='alspg-tabline.out')

      write(30,3000) pnum,nc,np,totnvs,n,nc+np,nc+totnvs,outiter,
     +               totiter,totfcnt,totgcnt,time,f,cnorm,nalpsupn,
     +               inform

      close(30)

C     WRITE SOLUTION

C     call drawsol(pnum,n,x,nx,ny,xstep,ystep,rmax) 

      stop

C     NON-EXECUTABLE STATEMENTS

 3000 format(I2,1X,I8,1X,I8,1X,I8,1X,I8,1X,I8,1X,I8,1X,I3,1X,I7,1X,I7,
     +       1X,I7,1X,F12.2,1X,1P,D10.4,1X,1P,D7.1,1X,1P,D7.1,1X,I3)
 9000 format (//'---------------------------------------'//)
 9010 format (/' F= ',8X,1P,D17.10,/' CNORM= ',4X,1P,D16.10,
     +        /' PGSUPNORM= ',1X,1P,D16.10,/' INFORM= ',6X,I1)
 9020 format (/' OUTITER= ',1X,I10,/' TOTITER= ',1X,I10,
     +        /' TOTFCNT= ',1X,I10,/' TOTGCNT= ',1X,I10)
 9030 format (/' TIME= ',F12.2,' SECONDS')

      end

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

      subroutine evalfc(n,x,f,m,c,inform)

C     This subroutine computes the objective function and 
C     the constraints.
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
C              This output parameter tells what happened in this 
C              subroutine, according to the following conventions:
C
C              0 = objective function and constraints were successfuly 
C                  evaluated,
C
C            < 0 = some error occured during the evaluation.


C     COMMON SCALARS
      double precision xdisp,xscal,ydisp,yscal

C     SCALAR ARGUMENTS
      double precision f
      integer n,m,inform

C     ARRAY ARGUMENTS
      double precision x(n),c(m)

C     LOCAL SCALARS
      integer i,ndist
      double precision diff1,diff2,dist

C     LOCAL ARRAYS
      double precision p(2)

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

C     COMMON BLOCKS
      common /ellip/ xdisp,xscal,ydisp,yscal

      inform = 0   

C     OBJECTIVE FUNCTION

      f = 0.0d0

      ndist = n / 2 - 1

      do i = 1,ndist
          diff1 = x(1) - x(2 * i + 1)
          diff2 = x(2) - x(2 * i + 2)
          dist = sqrt( diff1 ** 2 + diff2 ** 2)
          if ( dist .le. 1.0d-4 ) then
              write (*,fmt=*)
     +          'ERROR IN PROBLEM DEFINITION (DIST TOO SMALL)'
              inform = 1
          end if
          f = f + dist
      end do

      f = f / ndist

C     CONSTRAINTS

      do i = 1,m

         p(1) = ( x(2 * i - 1) - xdisp ) / xscal
         p(2) = ( x(2 * i)     - ydisp ) / yscal

         c(i) = p(1) ** 2 + p(2) ** 2 - 1.0d0

         if ( i .ne. 1 ) then
             c(i) = - c(i)
         end if

      end do

      end

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

      subroutine evalgjc(n,x,m,jcvar,jcfun,jcval,jcnnz,inform)

C     This subroutine computes the gradient of the objective function 
C     and the Jacobian of the constraints.
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     jcnnz    integer
C              number of non-null Jacobian elements,
C
C     jcvar,jcfun integers jcvar(jcnnz),jcfun(jcnnz)
C              see below,
C
C     jcval    double precision jcval(jcnnz)
C              the value of the non-zero Jacobian position [jcfun(i),
C              jcvar(i)] is jcval(i).
C
C              NOTE: jcfun(i) equal zero means objective function.
C
C     inform   integer
C
C              This output parameter tells what happened in this 
C              subroutine, according to the following conventions:
C
C              0 = gradient of the objective function and Jacobian of 
C                  the constraints were successfuly evaluated,
C
C            < 0 = some error occured during the evaluation.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax    =                630000 )
      integer npmax
      parameter ( npmax    =               1000000 )
      integer nmax
      parameter ( nmax     = 2 * ( npmax + ncmax ) )
      integer mmax
      parameter ( mmax     =         npmax + ncmax )
      integer jcnnzmax
      parameter ( jcnnzmax =       nmax + 2 * mmax )

C     COMMON SCALARS
      double precision xdisp,xscal,ydisp,yscal

C     SCALAR ARGUMENTS
      integer n,m,inform,jcnnz

C     ARRAY ARGUMENTS
      integer jcvar(jcnnzmax),jcfun(jcnnzmax)
      double precision x(n),jcval(jcnnzmax)

C     LOCAL SCALARS
      integer i,ndist
      double precision dcdx,dcdy,diff1,diff2,dist

C     LOCAL ARRAYS
      double precision g(nmax),p(2)

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

C     COMMON BLOCKS
      common /ellip/ xdisp,xscal,ydisp,yscal

      inform = 0

C     GRADIENT OF THE OBJECTIVE FUNCTION

      g(1) = 0.0d0
      g(2) = 0.0d0
      ndist = n / 2 - 1
      do i = 1,ndist
          diff1 = x(1) - x(2 * i + 1)
          diff2 = x(2) - x(2 * i + 2)
          dist = sqrt( diff1 ** 2 + diff2 ** 2)
          g(2*i+1) = - ( diff1 / dist ) / ndist
          g(2*i+2) = - ( diff2 / dist ) / ndist
          g(1) = g(1) - g(2 * i + 1)
          g(2) = g(2) - g(2 * i + 2)
      end do

      do i = 1,n
          jcfun(i) = 0
          jcvar(i) = i
          jcval(i) = g(i)
      end do

      jcnnz = n

C     CONSTRAINTS

      do i = 1,m

         p(1) = ( x(2 * i - 1) - xdisp ) / xscal
         p(2) = ( x(2 * i)     - ydisp ) / yscal

         dcdx = 2.0d0 * p(1) / xscal
         dcdy = 2.0d0 * p(2) / yscal

         if ( i .ne. 1 ) then
             dcdx = - dcdx
             dcdy = - dcdy
         end if

         jcnnz = jcnnz + 1

         jcfun(jcnnz) = i
         jcvar(jcnnz) = 2 * i - 1
         jcval(jcnnz) = dcdx

         jcnnz = jcnnz + 1

         jcfun(jcnnz) = i
         jcvar(jcnnz) = 2 * i
         jcval(jcnnz) = dcdy

      end do

      end

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

      subroutine evalc(n,x,m,c,inform)

C     This subroutine computes the constraints.
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     c        double precision c(m),
C              constraints at x,
C
C     inform   integer
C
C              This output parameter tells what happened in this 
C              subroutine, according to the following conventions:
C
C              0 = objective function and constraints were successfuly 
C                  evaluated,
C
C            < 0 = some error occured during the evaluation.

C     COMMON SCALARS
      double precision xdisp,xscal,ydisp,yscal

C     SCALAR ARGUMENTS
      integer n,m,inform

C     ARRAY ARGUMENTS
      double precision x(n),c(m)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS
      double precision p(2)

C     COMMON BLOCKS
      common /ellip/ xdisp,xscal,ydisp,yscal

      inform = 0

C     CONSTRAINTS

      do i = 1,m

         p(1) = ( x(2 * i - 1) - xdisp ) / xscal
         p(2) = ( x(2 * i)     - ydisp ) / yscal

         c(i) = p(1) ** 2 + p(2) ** 2 - 1.0d0

         if ( i .ne. 1 ) then
             c(i) = - c(i)
         end if

      end do

      end

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

      subroutine proj(n,x,inform)

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     On Entry:
C
C     n        integer,
C              size of the problem,
C
C     x        double precision x(n),
C              point that will be projected.
C
C     On Return
C
C     x        double precision x(n),
C              projected point,
C
C     inform   integer,
C
C              This output parameter tells what happened in this 
C              subroutine, according to the following conventions:
C
C              0 = the projection was successfully done,
C
C              1 = some error occurs in the projection.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax  =                630000 )
      integer npmax
      parameter ( npmax  =               1000000 )
      integer nvsmax
      parameter ( nvsmax =            13 * npmax )

C     SCALAR ARGUMENTS
      integer n,inform

C     ARRAY ARGUMENTS
      double precision x(n)

C     COMMON SCALARS
      integer nc,np,totnvs

C     COMMON ARRAYS 
      integer nvs(npmax)
      double precision center(2*ncmax),edges(3*nvsmax),ratios(ncmax),
     +        vert(2*nvsmax)

C     LOCAL SCALARS
      double precision xproj,yproj
      integer base,i

C     EXTERNAL SUBROUTINES
      external projc,projp

C     COMMON BLOCKS
      common /circl/center,ratios,nc
      common /poly1/ nvs,np,totnvs
      common /poly2/ vert
      common /poly3/ edges

      inform = 0

C     PROJECT ONTO THE POLYGONS

      base = 0

      do i = 1,np

        ! PROJECT z_i (i=1,np) ONTO POLYGON P_i

          call projp(i,base,x(2*i-1),x(2*i),xproj,yproj,inform)

          x(2*i-1) = xproj
          x(2*i)   = yproj

          base = base + nvs(i)

      end do

C     PROJECT ONTO THE CIRCLES

      do i = 1,nc

        ! PROJECT z_{np+i} (i=1,nc) ONTO CIRCLE C_i 

          call projc(i,x(2*(np+i)-1),x(2*(np+i)),xproj,yproj,inform)

          x(2 * ( np + i ) - 1) = xproj
          x(2 * ( np + i )    ) = yproj

      end do

      end

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

      subroutine projp(p,base,x,y,xproj,yproj,inform)

C     This subroutine complements proj subroutine projecting a point 
C     (x,y) in R^2 onto its corresponding polygon P_p.
C
C     On Entry:
C
C     p        integer,
C              index of the polygon onto which (x,y) will be projected,
C
C     base     integer,
C              index which indicates the positions where the edges of 
C              the polygon  P_p are stored inside the "edges" vector 
C              which is part of the polygons description,
C
C     x        double precision,
C              first component of the point to be projected,
C
C     y        double precision,
C              second component of the point to be projected,
C
C     On Return
C
C     xproj    double precision,
C              first component of the projected point,
C
C     yproj    double precision,
C              second component of the projected point,
C
C     inform   integer,
C
C              This output parameter tells what happened in this 
C              subroutine, according to the following conventions:
C
C              0 = the projection was successfully done,
C
C              1 = some error occurs in the projection.

C     PARAMETERS
      integer npmax
      parameter ( npmax  =               1000000 )
      integer nvsmax
      parameter ( nvsmax =            13 * npmax )

C     SCALAR ARGUMENTS
      integer base,inform,p
      double precision x,xproj,y,yproj

C     COMMON SCALARS
      integer np,totnvs

C     COMMON ARRAYS
      integer nvs(npmax)
      double precision edges(3*nvsmax),vert(2*nvsmax)

C     LOCAL SCALARS
      integer i
      logical fact
      double precision a,b,c,dist,mindis,nearvx,nearvy,px,py,v1x,v1y,
     +                 v2x,v2y,vx,vy

C     LOCAL ARRAYS
      logical satisf(nvsmax)

C     INTRINSIC FUNCTIONS
      intrinsic max,min,sqrt

C     COMMON BLOCKS
      common /poly1/ nvs,np,totnvs
      common /poly2/ vert
      common /poly3/ edges

      inform = 0

C     TEST ALL EDGES OF THE POLYGON FOR SATISFIABILITY

      fact = .true.
      do i = base + 1,base + nvs(p)
          a = edges(3*i-2)
          b = edges(3*i-1)
          c = edges(3*i)
          if (a*x+b*y+c.le.0.0d0) then
              satisf(i-base) = .true.

          else
              fact = .false.
              satisf(i-base) = .false.
          end if
      end do

      if (fact) then
          xproj = x
          yproj = y
          return

      end if

C     SEARCH FOR THE CLOSEST VERTEX

      mindis = 1.0d+99
      do i = base + 1,base + nvs(p)
          vx = vert(2*i-1)
          vy = vert(2*i)
          dist = sqrt((vx-x)**2+ (vy-y)**2)
          if (dist.le.mindis) then
              mindis = dist
              nearvx = vx
              nearvy = vy
          end if

      end do

      xproj = nearvx
      yproj = nearvy

C     PROJECT ONTO THE VIOLATED CONSTRAINTS

      do i = base + 1,base + nvs(p)
          if (.not.satisf(i-base)) then
              a = edges(3*i-2)
              b = edges(3*i-1)
              c = edges(3*i)
              if (a.ne.0.0d0) then
                  py = (y-b* (c/a+x)/a)/ (b**2/a**2+1)
                  px = - (c+b*py)/a

              else
                  if (b.ne.0.0d0) then
                      px = (x-a* (c/b+y)/b)/ (a**2/b**2+1)
                      py = - (c+a*px)/b

                  else
                      write (*,fmt=*)
     +                  'ERROR IN PROBLEM DEFINITION (a=b=0 for a ',
     +                  'constraint of the type a x + b y + c = 0)'
                      inform = 1
                  end if

              end if

              v1x = vert(2*i-1)
              v1y = vert(2*i)
              if (i.ne.base+nvs(p)) then
                  v2x = vert(2*i+1)
                  v2y = vert(2*i+2)

              else
                  v2x = vert(2*base+1)
                  v2y = vert(2*base+2)
              end if

              if (min(v1x,v2x).le.px .and. px.le.max(v1x,v2x) .and.
     +            min(v1y,v2y).le.py .and. py.le.max(v1y,v2y)) then
                  dist = sqrt((px-x)**2+ (py-y)**2)
                  if (dist.le.mindis) then
                      mindis = dist
                      xproj = px
                      yproj = py
                  end if

              end if

          end if

      end do

      end

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

      subroutine projc(c,x,y,xproj,yproj,inform)

C     This subroutine complements proj subroutine projecting a point 
C     (x,y) in R^2 onto its corresponding circle C_c.
C
C     On Entry:
C
C     c        integer,
C              index of the circle onto which (x,y) will be projected,
C
C     x        double precision,
C              first component of the point to be projected,
C
C     y        double precision,
C              second component of the point to be projected,
C
C     On Return
C
C     xproj    double precision,
C              first component of the projected point,
C
C     yproj    double precision,
C              second component of the projected point,
C
C     inform   integer,
C
C              This output parameter tells what happened in this 
C              subroutine, according to the following conventions:
C
C              0 = the projection was successfully done,
C
C              1 = some error occurs in the projection.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax  =                630000 )

C     SCALAR ARGUMENTS
      integer c,inform
      double precision x,xproj,y,yproj

C     COMMON SCALARS
      integer nc

C     COMMON ARRAYS
      double precision center(2*ncmax),ratios(ncmax)

C     LOCAL SCALARS
      double precision cx,cy,dist,r

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

C     COMMON BLOCKS
      common /circl/center,ratios,nc

      inform = 0

      cx = center(2 * c - 1)
      cy = center(2 * c)
      r  = ratios(c)

      dist = sqrt( ( x - cx ) ** 2 + ( y - cy ) ** 2 )

      if ( dist .le. r ) then

          xproj = x
          yproj = y

      else

          xproj = cx + ( x - cx ) * r / dist
          yproj = cy + ( y - cy ) * r / dist

      end if

      end

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

      subroutine drawsol(pnum,n,x,nx,ny,xstep,ystep,rmax)

C     This subroutine generate a metapost file with the
C     graphical representation of the problem.

C     On Entry:
C
C     pnum  integer,
C           problem identifier,
C
C     x     double precision x(n)
C           current approximation to the solution.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax  =                630000 )
      integer npmax
      parameter ( npmax  =               1000000 )
      integer nvsmax
      parameter ( nvsmax =            13 * npmax )

C     COMMON SCALARS
      integer nc,np,totnvs
      double precision xdisp,xscal,ydisp,yscal

C     COMMON ARRAYS
      integer nvs(npmax)
      double precision center(2*ncmax),edges(3*nvsmax),ratios(ncmax),
     +        vert(2*nvsmax)

C     SCALAR ARGUMENTS
      integer n,nx,ny,pnum
      double precision rmax,xstep,ystep

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      integer base,i,j
      double precision scale,vx,vy

C     INTRINSIC FUNCTIONS
      intrinsic atan,cos,sin,sqrt

C     COMMON BLOCKS
      common /circl/ center,ratios,nc
      common /poly1/ nvs,np,totnvs
      common /poly2/ vert
      common /poly3/ edges
      common /ellip/ xdisp,xscal,ydisp,yscal

      scale = min( 21.6d0 / ( ( nx - 1 ) * xstep + 2.0d0 * rmax ), 
     +             27.9d0 / ( ( ny - 1 ) * ystep + 2.0d0 * rmax ) )

C     DRAW THE PROBLEM

      open(unit=10,file='sol.mp')

      write(10,10) pnum,scale

C     ELLIPSE

      write(10,50) 2.0d0 * xscal, 2.0d0 * yscal,xdisp,ydisp

C     CIRCLES

      do i = 1,nc
          write(10,50) 2.0d0 * ratios(i),2.0d0 * ratios(i),
     +                 center(2 * i - 1),center(2 * i)
      end do

C     POLYGONS

      base = 0

      do i = 1,np

          if ( nvs(i) .ne. 0 ) then

              write(10,20)

              do j = base + 1,base + nvs(i)
                  vx = vert(2 * j - 1)
                  vy = vert(2 * j)
                  write(10,60) vx,vy
              end do

              write(10,40)

              base = base + nvs(i)

          end if

      end do

C     LINES

      do i = 1,np + nc
          write(10,70) x(1),x(2),x(2 * i - 1),x(2 * i)
      end do

C     DOTS WITH LABELS

C     write(10,80)

C     do i = 1,np + nc
C         write(10,90) i,x(2 * i - 1),x(2 * i)
C     end do

C     DOTS WITHOUT LABELS

      write(10,80)

      do i = 1,np + nc
          write(10,95) x(2 * i - 1),x(2 * i)
      end do

C     END

      write(10,100)

      close(10)

C     NON-EXECUTABLE STATEMENTS

  10  format('beginfig(',I2,');'/,'u = ',f20.10,' cm;') 
  20  format('draw')
  30  format(5X,'(',f20.10,'u,',f20.10,'u)..')
  40  format(5X,'cycle;')
  50  format('draw fullcircle',
     +       /,5X,'xscaled   ',f20.10,'u',
     +       /,5X,'yscaled   ',f20.10,'u',
     +       /,5X,'shifted (',f20.10,'u,',f20.10,'u);')
  60  format(5X,'(',f20.10,'u,',f20.10,'u)--')
  70  format('draw (',f20.10,'u,',f20.10,'u)--',
     +       /,5X,'(',f20.10,'u,',f20.10,'u) dashed evenly;')
  80  format('pickup pencircle scaled 4pt;')
  90  format('dotlabel.top (btex $',i4,'$ etex, (',f20.10,'u,',
     +        f20.10,'u) );')
  95  format('draw (',f20.10,'u,',f20.10,'u);')
 100  format('endfig;'/,'end;') 

      end

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

C     Last update of EASYALSPG: December 7th, 2004.

      subroutine easyalspg(n,x,m,lambda,rho,epsfc,epsop,maxoutit,
     +maxtotit,maxtotfc,iprint,ncomp,f,c,cnorm,nal,nalpsupn,outiter,
     +totiter,totfcnt,totgcnt,inform,wd)

C     SCALAR ARGUMENTS
      integer inform,iprint,m,maxoutit,maxtotfc,maxtotit,n,ncomp,
     +        outiter,totfcnt,totiter,totgcnt
      double precision cnorm,epsfc,epsop,f,nalpsupn

C     ARRAY ARGUMENTS
      double precision c(m),lambda(m),nal(n),rho(m),wd(5*n),x(n)

C     This subroutine aims to simplify the use of ALSPG. For this 
C     purpose it gives values to most of the ALSPG and SPG arguments and 
C     leaves to the user those arguments which he/she may would like to 
C     set by him/herself.
C
C     The arguments of EASYALSPG are the input and output arguments of 
C     ALSPG and SPG that are supposed to be useful for a common user. 
C     The input arguments are mostly related to basic problem 
C     information, like dimension and bounds, and the initial point. 
C     There are also input arguments related to simple stopping criteria 
C     (like norm of the projected gradient, and maximum number of 
C     iterations and functional evaluations). There are also two input 
C     arguments related to control the amount of information written 
C     into the screen. The output arguments are related to information 
C     of the solution and some few performance measurements. Basically, 
C     on return, EASYALSPG gives to the user the solution, the objective 
C     functional value and its gradient at the solution, Euclidian and 
C     sup-norm of the projected gradient at the solution, and the number 
C     of iterations, functional and gradient evaluations, used to reach 
C     the solution, and, finally, a flag that indicates the stopping 
C     criterion that was satisfied.
C
C     All the other arguments of ALSPG and SPG are setted with its 
C     default values by EASYALSPG. EASYALSPG divides the arguments of 
C     ALSPG and SPG in two sets. Those that are related to the behaviour
C     of ALSPG and SPG are declared as Fortran parameters (constants). 
C     The other arguments of ALSPG and SPG, most of them related to 
C     alternative stopping criteria, and that may depend of, for example, 
C     maxit, are declared as local variables of EASYALSPG.
C
C     SPG arguments that are defined as Fortran parameters in this 
C     subroutine are SPG arguments that should not be modified by a 
C     common user. These are arguments that modify the behaviour of 
C     SPG and whos values were selected because they are classical 
C     values in some cases or because some numerical experiments seemed 
C     to indicate that they are the best choices.
C
C     SPG arguments that are declared as local variables in this 
C     subroutine are SPG arguments that may be modified if, with 
C     their suggested values, SPG does not give the desired result. 
C     Most of them are related to disabled stopping criteria that 
C     may be useful in bad-scaled problems or problems with not 
C     trustable derivatives.
C
C     We describe below the meaning of the arguments of the EASYSPG
C     subroutine. More detailed descriptions as well as the 
C     descriptions of all the other SPG arguments that are not arguments 
C     of EASYSPG are also described at the begining of the SPG 
C     subroutine.
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 bound-constrainted subproblems, m 
C              would be the number of constraints, lambda the Lagrange 
C              multipliers approximation and rho the penalty parameters
C
C     epsfc    double precision
C              ALSPG considers that a point satisfies the feasibility
C              conditions if the sup-norm of the constraints is smaller 
C              than or equal to epsfc
C
C     epsop    double precision
C              ALSPG considers that a point satisfies the optimality 
C              conditions if its Augmented Lagrangian projected gradient 
C              sup-norm is smaller than or equal to epsop
C
C     maxoutit integer
C              ALSPG stops declaring ''maximum number of outer iteration 
C              achieved'' if the number of outer iterations exceeds 
C              maxoutit
C
C     maxtotit integer
C              ALSPG stops declaring ''maximum number of inner iteration 
C              achieved'' if the number of inner iterations exceeds 
C              maxtotit
C
C     maxtotfc integer 
C              ALSPG stops declaring ''maximum number of functional 
C              evaluations achieved'' if the number of inner iterations 
C              exceeds maxtotfc
C
C     iprint   integer
C              indicates the degree of details of the output generated 
C              by ALSPG. Setting iprint to a negative value will ALSPG
C              to generate no output at all. iprint equalt to zero will 
C              make ALSPG to print information related to the first and 
C              the last iterations. With iprint equal to 1, ALSPG will
C              show basic information of every outer iterations.  
C              An iprint value greater than or equal to 2 will generate 
C              information of every inner SPG iteration. An iprint value 
C              greater than or equal to 3 will also show information 
C              related to the line search procedure.
C
C     ncomp    integer
C              Sometimes, vectors like the current point x, the gradient 
C              of the objective function g, or the search directions,
C              among other vector, are showed in the screen. In such 
C              cases, if the problem dimension is large, to show just a 
C              few elements of these vectors may be preferable. Argument 
C              ncomp can be used to indicate how many array elements must 
C              be displayed.
C
C     wd       double precision wd(5*n)
C              double precision working space
C
C     On return:
C
C     x        double precision x(n)
C              estimation of the solution
C
C     f        double precision
C              objective function value at the solution
C
C     c        double precision c(m)
C              constraints at the solution
C
C     cnorm    double precision
C              sup-norm of the constraints at the solution
C
C     nal      double precision nal(n)
C              continuous projected gradient of the augmented Lagrangian 
C              at the solution
C
C     nalpsupn double precision
C              sup-norm of the continuous projected gradient of the 
C              augmented Lagrangian at the solution
C
C     totoutit integer
C              number of outer iterations used to reach the solution
C
C     totiter  integer
C              number of inner iterations used to reach the solution
C
C     totfcnt  integer
C              number of functional evaluations
C
C     totgcnt  integer
C              number of gradient evaluations
C
C     inform   integer

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     CONSTANTS FOR CLASSICAL LINE-SEARCH CONDITIONS

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

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

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

C     CONSTANTS FOR SPECIFIC PROCEDURES (NOT SO CLASSICAL)

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

C     In the line searches, when we need to interpolate and the result 
C     of the quadratic interpolation is rejected, the new step is 
C     computed as anew = aprev / nint. 

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

      integer mininterp
      parameter ( mininterp =         4 )

      double precision nint,lspgma,lspgmi
      parameter ( nint      =     2.0d0 )
      parameter ( lspgma    =   1.0d+10 )
      parameter ( lspgmi    =   1.0d-10 )

C     In the Augmented Lagrangian framework, it is assumed that each
C     Lagrange multiplier [lambda]_i belongs to [lammin,lammax].

      double precision lammax,lammin
      parameter ( lammax    =   1.0d+20 )
      parameter ( lammin    = - 1.0d+20 )

C     Also in the Augmented Lagrangian framework, every time the 
C     feasibility or the complementarity are not sufficiently improved
C     the penalty parameters must be increased. rhofrac is a fraction
C     between 0 and 1 used to determine whether there was enough
C     improvement in the feasibility and complementarity measurement
C     and rhomult is the constant (greather than 1) used to increase
C     the penalty parameters. 

      double precision rhofrac,rhomult
      parameter ( rhofrac   =     0.5d0 )
      parameter ( rhomult   =    10.0d0 )

C     Finally, there is a stopping criterion related to the non
C     improvement in feasibility. If, during a fixed number (maxitncp) 
C     of Augmented Lagrangian iterations the feasibility is not improved,
C     it may be the case that the problem is infeasible.

      integer maxitncp
      parameter ( maxitncp  =         9 )

C     HERE STARTS THE DESCRIPTION OF SOME SPG ARGUMENTS THAT ARE 
C     BEING SETTED INSIDE EASYSPG. THE FIRST SET OF ARGUMENTS ARE 
C     THOSE ARGUMENTS THAT WE WILL CALL ''CONSTANTS'' AND THAT, AS THEIR 
C     VALUES ALTER THE BEHAVIOUR OF SPG, SHOULD NOT BE MODIFIED BY A 
C     COMMON USER.

C     DIMENSIONS FOR SOME WORKING SPACES

C     In non-monotone line searches, given p, the last p objective 
C     functional values must be stored. For this reason we declare a 
C     vector with pmax double precision elements. So p must be less than
C     or equal to pmax.

      integer pmax
      parameter ( pmax =   100 )

C     Sometimes, is the problem is bad scaled, to request a small 
C     gradient norm at the solution may be inadequate. For this reason, 
C     a test to verify if this norm is not decreasing during maxitngp 
C     (MAXimum of ITerations with No Gradient Progress) consecutive 
C     iterations then we stop the method with a warning. As it is not 
C     expected a monotone decreasing of the gradient norm, again, the 
C     norm of the last maxitngp iterations must be saved. For this 
C     purpose, we declare a vector of tmax elements. So maxitngp must 
C     be less than or equal to tmax.
     
      integer tmax
      parameter ( tmax = 10000 )

C     ALSPG AND SPG INPUT ARGUMENTS THAT WILL BE SETTED BELOW

      integer gtype,maxfc,maxit,maxitnfp,maxitngp,p

      double precision epsnfp,epsopcon,epsopfin,epsopini,fmin

C     SPG WORKING VECTORS (WHICH DIMENSION IS NOT RELATED TO THE 
C     PROBLEM DIMENSION)

      double precision lastfv(pmax),lastgpns(tmax)

C     ARGUMENTS RELATED TO DERIVATIVES CALCULATIONS

C     gtype indicates in which way the gradient of the objective 
C     function will be computed. If the user have been implemented the 
C     user-supplied evalnal subroutine to compute the gradient of the 
C     objective function then gtype argument must be set to 0 (ZERO) and 
C     the user-supplied evalnal subroutine will be called by SPG any 
C     time the gradient would be required.
C
C     The prototype of the evalnal subroutine must be:
C
C         subroutine evalnal(n,x,m,lambda,rho,nal,flag)
C
C         SCALAR ARGUMENTS
C         integer n,m,flag
C
C         ARRAY ARGUMENTS
C         double precision x(n),lambda(m),rho(m),nal(n)
C
C         ''Here must be written the subroutine body that calculates the 
C         n-dimensional gradient vector of the objective function 
C         evaluated at x and saves it in nal. It also must set flag to 0 
C         (ZERO) if the gradient was successfully computed and to any 
C         other value if the gradient vector is not well defined at the 
C         required point x. If SPG is been used stand-alone to solve 
C         a unique convex-constrained problem then m, lambda and rho are 
C         dummy arguments. On the other hand, if SPG is been used in 
C         an Augmented Lagrangian framework then these arguments should 
C         be used for the number of constraints, the Lagrange 
C         multipliers approximation and the penalty parameters, 
C         respectively.''
C
C         end
C
C     If, on the other hand, the user is not able to provide evalnal 
C     subroutine, gtype argument must be set to 1 (ONE). In this case, 
C     every time SPG needs to compute the gradient of the objective 
C     function, an internal subroutine that approximates it by finite-
C     differences will be used (be aware that it maybe very time 
C     consuming). Moreover, note that the evalnal subroutine must still 
C     be present (with an empty body).

      gtype     =          0

C     ARGUMENTS RELATED TO STOPPING CRITERIA

C     For an explanation of maxitngp see above the explanation of tmax 
C     in ''DIMENSIONS FOR SOME WORKING SPACES''. Just note that the 
C     value of maxitngp must be less than or equal to tmax.

      maxitngp  =       tmax

C     maxit is the maximum allowed number of inner iterations to solve
C     a subproblem and maxfc is maximum allowed number of functional
C     evaluations to solve a subproblem

      maxit     =      20000
      maxfc     = 10 * maxit

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

      maxitnfp  =      maxit
      epsnfp    =      0.0d0

C     There is a stopping criterion that stops the method if a point 
C     with a functional value smaller than fmin is found. The idea 
C     behind this stopping criterion is to stop the method if the 
C     objective function is not bounded from bellow.

      fmin      =   - infrel

C     In a Grippo-Lampariello-Lucidi nonmontone line search, the 
C     Armijo-like condition is required to be satisfied using the
C     maximum functional value of the last p iterations instead of
C     the functional value of the previous iteration. Note that p
C     must be less than or equal to pmax (parameter declared above).

      p         =      10

C     The tolerance for solving the Augmented lagrangian subproblems
C     may be the same for all the subproblems or may be improved as
C     the method approximates to the solution. Numerical experiments
C     with Newtonian solvers for the subproblems show that to solve
C     each subproblem exactly is preferable. We have no experiments
C     when solving the subproblems using SPG. While we do not do the
C     experiments, we prefer to solve the subproblems exactly.

      epsopini =     epsop
      epsopfin =     epsop
      epsopcon =   1.0d-01

C     FINALLY, CALL ALSPG

      call alspg(n,x,m,lambda,rho,rhomult,rhofrac,epsfc,epsopini,
     +epsopfin,epsopcon,maxitncp,maxoutit,maxtotfc,maxtotit,gtype,p,
     +maxit,maxfc,maxitnfp,epsnfp,maxitngp,fmin,nint,mininterp,lastfv,
     +lastgpns,iprint,ncomp,f,c,cnorm,nal,nalpsupn,outiter,totiter,
     +totfcnt,totgcnt,inform,lammin,lammax,lspgma,lspgmi,gamma,sigma1,
     +sigma2,wd,sterel,steabs,epsrel,epsabs,infrel,infabs)

      end

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

      subroutine alspg(n,x,m,lambda,rho,rhomult,rhofrac,epsfc,epsopini,
     +epsopfin,epsopcon,maxitncp,maxoutit,maxtotfc,maxtotit,gtype,p,
     +maxit,maxfc,maxitnfp,epsnfp,maxitngp,fmin,nint,mininterp,lastfv,
     +lastgpns,iprint,ncomp,f,c,cnorm,nal,nalpsupn,outiter,totiter,
     +totfcnt,totgcnt,inform,lammin,lammax,lspgma,lspgmi,gamma,sigma1,
     +sigma2,wd,sterel,steabs,epsrel,epsabs,infrel,infabs)

C     SCALAR ARGUMENTS
      integer inform,iprint,gtype,m,maxfc,maxit,maxitncp,maxitnfp,
     +        maxitngp,maxoutit,maxtotfc,maxtotit,mininterp,n,ncomp,
     +        outiter,p,totfcnt,totgcnt,totiter
      double precision cnorm,epsabs,epsfc,epsnfp,epsopcon,epsopfin,
     +        epsopini,epsrel,f,fmin,gamma,infabs,infrel,lammax,lammin,
     +        lspgma,lspgmi,nalpsupn,nint,rhomult,rhofrac,sigma1,sigma2,
     +        steabs,sterel

C     ARRAY ARGUMENTS
      double precision c(m),lambda(m),lastfv(0:p-1),
     +        lastgpns(0:maxitngp-1),nal(n),rho(m),x(n),wd(5*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 a working paper by 
C     R. Andreani, E. G. Birgin, J. M. Martinez, and L. Schuverd.
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     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     lammin  double precision
C
C     lammax  double precision
C
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 values at the solution
C
C     cnorm   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]. cnorm 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     totfcnt integer 
C             total number of Conjugate Gradients iterations
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 proj subroutine.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax    =                630000 )
      integer npmax
      parameter ( npmax    =               1000000 )
      integer mmax
      parameter ( mmax     =         npmax + ncmax )

C     COMMON ARRAYS
      logical equatn(mmax)

C     LOCAL SCALARS
      integer fcnt,flag,gcnt,i,icrit,iter,mprint,nprint,maxsoutit,soutit
      double precision al,cnormprev,desired,epsgpen,epsgpsn,epsopk,
     +        epssub,nalpeucn2,nalpi

C     LOCAL ARRAYS
      double precision lambar(mmax)

C     COMMON BLOCKS
      common /pdata/ equatn

C     CONSTANTS

C     for the secondary loop
      maxsoutit = 10 * maxoutit

C     just for printing
      nprint   = min0(n,ncomp)
      mprint   = min0(m,ncomp)

C     INTIALIZATION

      icrit    = 0

      outiter  = 0
      soutit   = 0
      totiter  = 0
      totfcnt  = 0
      totgcnt  = 0

      cnormprev = 1.0d+99

      if ( m .ne. 0 ) then
          epsopk = epsopini
      else
          epsopk = epsopfin
      end if

      do i = 1,m
          lambar(i) = max( lammin, min( lammax, lambda(i) ) )
      end do

C     PROJECT INITIAL GUESS

      call proj(n,x,flag)

      if ( flag .ne. 0 ) then
        ! ERROR IN THE PROJECTION SUBROUTINE, STOP
          inform = - 92

          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,flag)

      if ( flag .ne. 0 ) then
        ! ERROR IN THE OBJECTIVE FUNCTION OR THE CONSTRAINTS 
        ! EVALUATION, STOP
          inform = - 90

          if ( iprint .ge. 1 ) then
              write(*, 2000) inform
              write(10,2000) inform
          end if

          go to 500
      end if

C     COMPUTE FEASIBILITY VIOLATION
 
      cnorm = 0.0d0
      do i = 1,m
          if ( equatn(i) ) then
              cnorm = max(cnorm, abs(c(i)))
          else
              cnorm = max(cnorm, abs(max(c(i), - lambar(i) / rho(i))))
          end if
      end do          

C     COMPUTE CONTINUOUS PROJECTED GRADIENT NORM
C     (using lambda instead of lambar, as is given by spg)

      if ( gtype .eq. 0 ) then
          call evalnal(n,x,m,lambda,rho,nal,inform)
      else ! if ( gtype .eq. 1 ) then
          call evalnaldiff(n,x,m,lambda,rho,nal,sterel,steabs,inform)
      end if
 
      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(*,  980) nprint,(x(i),i=1,nprint)
          write(*,  990) mprint,(lambda(i),i=1,mprint)
          write(*, 1020) cnorm,nalpsupn

          write(10, 900) n,m
          write(10, 980) nprint,(x(i),i=1,nprint)
          write(10, 990) mprint,(lambda(i),i=1,mprint)
          write(10,1020) cnorm,nalpsupn
      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) rho(1)
          write(*, 1020) cnorm,nalpsupn

          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) rho(1)
          write(10,1020) cnorm,nalpsupn
      end if

C     SAVING INTERMEDIATE DATA FOR CRASH REPORT

      open(20, file='alspg.output')

      write(20,3000) n,m,outiter,totiter,totfcnt,totgcnt,f,cnorm,
     +               nalpsupn

      close(20)

C     TEST STOPPING CRITERIA

      if ( cnorm .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 ( soutit .ge. maxsoutit ) then
        ! MAXIMUM NUMBER OF SECONDARY OUTER ITERATIONS EXCEEDED, STOP
          inform = 1

          if ( iprint .ge. 1 ) then
              write(*, 1075) inform,maxsoutit
              write(10,1075) inform,maxsoutit
          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 ( cnorm .gt. epsfc .and. cnorm .ge. cnormprev ) 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     SAVE PREVIOUS VIOLATION

      cnormprev = cnorm

C     CALL THE INNER-SOLVER

      epssub = epsopk

C     SECONDARY LOOP 

 200  continue

      epsgpsn  =     epssub
      epsgpen  =      0.0d0

      soutit = soutit + 1

      call spg(n,x,m,lambar,rho,p,epsgpen,epsgpsn,maxitnfp,epsnfp,
     +maxitngp,fmin,maxit,maxfc,nint,mininterp,gtype,iprint,ncomp,al,
     +nal,nalpeucn2,nalpsupn,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)

      write(*,*) ' SPG termination flag =  ',inform,nalpsupn

      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,flag)

      if ( flag .ne. 0 ) then
        ! ERROR IN THE OBJECTIVE FUNCTION OR THE CONSTRAINTS 
        ! EVALUATION, STOP
          inform = - 90

          if ( iprint .ge. 1 ) then
              write(*, 2000) inform
              write(10,2000) inform
          end if

          go to 500
      end if

C     COMPUTE FEASIBILITY VIOLATION
 
      cnorm = 0.0d0
      do i = 1,m
          if ( equatn(i) ) then
              cnorm = max(cnorm, abs(c(i)))
          else
              cnorm = max(cnorm, abs(max(c(i), - lambar(i) / rho(i))))
          end if
      end do          

C     COMPUTE CONTINUOUS PROJECTED GRADIENT NORM
C     (using lambda instead of lambar, as is given by spg)

      call evalnal(n,x,m,lambda,rho,nal,inform)

      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     TEST IF THE SUBPROBLEM WAS ABANDONED PREMATURELY

      if ( nalpsupn .gt. max( epsopfin, min( epsopk, cnorm ) ) .and.
     +     epssub .gt. epsopfin .and. soutit .lt. maxsoutit ) then
          epssub = max( epsopfin, epsopcon * epssub )
          go to 200
      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     UPDATE PENALTY PARAMETERS

      desired = max( rhofrac * cnormprev, epsfc )

      if ( cnorm .gt. desired ) then

          do i = 1,m
              rho(i) = rhomult * rho(i)
          end do

          if ( iprint .ge. 1 ) then
              write(*, 1030) desired,rhomult
              write(10,1030) desired,rhomult
          end if

      else

          if ( iprint .ge. 1 ) then
              write(*, 1050)
              write(10,1050)
          end if

      end if

C     COMPUTE OPTIMALITY REQUERIMENT FOR THE NEXT OUTER ITERATION

      epsopk = max( epsopfin, epsopk * epsopcon )

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( /,' Current Lagrange multipliers (first ',I7,
     +          ' components): ',
     +        /,6(1X,1PD11.4))
 1000 format( /,' Current penalty parameter: ',1PD11.4)
 1020 format( /,' Constraints norm              = ',1PD11.4,
     +        /,' AL projected gradient norm    = ',1PD11.4)
 1030 format( /,' The desired infeasibility improvement was not',
     +          ' achieved (',1PD11.4,').',
     +        /,' Penalty parameter will be increased multiplying',
     +          ' it by rhofrac = ',1PD11.4,'.')
 1050 format( /,' Desired infeasibility improvement was achieved.',
     +        /,' The 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 larger',
     +          ' complementarity',
     +        /,' violation smaller than ',1PD11.4,')')
 1070 format( /,' Flag of ALSPG = ',I2,
     +          ' (It was exceeded the maximum allowed number of outer',
     +        /,' iterations (maxoutit =',I7,').)')
 1075 format( /,' Flag of ALSPG = ',I2,
     +          ' (It was exceeded the maximum allowed number of outer',
     +        /,' secondary loop iterations (maxsoutit =',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 evalfc  subroutine', 
     +        /,' -91 : error in evalgjc subroutine', 
     +        /,' -92 : error in proj    subroutine',/) 
 3000 format(I8,1X,I8,1X,I3,1X,I7,1X,I7,1X,I7,1X,1P,D24.16,1X,1PD7.1,1X,
     +       1PD7.1)

      end

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

      subroutine evalal(n,x,m,lambda,rho,al,inform)

C     This subroutine computes the Augmented Lagrangian function.
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     On Return
C
C     al      double precision al,
C             value of the Augmented Lagrangian function,
C
C     inform  integer,
C
C             This output parameter tells what happened in this 
C             subroutine, according to the following conventions:
C 
C             0 = the Augmented Lagrangian function was succesfully 
C                 evaluated,
C
C             1 = some error occurs in the evaluation.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax    =                630000 )
      integer npmax
      parameter ( npmax    =               1000000 )
      integer mmax
      parameter ( mmax     =         npmax + ncmax )

C     COMMON ARRAYS
      logical equatn(mmax)

C     SCALAR ARGUMENTS
      double precision al
      integer n,m,inform

C     ARRAY ARGUMENTS
      double precision x(n),lambda(m),rho(m)

C     LOCAL SCALARS
      integer i,flag
      double precision f,p

C     LOCAL ARRAYS
      double precision c(mmax)

C     COMMON BLOCKS
      common /pdata/ equatn

C     INITIALIZE

      inform = 0

C     COMPUTE OBJECTIVE FUNTION AND CONSTRAINTS

      call evalfc(n,x,f,m,c,flag)
      if ( flag .ne. 0 ) then
          inform = - 90
          return
      end if

C     COMPUTES AL(x) = f(x) + \sum_i=1^m P(c(i),rho(i),lambda(i))

      al = f

      do i = 1,m
          call evalp(c(i),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 y,rho,lambda,p

      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,nal,inform)

C     This subroutine computes the gradient of the Augmented Lagrangian 
C     function.
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     On Return
C
C     nal   double precision nal(n),
C           gradient of the Augmented Lagrangian function,
C
C     inform integer,
C           termination parameter:
C           0 = the gradient was succesfully evaluated,
C           1 = some error occurs in the evaluation.

C     PARAMETERS
      integer ncmax
      parameter ( ncmax    =                630000 )
      integer npmax
      parameter ( npmax    =               1000000 )
      integer nmax
      parameter ( nmax     = 2 * ( npmax + ncmax ) )
      integer mmax
      parameter ( mmax     =         npmax + ncmax )
      integer jcnnzmax
      parameter ( jcnnzmax =       nmax + 2 * mmax )

C     COMMON ARRAYS
      logical equatn(mmax)

C     SCALAR ARGUMENTS
      integer n,m,inform

C     ARRAY ARGUMENTS
      double precision x(n),lambda(m),rho(m),nal(n)

C     LOCAL SCALARS
      integer flag,i,jcnnz

C     LOCAL ARRAYS
      integer jcfun(jcnnzmax),jcvar(jcnnzmax)
      double precision c(mmax),jcval(jcnnzmax),dpdc(mmax)
 
C     COMMON BLOCKS
      common /pdata/ equatn

C     INITIALIZE

      inform = 0

C     COMPUTE CONSTRAINTS

      call evalc(n,x,m,c,flag)
      if ( flag .ne. 0 ) then
          inform = - 90
          return
      end if

C     COMPUTE dP/dc

      do i = 1,m
          call evaldpdy(c(i),rho(i),lambda(i),equatn(i),dpdc(i))
      end do

C     COMPUTE \nabla L(x) = \nabla f(x) + \sum_{i=1}^m dPdc * dcdx

      do i = 1,n
          nal(i) = 0.0d0
      end do

      call evalgjc(n,x,m,jcvar,jcfun,jcval,jcnnz,flag)
      if ( flag .ne. 0 ) then
          inform = - 91
          return
      end if

      do i = 1,jcnnz
          if ( jcfun(i) .eq. 0 ) then
              nal(jcvar(i)) = nal(jcvar(i)) + jcval(i)
          else 
              nal(jcvar(i)) = nal(jcvar(i)) + dpdc(jcfun(i)) * jcval(i)
          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     *****************************************************************

      subroutine spg(n,x,m,lambda,rho,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
      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,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,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,g,inform)
      else ! if ( gtype .eq. 1 ) then
          call evalnaldiff(n,x,m,lambda,rho,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,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,g,inform)
      else ! if ( gtype .eq. 1 ) then
          call evalnaldiff(n,x,m,lambda,rho,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,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
      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,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,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,g,sterel,steabs,inform)

C     SCALAR ARGUMENTS
      integer n,m,inform
      double precision sterel,steabs

C     ARRAY ARGUMENTS
      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 is prepared to replace 
C     the subroutine evalnal (to evaluate the gradient vector) in the 
C     case of the lastest have not being provided by the user.
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              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     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,fplus,inform)
          if ( inform .lt. 0 ) then
              return
          end if

          x(j) = tmp - step
          call evalal(n,x,m,lambda,rho,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)
      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
