program bmbunchma

  implicit none

  ! LOCAL SCALARS
  logical :: pnewt
  integer :: allocstatus,fcnt,i,inform,iter,maxfcnt,maxit,n
  real(kind=8) :: eps,f,finish,gnorm,start,ftarget
  
  ! LOCAL ARRAYS
  real(kind=8), dimension(:), allocatable :: x,g

  eps = 1.0d-08
  ftarget = - 1.0d+99
  maxit = 1000
  maxfcnt = 10 * maxit
  
  n = 3
  
  write(*,*) 
  write(*,*) 
  write(*,*) '=========================================='
  write(*,*) '=========================================='
  write(*,*) 'Problem name: x^2 + 2 y^3 + x^2 y^2 + z^4'
  write(*,*) 'Number of variables n = ',n
  write(*,*) '=========================================='
  write(*,*) '=========================================='
  write(*,*) 
  write(*,*) 
     
  allocate(x(n),g(n),stat=allocstatus )
  if ( allocstatus .ne. 0 ) then
     write(*,*) 'Error in allocate!'
     stop
  end if

  x(1:n) = (/ (i,i=1,n) /)
  
  call cpu_time(start)
  call CRNewton(n,x,f,g,gnorm,eps,ftarget,maxit,maxfcnt,iter,fcnt,pnewt,inform)
  call cpu_time(finish)
        
  write(*,*) 'Inform                                 = ',inform
  write(*,*) 'There were pure Newton iterations only = ',pnewt
  write(*,*) 'Final objective functional value       = ',f
  write(*,*) 'Final sup-norm of the gradient         = ',gnorm
  write(*,*) 'Total number of iterations             = ',iter
  write(*,*) 'total number of functional evaluations = ',fcnt
  write(*,*) 'Elapsed CPU time in seconds            = ',finish-start
  
  deallocate(x,g,stat=allocstatus )
  if ( allocstatus .ne. 0 ) then
     write(*,*) 'Error in deallocate!'
     stop
  end if
  
  stop

end program bmbunchma

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

subroutine evalf(n,x,f,flag)

  implicit none

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

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

  flag = 0
  
  f = x(1) ** 2 + 2.0d0 * x(2) ** 3 + x(1) ** 2 * x(2) ** 2 + x(3) ** 4
  
end subroutine evalf

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

subroutine evalg(n,x,g,flag)

  implicit none

  !SCALAR ARGUMENTS
  integer, intent(in) :: n
  integer, intent(out) :: flag

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

  flag = 0
  
  g(1) = 2.0d0 * x(1) + 2.0d0 * x(1) * x(2) ** 2
  g(2) = 6.0d0 * x(2) ** 2 + 2.0d0 * x(2) * x(1) ** 2
  g(3) = 4.0d0 * x(3) ** 3
  
end subroutine evalg

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

subroutine evalh(n,x,H,flag)

  implicit none

  !SCALAR ARGUMENTS
  integer, intent(in) :: n
  integer, intent(out) :: flag

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

  ! Sparse Hessian's lower triangle
  
  flag = 0
  
  H(1:n,1:n) = 0.0d0

  H(1,1) = 2.0d0 + 2.0d0 * x(2) ** 2
  H(2,1) = 4.0d0 * x(1) * x(2)
  H(2,2) = 12.0d0 * x(2) + 2.0d0 * x(1) ** 2
  H(3,3) = 12.0d0 * x(3) ** 2
  
end subroutine evalh