program bmbunchma

  use mgh_dp
  
  implicit none

  ! LOCAL SCALARS
  integer :: allocstatus,fcnt,flag,inform,iter,maxfcnt,maxit,n,nprob,m,nfact,nprobini,ntimes
  real(kind=8) :: elapsed,eps,f,finish,gnorm,start,factor,ftarget
  
  ! LOCAL ARRAYS
  character(len=60) :: probname
  real(kind=8), dimension(:), allocatable :: x,g

  eps = 1.0d-08
  maxit = 500000
  maxfcnt = 10 * maxit
  
  open(unit=10,file='bmbunch-mgh-table.txt',access='append')
  write(10,*) 'Parameter: eps = ',eps
  write(10,*) 
  write(10,*) 'STOP =   0: xk is such that |g(xk)|<eps.'
  write(10,*) 'STOP =   1: xk is such that |g(xk)|<sqrt(eps) during  100 consecutive iterations.'
  write(10,*) 'STOP =   2: xk is such that |g(xk)|<eps^1/4   during 1000 consecutive iterations.'
  write(10,*) 'STOP =   3: xk is such that |g(xk)|<eps^1/8   during 5000 consecutive iterations.'
  write(10,*) 'STOP =   4: xk+s is such that s is the Newton step, xk+s did not give suff. decr., |s|<eps, and |g(xk+s)|<eps.'
  write(10,*) 'STOP =   5: xk+s is such that s is the Newton step, xk+s did not give sufficient decrease, and |s|<eps.'
  write(10,*) 'STOP =   6: xk is such that f(xk)<ftarget.'
  write(10,*) 'STOP =   7: xk+s is such that xk+s did not give sufficient decrease, but f(xk+s)<ftarget.'
  write(10,*) 'STOP =   8: xk=xk+1 and f(xk +/- hi ei) with hi = macheps max(1, |[xk]_i|) for i=1,n.'
  write(10,*) 'STOP =   9: xk and xk+1 are such that f(xk+1)=f(xk) during 10 consecutive iterations.'
  write(10,*) 'STOP =  10: xk is such that |x|>xlarge, with xlarge=1.0d+10.'
  write(10,*) 'STOP =  11: the maximum allowed number of iterations was exhausted (maxit = ',maxit,').'
  write(10,*) 'STOP =  12: the maximum allowed number of functional evaluations was exhausted (maxfcnt = ',maxfcnt,').'
  write(10,*) 'STOP = -90: xk is such that there was an error in the evaluation of the objective function.'
  write(10,*) 'STOP = -91: xk is such that there was an error in the evaluation of the gradient.'
  write(10,*) 'STOP = -92: xk is such that there was an error in the evaluation of the Hessian.'
  write(10,*)
  write(10,*) 'Each problem was run ntimes times, where ntimes is such that the CPU time for the ntimes runs is larger '
  write(10,*) 'not smaller than a second. Assuming that the time is measured with a precision of at least 0.01 seconds, the '
  write(10,*) 'time for the ntimes runs has an error smaller than 1%. The reported time is the average of the ntimes runs.'
  close(10)
  
  do nfact = 1,4
     factor = 10.0d0 ** ( nfact - 1 )
     if ( nfact .eq. 4 ) then
        factor = 1.0d0
     end if
     
     open(unit=10,file='bmbunch-mgh-table.txt',access='append')

     write(10,*)

     if ( nfact .eq. 1 ) then
        write(10,*) '****************'
        write(10,*) 'STARTING FROM x0'
        write(10,*) '****************'
     else if ( nfact .eq. 2 ) then
        write(10,*) '*******************'
        write(10,*) 'STARTING FROM 10 x0'
        write(10,*) '*******************'
     else if ( nfact .eq. 3 ) then
        write(10,*) '********************'
        write(10,*) 'STARTING FROM 100 x0'
        write(10,*) '********************'
     else if ( nfact .eq. 4 ) then
        write(10,*) '*************************************************'
        write(10,*) 'PROBLEMS WITH VARIABLE DIMENSION STARTING FROM x0'
        write(10,*) '*************************************************'
     end if

     write(10,*) '---------------------------------------------------------------------------------------------------------'
     write(10,*) 'Problem                             n    m          f          |g|        iter     fcnt         time STOP'
     write(10,*) '---------------------------------------------------------------------------------------------------------'

     close(10)

     nprobini = 1
     if ( nfact .eq. 4 ) then
        nprobini = 21
     end if

     do nprob = nprobini,35

        call mgh_set_problem(nprob,flag)
        if ( flag .ne. 0 ) then
           write(*,*) 'Error in mgh_set_problem'
           stop
        end if

        call mgh_get_name(probname)

        call mgh_get_dims(n,m)

        if ( nfact .eq. 4 ) then
           if ( nprob .eq. 21 ) then
              ! n even, m=n
              n = 1000
              m = n
           else if ( nprob .eq. 22 ) then
              ! n multiple of 4, m=n
              n = 1000
              m = n
           else if ( nprob .eq. 23 ) then
              ! n variable, m=n+1
              n = 1000
              m = n + 1
           else if ( nprob .eq. 24 ) then
              ! n variable, m=2n
              n = 1000
              m = 2 * n
           else if ( nprob .eq. 25 ) then
              ! n variable, m=n+2
              n = 1000
              m = n + 2
           else if ( nprob .ge. 26 .and. nprob .le. 31 ) then
              ! n variable, m=n
              n = 1000
              m = n
           else if ( nprob .ge. 32 .and. nprob .le. 34 ) then
              ! n variable, m>=n
              n = 1000
              m = n
           else if ( nprob .eq. 35 ) then
              ! n variable, m>=n
              n = 400
              m = n
           end if
           
           if ( nprob .ge. 20 .and. nprob .le. 35 ) then
              call mgh_set_dims(n,m,flag)
           end if
        end if
        
        ftarget = - 1.0d+10

        if ( nprob .eq. 10 ) then
           ftarget = 87.9459d0
           write(*,*) 'Valor de ftarget = ',ftarget
        
        else if ( nprob .eq. 33 .and. n .eq. 1000 ) then
           ftarget = real( m * ( m - 1 ), kind=8 ) / ( 2 * ( 2 * m + 1 ) )
           ftarget = ftarget + 1.0d-15 * abs( ftarget )
           write(*,*) 'Valor de ftarget = ',ftarget
        
        else if ( nprob .eq. 34 .and. n .eq. 1000 ) then
           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
     
        write(*,*) 
        write(*,*) 
        write(*,*) '=========================================='
        write(*,*) '=========================================='
        write(*,*) 'Problem name: ',probname
        write(*,*) 'Number of variables n = ',n,' ( m = ',m,')'
        write(*,*) '=========================================='
        write(*,*) '=========================================='
        write(*,*) 
        write(*,*) 
     
        allocate(x(n),g(n),stat=allocstatus )
        if ( allocstatus .ne. 0 ) then
           write(*,*) 'Error in allocate!'
           stop
        end if

        ntimes = 1
        
        call cpu_time(start)

10      continue
        
        call mgh_get_x0(x,factor)

        call CRNewton(n,x,f,g,gnorm,eps,ftarget,maxit,maxfcnt,iter,fcnt,inform)
           
        call cpu_time(finish)
           
        elapsed = finish - start

        if ( elapsed .le. 1.0d0 ) then
           ntimes = ntimes + 1
           go to 10
        end if
        
        elapsed = elapsed / ntimes 
        
        deallocate(x,g,stat=allocstatus )
        if ( allocstatus .ne. 0 ) then
           write(*,*) 'Error in deallocate!'
           stop
        end if
        
        open(unit=10,file='bmbunch-mgh-table.txt',access='append')
        write(10,100) nprob,probname,n,m,f,gnorm,iter,fcnt,elapsed,inform
        close(10)
        
        write(*,*) 'problem = ',nprob
        write(*,*) 'Inform                                    = ',inform
        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 (average of ',ntimes,' runs) = ',elapsed
        
     end do
  end do

  stop

100 format(I2,1X,A30,1X,I4,1X,I4,3X,1PD17.10,1X,1PD7.1,I8,1X,I8,1X,0P,F12.6,1X,I3)
  
end program bmbunchma

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

subroutine evalf(n,x,f,flag)

  use mgh_dp
  
  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 )

  flag = 0
  
  call mgh_evalf(x,f,flag)

  ! This is true if f is NaN
  if ( .not. ( f .eq. f ) ) then
     write(*,*) 'Point at which f was being evaluated, x = ',x(1:n)
     write(*,*) 'Computed value of f = ',f
     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)

  use mgh_dp
  
  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 )

  flag = 0
  
  call mgh_evalg(x,g,flag)

  ! This is true if there is a NaN
  if ( any( .not. ( g(1:n) .eq. g(1:n) ) ) ) then
     write(*,*) 'Point at which g was being evaluated, x = ',x(1:n)
     write(*,*) 'Computed value of g = ',g(1:n)
     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)

  use mgh_dp
  
  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)

  ! Lower triangle only

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

  ! LOCAL SCALARS
  integer :: j
  
  ! Lower triangle of Hessian, upper triangle with zeros.
  
  flag = 0
  
  call mgh_evalh(x,H,flag)

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

