!======================================================================
!
!     test program to see performance of IPOPT in case in using it 
!     as solver for non-linear equations.
!
!     this is achieved by setting min_x f(x) = 0 and setting the 
!     nonlinear equations as equality constraints.
!
!     min   0
!     s.t.  (x1-1)^2                   =  0
!
!
!  =====================================================================

      program example

 
      implicit none

!     include the Ipopt return codes

      include 'IpReturnCodes_SimonS.inc'

!     Size of the problem (number of variables and equality constraints)

      integer :: N, M, NELE_JAC, NELE_HESS, IDX_STY
      parameter(N = 1, M = 1, NELE_JAC = 1, NELE_HESS = 1)
      parameter(IDX_STY = 1) !index style = fortran style

!     Space for multipliers and constraints

      double precision LAM(M)
      double precision G(M)

!     Vector of variables

      double precision X(N)

!     Vector of lower and upper bounds

      double precision X_L(N), X_U(N), Z_L(N), Z_U(N)
      double precision G_L(M), G_U(M)

!      Private data for evaluation routines
!      This could be used to pass double precision and integer arrays untouched
!      to the evaluation subroutines EVAL_*

      double precision DAT(2)
      integer IDAT(1)

!     Place for storing the Ipopt Problem Handle

!C     for 32 bit platforms
!      integer IPROBLEM
!      integer IPCREATE
!     for 64 bit platforms:
      integer*8 IPROBLEM
      integer*8 IPCREATE

      integer IERR
      integer IPSOLVE, IPADDSTROPTION
      integer IPADDNUMOPTION, IPADDINTOPTION
      integer IPOPENOUTPUTFILE

      double precision F
      integer i

!     The following are the Fortran routines for computing the model
!     functions and their derivatives - their code can be found further
!     down in this file.

      external EV_F, EV_G, EV_GRAD_F, EV_JAC_G, EV_HESS
!C
!C     The next is an optional callback method.  It is called once per
!CC     iteration.
! 
!       external ITER_CB
! 
!      Set initial point and bounds:

      data X   / 0.d0/
      data X_L / -2d0 /
      data X_U / 2d0 /

!     Set bounds for the constraints

      data G_L / 0.0d0/
      data G_U / 0.0d0 /

!     First create a handle for the Ipopt problem (and read the options
!     file)

      IPROBLEM = IPCREATE(N, X_L, X_U, M, G_L, G_U, NELE_JAC, NELE_HESS, &
     &     IDX_STY, EV_F, EV_G, EV_GRAD_F, EV_JAC_G, EV_HESS)
      if (IPROBLEM.eq.0) then
         write(*,*) 'Error creating an Ipopt Problem handle.'
         stop
      endif

!     Open an output file

      IERR = IPOPENOUTPUTFILE(IPROBLEM, 'IPOPT.OUT', 5)
      if (IERR.ne.0 ) then
         write(*,*) 'Error opening the Ipopt output file.'
         goto 9000
      endif

!      Note: The following options are only examples, they might not be
!            suitable for your optimization problem.
! 
!      Set a string option

      IERR = IPADDSTROPTION(IPROBLEM, 'mu_strategy', 'adaptive')
      if (IERR.ne.0 ) goto 9990

!     Set an integer option

      IERR = IPADDINTOPTION(IPROBLEM, 'max_iter', 3000)
      if (IERR.ne.0 ) goto 9990

!     Set a double precision option

      IERR = IPADDNUMOPTION(IPROBLEM, 'tol', 1.d-7)
      if (IERR.ne.0 ) goto 9990

! C
! C     Set a callback function to give you control once per iteration.
! C     You can use it if you want to generate some output, or to stop
! C     the optimization early.
! C
!       call IPSETCALLBACK(IPROBLEM, ITER_CB)
! 
! 
!      As a simple example, we pass the constants in the constraints to
!      the EVAL_C routine via the "private" DAT array.

      DAT(1) = 0.d0
      DAT(2) = 0.d0

!     Call optimization routine

      IERR = IPSOLVE(IPROBLEM, X, G, F, LAM, Z_L, Z_U, IDAT, DAT)

!     Output:

      if( IERR.eq.IP_SOLVE_SUCCEEDED ) then
         write(*,*)
         write(*,*) 'The solution was found.'
         write(*,*)
         write(*,*) 'The final value of the objective function is ',F
         write(*,*)
         write(*,*) 'The optimal values of X are:'
         write(*,*)
         do i = 1, N
            write(*,*) 'X  (',i,') = ',X(i)
         enddo
         write(*,*)
         write(*,*) 'The multipliers for the lower bounds are:'
         write(*,*)
         do i = 1, N
            write(*,*) 'Z_L(',i,') = ',Z_L(i)
         enddo
         write(*,*)
         write(*,*) 'The multipliers for the upper bounds are:'
         write(*,*)
         do i = 1, N
            write(*,*) 'Z_U(',i,') = ',Z_U(i)
         enddo
         write(*,*)
         write(*,*) 'The multipliers for the equality constraints are:'
         write(*,*)
         do i = 1, M
            write(*,*) 'LAM(',i,') = ',LAM(i)
         enddo
         write(*,*)
      else
         write(*,*)
         write(*,*) 'An error occoured.'
         write(*,*) 'The error code is ',IERR
         write(*,*)
      endif

 9000 continue

!     Clean up

      call IPFREE(IPROBLEM)
      stop

 9990 continue
      write(*,*) 'Error setting an option'
      goto 9000
      end

!=============================================================================
! 
!                     Computation of objective function
! 
!=============================================================================

      subroutine EV_F(N, X, NEW_X, F, IDAT, DAT, IERR)
      implicit none
      integer N, NEW_X
      double precision F, X(N)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR
      F = 0.0
      IERR = 0
      return
      end

!  =============================================================================
! 
!                 Computation of gradient of objective function
! 
!  =============================================================================

      subroutine EV_GRAD_F(N, X, NEW_X, GRAD, IDAT, DAT, IERR)
      implicit none
      integer N, NEW_X
      double precision GRAD(N), X(N)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR
      GRAD(1) = 0.0
      IERR = 0
      return
      end

!  =============================================================================
! 
!                      Computation of equality constraints
! 
!  =============================================================================

      subroutine EV_G(N, X, NEW_X, M, G, IDAT, DAT, IERR)
      implicit none
      integer N, NEW_X, M
      double precision G(M), X(N)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR
      G(1) = (X(1)-1.)**2.    !
      IERR = 0
      return
      end

!  =============================================================================
! 
!                 Computation of Jacobian of equality constraints
! 
!  =============================================================================

      subroutine EV_JAC_G(TASK, N, X, NEW_X, M, NZ, ACON, AVAR, A, &
     &     IDAT, DAT, IERR)
      integer TASK, N, NEW_X, M, NZ
      double precision X(N), A(NZ)
      integer ACON(NZ), AVAR(NZ), I
      double precision DAT(*)
      integer IDAT(*)
      integer IERR

!     structure of Jacobian:
!
      integer AVAR1(1), ACON1(1)
      data  AVAR1 /1/
      data  ACON1 /1/
      save  AVAR1, ACON1

      if( TASK.eq.0 ) then
        do I = 1, 1!8
          AVAR(I) = AVAR1(I)
          ACON(I) = ACON1(I)
        enddo
      else
        A(1) = 2.*X(1)   
      endif
      IERR = 0
      return
      end

!  =============================================================================
! 
!                 Computation of Hessian of Lagrangian
! 
!  =============================================================================

      subroutine EV_HESS(TASK, N, X, NEW_X, OBJFACT, M, LAM, NEW_LAM, &
     &     NNZH, IRNH, ICNH, HESS, IDAT, DAT, IERR)
      implicit none
      integer TASK, N, NEW_X, M, NEW_LAM, NNZH, i
      double precision X(N), OBJFACT, LAM(M), HESS(NNZH)
      integer IRNH(NNZH), ICNH(NNZH)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR

!     structure of Hessian:

      integer IRNH1(1), ICNH1(1)
      data  IRNH1 /1/
      data  ICNH1 /1/
      save  IRNH1, ICNH1

      if( TASK.eq.0 ) then
         do i = 1, 1
            IRNH(i) = IRNH1(i)
            ICNH(i) = ICNH1(i)
         enddo
      else
         do i = 1, 1
            HESS(i) = LAM(1)*2d0
         enddo

      endif
      IERR = 0
      return
      end

!  =============================================================================
! 
!                    Callback method called once per iteration
! 
!  =============================================================================

      subroutine ITER_CB(ALG_MODE, ITER_COUNT,OBJVAL, INF_PR, INF_DU, &
     &     MU, DNORM, REGU_SIZE, ALPHA_DU, ALPHA_PR, LS_TRIAL, IDAT, &
     &     DAT, ISTOP)
      implicit none
      integer ALG_MODE, ITER_COUNT, LS_TRIAL
      double precision OBJVAL, INF_PR, INF_DU, MU, DNORM, REGU_SIZE
      double precision ALPHA_DU, ALPHA_PR
      double precision DAT(*)
      integer IDAT(*)
      integer ISTOP

!     You can put some output here

      write(*,*) 'Testing callback function in iteration ', ITER_COUNT

!     And set ISTOP to 1 if you want Ipopt to stop now.  Below is just a
!     simple example.

      if (INF_PR.le.1D-04) ISTOP = 1

      return
      end