program bmbunchma

  implicit none

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

  eps = 1.0d-08
  ftarget = - 1.0d+99
  maxit = 1000
  maxfcnt = 10 * maxit
  
  call getn(n)
  call gethnnz(hnnz)
  
  write(*,*) 
  write(*,*) 
  write(*,*) '=========================================='
  write(*,*) '=========================================='
  write(*,*) 'Problem name: toyprob'
  write(*,*) 'Number of variables n = ',n
  write(*,*) '=========================================='
  write(*,*) '=========================================='
  write(*,*) 
  write(*,*) 
     
  allocate(x(n),g(n),hrow(hnnz),hcol(hnnz),stat=allocstatus )
  if ( allocstatus .ne. 0 ) then
     write(*,*) 'Error in allocate!'
     stop
  end if

  call getx0(n,x)
  call gethpat(hnnz,hrow,hcol)
  
  call cpu_time(start)
  call CRNewton(n,x,f,g,gnorm,eps,ftarget,maxit,maxfcnt,hnnz,hrow,hcol,iter,fcnt,inform)
  call cpu_time(finish)
        
  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            = ',finish-start
  
  deallocate(x,g,hrow,hcol,stat=allocstatus )
  if ( allocstatus .ne. 0 ) then
     write(*,*) 'Error in deallocate!'
     stop
  end if
  
  stop

end program bmbunchma

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

subroutine getn(n)

  implicit none

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

  integer :: n_saved
  common /pdata/ n_saved
  
  write(*,*) 'Insert n (even): '
  read(*,*) n

  n_saved = n
  
end subroutine getn

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

subroutine getx0(n,x)

  implicit none

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

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

  ! LOCAL SCALARS
  integer :: i

  do i = 1,n,2
     x(i) = - 1.2d0
     x(i+1) = 1.0d0
  end do

end subroutine getx0

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

subroutine gethnnz(hnnz)

  implicit none

  !SCALAR ARGUMENTS
  integer, intent(out) :: hnnz

  integer :: n_saved
  common /pdata/ n_saved
  
  hnnz = 3 * n_saved / 2

end subroutine gethnnz

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

subroutine gethpat(hnnz,hrow,hcol)

  implicit none

  !SCALAR ARGUMENTS
  integer, intent(in) :: hnnz

  ! ARRAY ARGUMENTS
  integer, intent(out) :: hrow(hnnz),hcol(hnnz)

  integer :: n_saved
  common /pdata/ n_saved
  
  ! LOCAL SCALARS
  integer :: hnnz_local,i
  
  ! Sparsity pattern of the Hessian's lower triangle

  hnnz_local = 0
  
  do i = 1,n_saved,2
     hnnz_local = hnnz_local + 1

     hrow(hnnz_local) = i
     hcol(hnnz_local) = i
     
     hnnz_local = hnnz_local + 1

     hrow(hnnz_local) = i + 1
     hcol(hnnz_local) = i
     
     hnnz_local = hnnz_local + 1

     hrow(hnnz_local) = i + 1
     hcol(hnnz_local) = i + 1
  end do

  if ( hnnz_local .ne. hnnz ) then
     write(*,*) 'There is an error in gethpat!'
     stop
  end if
  
end subroutine gethpat

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

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)

  ! LOCAL SCALARS
  integer :: i
  
  flag = 0
  
  f = 0.0d0
  do i = 1,n,2
     f = f + ( 1.0d0 - x(i) ) ** 2 + 100.0d0  * ( x(i+1) - x(i) ** 2 ) ** 2
  end do
  
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)

  ! LOCAL SCALARS
  integer :: i
  
  flag = 0
  
  do i = 1,n,2
     g(i) = - 2.0d0 * ( 1.0d0 - x(i) ) - 400.0d0 * x(i) * ( x(i+1) - x(i) ** 2 )
     g(i+1) = 200.0d0 * ( x(i+1) - x(i) ** 2 )
  end do

end subroutine evalg

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

subroutine evalh(n,x,hnnz,hrow,hcol,hval,flag)

  implicit none

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

  ! ARRAY ARGUMENTS
  integer, intent(in) :: hcol(hnnz),hrow(hnnz)
  real(kind=8), intent(in) :: x(n)
  real(kind=8), intent(out) :: hval(hnnz)

  ! LOCAL SCALARS
  integer :: hnnz_local,i
  
  ! Sparse Hessian's lower triangle
  
  flag = 0

  hnnz_local = 0
  
  do i = 1,n,2
     hnnz_local = hnnz_local + 1
     
     if ( hrow(hnnz_local) .ne. i .or. hcol(hnnz_local) .ne. i ) then
        write(*,*) 'There is an error in evalh!'
        stop
     end if

     hval(hnnz_local) = 1200.0d0 * x(i) ** 2 - 400.0d0 * x(i+1) + 2.0d0
     
     hnnz_local = hnnz_local + 1

     if ( hrow(hnnz_local) .ne. i + 1 .or. hcol(hnnz_local) .ne. i ) then
        write(*,*) 'There is an error in evalh!'
        stop
     end if

     hval(hnnz_local) = - 400.0d0 * x(i)
     
     hnnz_local = hnnz_local + 1

     if ( hrow(hnnz_local) .ne. i + 1 .or. hcol(hnnz_local) .ne. i + 1 ) then
        write(*,*) 'There is an error in evalh!'
        stop
     end if

     hval(hnnz_local) = 200.0d0
  end do

  if ( hnnz_local .ne. hnnz ) then
     write(*,*) 'There is an error in evalh!'
     stop
  end if

end subroutine evalh

