program quadrule_prb
!
!*******************************************************************************
!
!! QUADRULE_PRB calls a set of tests for the QUADRULE library.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QUADRULE_PRB'
  write ( *, '(a)' ) '  Test problems for the QUADRULE library.'

  call test01
  call test02
  call test03
  call test04
  call test05
  call test06
  call test07
  call test08
  call test09
  call test10

  call test11
  call test12
  call test13
  call test14
  call test15
  call test16
  call test17
  call test18
  call test19
  call test20

  call test21
  call test22
  call test23
  call test235
  call test236
  call test24
  call test25
  call test26
  call test27
  call test28
  call test29
  call test30

  call test31
  call test32
  call test33
  call test34
  call test35
  call test36
  call test37
  call test38
  call test39

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QUADRULE_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests BASHFORTH_SET.
!! TEST01 tests SUMMER.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 8
!
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  BASHFORTH_SET sets up Adams-Bashforth quadrature;'
  write ( *, '(a)' ) '  SUMMER carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integration interval is [0,1].'
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '
  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call bashforth_set ( norder, xtab, weight )
 
        call summer ( func, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests BDFC_SET.
!! TEST02 tests BDF_SUM.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 10
!
  double precision diftab(maxorder)
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) &
    '  BDFC_SET sets up Backward Difference Corrector quadrature;'
  write ( *, '(a)' ) '  BDF_SUM carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integration interval is [0,1].'
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call bdfc_set ( norder, weight, xtab )
 
        call bdf_sum ( func, norder, weight, xtab, diftab, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests BDFP_SET.
!! TEST03 tests BDF_SUM.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 10
!
  double precision diftab(maxorder)
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) &
    '  BDFP_SET sets up Backward Difference Predictor quadrature;'
  write ( *, '(a)' ) '  BDF_SUM carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integration interval is [0,1].'
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call bdfp_set ( norder, weight, xtab )
 
        call bdf_sum ( func, norder, weight, xtab, diftab, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests CHEB_SET.
!! TEST04 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 9
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1

  xlo = -1.0D+00
  xhi =  1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  CHEB_SET sets up Chebyshev quadrature;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      if ( norder /= 8 ) then

        do i = ilo, ihi

          call func_set ( 'SET', i )

          call cheb_set ( norder, xtab, weight )

          call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
            result(i) )
   
        end do

        write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

      end if

    end do
 
  end do

  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 tests CHEB_TC_SET.
!! TEST05 tests SUMMER.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 6
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  a = - 1.0D+00
  b =   1.0D+00
  nsub = 1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  CHEB_TC_SET sets up Gauss-Chebyshev quadrature,'
  write ( *, '(a)' ) '  ( closed, first kind );'
  write ( *, '(a)' ) '  SUMMER carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is 1 / sqrt ( 1 - X**2 )'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 2, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call cheb_tc_set ( norder, xtab, weight )

        call summer ( func, norder, xtab, weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test06
!
!*******************************************************************************
!
!! TEST06 tests CHEB_TO_SET.
!! TEST06 tests SUMMER.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 6
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  a = -1.0D+00
  b =  1.0D+00
  nsub = 1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST06'
  write ( *, '(a)' ) '  CHEB_TO_SET sets up Gauss-Chebyshev quadrature,'
  write ( *, '(a)' ) '  ( open, first kind );'
  write ( *, '(a)' ) '  SUMMER carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is 1 / sqrt ( 1 - X**2 )'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '
    
    do norder = 1, maxorder
 
      do i = ilo, ihi
 
        call func_set ( 'SET', i )

        call cheb_to_set ( norder, xtab, weight )

        call summer ( func, norder, xtab, weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test07
!
!*******************************************************************************
!
!! TEST07 tests CHEB_U_SET.
!! TEST07 tests SUMMER.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 4
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  a = - 1.0D+00
  b =   1.0D+00
  nsub = 1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST07'
  write ( *, '(a)' ) '  CHEB_U_SET sets up Gauss-Chebyshev quadrature;'
  write ( *, '(a)' ) '  SUMMER carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is sqrt ( 1 - X**2 )'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '
 
    do norder = 1, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call cheb_u_set ( norder, xtab, weight )

        call summer ( func, norder, xtab, weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do
 
  end do

  return
end
subroutine test08
!
!*******************************************************************************
!
!! TEST08 tests HERMITE_COM.
!! TEST08 tests SUMMER.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST08'
  write ( *, '(a)' ) '  HERMITE_COM computes a Gauss-Hermite rule;'
  write ( *, '(a)' ) '  SUMMER carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The integration interval is ( -Infinity, +Infinity ).'

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '
 
    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call hermite_com ( norder, xtab, weight )
 
        call summer ( func, norder, xtab, weight, result(i) )
 
      end do
 
      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test09
!
!*******************************************************************************
!
!! TEST09 tests HERMITE_SET.
!! TEST09 tests SUMMER.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST09'
  write ( *, '(a)' ) '  HERMITE_SET sets up Gauss-Hermite quadrature;'
  write ( *, '(a)' ) '  SUMMER carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The integration interval is ( -Infinity, +Infinity ).'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '
 
    do norder = 1, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call hermite_set ( norder, xtab, weight )
 
        call summer ( func, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' )  norder, result(ilo:ihi)
 
    end do
 
  end do

  return
end
subroutine test10
!
!*******************************************************************************
!
!! TEST10 tests JACOBI_COM.
!! TEST10 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 10
!
  double precision a
  double precision alpha
  double precision b
  double precision beta
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer k
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  a = -1.0D+00
  b =  1.0D+00

  nsub = 1
  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST10'
  write ( *, '(a)' ) '  JACOBI_COM sets up Gauss-Jacobi quadrature;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '
  do k = 1, 2

    if ( k == 1 ) then
      alpha = 0.0D+00
      beta = 0.0D+00
    else if ( k == 2 ) then
      alpha = 1.0D+00
      beta = 0.0D+00
    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6)' ) '  ALPHA = ', alpha
    write ( *, '(a,g14.6)' ) '  BETA =  ', beta

    do ilo = 1, nfunc, 5

      ihi = min ( ilo + 4, nfunc )

      write ( *, '(a)' ) ' '
      write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
      write ( *, '(a)' ) ' '

      do norder = 1, maxorder

        do i = ilo, ihi

          call func_set ( 'SET', i )

          call jacobi_com ( norder, xtab, weight, alpha, beta )

          call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
            result(i) )

        end do

        write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

      end do

    end do

  end do 

  return
end
subroutine test11
!
!*******************************************************************************
!
!! TEST11 tests KRONROD_SET.
!! TEST11 tests LEGENDRE_SET.
!! TEST11 tests SUMMER_GK.
!
  implicit none
!
  integer, parameter :: norderg = 10
  integer, parameter :: norderk = 2 * norderg + 1
!
  double precision, external :: fx2sd1
  double precision resultg
  double precision resultk
  double precision weightg(norderg)
  double precision weightk(norderk)
  double precision xtabg(norderg)
  double precision xtabk(norderk)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST11'
  write ( *, '(a)' ) '  KRONROD_SET sets up Kronrod quadrature;'
  write ( *, '(a)' ) '  LEGENDRE_SET sets up Gauss-Legendre quadrature;'
  write ( *, '(a)' ) '  SUMMER_GK carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integration interval is [-1, 1].'
  write ( *, '(a)' ) '  Integrand is X**2 / SQRT ( 1.1 - X**2 ).'
  write ( *, '(a)' ) ' '

  call legendre_set ( norderg, xtabg, weightg )

  call kronrod_set ( norderk, xtabk, weightk )

  call summer_gk ( fx2sd1, norderg, weightg, resultg, norderk, xtabk, &
    weightk, resultk )

  write ( *, '(i2,2x,g16.10)' ) norderg, resultg
  write ( *, '(i2,2x,g16.10)' ) norderk, resultk
  write ( *, '(2x,2x,g16.10)' )          resultg-resultk

  return
end
subroutine test12
!
!*******************************************************************************
!
!! TEST12 tests KRONROD_SET.
!! TEST12 tests LEGENDRE_SET.
!! TEST12 tests SUM_SUB_GK.
!
  implicit none
!
  integer, parameter :: norderg = 7
  integer, parameter :: norderk = 2 * norderg + 1
!
  double precision a
  double precision b
  double precision error
  double precision, external :: fx2sd1
  integer nsub
  double precision resultg
  double precision resultk
  double precision weightg(norderg)
  double precision weightk(norderk)
  double precision xtabg(norderg)
  double precision xtabk(norderk)
!
  a = - 1.0D+00
  b =   1.0D+00
  nsub = 5

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST12'
  write ( *, '(a)' ) '  KRONROD_SET sets up Kronrod quadrature;'
  write ( *, '(a)' ) '  LEGENDRE_SET sets up Gauss-Legendre quadrature;'
  write ( *, '(a)' ) '  SUM_SUB_GK carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Integrand is X**2 / SQRT ( 1.1 - X**2 ).'
  write ( *, '(a)' ) ' '

  call legendre_set ( norderg, xtabg, weightg )

  call kronrod_set ( norderk, xtabk, weightk )

  call sum_sub_gk ( fx2sd1, a, b, nsub, norderg, weightg, resultg, &
    norderk, xtabk, weightk, resultk, error )

  write ( *, '(i2,2x,g16.10)' ) norderg, resultg
  write ( *, '(i2,2x,g16.10)' ) norderk, resultk
  write ( *, '(2x,2x,g16.10)' )          error

  return
end
subroutine test13
!
!*******************************************************************************
!
!! TEST13 tests LAGUERRE_COM.
!! TEST13 tests LAGUERRE_SUM.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision alpha
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ii
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST13'
  write ( *, '(a)' ) '  LAGUERRE_COM computes a Gauss-Laguerre rule;'
  write ( *, '(a)' ) '  LAGUERRE_SUM carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is EXP ( - X ).'
  write ( *, '(a)' ) ' '
!
  a = 1.0D+00
  alpha = 0.0D+00

  write ( *, '(a,g14.6,a)' ) '  The integration interval is [ ', &
    a, ', +Infinity ).'
  write ( *, '(a,g14.6)' ) '  ALPHA = ', alpha
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call laguerre_com ( norder, xtab, weight, alpha )
 
        call laguerre_sum ( func, a, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)
 
    end do

  end do

  return
end
subroutine test14
!
!*******************************************************************************
!
!! TEST14 tests LAGUERRE_COM.
!! TEST14 tests LAGUERRE_SUM.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision alpha
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ii
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST14'
  write ( *, '(a)' ) '  LAGUERRE_COM computes a Gauss-Laguerre rule;'
  write ( *, '(a)' ) '  LAGUERRE_SUM carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is EXP ( - X ).'
  write ( *, '(a)' ) ' '
!
  a = 0.0D+00
  alpha = 0.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a)' ) '  The integration interval is [ ',  &
    a, ', +Infinity ).'
  write ( *, '(a,g14.6)' ) '  ALPHA = ', alpha
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, 10

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call laguerre_com ( norder, xtab, weight, alpha )
 
        call laguerre_sum ( func, a, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)
 
    end do

  end do

  return
end
subroutine test15
!
!*******************************************************************************
!
!! TEST15 tests LAGUERRE_COM.
!! TEST15 tests LAGUERRE_SUM.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision alpha
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ii
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST15'
  write ( *, '(a)' ) '  LAGUERRE_COM computes a Gauss-Laguerre rule;'
  write ( *, '(a)' ) '  LAGUERRE_SUM carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is EXP ( - X ).'
  write ( *, '(a)' ) ' '
!
  a = 0.0D+00
  alpha = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a)' ) '  The integration interval is [ ', &
    a, ', +Infinity ).'
  write ( *, '(a,g14.6)' ) '  ALPHA = ', alpha
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5
  
    ihi = min ( ilo + 4, nfunc )
    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, 10

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call laguerre_com ( norder, xtab, weight, alpha )
 
        call laguerre_sum ( func, a, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5g14.6)' ) norder, result(ilo:ihi)
 
    end do

  end do

  return
end
subroutine test16
!
!*******************************************************************************
!
!! TEST16 tests LAGUERRE_COM.
!! TEST16 tests LAGUERRE_SUM.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision alpha
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ii
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST16'
  write ( *, '(a)' ) '  LAGUERRE_COM computes a Gauss-Laguerre rule;'
  write ( *, '(a)' ) '  LAGUERRE_SUM carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is EXP ( - X ).'
  write ( *, '(a)' ) ' '
!
  a = 0.0D+00
  alpha = 2.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a)' ) '  The integration interval is [ ', &
    a, ', +Infinity ).'
  write ( *, '(a,g14.6)' ) '  ALPHA = ', alpha
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, 10

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call laguerre_com ( norder, xtab, weight, alpha )
 
        call laguerre_sum ( func, a, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5g14.6)' ) norder, result(ilo:ihi)
 
    end do

  end do

  return
end
subroutine test17
!
!*******************************************************************************
!
!! TEST17 tests LAGUERRE_SET.
!! TEST17 tests LAGUERRE_SUM.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  a = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST17'
  write ( *, '(a)' ) '  LAGUERRE_SET sets up Gauss-Laguerre quadrature;'
  write ( *, '(a)' ) '  LAGUERRE_SUM carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6,a)' ) '  The integration interval is [ ', &
    a, ', +Infinity ).'
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) '  The weight function is EXP ( - X ).'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call laguerre_set ( norder, xtab, weight )
 
        call laguerre_sum ( func, a, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)
 
    end do

  end do
 
  return
end
subroutine test18
!
!*******************************************************************************
!
!! TEST18 compares LEGENDRE_COM.
!! TEST18 compares LEGENDRE_SET.
!
  implicit none
!
  integer, parameter :: norder = 19
!
  integer i
  integer iwdifmax
  integer ixdifmax
  double precision wdifmax
  double precision weight1(norder)
  double precision weight2(norder)
  double precision xdifmax
  double precision xtab1(norder)
  double precision xtab2(norder)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST16'
  write ( *, '(a)' ) '  LEGENDRE_COM computes Gauss-Legendre data;'
  write ( *, '(a)' ) '  LEGENDRE_SET looks up the same data.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Compare the data for NORDER = ', norder

  call legendre_com ( norder, xtab1, weight1 )
  call legendre_set ( norder, xtab2, weight2 )
 
  xdifmax = 0.0D+00
  ixdifmax = 0

  wdifmax = 0.0D+00
  iwdifmax = 0

  do i = 1, norder

    if ( abs ( xtab1(i) - xtab2(i) ) > xdifmax ) then
      xdifmax = abs ( xtab1(i) - xtab2(i) )
      ixdifmax = i
    end if

    if ( abs ( weight1(i) - weight2(i) ) > wdifmax ) then
      wdifmax = abs ( weight1(i) - weight2(i) )
      iwdifmax = i
    end if

  end do

  if ( ixdifmax /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,g25.18)' ) '  Maximum abscissa difference is ', xdifmax
    write ( *, '(a,i6)' ) '  for index I = ', ixdifmax
    write ( *, '(a,g25.18)' ) 'Computed:', xtab1(ixdifmax)
    write ( *, '(a,g25.18)' ) 'Stored:  ', xtab2(ixdifmax)
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  The computed and stored abscissas are identical.'
  end if

  if ( iwdifmax /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,g25.18)' ) '  Maximum weight difference is   ', wdifmax
    write ( *, '(a,i6)' ) '  for index I = ', iwdifmax
    write ( *, '(a,g25.18)' ) 'Computed:', weight1(iwdifmax)
    write ( *, '(a,g25.18)' ) 'Stored:  ', weight2(iwdifmax)
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  The computed and stored weights are identical.'
  end if
 
  return
end
subroutine test19
!
!*******************************************************************************
!
!! TEST19 tests LEGENDRE_COM.
!! TEST19 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: norder = 2
!
  double precision a
  double precision b
  double precision, external :: fx2sd1
  integer iexp
  integer nsub
  double precision result
  double precision weight(norder)
  double precision xhi
  double precision xlo
  double precision xtab(norder)
!
  a = 0.0D+00
  b = 1.0D+00

  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST19'
  write ( *, '(a)' ) '  LEGENDRE_COM computes a Gauss-Legendre rule;'
  write ( *, '(a)' ) '  SUM_SUB carries it out over subintervals.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Here, we use a fixed order NORDER = ', norder
  write ( *, '(a)' ) '  and use more and more subintervals.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  NSUB     Integral'
  write ( *, '(a)' ) ' '
 
  call legendre_com ( norder, xtab, weight )
 
  do iexp = 0, 9

    nsub = 2**iexp

    call sum_sub ( fx2sd1, a, b, nsub, norder, xlo, xhi, xtab, weight, result )

    write ( *, '(i4,g16.8)' ) nsub, result

  end do
 
  return
end
subroutine test20
!
!*******************************************************************************
!
!! TEST20 tests LEGENDRE_COM.
!! TEST20 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 10
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1

  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST20'
  write ( *, '(a)' ) '  LEGENDRE_COM sets up Gauss-Legendre quadrature;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_com ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test21
!
!*******************************************************************************
!
!! TEST21 tests LEGENDRE_SET.
!! TEST21 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1

  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST21'
  write ( *, '(a)' ) '  LEGENDRE_SET sets up Gauss-Legendre quadrature;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, &
          weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do
 
  end do

  return
end
subroutine test22
!
!*******************************************************************************
!
!! TEST22 tests LEGENDRE_SET.
!! TEST22 tests LEGENDRE_SET_X0_01
!! TEST22 tests RULE_ADJUST.
!
  implicit none
!
  integer, parameter :: norder = 5
!
  double precision a
  double precision b
  double precision c
  double precision d
  integer i
  double precision weight1(norder)
  double precision weight2(norder)
  double precision weight3(norder)
  double precision xtab1(norder)
  double precision xtab2(norder)
  double precision xtab3(norder)
!
  a = -1.0D+00
  b = +1.0D+00
  c =  0.0D+00
  d =  1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST22'
  write ( *, '(a)' ) '  LEGENDRE_SET sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating F(X) over [-1,1];'
  write ( *, '(a)' ) '  RULE_ADJUST adjusts a rule for a new interval.'
  write ( *, '(a)' ) '  LEGENDRE_SET_X0_01 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating F(X) over [0,1];'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We will use LEGENDRE_SET to get a rule for [-1,1],'
  write ( *, '(a)' ) '  adjust it to [0,1] using RULE_ADJUST,'
  write ( *, '(a)' ) '  and compare the results of LEGENDRE_SET_X0_01.'
  write ( *, '(a)' ) ' '

  call legendre_set ( norder, xtab1, weight1 )

  xtab2(1:norder) = xtab1(1:norder)
  weight2(1:norder) = weight1(1:norder)

  call rule_adjust ( a, b, c, d, norder, xtab2, weight2 )

  call legendre_set_x0_01 ( norder, xtab3, weight3 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Abscissas:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Original  Adjusted Stored'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i2,3f16.12)' ) i, xtab1(i), xtab2(i), xtab3(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Weights:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Original  Adjusted Stored'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i2,3f16.12)' ) i, weight1(i), weight2(i), weight3(i)
  end do

  return
end
subroutine test23
!
!*******************************************************************************
!
!! TEST23 tests LEGENDRE_SET_COS.
!! TEST23 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer iexp
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision d_pi
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  a = -0.5D+00 * d_pi ( )
  b = +0.5D+00 * d_pi ( )

  nsub = 1

  xlo = -0.5D+00 * d_pi ( )
  xhi = +0.5D+00 * d_pi ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST23'
  write ( *, '(a)' ) '  LEGENDRE_SET_COS sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    over [-PI/2,PI/2] with weight function COS(X);'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'

  call func_set ( 'COUNT', nfunc )

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do iexp = 0, 4

      norder = 2**iexp

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_cos ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, &
          weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test235
!
!*******************************************************************************
!
!! TEST235 tests LEGENDRE_SET_SQRTX_01.
!! TEST235 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer iexp
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  a = 0.0D+00
  b = 1.0D+00

  nsub = 1

  xlo = 0.0D+00
  xhi = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST235'
  write ( *, '(a)' ) '  LEGENDRE_SET_SQRTX_01 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    over [0,1] with weight function SQRT(X);'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'

  call func_set ( 'COUNT', nfunc )

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do iexp = 0, 3

      norder = 2**iexp

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_sqrtx_01 ( norder, xtab, weight )

        call summer ( func, norder, xtab, weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test236
!
!*******************************************************************************
!
!! TEST236 tests LEGENDRE_SET_SQRTX2_01.
!! TEST236 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer iexp
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  a = 0.0D+00
  b = 1.0D+00

  nsub = 1

  xlo = 0.0D+00
  xhi = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST236'
  write ( *, '(a)' ) '  LEGENDRE_SET_SQRTX2_01 sets up Gauss-Legendre'
  write ( *, '(a)' ) '    quadrature over [0,1] with weight function 1/SQRT(X);'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'

  call func_set ( 'COUNT', nfunc )

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do iexp = 0, 3

      norder = 2**iexp

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_sqrtx2_01 ( norder, xtab, weight )

        call summer ( func, norder, xtab, weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test24
!
!*******************************************************************************
!
!! TEST24 tests LEGENDRE_SET_COS2.
!! TEST24 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer iexp
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision d_pi
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  a = -0.5D+00 * d_pi ( )
  b = +0.5D+00 * d_pi ( )

  nsub = 1
  xlo = -0.5D+00 * d_pi ( )
  xhi = +0.5D+00 * d_pi ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST24'
  write ( *, '(a)' ) '  LEGENDRE_SET_COS2 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    over [0,PI/2] with weight function COS(X);'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'

  call func_set ( 'COUNT', nfunc )

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do iexp = 1, 4

      norder = 2**iexp

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_cos2 ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, &
          weight, result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test25
!
!*******************************************************************************
!
!! TEST25 tests LEGENDRE_SET_LOG.
!! TEST25 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
  integer, parameter :: ntest = 9
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer iexp
  integer ihi
  integer ilo
  integer itest
  integer nfunc
  integer norder
  integer, dimension ( ntest ) :: norder_test = (/ 1, 2, 3, 4, 5, 6, 7, 8, 16 /)
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  a = 0.0D+00
  b = 1.0D+00

  nsub = 1

  xlo = 0.0D+00
  xhi = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST25'
  write ( *, '(a)' ) '  LEGENDRE_SET_LOG sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    over [0,1] with weight function -LOG(X);'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'

  call func_set ( 'COUNT', nfunc )

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do itest = 1, ntest

      norder = norder_test(itest)
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_log ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder,  xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do
 
  end do

  return
end
subroutine test26
!
!*******************************************************************************
!
!! TEST26 tests LEGENDRE_SET_X0_01.
!! TEST26 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 8
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1

  xlo =  0.0D+00
  xhi =  1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST26'
  write ( *, '(a)' ) '  LEGENDRE_SET_X0_01 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating F(X) over [0,1];'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_x0_01 ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test27
!
!*******************************************************************************
!
!! TEST27 tests LEGENDRE_SET_X1.
!! TEST27 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 9
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1

  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST27'
  write ( *, '(a)' ) '  LEGENDRE_SET_X1 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating ( 1 + X ) * F(X) over [-1,1];'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_x1 ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test28
!
!*******************************************************************************
!
!! TEST28 tests LEGENDRE_SET_X1.
!! TEST28 tests LEGENDRE_SET_X1_01
!! TEST28 tests RULE_ADJUST.
!
  implicit none
!
  integer, parameter :: norder = 5
!
  double precision a
  double precision b
  double precision c
  double precision d
  integer i
  double precision weight1(norder)
  double precision weight2(norder)
  double precision weight3(norder)
  double precision xtab1(norder)
  double precision xtab2(norder)
  double precision xtab3(norder)
!
  a = -1.0D+00
  b = +1.0D+00
  c =  0.0D+00
  d =  1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST28'
  write ( *, '(a)' ) '  LEGENDRE_SET_X1 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating ( 1 + X ) * F(X) over [-1,1];'
  write ( *, '(a)' ) '  RULE_ADJUST adjusts a rule for a new interval.'
  write ( *, '(a)' ) '  LEGENDRE_SET_X1_01 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating X * F(X) over [0,1];'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We will use LEGENDRE_SET_X1 to get a rule for [-1,1],'
  write ( *, '(a)' ) '  adjust it to [0,1] using RULE_ADJUST,'
  write ( *, '(a)' ) '  make further adjustments because the weight function'
  write ( *, '(a)' ) '  is not 1, '
  write ( *, '(a)' ) '  and compare the results of LEGENDRE_SET_X1_01.'
  write ( *, '(a)' ) ' '

  call legendre_set_x1 ( norder, xtab1, weight1 )

  xtab2(1:norder) = xtab1(1:norder)
  weight2(1:norder) = weight1(1:norder)

  call rule_adjust ( a, b, c, d, norder, xtab2, weight2 )
!
!  Because the weight function W(X) is not 1, we need to do more
!  adjustments to the weight vector.
!
  weight2(1:norder) = weight2(1:norder) / 2.0D+00

  call legendre_set_x1_01 ( norder, xtab3, weight3 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Abscissas:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Original        Adjusted         Stored'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i2,3f16.12)' ) i, xtab1(i), xtab2(i), xtab3(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)') '  Weights:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Original  Adjusted Stored'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i2,3f16.12)' ) i, weight1(i), weight2(i), weight3(i)
  end do

  return
end
subroutine test29
!
!*******************************************************************************
!
!! TEST29 tests LEGENDRE_SET_X1_01.
!! TEST29 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 8
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1

  xlo = 0.0D+00
  xhi = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST29'
  write ( *, '(a)' ) '  LEGENDRE_SET_X1_01 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating X * F(X) over [0,1];'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_x1_01 ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test30
!
!*******************************************************************************
!
!! TEST30 tests LEGENDRE_SET_X2.
!! TEST30 tests SUM_SUB
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 9
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1
  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST30'
  write ( *, '(a)' ) '  LEGENDRE_SET_X2 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating (1+X)**2 * F(X) over [-1,1];'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_x2 ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test31
!
!*******************************************************************************
!
!! TEST31 tests LEGENDRE_SET_X2.
!! TEST31 tests LEGENDRE_SET_X2_01
!! TEST31 tests RULE_ADJUST.
!
  implicit none
!
  integer, parameter :: norder = 5
!
  double precision a
  double precision b
  double precision c
  double precision d
  integer i
  double precision weight1(norder)
  double precision weight2(norder)
  double precision weight3(norder)
  double precision xtab1(norder)
  double precision xtab2(norder)
  double precision xtab3(norder)
!
  a = -1.0D+00
  b = +1.0D+00
  c =  0.0D+00
  d =  1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST31'
  write ( *, '(a)' ) '  LEGENDRE_SET_X2 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating ( 1 + X )^2 * F(X) over [-1,1];'
  write ( *, '(a)' ) '  RULE_ADJUST adjusts a rule for a new interval.'
  write ( *, '(a)' ) '  LEGENDRE_SET_X2_01 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating X^2 * F(X) over [0,1];'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  We will use LEGENDRE_SET_X2 to get a rule for [-1,1],'
  write ( *, '(a)' ) '  adjust it to [0,1] using RULE_ADJUST,'
  write ( *, '(a)' ) '  make further adjustments because the weight function'
  write ( *, '(a)' ) '  is not 1, '
  write ( *, '(a)' ) '  and compare the results of LEGENDRE_SET_X2_01.'
  write ( *, '(a)' ) ' '

  call legendre_set_x2 ( norder, xtab1, weight1 )

  xtab2(1:norder) = xtab1(1:norder)
  weight2(1:norder) = weight1(1:norder)

  call rule_adjust ( a, b, c, d, norder, xtab2, weight2 )
!
!  Because the weight function W(X) is not 1, we need to do more
!  adjustments to the weight vector.
!
  weight2(1:norder) = weight2(1:norder) / 4.0D+00

  call legendre_set_x2_01 ( norder, xtab3, weight3 )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Abscissas:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Original  Adjusted Stored'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i2,3f16.12)' ) i, xtab1(i), xtab2(i), xtab3(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Weights:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Original  Adjusted Stored'
  write ( *, '(a)' ) ' '

  do i = 1, norder
    write ( *, '(i2,3f16.12)' ) i, weight1(i), weight2(i), weight3(i)
  end do

  return
end
subroutine test32
!
!*******************************************************************************
!
!! TEST32 tests LEGENDRE_SET_X2_01.
!! TEST32 tests SUM_SUB
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 8
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1
  xlo = 0.0D+00
  xhi = 1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST32'
  write ( *, '(a)' ) '  LEGENDRE_SET_X2_01 sets up Gauss-Legendre quadrature'
  write ( *, '(a)' ) '    for integrating X*X * F(X) over [0,1];'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '
  
    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call legendre_set_x2_01 ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test33
!
!*******************************************************************************
!
!! TEST33 tests LOBATTO_SET.
!! TEST33 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  a = -1.0D+00
  b =  1.0D+00

  nsub = 1
  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST33'
  write ( *, '(a)' ) '  LOBATTO_SET sets up Lobatto quadrature;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '
  
    do norder = 2, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call lobatto_set ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test34
!
!*******************************************************************************
!
!! TEST34 tests MOULTON_SET.
!! TEST34 tests SUMMER.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 10
!
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST34'
  write ( *, '(a)' ) '  MOULTON_SET sets up Adams-Moulton quadrature;'
  write ( *, '(a)' ) '  SUMMER carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Integration interval is [0,1].'
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call moulton_set ( norder, xtab, weight )
 
        call summer ( func, norder, xtab, weight, result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do

  return
end
subroutine test35
!
!*******************************************************************************
!
!! TEST35 tests NCC_SET.
!! TEST35 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1

  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST35'
  write ( *, '(a)' ) '  NCC_SET sets up a closed Newton Cotes rule;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo,ihi )
    write ( *, '(a)' ) ' '

    do norder = 2, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call ncc_set ( norder, xtab, weight )
 
        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do
 
  end do

  return
end
subroutine test36
!
!*******************************************************************************
!
!! TEST36 tests NCC_COM.
!! TEST36 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 20
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1
  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST36'
  write ( *, '(a)' ) '  NCC_COM computes a closed Newton Cotes rule;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '
  
    do norder = 2, maxorder
 
      do i = ilo, ihi

        call func_set ( 'SET', i )

        call ncc_com ( norder, xtab, weight )
 
        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )
 
      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do

  end do
 
  return
end
subroutine test37
!
!*******************************************************************************
!
!! TEST37 tests NCO_SET.
!! TEST37 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 9
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1
  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST37'
  write ( *, '(a)' ) '  NCO_SET sets up an open Newton-Cotes rule;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)') '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder
 
      if ( norder <= 7 .or. norder == 9 ) then

        do i = ilo, ihi

          call func_set ( 'SET', i )

          call nco_set ( norder, xtab, weight )
 
          call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
            result(i) )
 
        end do

        write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

      end if

    end do

  end do
 
  return
end
subroutine test38
!
!*******************************************************************************
!
!! TEST38 tests NCO_COM.
!! TEST38 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 9
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1
  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST38'
  write ( *, '(a)' ) '  NCO_COM sets up an open Newton-Cotes rule;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  Integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder
 
      if ( norder <= 7 .or. norder == 9 ) then

        do i = ilo, ihi

          call func_set ( 'SET', i )

          call nco_com ( norder, xtab, weight )
 
          call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, & 
            result(i) )
 
        end do

        write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

      end if

    end do

  end do
 
  return
end
subroutine test39
!
!*******************************************************************************
!
!! TEST39 tests RADAU_SET.
!! TEST39 tests SUM_SUB.
!
  implicit none
!
  integer, parameter :: maxfunc = 20
  integer, parameter :: maxorder = 15
!
  double precision a
  double precision b
  character ( len = 10 ) fname
  double precision, external :: func
  integer i
  integer ihi
  integer ilo
  integer nfunc
  integer norder
  integer nsub
  double precision result(maxfunc)
  double precision weight(maxorder)
  double precision xhi
  double precision xlo
  double precision xtab(maxorder)
!
  call func_set ( 'COUNT', nfunc )

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

  nsub = 1
  xlo = -1.0D+00
  xhi = +1.0D+00

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST39'
  write ( *, '(a)' ) '  RADAU_SET sets up Radau quadrature;'
  write ( *, '(a)' ) '  SUM_SUB carries it out.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,2g14.6)' ) '  The integration interval is ', a, b
  write ( *, '(a,i6)' ) '  Number of subintervals is ', nsub
  write ( *, '(a)' ) '  Quadrature order will vary.'
  write ( *, '(a)' ) '  Integrand will vary.'
  write ( *, '(a)' ) ' '

  do ilo = 1, nfunc, 5

    ihi = min ( ilo + 4, nfunc )

    write ( *, '(a)' ) ' '
    write ( *, '(4x, 5a14 )' ) ( fname(i), i = ilo, ihi )
    write ( *, '(a)' ) ' '

    do norder = 1, maxorder

      do i = ilo, ihi

        call func_set ( 'SET', i )

        call radau_set ( norder, xtab, weight )

        call sum_sub ( func, a, b, nsub, norder, xlo, xhi, xtab, weight, &
          result(i) )

      end do

      write ( *, '(i2,2x,5f14.8)' ) norder, result(ilo:ihi)

    end do
 
  end do

  return
end
function f1sd1 ( x )
!
!*******************************************************************************
!
!! F1SD1 evaluates the function 1.0D+00/ sqrt ( 1.1 - x**2 ).
!
!
!  Parameters:
!
!    Input, double precision X, the argument of the function.
!
!    Output, double precision F1SD1, the value of the function.
!
  implicit none
!
  double precision f1sd1
  double precision x
!
  f1sd1 = 1.0D+00 / sqrt ( 1.1D+00 - x**2 )
 
  return
end
function fxsd1 ( x )
!
!*******************************************************************************
!
!! FXSD1 evaluates the function x / sqrt ( 1.1 - x**2 ).
!
!
!  Parameters:
!
!    Input, double precision X, the argument of the function.
!
!    Output, double precision FXSD1, the value of the function.
!
  implicit none
!
  double precision fxsd1
  double precision x
!
  fxsd1 = x / sqrt ( 1.1D+00 - x**2 )
 
  return
end
function fx2sd1 ( x )
!
!*******************************************************************************
!
!! FX2SD1 evaluates the function x**2 / sqrt ( 1.1 - x**2 ).
!
!
!  Parameters:
!
!    Input, double precision X, the argument of the function.
!
!    Output, double precision FX2SD1, the value of the function.
!
  implicit none
!
  double precision fx2sd1
  double precision x
!
  fx2sd1 = x**2 / sqrt ( 1.1D+00 - x**2 )
 
  return
end
function func ( x )
!
!*******************************************************************************
!
!! FUNC evaluates a function of X, as chosen by the user.
!
!
!  Modified:
!
!    21 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision X, the argument of the function.
!
!    Output, double precision FUNC, the value of the function.
!
  implicit none
!
  double precision func
  integer ifunc
  double precision x
!
  call func_set ( 'GET', ifunc )

  if ( ifunc == 1 ) then
    func = 1.0D+00
  else if ( ifunc == 2 ) then
    func = x
  else if ( ifunc == 3 ) then
    func = x**2
  else if ( ifunc == 4 ) then
    func = x**3
  else if ( ifunc == 5 ) then
    func = x**4
  else if ( ifunc == 6 ) then
    func = x**5
  else if ( ifunc == 7 ) then
    func = x**6
  else if ( ifunc == 8 ) then
    func = x**7
  else if ( ifunc == 9 ) then
    func = sin ( x )
  else if ( ifunc == 10 ) then
    func = exp ( x )
  else if ( ifunc == 11 ) then
    func = sqrt ( abs ( x ) )
  else
    func = 0.0D+00
  end if

  return
end
subroutine func_set ( action, i )
!
!*******************************************************************************
!
!! FUNC_SET sets the function to be returned by FUNC.
!
!
!  Modified:
!
!    21 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION, the action to be carried out.
!    'COUNT' means the call is made to count the number of functions available.
!    'GET' means the call is made to find out the current function index.
!    'SET' means the call is made to set the current function index.
!
!    Input/output, integer I.
!    For "COUNT', I is the number of functions available;
!    For 'GET', I is the currently chosen function;
!    For 'SET', I is the user's new choice for the function.
!
  implicit none
!
  character ( len = * ) action
  integer i
  integer, save :: ival = 0
!
  if ( action == 'COUNT' ) then
    i = 11
  else if ( action == 'GET' ) then
    i = ival
  else if ( action == 'SET' ) then
    ival = i
  end if

  return
end
function fname ( ifunc )
!
!*******************************************************************************
!
!! FNAME returns the name of the function that will be evaluated in FUNC.
!
!
!  Modified:
!
!    21 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IFUNC, the index of the function.
!
!    Output, character ( len = 10 ) FNAME, the name of the function.
!
  implicit none
!
  character ( len = 10 ) fname
  integer ifunc
!
  if ( ifunc == 1 ) then
    fname = '         1'
  else if ( ifunc == 2 ) then
    fname = '         X'
  else if ( ifunc == 3 ) then
    fname = '      X**2'
  else if ( ifunc == 4 ) then
    fname = '      X**3'
  else if ( ifunc == 5 ) then
    fname = '      X**4'
  else if ( ifunc == 6 ) then
    fname = '      X**5'
  else if ( ifunc == 7 ) then
    fname = '      X**6'
  else if ( ifunc == 8 ) then
    fname = '      X**7'
  else if ( ifunc == 9 ) then
    fname = '    SIN(X)'
  else if ( ifunc == 10 ) then
    fname = '    EXP(X)'
  else if ( ifunc == 11 ) then
    fname = ' SQRT(|X|)'
  else
    fname = '??????'
  end if

  return
end