program main

!*****************************************************************************80
!
!! nintlib_test() tests nintlib().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    11 September 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer, parameter :: test_num = 3

  real ( kind = rk8 ) a
  real ( kind = rk8 ) b
  integer dim_num
  integer, dimension ( test_num ) :: dim_num_test = (/ 2, 3, 4 /)
  real ( kind = rk8 ), external :: f1dn
  real ( kind = rk8 ), external :: fbdn
  real ( kind = rk8 ), external :: fedn
  real ( kind = rk8 ), external :: fxdn
  real ( kind = rk8 ), external :: fx2dn
  real ( kind = rk8 ), external :: fx3dn
  integer test

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'NINTLIB_TEST'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Test the NINTLIB library.'

  a = 0.0D+00
  b = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TESTND'
  write ( *, '(a)' ) '  Test routines for estimating the integral of'
  write ( *, '(a)' ) '  of F(X) in the hypercube [A,B]^DIM_NUM.'
  write ( *, '(a)' ) ' '

  do test = 1, test_num

    dim_num = dim_num_test(test)

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  DIM_NUM = ', dim_num
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a,g20.12)' ) '  A(1:DIM_NUM) = ', a
    write ( *, '(a,g20.12)' ) '  B(1:DIM_NUM) = ', b

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  F(X(1:DIM_NUM)) = 1'
    write ( *, '(a)' ) ' '

    call testnd ( dim_num, f1dn )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  F(X(1:DIM_NUM)) = sum ( X(1:DIM_NUM) )'
    write ( *, '(a)' ) ' '

    call testnd ( dim_num, fxdn )

    write ( *, '(a)') ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  F(X(1:DIM_NUM)) = sum( X(1:DIM_NUM)^2 )'
    write ( *, '(a)' ) ' '

    call testnd ( dim_num, fx2dn )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  F(X(1:DIM_NUM)) = sum ( X(1:DIM_NUM)^3 )'
    write ( *, '(a)' ) ' '

    call testnd ( dim_num, fx3dn )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  F(X(1:DIM_NUM)) = exp(sum(X(1:DIM_NUM)))'
    write ( *, '(a)' ) ' '

    call testnd ( dim_num, fedn )

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  F(X(1:DIM_NUM)) = 1/(1+sum(X(1:DIM_NUM)^2))'
    write ( *, '(a)' ) ' '

    call testnd ( dim_num, fbdn )

  end do
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'NINTLIB_TEST'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop 0
end
subroutine testnd ( dim_num, func )

!*****************************************************************************80
!
!! TESTND tests the integrators on a particular function.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    11 September 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ), external FUNC, the name of the function
!    to be integrated.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ), external :: func
  integer dim_num

  call test01 ( dim_num, func )
  call test02 ( dim_num, func )
  call test03 ( dim_num, func )
  call test04 ( dim_num, func )
!
!  TEST05 is only set up for DIM_NUM = 2.
!
  if ( dim_num == 2 ) then
    call test05 ( dim_num, func )
  end if

  call test06 ( dim_num, func )

  return
end
subroutine test01 ( dim_num, func )

!*****************************************************************************80
!
!! TEST01 tests BOX_ND.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    11 September 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ), external FUNC, the name of the function
!    to be integrated.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num
  integer, parameter :: order = 5

  integer eval_num
  real ( kind = rk8 ), external :: func
  real ( kind = rk8 ) result
  real ( kind = rk8 ), dimension ( order ) :: wtab = (/ &
    0.236926885056189087514264040720D+00, &
    0.478628670499366468041291514836D+00, &
    0.568888888888888888888888888889D+00, &
    0.478628670499366468041291514836D+00, &
    0.236926885056189087514264040720D+00 /)
  real ( kind = rk8 ), dimension ( order ) :: wtab2
  real ( kind = rk8 ), dimension ( order ) :: xtab = (/ &
    -0.906179845938663992797626878299D+00, &
    -0.538469310105683091036314420700D+00, &
     0.0D+00, &
     0.538469310105683091036314420700D+00, &
     0.906179845938663992797626878299D+00 /)
  real ( kind = rk8 ), dimension ( order ) :: xtab2
!
!  Adjust the quadrature rule from [-1,1] to [0,1]:
!
  xtab2(1:order) = ( xtab(1:order) + 1.0D+00 ) / 2.0D+00
  wtab2(1:order) = 0.5D+00 * wtab(1:order)

  call box_nd ( func, dim_num, order, xtab2, wtab2, result, eval_num )

  write ( *, '(a,g20.12,2x,i8)' ) '  BOX_ND:         ', result, eval_num

  return
end
subroutine test02 ( dim_num, func )

!*****************************************************************************80
!
!! TEST02 tests P5_ND.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    11 September 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ), external FUNC, the name of the function
!    to be integrated.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) a(dim_num)
  real ( kind = rk8 ) b(dim_num)
  integer eval_num
  real ( kind = rk8 ), external :: func
  real ( kind = rk8 ) result
!
!  Set the integration limits.
!
  a(1:dim_num) = 0.0D+00
  b(1:dim_num) = 1.0D+00

  call p5_nd ( func, dim_num, a, b, result, eval_num )

  write ( *, '(a,g20.12,2x,i8)' ) '  P5_ND:          ', result, eval_num

  return
end
subroutine test03 ( dim_num, func )

!*****************************************************************************80
!
!! TEST03 tests ROMBERG_ND.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    11 September 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ), external FUNC, the name of the function
!    to be integrated.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) a(dim_num)
  real ( kind = rk8 ) b(dim_num)
  integer eval_num
  real ( kind = rk8 ), external :: func
  integer ind
  integer, parameter :: it_max = 3
  real ( kind = rk8 ) result
  integer sub_num(dim_num)
  real ( kind = rk8 ) tol
!
!  Set the integration limits.
!
  a(1:dim_num) = 0.0D+00
  b(1:dim_num) = 1.0D+00

  tol = 0.001D+00

  sub_num(1:dim_num) = 10

  call romberg_nd ( func, a, b, dim_num, sub_num, it_max, tol, &
    result, ind, eval_num )

  write ( *, '(a,g20.12,2x,i8)' ) '  ROMBERG_ND:     ', result, eval_num

  return
end
subroutine test04 ( dim_num, func )

!*****************************************************************************80
!
!! TEST04 tests SAMPLE_ND.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    11 September 2006
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ), external FUNC, the name of the function
!    to be integrated.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  integer, parameter :: k2 = 4

  real ( kind = rk8 ) dev1(k2)
  real ( kind = rk8 ) dev2(k2)
  real ( kind = rk8 ) err1(k2)
  real ( kind = rk8 ) est1(k2)
  real ( kind = rk8 ) est2(k2)
  real ( kind = rk8 ) err2(k2)
  integer eval_num
  real ( kind = rk8 ), external :: func
  integer k1

  k1 = 1

  call sample_nd ( func, k1, k2, dim_num, est1, err1, dev1, est2, err2, &
    dev2, eval_num )

  write ( *, '(a,g20.12,2x,i8)' ) '  SAMPLE_ND:      ', est2(k2), eval_num

  return
end
subroutine test05 ( dim_num, func )

!*****************************************************************************80
!
!! TEST05 demonstrates how to refine N-dimensional integration results.
!
!  Discussion:
!
!    This routine is only set up for DIM_NUM = 2 for now.
!
!    We are given a routine, NDP5, which will integrate over a
!    DIM_NUM dimensional hypercube using a fixed method.  In order to
!    improve the approximation to an integral, we can subdivide
!    the hypercube and call NDP5 to integrate again over each of
!    these regions.
!
!    The information that we gather can be used to tell us when
!    to expect that we have achieved a certain degree of accuracy.
!
!    With a little more work, we could make this code adaptive.
!    That is, it would only refine SOME of the subregions, where
!    the approximation to the integral was still not good enough.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    25 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ), external FUNC, the name of the function
!    to be integrated.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) a(dim_num)
  real ( kind = rk8 ) b(dim_num)
  integer eval_num
  integer eval_total
  real ( kind = rk8 ), external :: func
  integer i
  integer igrid
  integer j
  integer ngrid
  real ( kind = rk8 ) result
  real ( kind = rk8 ) result_total
  real ( kind = rk8 ) xlo(dim_num)
  real ( kind = rk8 ) xhi(dim_num)

  xlo(1:dim_num) = 0.0D+00
  xhi(1:dim_num) = 1.0D+00

  do igrid = 1, 6

    ngrid = 2**( igrid - 1 )

    result_total = 0.0D+00
    eval_total = 0

    do i = 1, ngrid

      a(1) = ( real ( ngrid - i + 1, kind = rk8 ) * xlo(1)   &
             + real (         i - 1, kind = rk8 ) * xhi(1) ) &
             / real ( ngrid,         kind = rk8 )

      b(1) = ( real ( ngrid - i, kind = rk8 ) * xlo(1)   &
             + real (         i, kind = rk8 ) * xhi(1) ) &
             / real ( ngrid,     kind = rk8 )

      do j = 1, ngrid

        a(2) = ( real ( ngrid - j + 1, kind = rk8 ) * xlo(2)   &
               + real (         j - 1, kind = rk8 ) * xhi(2) ) &
               / real ( ngrid,         kind = rk8 )

        b(2) = ( real ( ngrid - j, kind = rk8 ) * xlo(2)   &
               + real (         j, kind = rk8 ) * xhi(2) ) &
               / real ( ngrid,     kind = rk8 )

        call p5_nd ( func, dim_num, a, b, result, eval_num )

        result_total = result_total + result
        eval_total = eval_total + eval_num

      end do

    end do

    write ( *, '(a,g20.12,2x,i8)' ) &
      '  P5_ND+:         ', result_total, eval_total

  end do

  return
end
subroutine test06 ( dim_num, func )

!*****************************************************************************80
!
!! TEST06 tests MONTE_CARLO_ND.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    25 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ), external FUNC, the name of the function
!    to be integrated.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num
  integer, parameter :: test_num = 3

  real ( kind = rk8 ) a(dim_num)
  real ( kind = rk8 ) b(dim_num)
  integer eval_num
  real ( kind = rk8 ), external :: func
  real ( kind = rk8 ) result
  integer test
!
!  Set the integration limits.
!
  a(1:dim_num) = 0.0D+00
  b(1:dim_num) = 1.0D+00

  do test = 1, test_num

    eval_num = 8**test * 10000

    call monte_carlo_nd ( func, dim_num, a, b, eval_num, result )

    write ( *, '(a,g20.12,2x,i8)' ) '  MONTE_CARLO_ND: ', result, eval_num

  end do

  return
end
function f1dn ( dim_num, x )

!*****************************************************************************80
!
!! F1DN(X(1:DIM_NUM)) = 1.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ) X(DIM_NUM), the argument.
!
!    Output, real ( kind = rk8 ) F1DN, the value of the function at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) f1dn
  real ( kind = rk8 ) x(dim_num)

  call r8_fake_use ( x(1) )

  f1dn = 1.0D+00

  return
end
function fbdn ( dim_num, x )

!*****************************************************************************80
!
!! FBDN(X(1:DIM_NUM)) = 1 / ( 1 + sum ( X(1:DIM_NUM)**2 ) )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ) X(DIM_NUM), the argument.
!
!    Output, real ( kind = rk8 ) FBDN, the value of the function at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) fbdn
  real ( kind = rk8 ) x(dim_num)

  fbdn = 1.0D+00 / ( 1.0D+00 + sum ( x(1:dim_num)**2 ) )

  return
end
function fedn ( dim_num, x )

!*****************************************************************************80
!
!! FEDN(X(1:DIM_NUM)) = EXP ( sum ( X(1:DIM_NUM) ) )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ) X(DIM_NUM), the argument.
!
!    Output, real ( kind = rk8 ) FEDN, the value of the function at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) fedn
  real ( kind = rk8 ) x(dim_num)

  fedn = exp ( sum ( x(1:dim_num) ) )

  return
end
function fxdn ( dim_num, x )

!*****************************************************************************80
!
!! FXDN(X(1:DIM_NUM)) = sum ( X(1:DIM_NUM) )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ) X(DIM_NUM), the argument.
!
!    Output, real ( kind = rk8 ) FXDN, the value of the function at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) fxdn
  real ( kind = rk8 ) x(dim_num)

  fxdn = sum ( x(1:dim_num) )

  return
end
function fx2dn ( dim_num, x )

!*****************************************************************************80
!
!! FX2DN(X(1:DIM_NUM)) = sum ( X(1:DIM_NUM)**2 )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ) X(DIM_NUM), the argument.
!
!    Output, real ( kind = rk8 ) FX2DN, the value of the function at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) fx2dn
  real ( kind = rk8 ) x(dim_num)

  fx2dn = sum ( x(1:dim_num)**2 )

  return
end
function fx3dn ( dim_num, x )

!*****************************************************************************80
!
!! FX3DN(X(1:DIM_NUM)) = sum ( X(1:DIM_NUM)**3 )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    24 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real ( kind = rk8 ) X(DIM_NUM), the argument.
!
!    Output, real ( kind = rk8 ) FX3DN, the value of the function at X.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk8 ) fx3dn
  real ( kind = rk8 ) x(dim_num)

  fx3dn = sum ( x(1:dim_num)**3 )

  return
end
subroutine r8_fake_use ( x )

!*****************************************************************************80
!
!! r8_fake_use() pretends to use an R8 variable.
!
!  Discussion:
!
!    Some compilers will issue a warning if a variable is unused.
!    Sometimes there's a good reason to include a variable in a program,
!    but not to use it.  Calling this function with that variable as
!    the argument will shut the compiler up.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 April 2020
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) X, the variable to be "used".
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) x

  if ( x /= x ) then
    write ( *, '(a)' ) '  r8_fake_use(): variable is NAN.'
  end if

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end