SANDIA_RULES
Quadrature Rules of Gaussian Type

SANDIA_RULES is a FORTRAN90 library which generates a variety of quadrature rules of various orders and types.

This library is used, in turn, by several other libraries, including SPARSE_GRID_MIXED and SGMGA. This means that a program that calls any of those libraries must have access to a compiled copy of SANDIA_RULES as well.

Name	Usual domain	Weight function
Gauss-Legendre	[-1,+1]	1
Clenshaw-Curtis	[-1,+1]	1
Fejer Type 2	[-1,+1]	1
Gauss-Chebyshev 1	[-1,+1]	1/sqrt(1-x2)
Gauss-Chebyshev 2	[-1,+1]	sqrt(1-x2)
Gauss-Gegenbauer	[-1,+1]	(1-x2)alpha
Gauss-Jacobi	[-1,+1]	(1-x)alpha (1+x)beta
Gauss-Laguerre	[0,+oo)	e-x
Generalized Gauss-Laguerre	[0,+oo)	xalpha e-x
Gauss-Hermite	(-oo,+oo)	e-x*x
Generalized Gauss-Hermite	(-oo,+oo)	|x|alpha e-x*x
Hermite Genz-Keister	(-oo,+oo)	e-x*x
Newton-Cotes-Closed	[-1,+1]	1
Newton-Cotes-Open	[-1,+1]	1
Newton-Cotes-Open-Half	[-1,+1]	1

For example, a Gauss-Gegenbauer quadrature rule is used to approximate:

        Integral ( -1 <= x <= +1 ) f(x) (1-x^2)^alpha dx
      

The approximation to the integral is formed by computing a weighted sum of function values at specific points:

        Sum ( 1 <= i <= n ) w(i) * f(x(i))
      

As a matter of convenience, most of the quadrature rules are available through three related functions:

• name_COMPUTE returns points X and weights W;
• name_COMPUTE_POINTS returns points X;
• name_COMPUTE_WEIGHTS returns weights W;

Some of these quadrature rules expect a parameter ALPHA, and perhaps also a parameter BETA, in order to fully define the rule. Therefore, the argument lists of these functions vary. They always include the input quantity ORDER, but may have one or two additional inputs. In order to offer a uniform interface, there is also a family of functions with a standard set of input arguments, ORDER, NP, and P. Here NP is parameter counter, and P is the parameter value vector P. Using this interface, it is possible to call all the quadrature functions with the same argument list. The uniform interface functions can be identified by the suffix _NP that appears in their names. Generally, these functions "unpack" the parameter vector where needed, and then call the corresponding basic function. Of course, for many rules NP is zero and P may be a null pointer.

• name_COMPUTE_NP ( ORDER, NP, P, X, W ) unpacks parameters, calls name_COMPUTE, returns points X and weights W;
• name_COMPUTE_POINTS_NP ( ORDER, NP, P, X ) unpacks parameters, calls name_COMPUTE_POINTS, returns points X;
• name_COMPUTE_WEIGHTS_NP ( ORDER, NP, P, W ) unpacks parameters, calls name_COMPUTE_WEIGHTS, returns weights W;

Licensing:

The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.

Languages:

SANDIA_RULES is available in a C version and a C++ version and a FORTRAN90 version and a MATLAB version.

Related Data and Programs:

QUADRULE, a FORTRAN90 library which defines quadrature rules for 1-dimensional domains.

R8LIB, a FORTRAN90 library which contains many utility routines, using "R8" or "double precision real" arithmetic.

SGMGA, a FORTRAN90 library which creates sparse grids based on a mixture of 1D quadrature rules, allowing anisotropic weights for each dimension.

Reference:

1. Milton Abramowitz, Irene Stegun,
   Handbook of Mathematical Functions,
   National Bureau of Standards, 1964,
   ISBN: 0-486-61272-4,
   LC: QA47.A34.
2. William Cody,
   An Overview of Software Development for Special Functions,
   in Numerical Analysis Dundee, 1975,
   edited by GA Watson,
   Lecture Notes in Mathematics 506,
   Springer, 1976.
3. Philip Davis, Philip Rabinowitz,
   Methods of Numerical Integration,
   Second Edition,
   Dover, 2007,
   ISBN: 0486453391,
   LC: QA299.3.D28.
4. Sylvan Elhay, Jaroslav Kautsky,
   Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of Interpolatory Quadrature,
   ACM Transactions on Mathematical Software,
   Volume 13, Number 4, December 1987, pages 399-415.
5. Alan Genz, Bradley Keister,
   Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight,
   Journal of Computational and Applied Mathematics,
   Volume 71, 1996, pages 299-309.
6. John Hart, Ward Cheney, Charles Lawson, Hans Maehly, Charles Mesztenyi, John Rice, Henry Thatcher, Christoph Witzgall,
   Computer Approximations,
   Wiley, 1968,
   LC: QA297.C64.
7. Knut Petras,
   Smolyak Cubature of Given Polynomial Degree with Few Nodes for Increasing Dimension,
   Numerische Mathematik,
   Volume 93, Number 4, February 2003, pages 729-753.
8. Arthur Stroud, Don Secrest,
   Gaussian Quadrature Formulas,
   Prentice Hall, 1966,
   LC: QA299.4G3S7.
9. Shanjie Zhang, Jianming Jin,
   Computation of Special Functions,
   Wiley, 1996,
   ISBN: 0-471-11963-6,
   LC: QA351.C45

Source Code:

• sandia_rules.f90, the source code.

Examples and Tests:

• sandia_rules_test.f90, a sample calling program.
• sandia_rules_test.txt, the output file.

List of Routines:

• BINARY_VECTOR_NEXT generates the next binary vector.
• CCN_COMPUTE computes a nested Clenshaw Curtis quadrature rule.
• CCN_COMPUTE_NP computes a nested Clenshaw Curtis rule.
• CCN_COMPUTE_POINTS: compute nested Clenshaw Curtis points.
• CCN_COMPUTE_POINTS_NP: abscissas of a nested Clenshaw Curtis rule.
• CCN_COMPUTE_WEIGHTS: weights for nested Clenshaw Curtis rule.
• CCN_COMPUTE_WEIGHTS_NP computes nested Clenshaw Curtis weights.
• CHEBYSHEV1_COMPUTE computes a Chebyshev type 1 quadrature rule.
• CHEBYSHEV1_INTEGRAL evaluates a monomial Chebyshev type 1 integral.
• CHEBYSHEV2_COMPUTE computes a Chebyshev type 2 quadrature rule.
• CHEBYSHEV2_INTEGRAL evaluates a monomial Chebyshev type 2 integral.
• CLENSHAW_CURTIS_COMPUTE computes a Clenshaw Curtis quadrature rule.
• CLENSHAW_CURTIS_COMPUTE_POINTS: abscissas of a Clenshaw Curtis rule.
• CLENSHAW_CURTIS_COMPUTE_POINTS_NP: abscissas of a Clenshaw Curtis rule.
• CLENSHAW_CURTIS_COMPUTE_WEIGHTS computes Clenshaw Curtis weights.
• CLENSHAW_CURTIS_COMPUTE_WEIGHTS_NP computes Clenshaw Curtis weights.
• COMP_NEXT computes the compositions of the integer N into K parts.
• DIF_TO_R8POLY converts a divided difference table to a standard polynomial.
• DIF_VALS evaluates a divided difference polynomial at a set of points.
• FEJER2_COMPUTE computes a Fejer Type 2 quadrature rule.
• FEJER2_COMPUTE_POINTS returns the abscissas of a Fejer type 2 rule.
• FEJER2_COMPUTE_POINTS_NP returns the abscissas of a Fejer type 2 rule.
• FEJER2_COMPUTE_WEIGHTS computes weights for a Fejer type 2 quadrature rule.
• FEJER2_COMPUTE_WEIGHTS_NP: weights for a Fejer type 2 quadrature rule.
• GEGENBAUER_COMPUTE computes a Gegenbauer quadrature rule.
• GEGENBAUER_INTEGRAL integrates a monomial with Gegenbauer weight.
• GEGENBAUER_RECUR finds the value and derivative of a Gegenbauer polynomial.
• GEGENBAUER_ROOT improves an approximate root of a Gegenbauer polynomial.
• GEN_HERMITE_COMPUTE computes a generalized Gauss-Hermite quadrature rule.
• GEN_HERMITE_COMPUTE_POINTS: abscissas of a Generalized Hermite rule.
• GEN_HERMITE_COMPUTE_POINTS_NP: abscissas of a Generalized Hermite rule.
• GEN_HERMITE_COMPUTE_WEIGHTS: weights of a Generalized Hermite rule.
• GEN_HERMITE_COMPUTE_WEIGHTS_NP: weights of a Generalized Hermite rule.
• GEN_HERMITE_DR_COMPUTE computes a Generalized Hermite rule.
• GEN_HERMITE_INTEGRAL evaluates a monomial Generalized Hermite integral.
• GEN_LAGUERRE_COMPUTE: generalized Gauss-Laguerre quadrature rule.
• GEN_LAGUERRE_COMPUTE_POINTS: points of a Generalized Laguerre rule.
• GEN_LAGUERRE_COMPUTE_POINTS_NP: points of a Generalized Laguerre rule.
• GEN_LAGUERRE_COMPUTE_WEIGHTS: weights of a Generalized Laguerre rule.
• GEN_LAGUERRE_COMPUTE_WEIGHTS_NP: weights of a Generalized Laguerre rule.
• GEN_LAGUERRE_INTEGRAL evaluates a monomial genearlized Laguerre integral.
• GEN_LAGUERRE_SS_COMPUTE computes a Generalized Laguerre quadrature rule.
• GEN_LAGUERRE_SS_RECUR evaluates a Generalized Laguerre polynomial.
• GEN_LAGUERRE_SS_ROOT seeks roots of a Generalized Laguerre polynomial.
• GET_UNIT returns a free FORTRAN unit number.
• HC_COMPUTE_WEIGHTS_FROM_POINTS: Hermite-Cubic weights, user-supplied points.
• HCC_COMPUTE computes a Hermite-Cubic-Chebyshev-Spacing quadrature rule.
• HCC_COMPUTE_NP computes a Hermite-Cubic-Chebyshev-Spacing quadrature rule.
• HCC_COMPUTE_POINTS: abscissas of a Hermite-Cubic-Chebyshev-Spacing rule.
• HCC_COMPUTE_POINTS_NP: abscissas of a Hermite-Cubic-Chebyshev-Spacing rule.
• HCC_COMPUTE_WEIGHTS computes Hermite-Cubic-Chebyshev-Spacing weights.
• HCC_COMPUTE_WEIGHTS_NP computes Hermite-Cubic-Chebyshev-Spacing weights.
• HCE_COMPUTE computes a Hermite-Cubic-Equal-Spacing quadrature rule.
• HCE_COMPUTE_NP computes a Hermite-Cubic-Equal-Spacing quadrature rule.
• HCE_COMPUTE_POINTS: abscissas of a Hermite-Cubic-Equal-Spacing rule.
• HCE_COMPUTE_POINTS_NP: abscissas of a Hermite-Cubic-Equal-Spacing rule.
• HCE_COMPUTE_WEIGHTS computes Hermite-Cubic-Equal-Spacing weights.
• HCE_COMPUTE_WEIGHTS_NP computes Hermite-Cubic-Equal-Spacing weights.
• HDATA_TO_DIF sets up a divided difference table from Hermite data.
• HERMITE_COMPUTE computes a Gauss-Hermite quadrature rule.
• HERMITE_COMPUTE_POINTS computes points of a Hermite quadrature rule.
• HERMITE_COMPUTE_POINTS_NP computes points of a Hermite quadrature rule.
• HERMITE_COMPUTE_WEIGHTS computes weights of a Hermite quadrature rule.
• HERMITE_COMPUTE_WEIGHTS_NP computes weights of a Hermite quadrature rule.
• HERMITE_GENZ_KEISTER_LOOKUP returns a Genz-Keister rule for Hermite problems.
• HERMITE_GENZ_KEISTER_LOOKUP_POINTS: abscissas of a Genz-Keister Hermite rule.
• HERMITE_GENZ_KEISTER_LOOKUP_POINTS_NP: Genz-Keister Hermite abscissas.
• HERMITE_GENZ_KEISTER_LOOKUP_WEIGHTS: weights for Genz-Keister Hermite rule.
• HERMITE_GENZ_KEISTER_LOOKUP_WEIGHTS_NP sets weights for a Patterson rule.
• HERMITE_GK18_LOOKUP_POINTS: abscissas of a Hermite Genz-Keister 18 rule.
• HERMITE_GK22_LOOKUP_POINTS: abscissas of a Genz-Keister 22 Hermite rule.
• HERMITE_GK24_LOOKUP_POINTS: abscissas of a Genz-Keister 24 Hermite rule.
• HERMITE_INTEGRAL evaluates a monomial Hermite integral.
• HERMITE_INTERPOLANT_RULE: quadrature rule for a Hermite interpolant.
• HERMITE_LOOKUP_POINTS returns the abscissas of a Hermite rule.
• HERMITE_LOOKUP_WEIGHTS returns weights for Hermite quadrature rules.
• HERMITE_SS_COMPUTE computes a Hermite quadrature rule.
• HERMITE_SS_RECUR finds the value and derivative of a Hermite polynomial.
• HERMITE_SS_ROOT improves an approximate root of a Hermite polynomial.
• I4_CHOOSE computes the binomial coefficient C(N,K) as an I4.
• I4_LOG_2 returns the integer part of the logarithm base 2 of an I4.
• I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed.
• I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4MAT.
• I4MAT_WRITE writes an I4MAT file.
• I4VEC_MIN_MV determines U(1:N) /\ V for vectors U and a single vector V.
• I4VEC_PRINT prints an I4VEC.
• IMTQLX diagonalizes a symmetric tridiagonal matrix.
• JACOBI_COMPUTE: Elhay-Kautsky method for Gauss-Jacobi quadrature rule.
• JACOBI_COMPUTE_POINTS returns the points of a Jacobi rule.
• JACOBI_COMPUTE_POINTS_NP returns the points of a Jacobi rule.
• JACOBI_COMPUTE_WEIGHTS returns the weights of a Jacobi rule.
• JACOBI_COMPUTE_WEIGHTS_NP returns the weights of a Jacobi rule.
• JACOBI_INTEGRAL evaluates the integral of a monomial with Jacobi weight.
• JACOBI_SS_COMPUTE computes a Jacobi quadrature rule.
• JACOBI_SS_RECUR finds the value and derivative of a Jacobi polynomial.
• JACOBI_SS_ROOT improves an approximate root of a Jacobi polynomial.
• LAGUERRE_COMPUTE: Laguerre quadrature rule by the Elhay-Kautsky method.
• LAGUERRE_COMPUTE_POINTS computes points of a Laguerre quadrature rule.
• LAGUERRE_COMPUTE_POINTS_NP computes points of a Laguerre quadrature rule.
• LAGUERRE_COMPUTE_WEIGHTS computes weights of a Laguerre quadrature rule.
• LAGUERRE_COMPUTE_WEIGHTS_NP computes weights of a Laguerre quadrature rule.
• LAGUERRE_INTEGRAL evaluates a monomial Laguerre integral.
• LAGUERRE_LOOKUP_POINTS returns the abscissas of a Laguerre rule.
• LAGUERRE_LOOKUP_WEIGHTS returns weights for Laguerre quadrature rules.
• LAGUERRE_SS_COMPUTE computes a Laguerre quadrature rule.
• LAGUERRE_SS_RECUR finds the value and derivative of a Laguerre polynomial.
• LAGUERRE_SS_ROOT improves an approximate root of a Laguerre polynomial.
• LEGENDRE_COMPUTE: Legendre quadrature rule by the Elhay-Kautsky method.
• LEGENDRE_COMPUTE_POINTS computes abscissas of a Legendre quadrature rule.
• LEGENDRE_COMPUTE_POINTS_NP computes abscissas of a Legendre quadrature rule.
• LEGENDRE_COMPUTE_WEIGHTS computes weights of a Legendre quadrature rule.
• LEGENDRE_COMPUTE_WEIGHTS_NP computes weights of a Legendre quadrature rule.
• LEGENDRE_DR_COMPUTE computes a Legendre quadrature rule.
• LEGENDRE_INTEGRAL evaluates a monomial Legendre integral.
• LEGENDRE_ZEROS computes the zeros of the Legendre polynomial of degree N.
• LEVEL_GROWTH_TO_ORDER: convert Level and Growth to Order.
• LEVEL_TO_ORDER_DEFAULT: default growth.
• LEVEL_TO_ORDER_EXPONENTIAL: exponential growth.
• LEVEL_TO_ORDER_EXPONENTIAL_SLOW: slow exponential growth.
• LEVEL_TO_ORDER_LINEAR: linear growth.
• NC_COMPUTE computes a Newton-Cotes quadrature rule.
• NCC_COMPUTE_POINTS: Newton-Cotes Closed points
• NCC_COMPUTE_WEIGHTS: Newton-Cotes Closed weights.
• NCO_COMPUTE_POINTS: points for a Newton-Cotes Open quadrature rule.
• NCO_COMPUTE_WEIGHTS: weights for a Newton-Cotes Open quadrature rule.
• NCOH_COMPUTE_POINTS: points for a Newton-Cotes Open Half quadrature rule.
• NCOH_COMPUTE_WEIGHTS: weights for a Newton-Cotes Open Half quadrature rule.
• PATTERSON_LOOKUP returns the abscissas and weights of a Patterson rule.
• PATTERSON_LOOKUP_POINTS returns the abscissas of a Patterson rule.
• PATTERSON_LOOKUP_POINTS_NP returns the abscissas of a Patterson rule.
• PATTERSON_LOOKUP_WEIGHTS sets weights for a Patterson rule.
• PATTERSON_LOOKUP_WEIGHTS_NP sets weights for a Patterson rule.
• POINT_RADIAL_TOL_UNIQUE_COUNT counts the tolerably unique points.
• POINT_RADIAL_TOL_UNIQUE_COUNT counts the tolerably unique points.
• POINT_RADIAL_TOL_UNIQUE_COUNT_INC1 counts the tolerably unique points.
• POINT_RADIAL_TOL_UNIQUE_COUNT_INC2 counts the tolerably unique points.
• POINT_RADIAL_TOL_UNIQUE_INDEX indexes the tolerably unique points.
• POINT_RADIAL_TOL_UNIQUE_INDEX_OLD indexes the tolerably unique points.
• POINT_RADIAL_TOL_UNIQUE_INDEX_INC1 indexes the tolerably unique points.
• POINT_RADIAL_TOL_UNIQUE_INDEX_INC2 indexes unique temporary points.
• POINT_RADIAL_TOL_UNIQUE_INDEX_INC3 merges index data.
• POINT_UNIQUE_INDEX indexes unique points.
• PRODUCT_MIXED_WEIGHT computes the weights of a mixed product rule.
• R8_CEILING rounds an R8 "up" (towards +oo) to the next integer.
• R8_CHOOSE computes the binomial coefficient C(N,K) as an R8.
• R8_EPSILON returns the R8 roundoff unit.
• R8_FACTORIAL computes the factorial function.
• R8_FACTORIAL2 computes the double factorial function.
• R8_FLOOR rounds an R8 "down" (towards -infinity) to the next integer.
• R8_GAMMA evaluates Gamma(X) for a real argument.
• R8_HUGE returns a very large R8.
• R8_HYPER_2F1 evaluates the hypergeometric function F(A,B,C,X).
• R8_MOP returns the I-th power of -1 as an R8 value.
• R8_PSI evaluates the function Psi(X).
• R8COL_COMPARE compares columns in an R8COL.
• R8COL_SORT_HEAP_A ascending heapsorts an R8COL.
• R8COL_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R8COL.
• R8COL_SORTED_TOL_UNDEX indexes tolerably unique entries in a sorted R8COL.
• R8COL_SORTED_TOL_UNIQUE_COUNT: tolerably unique elements in a sorted R8COL.
• R8COL_SORTED_UNIQUE_COUNT counts unique elements in a sorted R8COL.
• R8COL_SWAP swaps columns I and J of an R8COL.
• R8COL_TOL_UNDEX indexes tolerably unique entries of an R8COL.
• R8COL_TOL_UNIQUE_COUNT counts tolerably unique entries in an R8COL.
• R8COL_UNDEX returns unique sorted indexes for an R8COL.
• R8COL_UNIQUE_INDEX indexes the unique occurrence of values in an R8COL.
• R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed.
• R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
• R8MAT_WRITE writes an R8MAT file.
• R8POLY_ANT_VAL evaluates the antiderivative of a polynomial in standard form.
• R8VEC_CHEBYSHEV creates a vector of Chebyshev spaced values.
• R8VEC_COMPARE compares two R8VEC's.
• R8VEC_DIRECT_PRODUCT2 creates a direct product of R8VEC's.
• R8VEC_INDEX_SORTED_RANGE: search index sorted vector for elements in a range.
• R8VEC_LEGENDRE creates a vector of Legendre-spaced values.
• R8VEC_LINSPACE creates a vector of linearly spaced values.
• R8VEC_MIN_POS returns the minimum positive value of an R8VEC.
• R8VEC_PRINT prints an R8VEC.
• R8VEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R8VEC.
• R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
• SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order.
• TIMESTAMP prints the current YMDHMS date as a time stamp.
• VEC_COLEX_NEXT3 generates vectors in colex order.

You can go up one level to the FORTRAN90 source codes.