      program ivmmaispgc

C     IVM Test driver. This master file uses IVM (Inexact Variable
C     Metric Method) to solve some CUTE linearly inequality constrained 
C     problems with interior initial point.
C
C     R. Andreani, E. G. Birgin, J. M. Martínez and J. Yuan, "Spectral 
C     projected gradient and variable metric methods for optimization 
C     with linear inequalities", IMA Journal of Numerical Analysis 25, 
C     pp. 221-252, 2005.
C
C     This version 25 SET 2003 by E.G.Birgin, J.M.Martinez and J.Yuan.
C     Final revision ... by E.G.Birgin, J.M.Martinez and J. Yuan.

C     PARAMETERS
      integer mmax
      parameter (mmax    =   100 000)
      integer nmax
      parameter (nmax    =     4 000)
      integer anzemax
      parameter (anzemax =  3 * mmax)

C     COMMON SCALARS
      integer anze,n,morig

C     COMMON ARRAYS
      integer arow(anzemax),acol(anzemax)
      double precision aval(anzemax)

C     LOCAL SCALARS
      logical extrap,hess,efirst,lfirst,nvfirst,interior
      integer fcnt,flag,gcnt,hcnt,i,innit,j,kprint,m,maxfc,maxinnit,
     +        maxoutit,outit,output,p
      double precision amax,dsupn,deucn2,epsdsn,epsden,epsmac,epsnlsn,
     +        eta,f,gamma,infty,lmin,lmax,steabs,sterel,tau,tau1,tau2,
     +        theta,tol,tmp
      character * 10 pname,xnames(nmax),gnames(mmax)
      real time,dtime

C     LOCAL ARRAYS
      logical equatn(mmax),linear(mmax)
      double precision Ax(mmax),b(mmax),l(nmax),u(nmax),x(nmax),
     +        y(mmax),cl(mmax),cu(mmax),czero(mmax),gc(nmax),xtemp(nmax)
      real dum(2)

C     EXTERNAL SUBROUTINES
      external ivmispg,ccifg,cdh,cfn,cofg

C     INTRINSIC FUNCTIONS
      intrinsic dtime

C     COMMON BLOCKS
      common /orignc/ morig
      common /ldata1/ arow,acol,aval,anze,n

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,x,l,u,nmax,equatn,linear,y,cl,cu,mmax,
     +            efirst,lfirst,nvfirst)

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

      close(10)

C     MOUNT LINEAR CONSTRAINTS IN THE FORMAT A x <= b 
C     (also including bound constraints)

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

      call cfn(n,morig,xtemp,f,mmax,czero)

      anze = 0

      m = morig

      do i = 1,morig

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

          else 

              call ccifg(n,i,x,tmp,gc,.true.)

              if ( cl(i) .gt. -1.0d+20 .and. cu(i) .ge. 1.0d+20 ) then

                  do j = 1,n
                      if ( gc(j) .ne. 0.0d0 ) then
                          anze = anze + 1
                          arow(anze) = i
                          acol(anze) = j
                          aval(anze) = - gc(j)
                      end if
c                     A(i,j) = - gc(j)
                  end do
                  b(i) = - ( cl(i) - czero(i) )

              else if ( cl(i).le.-1.0d+20 .and. cu(i).lt.1.0d+20 ) then

                  do j = 1,n
                      if ( gc(j) .ne. 0.0d0 ) then
                          anze = anze + 1
                          arow(anze) = i
                          acol(anze) = j
                          aval(anze) = gc(j)
                      end if
c                     A(i,j) = gc(j)
                  end do
                  b(i) = ( cu(i) - czero(i) )

              else

                  do j = 1,n
                      if ( gc(j) .ne. 0.0d0 ) then
                          anze = anze + 1
                          arow(anze) = i
                          acol(anze) = j
                          aval(anze) = - gc(j)
                      end if
c                     A(i,j) = - gc(j)
                  end do
                  b(i) = - ( cl(i) - czero(i) )

                  m = m + 1

                  if ( m .gt. mmax ) then
                      write(*,*) 'Increase mmax and run again'
                      stop
                  end if

                  do j = 1,n
                      if ( gc(j) .ne. 0.0d0 ) then
                          anze = anze + 1
                          arow(anze) = m
                          acol(anze) = j
                          aval(anze) = gc(j)
                      end if
c                     A(m,j) = gc(j)
                  end do
                  b(m) = ( cu(i) - czero(i) )

              end if

          end if

      end do

      do i = 1,n
          if ( l(i) .gt. -1.0d+20 ) then
              m = m + 1
c             do j = 1,n
c                 A(m,j) = 0.0d0
c             end do
c             A(m,i) = - 1.0d0
              anze = anze + 1
              arow(anze) = m
              acol(anze) = i
              aval(anze) = - 1.0d0
              b(m) = - l(i)
          end if
      end do

      do i = 1,n
          if ( u(i) .lt. 1.0d+20 ) then
              m = m + 1
c             do j = 1,n
c                 A(m,j) = 0.0d0
c             end do
c             A(m,i) = 1.0d0
              anze = anze + 1
              arow(anze) = m
              acol(anze) = i
              aval(anze) = 1.0d0
              b(m) = u(i)
          end if
      end do

C     SAVE A IN SPARSE FORMAT

c     anze = 0
c     do i = 1,m
c         do j = 1,n
c             if ( A(i,j) .ne. 0.0d0 ) then
c                 anze = anze + 1
c                 arow(anze) = i
c                 acol(anze) = j
c                 aval(anze) = A(i,j)
c             end if
c         end do
c     end do

C     TEST WETHER THE INITIAL POINT IS INTERIOR

      interior = .true.

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

      do i = 1,anze
          Ax(arow(i)) = Ax(arow(i)) + aval(i) * x(acol(i))
      end do

      do i = 1,m
c         tmp = 0.0d0
c         do j = 1,n
c             tmp = tmp + A(i,j) * x(j)
c         end do
c         if ( b(i) - tmp .le. 0.0d0 ) then
          if ( b(i) - Ax(i) .le. 0.0d0 ) then
              interior = .false.
          end if
      end do

      if ( .not. interior ) then
          write(*,*) 'The initial point is not interior'
          stop
      end if

C     SET THE INITIAL ESTIMATION OF THE LAGRANGE MULTIPLIERS

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

C     SETUP OPTIMIZER PARAMETERS 

      call iniivm(epsdsn,epsden,epsnlsn,maxoutit,maxinnit,maxfc,p,
     +            extrap,hess,output,kprint,eta,theta,tau,tol,
     +            lmin,lmax,gamma,sterel,steabs,epsmac,infty,tau1,
     +            tau2)

C     CALL THE OPTIMIZER

      time= dtime(dum)

      call ivmispg(n,m,x,y,arow,acol,aval,anze,b,epsdsn,epsden,epsnlsn,
     +             maxoutit,maxinnit,maxfc,p,extrap,hess,output,kprint,
     +             eta,theta,tau,tol,lmin,lmax,gamma,sterel,steabs,
     +             epsmac,infty,tau1,tau2,f,dsupn,deucn2,amax,outit,
     +             innit,fcnt,gcnt,hcnt,flag)

      time= dtime(dum)

      time= dum(1)

C     WRITE STATISTICS

      open(30, file='ivmline.txt')
      write(30,fmt=9000) pname,n,m,outit,fcnt,gcnt,hcnt,innit,time,f,
     +                   dsupn,sqrt(deucn2),amax,flag
      close(30)

      write(*,*)
      write(*,*) 'Solution: '
      do i = 1,n/3
          write(*,fmt=9010) (x(3*(i-1)+j),j=1,3)
      end do

      stop

C     NON EXECUTABLE STATEMENTS

 9000 format(A8,' &',I4,' &',I4,' &',I6,' &',I6,' &',I6,' &',I6,' &',
     +       I6,' &',F8.2,' &',1PD11.3,' &',1PD9.1,' &',1PD9.1,' &',
     +       1PD9.1,' &',I1)
 9010 format(3(1X,1P,D17.9))

      end

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

      subroutine evalf(n,x,f,flag)

C     This subroutine computes the objective function.
C
C     On Entry:
C
C     n     integer,
C           number of variables,
C
C     x     double precision x(n),
C           current point,
C
C     On Return
C
C     f     double precision,
C           objective function value at x,
C
C     flag integer
C           0 means "no errors",
C           1 means "some error occured in the objective function evaluation".

C     SCALAR ARGUMENTS
      double precision f
      integer n,flag

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      double precision g(1)

C     EXTERNAL SUBROUTINES
      external cofg
 
      flag = 0

      call cofg(n,x,f,g,.false.)

      end


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

      subroutine evalg(n,x,g,flag)

C     This subroutine computes the gradient of the objective 
C     function.
C
C     On Entry:
C
C     n     integer,
C           number of variables,
C
C     x     double precision x(n),
C           current point,
C
C     On Return
C
C     g     double precision g(n),
C           gradient of the objective function at x,
C
C     flag integer
C           0 means "no errors",
C           1 means "some error occured in the gradient evaluation".

C     SCALAR ARGUMENTS
      integer n,flag

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

C     LOCAL SCALARS
      double precision f

C     EXTERNAL SUBROUTINES
      external cofg

      flag = 0
 
      call cofg(n,x,f,g,.true.)

      end

c     ******************************************************
c     ******************************************************

      subroutine iniivm(epsdsn,epsden,epsnlsn,maxoutit,maxinnit,maxfc,p,
     +                  extrap,hess,output,kprint,eta,theta,tau,tol,
     +                  lmin,lmax,gamma,sterel,steabs,epsmac,infty,tau1,
     +                  tau2)

C     This subroutine set the IVM parameters with its default values

C     SCALAR ARGUMENTS
      logical hess,extrap
      integer kprint,maxfc,maxinnit,maxoutit,output,p
      double precision epsdsn,epsden,epsnlsn,eta,theta,tau,tol,lmin,
     +        lmax,gamma,tau1,tau2,sterel,steabs,epsmac,infty

      epsdsn   =  1.0d-06
      epsden   =  1.0d-06
c     epsnlsn  =  1.0d-08 <-- For all problems
c     and the line below just for KSIP and SIPOW4G
      epsnlsn  =  1.0d-05
      maxoutit =      100 
      maxinnit =   100000   
      maxfc    =     2000
      p        =       10
      extrap   =    .true.
      hess     =    .true.
      output   =        3
      kprint   =        5
      eta      =   0.50d0
      theta    =   0.75d0
      tau      =  1.0d-01
      tol      =   0.90d0
      lmin     =  1.0d-03
      lmax     =  1.0d+03
      gamma    =  1.0d-04
      sterel   =  1.0d-07
      steabs   =  1.0d-10
      infty    =  1.0d+20

      call mcheps(epsmac)

      tau1 = epsmac ** ( 1.0d0 / 3.0d0 )
      tau2 = epsmac ** ( 1.0d0 / 3.0d0 )

      end

c     ******************************************************
c     ******************************************************

      subroutine mcheps(epsmac)

C     This subroutine computes an approximation of machine epsilon.

C     SCALAR ARGUMENTS
      double precision epsmac

      epsmac = 1.0d0

 20   if ( ( 1.0d0 + epsmac ) .ne. 1.0d0 ) then
          epsmac = epsmac / 2.0d0
          go to 20
      end if

      epsmac = 2.0d0 * epsmac

      end
