program nintlib_prb
!
!***********************************************************************
!
!! NINTLIB_PRB runs the NINTLIB tests.
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'NINTLIB_PRB'
  write ( *, '(a)' ) '  Sample problems for the NINTLIB library.'
 
  call testnd
  
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'NINTLIB_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine testnd
!
!***********************************************************************
!
  real a
  real b
  real, external :: f1dn
  real, external :: fbdn
  real, external :: fedn
  real, external :: fxdn
  real, external :: fx2dn
  real, external :: fx3dn
!
  a = 0.0E+00
  b = 1.0E+00
  write ( *, * ) ' '
  write ( *, * ) 'TESTND'
  write ( *, * ) '  Test N-D quadrature codes'
  write ( *, * ) '  for integral of F(X) on [A,B]**N.'
  write ( *, * ) ' '
  write ( *, * ) '  A=', a
  write ( *, * ) '  B=', b
 
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  F(X)=1'
  write ( *, * ) ' '

  call tstnd ( f1dn )
 
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  F(X)=X'
  write ( *, * ) ' '

  call tstnd ( fxdn )
 
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  F(X)=X*2'
  write ( *, * ) ' '

  call tstnd ( fx2dn )
 
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  F(X)=X*3'
  write ( *, * ) ' '

  call tstnd ( fx3dn )
 
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  F(X)=EXP(X)'
  write ( *, * ) ' '

  call tstnd ( fedn )
 
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  F(X)=1/(1+X*X)'
  write ( *, * ) ' '

  call tstnd ( fbdn )

  return
end
subroutine tstnd ( func )
!
!***********************************************************************
!
!! TESTND tests the ND integrators.
!
  real, external :: func
!
  call test32 ( func )
  call test325 ( func )
  call test33 ( func )
  call test34 ( func )
  call test35 ( func )
 
  return
end
subroutine test32 ( func )
!
!***********************************************************************
!
!! TEST32 tests NDP5.
!
  integer, parameter :: ndim = 2
!
  real aval(ndim)
  real bval(ndim)
  real, external :: func
  integer i
  real result
  real work(ndim)
!
!  Set the integration limits.
!
  aval(1:ndim) = 0.0E+00
  bval(1:ndim) = 1.0E+00
 
  call ndp5 ( func, ndim, aval, bval, work, result )
 
  write ( *, '(a6,g13.6)' ) 'NDP5  ', result
 
  return
end
subroutine test325 ( func )
!
!***********************************************************************
!
!! TEST325 tests NDBOX.
!
  integer, parameter :: ndim = 2
  integer, parameter :: norder = 5
!
  real, external :: func
  real result
  real wtab(norder)
  real xtab(norder)
!
  xtab(1) = - 0.906179845938663992797626878299E+00
  xtab(2) = - 0.538469310105683091036314420700E+00
  xtab(3) =   0.0E+00
  xtab(4) =   0.538469310105683091036314420700E+00
  xtab(5) =   0.906179845938663992797626878299E+00
 
  wtab(1) = 0.236926885056189087514264040720E+00
  wtab(2) = 0.478628670499366468041291514836E+00
  wtab(3) = 0.568888888888888888888888888889E+00
  wtab(4) = 0.478628670499366468041291514836E+00
  wtab(5) = 0.236926885056189087514264040720E+00
!
!  Adjust the quadrature rule from [-1,1] to [0,1]:
!
  xtab(1:norder) = ( xtab(1:norder) + 1.0E+00 ) / 2.0E+00
  wtab(1:norder) = 0.5E+00 * wtab(1:norder)

  call ndbox ( func, ndim, norder, xtab, wtab, result )

  write ( *, '(a6,g13.6)' ) 'NDBOX ', result

  return
end
subroutine test33 ( func )
!
!***********************************************************************
!
!! TEST33 tests NDSAMP.
!
  integer, parameter :: k2 = 4
!
  real dev1(k2)
  real dev2(k2)
  real err1(k2)
  real est1(k2)
  real est2(k2)
  real err2(k2)
  real, external :: func
  integer k1
  integer ndim
!
  k1 = 1
  ndim = 2
  call ndsamp ( func, k1, k2, ndim, est1, err1, dev1, est2, err2, dev2 )
 
  write ( *, '(a6,g13.6)' ) 'NDSAMP', est2(k2)
 
  return
end
subroutine test34 ( func )
!
!***********************************************************************
!
!! TEST34 tests NDROMB.
!
  integer, parameter :: maxit = 3
  integer, parameter :: ndim = 2

  integer, parameter :: nwork = maxit + ndim
!
  real aval(ndim)
  real bval(ndim)
  real eps
  real, external :: func
  integer i
  integer ind
  integer iwork(nwork)
  integer nsub(ndim)
  real result
  real work(nwork)
!
!  Set the integration limits.
!
  aval(1:ndim) = 0.0E+00
  bval(1:ndim) = 1.0E+00
 
  eps = 0.001E+00
 
  nsub(1:2) = (/ 10, 10 /)

  call ndromb ( func, aval, bval, ndim, nsub, maxit, eps, iwork, &
    work, nwork, result, ind )
 
  write ( *, '(a6,g13.6)' ) 'NDROMB', result
 
  return
end
subroutine test35 ( func )
!
!***********************************************************************
!
!! TEST35 demonstrates how to refine N-dimensional integration results.
!
!  We are given a routine, NDP5, which will integrate over an
!  N 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.
!
  integer, parameter :: ndim = 2
!
  real aval(ndim)
  real bval(ndim)
  real, external :: func
  integer i
  integer igrid
  integer j
  integer ngrid
  real result
  real total
  real work(ndim)
  real xlo(ndim)
  real xhi(ndim)
!
  xlo(1:2) = 0.0E+00
  xhi(1:2) = 1.0E+00
 
  do igrid = 1, 6
 
    ngrid = 2**( igrid - 1 )
    total = 0.0E+00
 
    do i = 1, ngrid
 
      aval(1) = ( real ( ngrid + 1 - i ) * xlo(1) + real ( i - 1 ) * xhi(1) ) / real ( ngrid )
      bval(1) = ( real ( ngrid - i ) * xlo(1) + real ( i ) * xhi(1) ) / real ( ngrid )
  
      do j = 1, ngrid

        aval(2) = ( real ( ngrid + 1 - j ) * xlo(2)+ real ( j - 1 ) * xhi(2) ) / real ( ngrid )
        bval(2) = ( real ( ngrid - j ) * xlo(2) + real ( j ) * xhi(2) ) / real ( ngrid )

        call ndp5 ( func, ndim, aval, bval, work, result )

        total = total + result

      end do
 
    end do
 
    write ( *, '(a6,g13.6)' ) 'NDP5-R', total
 
  end do
 
  return
end
function fbdn ( n, x )
!
!***********************************************************************
!
!! FBDN(X(1:N)) = 1 / ( 1 + SUM ( X(1:N)**2 ) )
!
  integer n
!
  real fbdn
  real x(n)
!
  fbdn = 1.0E+00 / ( 1.0E+00 + sum ( x(1:n)**2 ) )
 
  return
end
function f1dn ( n, x )
!
!***********************************************************************
!
!! F1DN(X(1:N)) = 1.
!
  integer n
!
  real f1dn
  real x(n)
!
  f1dn = 1.0E+00
 
  return
end
function fxdn ( n, x )
!
!***********************************************************************
!
!! FXDN(X(1:N)) = SUM ( X(1:N) )
!
  integer n
!
  real fxdn
  real x(n)
!
  fxdn = sum ( x(1:n) )
 
  return
end
function fx2dn ( n, x )
!
!***********************************************************************
!
!! FX2DN(X(1:N)) = SUM ( X(1:N)**2 )
!
  integer n
!
  real fx2dn
  real x(n)
!
  fx2dn = sum ( x(1:n)**2 )

  return
end
function fx3dn ( n, x )
!
!***********************************************************************
!
!! FX3DN(X(1:N)) = SUM ( X(1:N)**3 )
!
  integer n
!
  real fx3dn
  real x(n)
!
  fx3dn = sum ( x(1:n)**3 )
 
  return
end
function fedn ( n, x )
!
!***********************************************************************
!
!! FEDN(X(1:N)) = EXP ( SUM ( X(1:N) ) )
!
  integer n
!
  real fedn
  real x(n)
!
  fedn = exp ( sum ( x(1:n) ) )
 
  return
end