program intlib_prb
!
!***********************************************************************
!
!! INTLIB_PRB runs the INTLIB tests.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INTLIB_PRB'
  write ( *, '(a)' ) '  Sample problems for the INTLIB library.'
 
  call test1d
  
  call test27
  call test28
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INTLIB_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test1d
!
!***********************************************************************
!
!! TEST1D tests all the 1D integration codes.
!
  implicit none
!
  real a
  real b
  real, external :: f1d1
  real, external :: fbd1
  real, external :: fed1
  real, external :: fqd1
  real, external :: fxd1
  real, external :: fx2d1
  real, external :: fx3d1
!
  a = 0.0E+00
  b = 1.0E+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST1D'
  write ( *, '(a)' ) '  Test 1D quadrature codes'
  write ( *, '(a)' ) '  for integral of F(X) on [A,B].'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  A = ', a
  write ( *, '(a,g14.6)' ) '  B = ', b
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)=1'
  write ( *, '(a)' ) ' '

  call tst1d ( f1d1 )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)=X'
  write ( *, '(a)' ) ' '
 
  call tst1d ( fxd1 )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)=X*2'
  write ( *, '(a)' ) ' '

  call tst1d ( fx2d1 )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)=X*3'
  write ( *, '(a)' ) ' '

  call tst1d ( fx3d1 )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)=EXP(X)'
  write ( *, '(a)' ) ' '

  call tst1d ( fed1 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)=SQRT(X)'
  write ( *, '(a)' ) ' '

  call tst1d ( fqd1 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  F(X)=1/(1+X*X)'
  write ( *, '(a)') ' '

  call tst1d ( fbd1 )

  return
end
subroutine tst1d ( func )
!
!***********************************************************************
!
  implicit none
!
  real, external :: func
!
  call test01 ( func )
  call test02 ( func )
  call test03 ( func )
  call test04 ( func )
  call test05 ( func )
  call test06 ( func )
  call test07 ( func )
  call test08 ( func )
  call test09 ( func )
  call test10 ( func )
  call test11 ( func )
  call test12 ( func )
  call test13 ( func )
  call test14 ( func )
  call test15 ( func )
  call test16 ( func )
  call test17 ( func )

  return
end
subroutine test01 ( func )
!
!***********************************************************************
!
!! TEST01 tests PLINT.
!
  implicit none
!
  integer, parameter :: ntab = 11
!
  real a
  real b
  real ftab(ntab)
  real, external :: func
  integer i
  real result
  real xtab(ntab)
!
  a = 0.0E+00
  b = 1.0E+00
 
  call rvec_even ( a, b, ntab, xtab )

  do i = 1, ntab
    ftab(i) = func ( xtab(i) )
  end do
 
  call plint ( ftab, xtab, ntab, a, b, result )
 
  write ( *, '(a6,g13.6)' ) 'PLINT ', result
 
  return
end
subroutine test02 ( func )
!
!***********************************************************************
!
!! TEST02 tests AVINT.
!
  implicit none
!
  integer, parameter :: ntab = 11
!
  real a
  real b
  real ftab(ntab)
  real, external :: func
  integer i
  real result
  real xtab(ntab)
!
  a = 0.0E+00
  b = 1.0E+00
 
  call rvec_even ( a, b, ntab, xtab )

  do i = 1, ntab
    ftab(i) = func ( xtab(i) )
  end do
 
  call avint ( ftab, xtab, ntab, a, b, result )
 
  write ( *, '(a6,g13.6)' ) 'AVINT ', result
 
  return
end
subroutine test03 ( func )
!
!***********************************************************************
!
!! TEST03 tests CUBINT.
!
  implicit none
!
  integer, parameter :: ntab = 11
!
  real a
  real b
  real error
  real ftab(ntab)
  real, external :: func
  integer i
  integer ia
  integer ib
  real result
  real xtab(ntab)
!
  a = 0.0E+00
  b = 1.0E+00
 
  ia = 1
  ib = ntab

  call rvec_even ( a, b, ntab, xtab )

  do i = 1, ntab
    ftab(i) = func ( xtab(i) )
  end do
 
  call cubint ( ftab, xtab, ntab, ia, ib, result, error )
 
  write ( *, '(a6,g13.6)' ) 'CUBINT', result
 
  return
end
subroutine test04 ( func )
!
!***********************************************************************
!
!! TEST04 tests WEDINT.
!
  implicit none
!
  integer, parameter :: n = 2
  integer, parameter :: ntab = 6 * n + 1
!
  real a
  real b
  real ftab(ntab)
  real, external :: func
  real h
  integer i
  real result
  real xtab(ntab)
!
  a = 0.0E+00
  b = 1.0E+00
 
  h = ( b - a ) / real ( ntab - 1 )
 
  call rvec_even ( a, b, ntab, xtab )

  do i = 1, ntab
    ftab(i) = func ( xtab(i) )
  end do
 
  call wedint ( ftab, h, ntab, result )
 
  write ( *, '(a6,g13.6)' ) 'WEDINT', result
 
  return
end
subroutine test05 ( func )
!
!***********************************************************************
!
!! TEST05 tests CSPINT.
!
  implicit none
!
  integer, parameter :: ntab = 13
!
  real a
  real aleft
  real b
  real brite
  real e(ntab)
  real ftab(ntab)
  real, external :: func
  integer i
  real result
  real xtab(ntab)
  real work(ntab)
  real y(3,ntab)
!
  a = 0.0E+00
  b = 1.0E+00
!
!  Note that, for accuracy, it is useful to have two data points
!  outside the interval (A,B).
!
  aleft = a - ( b - a ) / real ( ntab - 3 )
  brite = b + ( b - a ) / real ( ntab - 3 )
 
  call rvec_even ( aleft, brite, ntab, xtab )

  do i = 1, ntab
    ftab(i) = func ( xtab(i) )
  end do
 
  call cspint ( ftab, xtab, ntab, a, b, y, e, work, result )
 
  write ( *, '(a6,g13.6)' ) 'CSPINT', result
 
  return
end
subroutine test06 ( func )
!
!***********************************************************************
!
!! TEST06 tests GAUS8.
!
  implicit none
!
  real a
  real b
  real err
  real, external :: func
  integer ier
  real result
!
  a = 0.0E+00
  b = 1.0E+00
  err = 0.001E+00
 
  call gaus8 ( func, a, b, err, result, ier )
 
  write ( *, '(a6,g13.6)' ) 'GAUS8 ', result
 
  return
end
subroutine test07 ( func )
!
!***********************************************************************
!
!! TEST07 tests QNC79.
!
  implicit none
!
  real a
  real b
  real err
  real, external :: func
  integer ier
  integer k
  real result
 !
  a = 0.0E+00
  b = 1.0E+00
  err = 0.001E+00
 
  call qnc79 ( func, a, b, err, result, ier, k )
 
  write ( *, '(a6,g13.6)' ) 'QNC79 ', result
 
  return
end
subroutine test08 ( func )
!
!***********************************************************************
!
!! TEST08 tests QUAD.
!
  implicit none
!
  integer, parameter :: nleast = 3
  integer, parameter :: nmost = nleast + 2
!
  real a
  real abserr
  real b
  real, external :: func
  real relerr
  real result
  real work(nmost+1)
 !
  a = 0.0E+00
  b = 1.0E+00
  abserr = 0.001E+00
  relerr = 0.001E+00
 
  call quad ( func, a, b, abserr, relerr, nleast, nmost, work, result )
 
  write ( *, '(a6,g13.6)' ) 'QUAD  ', result
 
  return
end
subroutine test09 ( func )
!
!***********************************************************************
!
!! TEST09 tests RMINSP.
!
  implicit none
!
  real a
  real b
  real epsin
  real epsout
  real, external :: func
  integer iop
  real result
!
  a = 0.0E+00
  b = 1.0E+00
  epsin = 0.001E+00
  iop = 1
 
  call rminsp ( func, a, b, epsin, epsout, iop, result )
 
  write ( *, '(a6,g13.6)' ) 'RMINSP', result
 
  return
end
subroutine test10 ( func )
!
!***********************************************************************
!
!! TEST10 tests RMINSP with cosine transformation.
!
  implicit none
!
  real a
  real b
  real epsin
  real epsout
  real, external :: func
  integer iop
  real result
!
  a = 0.0E+00
  b = 1.0E+00
  epsin = 0.001E+00
  iop = 2
 
  call rminsp ( func, a, b, epsin, epsout, iop, result )
 
  write ( *, '(a6,g13.6)' ) 'RMINSP', result
 
  return
end
subroutine test11 ( func )
!
!***********************************************************************
!
!! TEST11 tests IRATEX.
!
  implicit none
!
  real a
  real b
  real epsin
  real epsout
  real, external :: func
  integer ind
  real result
!
  a = 0.0E+00
  b = 1.0E+00
  epsin = 0.001E+00
 
  call iratex ( func, a, b, epsin, epsout, result, ind )
 
  write ( *, '(a6,g13.6)' ) 'IRATEX', result
 
  return
end
subroutine test12 ( func )
!
!***********************************************************************
!
!! TEST12 tests CADRE.
!
  implicit none
!
  real a
  real abserr
  real b
  real error
  real, external :: func
  integer ind
  real relerr
  real result
!
  a = 0.0E+00
  b = 1.0E+00
  abserr = 0.001E+00
  relerr = 0.001E+00
 
  call cadre ( func, a, b, abserr, relerr, error, result, ind )
 
  write ( *, '(a6,g13.6)' ) 'CADRE ', result
 
  return
end
subroutine test13 ( func )
!
!***********************************************************************
!
!! TEST13 tests CHINSP.
!
  implicit none
!
  real a
  real b
  real epsin
  real epsout
  real, external :: func
  real result
!
  a = 0.0E+00
  b = 1.0E+00
  epsin = 0.001E+00
 
  call chinsp ( func, a, b, epsin, epsout, result )
 
  write ( *, '(a6,g13.6)' ) 'CHINSP', result
 
  return
end
subroutine test14 ( func )
!
!***********************************************************************
!
!! TEST14 tests SIMP.
!
  implicit none
!
  real a
  real b
  real eps
  real, external :: func
  real result
!
  a = 0.0E+00
  b = 1.0E+00
  eps = 0.001E+00
 
  call simp ( func, a, b, eps, result )
 
  write ( *, '(a6,g13.6)' ) 'SIMP  ', result
 
  return
end
subroutine test15 ( func )
!
!***********************************************************************
!
!! TEST15 tests HIORDQ.
!
  implicit none
!
  integer, parameter :: n = 11
  integer, parameter :: nwork = 2 * ( n - 1 )
!
  real a
  real b
  real, external :: func
  real h
  integer i
  real result
  real work(nwork)
  real x(n)
  real y(n)
!
  a = 0.0E+00
  b = 1.0E+00
 
  call rvec_even ( a, b, n, x )

  do i = 1, n
    y(i) = func ( x(i) )
  end do
 
  h = ( b - a ) / real ( n - 1 )
 
  call hiordq ( n, y, h, work, result )
 
  write ( *, '(a6,g13.6)' ) 'HIORDQ', result
 
  return
end
subroutine test16 ( func )
!
!***********************************************************************
!
!! TEST16 tests SIMPSN.
!
  implicit none
!
  integer, parameter :: n = 11
!
  real a
  real b
  real, external :: func
  real h
  integer i
  real result
  real x(n)
  real y(n)
!
  a = 0.0E+00
  b = 1.0E+00
 
  call rvec_even ( a, b, n, x )

  do i = 1, n
    y(i) = func ( x(i) )
  end do
 
  h = ( b - a ) / real ( n - 1 )
 
  call simpsn ( h, y, n, result )
 
  write ( *, '(a6,g13.6)' ) 'SIMPSN', result
 
  return
end
subroutine test17 ( func )
!
!***********************************************************************
!
!! TEST17 tests SIMPNE.
!
  implicit none
!
  integer, parameter :: n = 11
!
  real a
  real b
  real, external :: func
  integer i
  real result
  real x(n)
  real y(n)
!
  a = 0.0E+00
  b = 1.0E+00
 
  call rvec_even ( a, b, n, x )

  do i = 1, n
    y(i) = func ( x(i) )
  end do
 
  call simpne ( x, y, n, result )
 
  write ( *, '(a6,g13.6)' ) 'SIMPNE', result
 
  return
end
subroutine test27
!
!***********************************************************************
!
!! TEST27 tests FILON_COS.
!
  implicit none
!
  integer, parameter :: ntab = 11
!
  real a
  real b
  real exact
  real ftab(ntab)
  integer i
  integer ii
  integer j
  real result
  real t
  real xtab(ntab)
!
  a = 0.0E+00
  b = 1.0E+00

  call rvec_even ( a, b, ntab, xtab )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST27'
  write ( *, '(a)') '  FILON_COS estimates the integral of.'
  write ( *, '(a)' ) '  F(X) * COS ( T * X )'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integrate F(X)*COS(T*X):'
  write ( *, '(a)' ) '  with F(X)=1, X, X**2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  A = ', a
  write ( *, '(a,g14.6)' ) '  B = ', b
  write ( *, '(a,i6)' ) '  NTAB = ', ntab
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       T       Approximate      Exact'
  write ( *, '(a)' ) ' '
 
  do i = 1, 9
 
    ii = mod ( i - 1, 3 ) + 1

    do j = 1, ntab
 
      if ( i == 1 )then
        ftab(j) = 1.0E+00
      else
        ftab(j) = xtab(j)**(ii-1)
      end if
 
    end do
 
    if ( i <= 3 ) then
      t = 1.0E+00
    else if ( i <= 6 ) then
      t = 2.0E+00
    else
      t = 10.0E+00
    end if

    call filon_cos ( ftab, a, b, ntab, t, result )

    if ( ii == 1 ) then
      exact = ( sin ( t * b ) - sin ( t * a ) ) / t
    else if ( ii == 2 ) then
      exact = ( ( cos ( t * b ) + t * b * sin ( t * b ) ) &
              - ( cos ( t * a ) + t * a * sin ( t * a ) ) ) / t**2
    else if ( ii == 3 ) then
      exact = ( ( 2.0E+00 * t * b * cos ( t * b ) &
            + ( t**2 * b**2 - 2.0E+00 ) * sin ( t * b ) ) &
              - ( 2.0E+00 * t * a * cos ( t * a ) &
            + ( t**2 * a**2 - 2.0E+00 ) * sin ( t * a ) ) ) / t**3
    end if

    if ( ii == 1 ) then
      write ( *, '(a)' ) ' '
    end if

    write ( *, '(3g14.6)' ) t, result, exact
 
  end do
 
  return
end
subroutine test28
!
!***********************************************************************
!
!! TEST28 tests FILON_SIN.
!
  implicit none
!
  integer, parameter :: ntab = 11
!
  real a
  real b
  real exact
  real ftab(ntab)
  integer i
  integer ii
  integer j
  real result
  real t
  real x(ntab)
!
  a = 0.0E+00
  b = 1.0E+00

  call rvec_even ( a, b, ntab, x )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST28'
  write ( *, '(a)' ) '  FILON_SIN estimates the integral of.'
  write ( *, '(a)' ) '  F(X) * SIN ( T * X )'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  A = ', a
  write ( *, '(a,g14.6)' ) '  B = ', b
  write ( *, '(a,i6)' ) '  NTAB = ', ntab
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integrate F(X)*SIN(T*X)'
  write ( *, '(a)' ) '  with F(X)=1, X, X**2.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       T       Approximate      Exact'
  write ( *, '(a)' ) ' '

  do i = 1, 9
 
    ii = mod ( i - 1, 3 ) + 1

    do j = 1, ntab
 
      if ( ii == 1 ) then
        ftab(j) = 1.0E+00
      else
        ftab(j) = x(j)**(ii-1)
      end if
 
    end do
 
    if ( i <= 3 ) then
      t = 1.0E+00
    else if ( i <= 6 ) then
      t = 2.0E+00
    else
      t = 10.0E+00
    end if

    call filon_sin ( ftab, a, b, ntab, t, result )

    if ( ii == 1 ) then
      exact = ( - cos ( t * b ) + cos ( t * a ) ) / t
    else if ( ii == 2 ) then
      exact = ( ( sin ( t * b ) - t * b * cos ( t * b ) ) &
              - ( sin ( t * a ) - t * a * cos ( t * a ) ) ) / t**2
    else if ( ii == 3 ) then
      exact = ( ( 2.0E+00 * t * b * sin ( t * b ) &
              + ( 2.0E+00 - t**2 * b**2 ) * cos ( t * b ) ) &
              - ( 2.0E+00 * t * a * sin ( t * a ) &
              + ( 2.0E+00 - t**2 * a**2 ) * cos ( t * a ) ) ) / t**3
    end if
 
    if ( ii == 1 ) then
      write ( *, '(a)' ) ' '
    end if

    write ( *, '(3g14.6)' ) t, result, exact
 
  end do
 
  return
end
function f1d1 ( x )
!
!***********************************************************************
!
!! F1D1(X) = 1.
!
  implicit none
!
  real f1d1
  real x
!
  f1d1 = 1.0E+00
 
  return
end
function f1d2 ( x, y )
!
!***********************************************************************
!
!! F1D2(X,Y) = 1.
!
  implicit none
!
  real f1d2
  real x
  real y
!
  f1d2 = 1.0E+00
 
  return
end
function fbd1 ( x )
!
!***********************************************************************
!
!! FBD1(X) = 1 / ( 1 + X**2 )
!
  implicit none
!
  real fbd1
  real x
!
  fbd1 = 1.0E+00 / ( 1.0E+00 + x**2 )
 
  return
end
function fchby ( x )
!
!***********************************************************************
!
!! FCHBY ...
!
  implicit none
!
  real fchby
  real g
  integer iprob
  real x
!
  common /problm/ iprob
!
  if ( iprob == 1 )then
    g = 1.0E+00
  else
    g = x**(iprob-1)
  end if
 
  fchby = g / sqrt ( 1.0E+00 - x*x )
 
  return
end
function fed1 ( x )
!
!***********************************************************************
!
!! FED1(X) = EXP(X).
!
  implicit none
!
  real fed1
  real x
!
  fed1 = exp ( x )
 
  return
end
function fqd1 ( x )
!
!***********************************************************************
!
!! FQD1(X) = SQRT ( X ).
!
  implicit none
!
  real fqd1
  real x
!
  fqd1 = sqrt ( abs ( x ) )

  return
end
function fxd1 ( x )
!
!***********************************************************************
!
!! FXD1(X) = X.
!
  implicit none
!
  real fxd1
  real x
!
  fxd1 = x
 
  return
end
function fx2d1 ( x )
!
!***********************************************************************
!
!! FX2D1(X) = X**2.
!
  implicit none
!
  real fx2d1
  real x
!
  fx2d1 = x**2
 
  return
end
function fx3d1 ( x )
!
!***********************************************************************
!
!! FX3D1(X) = X**3.
!
  implicit none
!
  real fx3d1
  real x
!
  fx3d1 = x**3
 
  return
end
function f1sd1 ( x, y )
!
!***********************************************************************
!
!! F1SD1(X,Y) = 1 / SQRT ( 1 - X**2 )
!
  implicit none
!
  real f1sd1
  real x
  real y
!
  f1sd1 = 1.0E+00 / sqrt ( 1.0E+00 - x**2)
 
  return
end
function fxd2 ( x, y )
!
!***********************************************************************
!
!! FXD2(X,Y) = X + Y.
!
  implicit none
!
  real fxd2
  real x
  real y
!
  fxd2 = x + y
 
  return
end
function fx2d2 ( x, y )
!
!***********************************************************************
!
!! FX2D2(X,Y) = X**2 + Y**2
!
  implicit none
!
  real fx2d2
  real x
  real y
!
  fx2d2 = x**2 + y**2
 
  return
end
function fx3d2 ( x, y )
!
!***********************************************************************
!
!! FX3D2(X,Y) = X**3 + Y**3
!
  implicit none
!
  real fx3d2
  real x
  real y
!
  fx3d2 = x**3 + y**3
 
  return
end
function fxsd1 ( x )
!
!***********************************************************************
!
!! FXSD1(X) = X / SQRT ( 1 - X**2 )
!
  implicit none
!
  real fxsd1
  real x
!
  fxsd1 = x / sqrt ( 1.0E+00 - x**2 )
 
  return
end
function fx2sd1 ( x, y )
!
!***********************************************************************
!
!! FX2SD1(X,Y) = X**2 / SQRT ( 1 - X**2 )
!
  implicit none
!
  real fx2sd1
  real x
  real y
!
  fx2sd1 = x**2 / sqrt ( 1.0E+00 - x**2 )
 
  return
end
function fed2 ( x, y )
!
!***********************************************************************
!
!! FED2(X,Y) = EXP(X+Y).
!
  implicit none
!
  real fed2
  real x
  real y
!
  fed2 = exp ( x + y )
 
  return
end
function fbd2 ( x, y )
!
!***********************************************************************
!
!! FBD2(X,Y) = 1 / ( 1 + X**2 + Y**2 )
!
  implicit none
!
  real fbd2
  real x
  real y
!
  fbd2 = 1.0E+00 / ( 1.0E+00 + x**2 + y**2 )
 
  return
end