program algpma

    use bgmalgp
    use iso_c_binding, only: c_ptr, c_loc, c_f_pointer
  
    implicit none
  
    type pdata_type
       integer :: counters(3) = 0
    end type pdata_type
  
    ! LOCAL SCALARS
    logical :: extallowed
    integer :: hnnzmax,ierr,istop,iter,maxit,n,nbds
    real(kind=8) :: eps,f,finish,ftarget,gpsupn,start
    type(pdata_type), target :: pdata
  
    ! LOCAL ARRAYS
    logical :: lind(2),uind(2)
    real(kind=8) :: g(2),lbnd(2),ubnd(2),x(2)

    n = 2 
    lind(1) = .true. 
    lbnd(1) = 0.0d0 
    uind(1) = .true. 
    ubnd(1) = 10.0d0 

    lind(2) = .false. 
    uind(2) = .false. 
  
    x(1:2) = 10.0d0 
    hnnzmax = 100000
      
    ftarget     =  - 1.0d+12
    eps         =    1.0d-08
    maxit       =      huge(1)
    extallowed  =     .true.
  
    nbds = count( lind(1:n) ) + count( uind(1:n) )

    call cpu_time(start)
      
    if ( nbds .eq. 0 ) then
       write(*,*) 'The problem is unconstrained.'
       call algpunc(evalf,evalg,evalh,hnnzmax,n,x,f,g,gpsupn, &
            ftarget,eps,maxit,extallowed,iter,ierr,istop,pdata=c_loc(pdata))
    else
       write(*,*) 'The problem is a bound constrained problem.'
       call algp(evalf,evalg,evalh,hnnzmax,n,x,lind, &
              lbnd,uind,ubnd,f,g,gpsupn,ftarget,eps,maxit,extallowed,iter, &
              ierr,istop,pdata=c_loc(pdata))
    end if
    
    call cpu_time(finish)
  
    write(*,*)
    write(*,*) 'Number of variables                                    = ',n
    write(*,*) 'Number of lower and upper bounds                       = ',count( lind(1:n) ) + count( uind(1:n) )
    write(*,*)
    write(*,*) '(REPORTED BY SOLVER) istop                             = ',istop
    write(*,*) '(REPORTED BY SOLVER) ierr                              = ',ierr
    write(*,*) '(REPORTED BY SOLVER) f                                 = ',f
    write(*,*) '(REPORTED BY SOLVER) gpsupn                            = ',gpsupn
    write(*,*) '(REPORTED BY SOLVER) Number of iterations              = ',iter
  
    write(*,*)
    write(*,*) '(COMPUTED BY CALLER) Number of function evaluations    = ',pdata%counters(1)
    write(*,*) '(COMPUTED BY CALLER) Number of gradient evaluations    = ',pdata%counters(2)
    write(*,*) '(COMPUTED BY CALLER) Number of Hessian  evaluations    = ',pdata%counters(3)
    write(*,*) '(COMPUTED BY CALLER) CPU time in seconds               = ',finish - start

  contains
    
  
  ! *****************************************************************
  ! *****************************************************************

  subroutine evalf(n,x,f,inform,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    integer, intent(inout) :: inform
    real(kind=8), intent(out) :: f
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)

    ! This routine must compute the objective function.
    
    ! LOCAL SCALARS
    type(pdata_type), pointer :: pdata
    
    call c_f_pointer(pdataptr,pdata)
    pdata%counters(1) = pdata%counters(1) + 1
    
    f = ( x(1) + 4.0d0 ) ** 4.0d0 + x(2) ** 2.0d0
    
    inform = 0
  end subroutine evalf

  ! *****************************************************************
  ! *****************************************************************

  subroutine evalg(n,x,g,inform,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    integer, intent(inout) :: inform
    type(c_ptr), optional, intent(in) :: pdataptr
  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: g(n)
    
    ! This routine must compute the gradient of the objective
    ! function.
    
    ! LOCAL SCALARS
    type(pdata_type), pointer :: pdata
    
    call c_f_pointer(pdataptr,pdata)
    pdata%counters(2) = pdata%counters(2) + 1
    
    g(1) = 4.0d0 * ( x(1) + 4.0d0 ) ** 3.0d0
    g(2) = 2.0d0 * x(2)

    inform = 0

  end subroutine evalg

  ! *****************************************************************
  ! *****************************************************************

  subroutine evalh(n,x,lim,hnnz,hrow,hcol,hval,inform,pdataptr)

    implicit none
  
      ! SCALAR ARGUMENTS
      integer, intent(in) :: lim,n
      integer, intent(inout) :: inform
      integer, intent(out) :: hnnz
      type(c_ptr), optional, intent(in) :: pdataptr
  
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
      integer, intent(out) :: hcol(lim),hrow(lim)
      real(kind=8), intent(out) :: hval(lim)
  
      ! LOCAL SCALARS
      type(pdata_type), pointer :: pdata
    
      call c_f_pointer(pdataptr, pdata)
      pdata%counters(3) = pdata%counters(3) + 1

      if ( lim .lt. 2 ) then 
        inform = -94
      end if

      hrow(1) = 1
      hcol(1) = 1
      hval(1) = 12*(x(1) + 4) ** 2

      hrow(2) = 2
      hcol(2) = 2
      hval(2) = 2
  
      hnnz = 2

      inform = 0
  end subroutine evalh

end program algpma
