      program algencanma

C     ALGENCAN Test driver. This master file uses ALGENCAN (an Augmented
C     Lagrangian method for inequality constraints and bounds that uses 
C     GENCAN to solve the subproblems) to solve a very easy problem. 67 
C     different Augmented Lagrangian methods are implemented.
C
C     E. G. Birgin, R. Castillo and J. M. Martínez, "Numerical comparison 
C     of Augmented Lagrangian algorithms for nonconvex problems",
C     Computational Optimization and Applications 31, pp. 31-56, 2005.

C     PARAMETERS

      logical nearlyq
      integer nmax,mmax,mininterp,maxextrap,ncomp,cgscre,gtype,trtype,
     +        cgmaxit,maxit,maxfc,maxitnfp,maxitngp,maxitnqmp,iprint
      double precision eta,delmin,lammax,lammin,theta,gamma,beta,sigma1,
     +        sigma2,nint,next,sterel,steabs,epsrel,epsabs,infrel,
     +        infabs,fmin,delta0,cggpnf,cgepsi,cgepsf,epsnqmp,epsgpen,
     +        epsgpsn,epsfc,epsnfp

      parameter (nmax      =  100000)
      parameter (mmax      =  100000)
      parameter (eta       =   0.9d0)
      parameter (delmin    =   0.1d0)
      parameter (lammax    = 1.0d+02)
      parameter (lammin    = 1.0d-02)
      parameter (theta     = 1.0d-06)
      parameter (gamma     = 1.0d-04)
      parameter (beta      =   0.5d0)
      parameter (sigma1    =   0.1d0)
      parameter (sigma2    =   0.9d0)
      parameter (nint      =   2.0d0)
      parameter (next      =   2.0d0)
      parameter (sterel    = 1.0d-07)
      parameter (steabs    = 1.0d-10)
      parameter (epsrel    = 1.0d-10)
      parameter (epsabs    = 1.0d-20)
      parameter (infrel    = 1.0d+20)
      parameter (infabs    = 1.0d+99)
      parameter (mininterp =       4)
      parameter (maxextrap =       3)
      parameter (ncomp     =      12)
      parameter (cgscre    =       2)
      parameter (gtype     =       0)
      parameter (trtype    =       0)
      parameter (nearlyq   = .false.)
      parameter (fmin      = -infrel)
      parameter (delta0    =  -1.0d0)
      parameter (cggpnf    =  1.0d-4)
      parameter (cgepsi    =  1.0d-1)
      parameter (cgepsf    =  1.0d-5)
      parameter (epsnqmp   =  1.0d-4)
      parameter (cgmaxit   =      -1)
      parameter (maxit     =    5000)
      parameter (maxfc     =  100000)
      parameter (maxitnfp  =   maxit)
      parameter (maxitngp  =   maxit)
      parameter (maxitnqmp =       5)
      parameter (iprint    =       3)
      parameter (epsgpen   =  1.0d-4)
      parameter (epsgpsn   =  1.0d-4)
      parameter (epsfc     =  1.0d-4)
      parameter (epsnfp    =   0.0d0)

C     LOCAL SCALARS

      integer i,n,outiter,totiter,totfcnt,totgcnt,totcgit,inform,
     +        totspgit,tottnit,htvtype,inform2,icomb,icomb1,
     +        icomb2
      double precision f,al,gpeucn2,gpsupn,csupn,mumenor
      real elatime

C     LOCAL ARRAYS

      integer w2(nmax)
      double precision x(nmax),l(nmax),u(nmax),c(mmax),mu(mmax),
     +        cprev(mmax),nal(nmax),w1(10*nmax+maxitngp),xini(nmax),
     +        lini(nmax),uini(nmax)
      character * 10 pname
      character * 1 fact,opti,crit,erro

C     COMMON SCALAR

      integer nc,m,ptype,ttype
      double precision tconst,umuini,urhoini,urhomul,urhofrac,umumin,
     +        umumax

C     COMMON ARRAYS

      double precision xs(nmax,2),cs(mmax,2),mubar(mmax),
     +        zeta(mmax),rho(mmax)

C     COMMON BLOCKS

      common /constr/xs,cs,nc
      common /auglag/mubar,zeta,rho,m
      common /p/ ptype,ttype,tconst
      common /user/umuini,urhoini,urhomul,urhofrac,umumin,umumax

C     EXTERNAL SUBROUTINES

      integer i1,i2,i3,i4,i5
      double precision muini(162),rhoini(162),rhomul(162),rhofrac(162),
     +        mumin(162),mumax(162),muiniop(2),rhoiniop(3),rhomultop(3),
     +        rhofracop(3),muminop(3),mumaxop(3)

      muiniop(1)   =   1.0d-6
      muiniop(2)   =   1.0d0

      rhoiniop(1)  =   1.0d-3
      rhoiniop(2)  =   1.0d0
      rhoiniop(3)  =  10.0d0

      rhomultop(1) =   2.0d0
      rhomultop(2) =  10.0d0
      rhomultop(3) = 100.0d0

      rhofracop(1) = 1.0d-2
      rhofracop(2) = 1.0d-1
      rhofracop(3) = 0.5d0

      muminop(1)   = 1.0d-6
      mumaxop(1)   = 1.0d+20
      muminop(2)   = 1.0d-6
      mumaxop(2)   = 1.0d+6
      muminop(3)   = 1.0d-3
      mumaxop(3)   = 1.0d+3

      i = 1
      do i1 = 1,2
          do i2 = 1,3
              do i3 = 1,3
                  do i4 = 1,3
                      do i5 = 1,3
                          muini  (i) = muiniop  (i1)
                          rhoini (i) = rhoiniop (i2)
                          rhomul (i) = rhomultop(i3)
                          rhofrac(i) = rhofracop(i4)
                          mumin  (i) = muminop  (i5)
                          mumax  (i) = mumaxop  (i5)
                          i = i + 1
                      end do
                  end do
              end do
          end do
      end do

C     SETUP PROBLEM DATA

      pname = 'ML1       '

      n = 2

      do i = 1,n
          lini(i) = -1.0d+99
          uini(i) =  1.0d+99
          xini(i) =  10.0d0
      end do

      m = 2

C     SOLVE THE PROBLEM WITH THE DIFFERENT METHODS

      open(30,file='phrline.txt')
      open(40,file='phrline2.txt')
      open(50,file='param.txt')

 10   continue

C     SET UP METHOD

      read(50,*) ptype,ttype,tconst,icomb1,icomb2

      if ( ptype .eq. 99 ) then
          go to 30
      end if
 
      if ( ptype .eq. 4 .and. ttype .eq. 1 ) then
          htvtype = 0
      else
          htvtype = 1
      end if

C     SOLVE IT

      do icomb = icomb1,icomb2

          umuini   = muini(icomb)
          urhoini  = rhoini(icomb)
          urhomul  = rhomul(icomb)
          urhofrac = rhofrac(icomb)
          umumin   = mumin(icomb)
          umumax   = mumax(icomb)

          nc = 0

          do i = 1,n
              x(i) = xini(i)
              l(i) = lini(i)
              u(i) = uini(i)
          end do

          call algencan(pname,n,x,l,u,m,c,mu,mubar,rho,zeta,epsfc,
     +    maxit,maxfc,epsnqmp,maxitnqmp,al,nal,gpeucn2,gpsupn,csupn,
     +    outiter,totiter,totfcnt,totgcnt,totcgit,totspgit,
     +    tottnit,inform,nearlyq,epsgpen,epsgpsn,epsnfp,fmin,delta0,
     +    cgmaxit,cgscre,cggpnf,cgepsi,cgepsf,maxitnfp,maxitngp,gtype,
     +    htvtype,trtype,iprint,w1,w2,cprev,eta,delmin,lammax,lammin,
     +    theta,gamma,beta,sigma1,sigma2,nint,next,mininterp,maxextrap,
     +    ncomp,sterel,steabs,epsrel,epsabs,infrel,infabs)

          elatime = 0.0d0

          if ( inform .ge. 0 ) then
              call evalfc(n,x,f,m,c,inform2)
          end if

          if ( inform .lt. 0 ) then
              erro = 'E'
              fact = 'N'
              opti = 'N'
              crit = 'N'

              csupn  = 1.0d+99
              gpsupn = 1.0d+99
              f      = 1.0d+99
          else
              erro = 'N'

              if ( csupn .le. epsfc ) then
                  fact = 'F'
              else
                  fact = 'N'
              end if

              if ( gpsupn .le. epsgpsn ) then
                  opti = 'O'
              else
                  opti = 'N'
              end if

              if ( fact .eq. 'F' .and. opti .eq. 'O' ) then
                  crit = 'C'
              else
                  crit = 'N'
              end if
          end if

          mumenor = 1.0d+99
          do i = 1,m
              mumenor = min( mumenor, abs( mu(i) ) )
          end do

          write(30,997) ptype,ttype,pname,n,m,outiter,totiter,totspgit,
     +                  tottnit,totfcnt,totgcnt,totcgit,elatime,f,csupn,
     +                  gpsupn,mumenor,inform
          write(*, 998) pname,n,m,ptype,ttype,elatime,f,csupn,
     +                  gpsupn,fact,opti,crit,erro
          write(40,998) pname,n,m,ptype,ttype,elatime,f,csupn,
     +                  gpsupn,fact,opti,crit,erro

      end do

C     ITERATE
 
      go to 10
 
C     CLOSE FILES AND STOP

 30   continue

      close(30)
      close(40)
      close(50)

      stop

 997  format(i2,'-',i2,'&',a8,' (',i6,' , ',i6,') &',i2,' &',i7,' &',i7,
     +       ' & ',i7,' &',i7,' &',i7,' &',i7,' &',f8.2,' &',1PD12.4,
     +       ' &',1PD11.3,' &',1PD11.3,' &',1PD11.3,' &',i3)

 998  format(a8,' ',I4,' ',I5,' ',I1,' ',I2,' ',F15.10,' ',
     +       1PD14.7,' ',1PD7.1,' ',1PD7.1,' ',a1,' ',a1,' ',a1,' ',a1)

 999  format(A8,1X,I2,1X,I2,1X,18(1X,I1),
     +       /,15X,18(1X,I1),
     +       /,15X,18(1X,I1))

      end


C     ******************************************************
C     ******************************************************
      subroutine evalfc(n,x,f,m,c,flag)

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     flag integer
C           0 means "no errors",
C           1 means "some error occured in the objective function evaluation",
C           2 means "some error occured in the constraints evaluation".

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

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

      flag = 0

      f = x(2)

      c(1) = - x(2)
      c(2) =   x(1) ** 2 - x(2)


c     f = x(2)

c     c(1) =   x(1) ** 2 - x(2)
c     c(2) = - x(1) ** 2 - x(2)


c     f = x(2)

c     c(1) =         x(1) ** 2 - x(2)
c     c(2) = 2.0d0 * x(1) ** 2 - x(2)

      end


C     ******************************************************
C     ******************************************************
      subroutine evalgjc(n,x,m,jcvar,jcfun,jcval,jcnze,flag)

C     This subroutine computes the gradient of the objective 
C     function 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     g     double precision g(n),
C           gradient of the objective function at x,
C
C     jcnze integer
C           number of non-null Jacobian elements,
C
C     jcvar,jcfun integer jcvar(jcnze),jcfun(jcnze)
C           see below,
C
C     jcval double precision jcval(jcnze)
C           the value of the non-zero Jacobian position 
C           [jcfun(i),jcvar(i)] is jcval(i).
C
C           NOTE: jcfun(i) equal zero means objective function.
C
C     flag integer
C           0 means "no errors",
C           1 means "some error occured in the gradient evaluation",
C           2 means "some error occured in the Jacobian evaluation".

C     PARAMETERS
      integer jcnzemax
      parameter (jcnzemax=100000)

C     SCALAR ARGUMENTS
      integer n,m,flag,jcnze

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

      flag = 0

      jcnze = 4

      jcfun(1) = 0
      jcvar(1) = 2
      jcval(1) = 1.0d0

      jcfun(2) =   1
      jcvar(2) =   2
      jcval(2) = - 1.0d0

      jcfun(3) =   2
      jcvar(3) =   1
      jcval(3) =   2.0d0 * x(1)

      jcfun(4) =   2
      jcvar(4) =   2
      jcval(4) = - 1.0d0


c     jcnze = 5

c     jcfun(1) = 0
c     jcvar(1) = 2
c     jcval(1) = 1.0d0

c     jcfun(2) =   1
c     jcvar(2) =   1
c     jcval(2) =   2.0d0 * x(1)

c     jcfun(3) =   1
c     jcvar(3) =   2
c     jcval(3) = - 1.0d0

c     jcfun(4) =   2
c     jcvar(4) =   1
c     jcval(4) = - 2.0d0 * x(1)

c     jcfun(5) =   2
c     jcvar(5) =   2
c     jcval(5) = - 1.0d0


c     jcnze = 5

c     jcfun(1) = 0
c     jcvar(1) = 2
c     jcval(1) = 1.0d0

c     jcfun(2) =   1
c     jcvar(2) =   1
c     jcval(2) =   2.0d0 * x(1)

c     jcfun(3) =   1
c     jcvar(3) =   2
c     jcval(3) = - 1.0d0

c     jcfun(4) =   2
c     jcvar(4) =   1
c     jcval(4) =   4.0d0 * x(1)

c     jcfun(5) =   2
c     jcvar(5) =   2
c     jcval(5) = - 1.0d0

      end
