      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 problems from the CUTE
C     collection. 67 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    =       0)
      parameter (epsgpen   =  1.0d-4)
      parameter (epsgpsn   =  1.0d-4)
      parameter (epsfc     =  1.0d-4)
      parameter (epsnfp    =   0.0d0)

C     LOCAL SCALARS

      logical efirst,lfirst,nvfirst
      integer i,j,n,outiter,totiter,totfcnt,totgcnt,totcgit,inform,
     +        totspgit,tottnit,ntimes,htvtype,inform2,icomb,icomb1,
     +        icomb2
      double precision f,al,gpeucn2,gpsupn,csupn,mumenor
      real elatime,dum,cputim

C     LOCAL ARRAYS

      logical equatn(mmax),linear(mmax)
      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 xnames(nmax),gnames(mmax),pname
      character * 1 fact,opti,crit,erro

C     COMMON SCALAR

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

C     COMMON ARRAYS

      integer rcl(mmax),rcu(mmax)
      double precision xs(nmax,2),cs(mmax,2),mubar(mmax),
     +        zeta(mmax),rho(mmax),cl(mmax),cu(mmax)

C     COMMON BLOCKS

      common /constr/xs,cs,nc
      common /auglag/mubar,zeta,rho,m
      common /cbds/cl,cu,rcl,rcu,morig
      common /p/ ptype,ttype,tconst
      common /user/umuini,urhoini,urhomul,urhofrac,umumin,umumax

C     EXTERNAL SUBROUTINES

      external csetup,cnames,cfn,cputim

      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     SET UP SIF DATA

      open(10,file='OUTSDIF.d',form='formatted',status='old')

      rewind 10

      efirst  = .true.
      lfirst  = .false.
      nvfirst = .false.

      call csetup(10,20,n,morig,xini,lini,uini,nmax,equatn,linear,mu,
     +            cl,cu,mmax,efirst,lfirst,nvfirst)

      do i = 1,n
          if ( lini(i) .gt. uini(i) ) then
              write(*,*) 'Wrong bounds. Some lb is greater than the '
              write(*,*) 'corresponding ub. Infeasible problem.     '
              stop
          end if
      end do

      do i = 1,morig
          if ( equatn(i) ) then
              write(*,*) 'This master file is not ready to handle ',
     +                   'equalities'
              stop
          end if
      end do

      do i = 1,morig
          if ( mu(i) .ne. 0.0d0 ) then
              write(*,*) 'Opa! Initial estimation for the Lagrangian ',
     +                   'multipliers given by CUTEr is not zero (as ',
     +                   'we assumed). Zero will be used instead of ',
     +                    mu(i)
          end if
      end do

      m = morig

      do i = 1,morig

          if ( cl(i) .le. - infrel .and. cu(i) .ge. infrel ) then
              write(*,*) 'This master file is not ready to handle ',
     +                   'unreal inequaltities of type ',
     +                   '- infrel <= c(x) <= infrel'
              stop

          else if ( cl(i) .gt. -infrel .and. cu(i) .ge. infrel ) then
              rcl(i) = i
              rcu(i) = 0

          else if ( cl(i) .le. -infrel .and. cu(i) .lt. infrel ) then
              rcl(i) = 0
              rcu(i) = i

          else
              m      = m + 1
              rcl(i) = i
              rcu(i) = m
          end if

      end do

      if ( m .gt. mmax ) then
          write(*,*) 'Increase mmax!',' m = ',m,' mmax = ',mmax
          stop
      end if

      call cnames(n,morig,pname,xnames,gnames)

      close(10)

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

          elatime = cputim(dum)

          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 = cputim(dum) - elatime

C         CHECK IF MORE THAN ONE RUN IS NECESSARY

          if ( inform .ge. 0 .and. csupn .le. epsfc ) then

              if ( elatime .lt. 0.1d0 ) then
                  ntimes =   1000
              else if ( elatime .lt. 1.0d0 ) then
                  ntimes =    100
              else if ( elatime .lt. 10.0d0 ) then
                  ntimes =     10
              else
                  ntimes =      1
              end if

          else
              ntimes  = 0
          end if

C         IF IT IS, RUN AS MANY TIMES AS REQUIRED

 20       continue

          if ( ntimes .gt. 1 ) then

              elatime = cputim(dum)

              do j = 1,ntimes

                  nc= 0

                  do i = 1,n
                      x(i) = xini(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)

              end do

              elatime = cputim(dum) - elatime

              if ( elatime .lt. 0.3d0 .and. ntimes .ne. 10000 ) then
                  ntimes = 10000
                  go to 20
              end if

              elatime = elatime / ntimes

          end if

          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     PARAMETERS
      integer mmax
      parameter (mmax=100000)

C     COMMON SCALARS
      integer morig

C     COMMON ARRAYS
      integer rcl(mmax),rcu(mmax)
      double precision cl(mmax),cu(mmax)

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

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

C     LOCAL SCALARS
      integer i
      double precision tmp

C     COMMON BLOCKS
      common /cbds/cl,cu,rcl,rcu,morig

C     EXTERNAL SUBROUTINES
      external cfn
 
      flag = 0

      call cfn(n,morig,x,f,mmax,c)

      do i = 1,morig

          tmp = c(i)

          if ( rcl(i) .ne. 0 ) then
              c(rcl(i)) = - tmp + cl(i)
          end if

          if ( rcu(i) .ne. 0 ) then
              c(rcu(i)) =   tmp - cu(i)
          end if

      end do

      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 mmax,jcnzemax
      parameter (mmax=100000,jcnzemax=100000)

C     COMMON SCALARS
      integer morig

C     COMMON ARRAYS
      integer rcl(mmax),rcu(mmax)
      double precision cl(mmax),cu(mmax)

C     SCALAR ARGUMENTS
      integer n,m,flag,jcnze

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

C     LOCAL SCALARS
      integer i,dummy1,jcnzeorig
      double precision dummy2(1)

C     EXTERNAL SUBROUTINES
      external csgr

C     COMMON BLOCKS
      common /cbds/cl,cu,rcl,rcu,morig

      flag = 0
 
      call csgr(n,morig,.false.,dummy1,dummy2,x,jcnze,jcnzemax,
     +jcval,jcvar,jcfun)

      if ( 2 * jcnze .gt. jcnzemax ) then
          write(*,*) 'jcnzemax must be increased up to ',2 * jcnze
          stop
      end if

      jcnzeorig = jcnze

      do i = 1,jcnzeorig
          if ( jcfun(i) .ne. 0 ) then
              if ( rcl(jcfun(i)) .ne. 0 ) then

                  if ( rcu(jcfun(i)) .ne. 0 ) then

                      jcnze = jcnze + 1
                      jcfun(jcnze) = rcu(jcfun(i))
                      jcvar(jcnze) = jcvar(i)
                      jcval(jcnze) = jcval(i)

                  end if

                  jcval(i) = - jcval(i)

              end if
          end if
      end do

      end
