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

program bmbunchma

  implicit none

  ! LOCAL SCALARS
  logical :: pnewt
  integer :: allocstatus,fcnt,inform,iter,m,maxfcnt,maxit,n,status
  real(kind=8) :: eps,f,finish,ftarget,gnorm,start

  ! LOCAL ARRAYS
  character(len=10) :: pname
  character(len=10), allocatable :: xnames(:)
  real(kind=8), allocatable :: g(:),l(:),u(:),x(:)

  ! COMMON BLOCK
  integer :: hnnzmax
  integer, pointer :: hcol(:),hrow(:)
  real(kind=8), pointer :: hval(:)
  common /sparsehess/ hcol,hrow,hval,hnnzmax
  
  open(10,file='OUTSDIF.d',form='formatted',status='old')

  rewind 10

  call cutest_udimen(status,10,n)
  if ( status .ne. 0 ) then
     write(*,*) 'cutest_udimen: status not equal to zero. status = ',status
     stop
  end if

  allocate(g(n),x(n),l(n),u(n),xnames(n),stat=allocstatus )
  if ( allocstatus .ne. 0 ) then
     write(*,*) 'Error in allocate!'
     stop
  end if

  call cutest_usetup(status,10,6,4,n,x,l,u)
  if ( status .ne. 0 ) then
     write(*,*) 'cutest_usetup: status not equal to zero. status = ',status
     stop
  end if

  call cutest_unames(status,n,pname,xnames)

  call cutest_udimsh(status,hnnzmax)
  if ( status .ne. 0 ) then
     write(*,*) 'cutest_udimsh: status not equal to zero. status = ',status
     stop
  end if

  allocate(hrow(hnnzmax),hcol(hnnzmax),hval(hnnzmax),stat=allocstatus )
  if ( allocstatus .ne. 0 ) then
     write(*,*) 'Error in allocate!'
     stop
  end if

  close(10)

  write(*,*) 
  write(*,*) 
  write(*,*) '=========================================='
  write(*,*) '=========================================='
  write(*,*) 'Problem name: ',pname
  write(*,*) 'Number of variables n = ',n
  write(*,*) '=========================================='
  write(*,*) '=========================================='
  write(*,*) 
  write(*,*) 
  
  eps = 1.0d-08
  ftarget = - 1.0d+10
  maxit = 50000
  maxfcnt = 10 * maxit
  
  if ( pname .eq. 'ARGLINC' ) then
     m = 1000
     ftarget = real( ( m ** 2 + 3 * m - 6 ), kind=8 ) / ( 2 * ( 2 * m - 3 ) )
     ftarget = ftarget + 1.0d-15 * abs( ftarget )
     write(*,*) 'Valor de ftarget = ',ftarget
  end if

  open(10,file='bmbunch-cutest-ran-out-of-time.txt')
  close(10)

  call cpu_time(start)
  call CRNewton(n,x,f,g,gnorm,eps,ftarget,maxit,maxfcnt,iter,fcnt,pnewt,inform)
  call cpu_time(finish)
        
  open(10,file='bmbunch-cutest-ran-out-of-time.txt')
  close(10,status='delete')

  open(10,file='bmbunch-cutest-table.txt',access='append')
  write(10,100) pname,n,f,gnorm,iter,fcnt,pnewt,finish-start,inform
  close(10)
        
  write(*,*) 'problem                                = ',pname
  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(g,x,l,u,xnames,hrow,hcol,hval,stat=allocstatus )
  if ( allocstatus .ne. 0 ) then
     write(*,*) 'Error in deallocate!'
     stop
  end if

  stop

100 format(A10,1X,I4,3X,1PD17.10,1X,1PD7.1,I6,1X,I6,1X,1L,1X,0P,F8.2,1X,I3)

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)

  ! SCALAR PARAMETERS
  real( kind=8), parameter :: lrgstreal = huge( 1.0d0 )

  ! LOCAL SCALARS
  integer :: status

  flag = 0

  call cutest_ufn(status,n,x,f)
  if ( status .ne. 0 ) flag = - 1

  ! This is true if f is NaN
  if ( .not. ( f .eq. f ) ) then
     flag = - 1
  end if
  
  ! This is true if f is a +/- Infty
  if ( .not. abs( f ) .le. lrgstreal ) then
     flag = 1
     return
  end if

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)

  ! SCALAR PARAMETERS
  real( kind=8), parameter :: lrgstreal = huge( 1.0d0 )

  ! LOCAL SCALARS
  integer :: status

  flag = 0

  call cutest_ugr(status,n,x,g)
  if ( status .ne. 0 ) flag = - 1

  ! This is true if there is a NaN
  if ( .not. all( g(1:n) .eq. g(1:n) ) ) then
     flag = - 1
  end if
  
  ! This is true if there is a +/- Infty
  if ( .not. all( abs( g(1:n) ) .le. lrgstreal ) ) then
     flag = 1
     return
  end if
  
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)

  ! COMMON BLOCK
  integer :: hnnzmax
  integer, pointer :: hcol(:),hrow(:)
  real(kind=8), pointer :: hval(:)
  common /sparsehess/ hcol,hrow,hval,hnnzmax
  
  ! SCALAR PARAMETERS
  real( kind=8), parameter :: lrgstreal = huge( 1.0d0 )

  ! LOCAL SCALARS
  integer :: hnnz,k,status,row,col

  ! Lower triangle of Hessian, upper triangle with zeros.
    
  flag = 0

  ! This routine computes the upper triangle and we want the lower one.
  call cutest_ush(status,n,x,hnnz,hnnzmax,hval,hrow,hcol)
  if ( status .ne. 0 ) then
     flag = - 1
     return
  end if

  H(1:n,1:n) = 0.0d0

  do k = 1,hnnz
     row = hcol(k)
     col = hrow(k)
     
     H(row,col) = H(row,col) + hval(k)
  end do

  ! This is true if there is a NaN
  if ( .not. all( H(1:n,1:n) .eq. H(1:n,1:n) ) ) then
     flag = - 1
  end if

  ! This is true if there is a +/- Infty
  if ( .not. all( abs( H(1:n,1:n) ) .le. lrgstreal ) ) then
     flag = 1
     return
  end if
  
end subroutine evalh
