program bobyqama

  use bmproblem, only: pdata_type,iniprob,endprob,evalls,evalvero,frompvtox2
  use iso_c_binding, only: c_ptr,c_loc

  implicit none

  ! LOCAL SCALARS
  integer :: allocstat,fcnt,i,iprint,k,ksig,maxfun,n,npt
  real(kind=8) :: bdsvio,fbobyqa,feps,finish,fvero,fvero2,start,rhobeg,rhoend,seed,sigini,sigend,efe
  type(pdata_type), target :: pdata

  ! LOCAL ARRAYS
  real(kind=8), allocatable :: l(:),u(:),x(:),w(:),p(:),v(:),xtest(:)

  ! EXTERNAL FUNCTIONS
  real(kind=8), external :: drand

  fcnt = 0
  
  call iniprob(c_loc(pdata))

  n = 2 * pdata%inodes + 2

  allocate(l(n),u(n),x(n+1),stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in allocate.'
     stop
  end if

  l(1:pdata%inodes) = 0.0d0
  u(1:pdata%inodes) = 1.0d0
  
  l(pdata%inodes+1:2*pdata%inodes+2) = 0.0d0
  u(pdata%inodes+1:2*pdata%inodes+2) = huge( 1.0d0 )

  seed = 123456.0d0

  do i = 1,pdata%inodes
     x(i) = drand(seed)
  end do
  
  x(pdata%inodes+1:2*pdata%inodes+2) = 0.0d0

  !================
  !================
  allocate(p(0:pdata%inodes+1),v(0:pdata%inodes+1),xtest(pdata%ndiv),stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in allocate.'
     stop
  end if

  efe = 3.3663306720355467d-002
  p(0:pdata%inodes+1) = (/ (real(i)/(pdata%inodes+1),i=0,pdata%inodes+1) /)
  v(0:pdata%inodes+1) = efe + 10.0d0 * efe * &
       exp( -1.0d-04 * ( pdata%dx * ( 1 + p(0:pdata%inodes+1) * ( pdata%ndiv - 1 ) - 0.5d0 * pdata%ndiv ) ) ** 2 )

  x(1:n) = (/ p(1:pdata%inodes), v(0:pdata%inodes+1) /)

  call frompvtox2(n,x,pdata%ndiv,xtest)

!!$  do i = 1,pdata%ndiv
!!$     write(*,*) i,p(i-1),v(i-1),xtest(i),pdata%xsol(i),abs(xtest(i)-pdata%xsol(i))
!!$  end do
  
  call evalls(n,x,fbobyqa,c_loc(pdata))
  write(*,*) 'Value of f at initial guess = ',fbobyqa
  !================
  !================
  
  iprint = 2

  npt = 2 * n + 1
! maxfun = 400 * ( n + 1 )
  maxfun = huge( 1 )

  allocate(w((npt+5)*(npt+n)+3*n*(n+5)/2),stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in allocate.'
     stop
  end if

  call cpu_time(start)

  rhobeg = 0.1d0
  rhoend = 1.0d-20

  feps = 0.0d0
  call bobyqa(n,npt,x,l,u,rhobeg,rhoend,iprint,maxfun,w,fbobyqa,feps,calfun)

  call cpu_time(finish)

  bdsvio = max( 0.0d0, maxval( l(1:n) - x(1:n) ), maxval( x(1:n) - u(1:n) ) )
  
  write(*,*) 'Solution: ',x(1:n)
  write(*,*) 'Solution (two first components): ',x(1:min(2,n))
  write(*,*) 'Box-violation sup-norm = ',bdsvio
  write(*,*) 'Functional evaluations = ',fcnt
  write(*,*) 'CPU time in seconds = ',finish - start

  sigini = 0.0d0
  sigend = 0.4d0
  ksig = 101
    
  do k = 1,ksig
     x(n+1) = sigini + ( ( real( k - 1 ) / ( ksig - 1 ) ) ) * ( sigend - sigini )
     call evalvero(n+1,x,fvero,fvero2,c_loc(pdata))
     write(*,*) k,x(n+1),fvero,fvero2
  end do
  
  deallocate(l,u,x,w,p,v,stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in deallocate.'
     stop
  end if

  call endprob(c_loc(pdata))
  
  stop
  
contains

  subroutine calfun(n,x,f)

    implicit none
    
    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    real(kind=8), intent(out) :: f

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

    call evalls(n,x,f,c_loc(pdata))
    fcnt = fcnt + 1

!!$    write(*,*) 'In calfun, x = ',x(1:n),' f = ',f
    
  end subroutine calfun
  
end program bobyqama
