C     =================================================================
C     File: ipoptma.f
C     =================================================================

C     =================================================================
C     Module: Main program
C     =================================================================

C     Last update of any of the component of this module: 
C
C     November 27, 2006.

C     Users are encouraged to download periodically updated versions of 
C     this code at LA CUMPARSITA home page:
C
C     www.ime.usp.br/~egbirgin/tango/collection/

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

      program ipoptma

      implicit none

C     PARAMETERS
      integer jcnnzmax,liw,lrw,mmax,nmax
      parameter ( jcnnzmax  = 5000000 )
      parameter ( liw       =       0 )
      parameter ( lrw       =       0 )
      parameter ( mmax      =  500000 )
      parameter ( nmax      =  500000 )

C     LOCAL SCALARS
      logical wstart
      integer i,inform,iter,itmp,iw,m,mtmp,ntmp,n,nz,nargs,ln,un
      double precision cnorm,dat,f,gnorm,rw,snorm,time

C     LOCAL ARRAYS
      character * 20 cargs(50)
      logical equatn(mmax),linear(mmax)
      integer acon(jcnnzmax),avar(jcnnzmax),idat(1),lvar(nmax),
     +        uvar(nmax)
      double precision args(50),a(jcnnzmax),c(mmax),g(nmax),
     +        lambda(mmax),l(nmax),lval(nmax),llambda(nmax),u(nmax),
     +        uval(nmax),ulambda(nmax),rho(mmax),x(nmax)

C     EXTERNAL FUNCTIONS
      double precision ipopt

C     EXTERNAL SUBROUTINES
      external ev_f,ev_c,ev_g,ev_a,ev_h,ev_hlv_dummy,ev_hov_dummy,
     +        ev_hcv_dummy

C     DECLARATIONS RELATED TO THE TIME MEASUREMENT USING DTIME
      real dtime
      external dtime

      real dum(2)
      data dum/0.0,0.0/

C     SET UP PROBLEM DATA
      call inipt(n,x,l,u,m,lambda,rho,equatn,linear)

C     Set bounds in sparse format and initialize Lagrange multipliers 
C     for bound constraints

      ln = 0 
      do i = 1,n
          if ( l(i) .gt. - 1.0d+20 ) then
              ln      = ln + 1
              lvar(ln)    = i
              lval(ln)    = l(i)
              llambda(ln) = 0.0d0
          end if
      end do

      un = 0 
      do i = 1,n
          if ( u(i) .lt. 1.0d+20 ) then
              un          = un + 1
              uvar(un)    = i
              uval(un)    = u(i)
              ulambda(un) = 0.0d0
          end if
      end do

C     Subroutine ev_a needs to know the number of constraints
      idat(1) = m

C     LOOK FOR A WARM-START INITIAL POINT

      wstart = .false.

      open(10,err=100,file='inipoint.txt',status='old')

      read(10,*) ntmp,mtmp
      if ( ntmp .ne. n .or. mtmp .ne. m ) then
          close(10)
          write(*,*) 'Wrong dimensions in inipoint.txt.'
          go to 200
      end if

      do i = 1,n
          read(10,*) itmp,x(i)
          if ( itmp .ne. i ) then
              close(10) 
              write(*,*) 'Wrong index for x in inipoint.txt.'
              go to 200
          end if
      end do

      do i = 1,m
          read(10,*) itmp,lambda(i)
          if ( itmp .ne. i ) then
              close(10) 
              write(*,*) 'Wrong index for lambda in inipoint.txt.'
              go to 200
          end if
      end do

      do i = 1,ln
          read(10,*) itmp,llambda(i)
          if ( itmp .ne. lvar(i) ) then
              close(10) 
              write(*,*) 'Wrong index for llambda in inipoint.txt.'
              go to 200
          end if
      end do

      do i = 1,un
          read(10,*) itmp,ulambda(i)
          if ( itmp .ne. uvar(i) ) then
              close(10)
              write(*,*) 'Wrong index for ulambda in inipoint.txt.'
              go to 200
          end if
      end do

      close(10)

      write(*,*) 'Using warm-start!'
      wstart = .true.
      go to 200

 100  continue
      write(*,*) 'inipoint.txt file was not found.'
      write(*,*) 'The initial point defined in inip will be used.'

 200  continue

C     SET SOME SOLVER ARGUMENTS
      if ( .not. wstart ) then
          nargs = 0

      else
          nargs = 1

          cargs(1) = 'IMUINIT'
          args(1)  = 0

C         cargs(2) = 'DMU0'
C         args(2)  = 1.0d-08

C         cargs(3) = 'DBNDFRAC'
C         args(3)  = 1.0d-08

C         cargs(4) = 'DBNDPUSH'
C         args(4)  = 1.0d-08

C         cargs(5) = 'DMOVEBOUNDS'
C         args(5)  = 1.0d-08

      end if

C     CALL OPTIMIZATION SOLVER
      time = dtime(dum)

      f = ipopt(n,x,m,ln,lvar,lval,un,uvar,uval,llambda,ulambda,lambda,
     +    c,lrw,rw,liw,iw,iter,inform,ev_f,ev_c,ev_g,ev_a,ev_h,
     +    ev_hlv_dummy,ev_hov_dummy,ev_hcv_dummy,dat,idat,nargs,args,
     +    cargs)

      time = dtime(dum)
      time = dum(1)

C     SAVE THE SOLUTION

      open(10,file='solution.txt')

C     Solution point
      write(10,9000)
      do i = 1,n
          write(10,9010) i,x(i)
      end do

C     Lagrange multipliers
      if ( m .gt. 0 ) then
          write(10,9020)
          do i = 1,m
              write(10,9010) i,lambda(i)
          end do
      end if

C     Lagrange multipliers for the bound constraints
      write(10,9030)
      do i = 1,ln
          write(10,9010) lvar(i),llambda(i)
      end do

      write(10,9040)
      do i = 1,un
          write(10,9010) uvar(i),ulambda(i)
      end do

      close(10)

C     COMPUTE CONSTRAINTS AND ITS NORM

      call ev_c(n,x,m,c,dat,idat)

      cnorm = 0.0d0
      do i = 1,m
          cnorm = max( cnorm, abs( c(i) ) )
      end do

C     COMPUTE THE GRADIENT OF THE LAGRANGIAN AND ITS NORM

      call ev_g(n,x,g,dat,idat)
      call ev_a(1,n,x,nz,a,acon,avar,dat,idat)

      do i = 1,nz
          g(avar(i)) = g(avar(i)) + a(i) * lambda(acon(i))
      end do

      do i = 1,ln
          g(lvar(i)) = g(lvar(i)) - llambda(i)
      end do

      do i = 1,un
          g(uvar(i)) = g(uvar(i)) + ulambda(i)
      end do

      gnorm = 0.0d0
      do i = 1,n
          g(i) = max( l(i), min( x(i) - g(i), u(i) ) ) - x(i)
          gnorm = max( gnorm, abs( g(i) ) )
      end do

C     COMPUTE COMPLEMENTARITY

      snorm = 0.0d0
      do i = 1,ln
          snorm = max(snorm, abs(( x(lvar(i)) - lval(i) ) * llambda(i)))
      end do
      do i = 1,un
          snorm = max(snorm, abs(( uval(i) - x(uvar(i)) ) * ulambda(i)))
      end do

C     WRITE STATISTICS

      open(10,file='ipopt-tabline.out')

      write(10,9050) time,inform,n,m,iter,f,cnorm,gnorm,snorm

      close(10)

C     NON-EXECUTABLE STATEMENTS

 9000 format(/,'FINAL POINT:',//,2X,'INDEX',16X,'X(INDEX)')
 9010 format(  I7,1P,D24.16)
 9020 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS',
     +       //,2X,'INDEX',11X,'LAMBDA(INDEX)')
 9030 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS FOR ',
     +         'THE LOWER BOUNDS: ',//,2X,'INDEX',5X,
     +         'LAMBDA-LOWER(INDEX)')
 9040 format(/,'FINAL ESTIMATION OF THE LAGRANGE MULTIPLIERS FOR ',
     +         'THE UPPER BOUNDS: ',//,2X,'INDEX',5X,
     +         'LAMBDA-UPPER(INDEX)')
 9050 format(F8.2,1X,I3,1X,I6,1X,I6,1X,I7,1P,D24.16,1X,1P,D7.1,1X,
     +       1P,D7.1,1X,1P,D7.1)

      stop

      end

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

      subroutine ev_f(n,x,f,dat,idat)

      implicit none

C     SCALAR ARGUMENTS
      integer n
      double precision f

C     ARRAY ARGUMENTS
      integer idat(*)
      double precision dat(*),x(n)

C     LOCAL SCALARS
      integer flag

      call evalfs(n,x,f,flag)

      if ( flag .ne. 0 ) then
          write(*,*) 'There was an error in evalfs subroutine.'
          stop
      end if

      end

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

      subroutine ev_g(n,x,g,dat,idat)

      implicit none

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      integer idat(*)
      double precision dat(*),g(n),x(n)

C     LOCAL SCALARS
      integer flag

      call evalgs(n,x,g,flag)

      if ( flag .ne. 0 ) then
          write(*,*) 'There was an error in evalgs subroutine.'
          stop
      end if

      end

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

      subroutine ev_c(n,x,m,c,dat,idat)

      implicit none

C     SCALAR ARGUMENTS
      integer m,n

C     ARRAY ARGUMENTS
      integer idat(*)
      double precision c(m),dat(*),x(n)

C     LOCAL SCALARS
      integer flag,i

      do i = 1,m
          call evalcs(n,x,i,c(i),flag)

          if ( flag .ne. 0 ) then
              write(*,*) 'There was an error in evalcs subroutine.'
              stop
          end if
      end do

      end

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

      subroutine ev_a(task,n,x,nz,a,acon,avar,dat,idat)

      implicit none

C     SCALAR ARGUMENTS
      integer n,nz,task

C     ARRAY ARGUMENTS
      integer acon(*),avar(*),idat(*)
      double precision a(*),dat(*),x(n)

C     LOCAL SCALARS
      integer flag,i,j,m,nztmp

      m = idat(1)

      if ( task .eq. 0 ) then
          nz = 0

          do j = 1,m
              call evaljact(n,x,j,task,avar,a,nztmp,flag)

              if ( flag .ne. 0 ) then
                  write(*,*) 
     +            'There was an error in evaljact subroutine.'
                  stop
              end if

              nz = nz + nztmp
          end do

      else
          nz = 0

          do j = 1,m
              call evaljact(n,x,j,task,avar(nz+1),a(nz+1),nztmp,flag)

              if ( flag .ne. 0 ) then
                  write(*,*) 
     +            'There was an error in evaljact subroutine.'
                  stop
              end if

              do i = 1,nztmp
                  acon(nz+i) = j
              end do

              nz = nz + nztmp
          end do
      end if

      end

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

      subroutine ev_h(task,n,x,m,lam,nnzh,hess,irnh,icnh,dat,idat)

      implicit none

C     SCALAR ARGUMENTS
      integer m,n,nnzh,task

C     ARRAY ARGUMENTS
      integer icnh(*),idat(*),irnh(*)
      double precision dat(*),hess(*),lam(m),x(n)

C     LOCAL SCALARS
      integer flag

      call evalhlt(n,x,m,lam,task,irnh,icnh,hess,nnzh,flag)

      if ( flag .ne. 0 ) then
          write(*,*) 'There was an error in evalhlt subroutine.'
          stop
      end if

      end

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

      subroutine inipt(n,x,l,u,m,lambda,rho,equatn,linear)

      implicit none

C     SCALAR ARGUMENTS
      integer m,n

C     ARRAY ARGUMENTS
      logical equatn(*),linear(*)
      double precision l(*),lambda(*),rho(*),u(*),x(*)

C     PARAMETERS
      integer hnnzmax,jcnnzmax,mmax
      parameter ( hnnzmax   = 5000000 )
      parameter ( jcnnzmax  = 5000000 )
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      logical cspace
      integer fhlnnz

C     COMMON ARRAYS
      integer fjcnnz(mmax)

C     LOCAL SCALARS
      integer flag,i,jcnnz,hlnnz

C     LOCAL ARRAYS
      integer jcvar(jcnnzmax),hllin(hnnzmax),hlcol(hnnzmax)
      double precision jcval(jcnnzmax),hlval(hnnzmax)

C     COMMON BLOCKS
      common /fixdim/ fjcnnz,fhlnnz,cspace
      save   /fixdim/

      call inips(n,x,l,u,m,lambda,rho,equatn,linear)

C     Some solvers requires the number of non null elements in the
C     Jacobian of the constraints and in the Hessian of the Lagrangian.
C     We compute those figures here and return them when required.

      cspace = .true.

      if ( cspace ) then

C         Jacobian of the constraints
          do i = 1,m
              call evaljacs(n,x,i,jcvar,jcval,jcnnz,flag)
              if ( flag .ne. 0 ) then
                  return
              end if

              fjcnnz(i) = jcnnz
          end do

C         Hessian of the Lagrangian
          call evalhls(n,x,m,lambda,hllin,hlcol,hlval,hlnnz,flag)
          if ( flag .ne. 0 ) then
              return
          end if

          fhlnnz = hlnnz

      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evaljact(n,x,ind,task,jcvar,jcval,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,jcnnz,task

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

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      logical cspace
      integer fhlnnz

C     COMMON ARRAYS
      integer fjcnnz(mmax)

C     COMMON BLOCKS
      common /fixdim/ fjcnnz,fhlnnz,cspace
      save   /fixdim/

      if ( task .eq. 0 ) then
          flag = 0

          if ( .not. cspace ) then
              write(*,100)
              flag = - 1
          end if

          jcnnz = fjcnnz(ind)
      else
          call evaljacs(n,x,ind,jcvar,jcval,jcnnz,flag)
      end if

C     Non-executable statements

 100  format(/,1X,'ERROR: Subroutine evaljact was called with task=0 ',
     +         'but the required dimension has not',/,1X,'been ',
     +         'computed. To fix it, set cspace equal to TRUE in ',
     +         'inipt subroutine compile and run again.')

      end

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

      subroutine evalhlt(n,x,m,lambda,task,hllin,hlcol,hlval,hlnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,hlnnz,m,n,task

C     ARRAY ARGUMENTS
      integer hllin(*),hlcol(*)
      double precision hlval(*),lambda(m),x(n)

C     PARAMETERS
      integer mmax
      parameter ( mmax      =  500000 )

C     COMMON SCALARS
      logical cspace
      integer fhlnnz

C     COMMON ARRAYS
      integer fjcnnz(mmax)

C     COMMON BLOCKS
      common /fixdim/ fjcnnz,fhlnnz,cspace
      save   /fixdim/

      if ( task .eq. 0 ) then
          flag = 0

          if ( .not. cspace ) then
              write(*,100)
              flag = - 1
          end if

          hlnnz = fhlnnz
      else
          call evalhls(n,x,m,lambda,hllin,hlcol,hlval,hlnnz,flag)
      end if

C     Non-executable statements

 100  format(/,1X,'ERROR: Subroutine evalhlt was called with task=0 ',
     +         'but the required dimension has not',/,1X,'been ',
     +         'computed. To fix it, set cspace equal to TRUE in ',
     +         'inipt subroutine compile and run again.')

      end

