Multidimensional integration ("cubature") code by Steven G. Johnson
<stevenj@alum.mit.edu>, based in part on code from the GNU Scientific
Library (GSL) by Brian Gough and others and from the HIntLib
numeric-integration library by Rudolf Schuerer.  Adaptive integration
of either one integrand or a vector of integrands is supported.

The latest version can be downloaded from:

    http://ab-initio.mit.edu/cubature/

 * Copyright (c) 2005-2009 Steven G. Johnson
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * Portions (see comments) based on HIntLib (also distributed under
 * the GNU GPL, v2 or later), copyright (c) 2002-2005 Rudolf Schuerer.
 *      (http://www.cosy.sbg.ac.at/~rschuer/hintlib/)
 * 
 * Portions (see comments) based on GNU GSL (also distributed under
 * the GNU GPL, v2 or later), copyright (c) 1996-2000 Brian Gough.
 *      (http://www.gnu.org/software/gsl/)

The basic algorithm is based on the adaptive cubature described in
 
     A. C. Genz and A. A. Malik, "An adaptive algorithm for numeric
     integration over an N-dimensional rectangular region,"
     J. Comput. Appl. Math. 6 (4), 295-302 (1980).

and subsequently extended to integrating a vector of integrands in

     J. Berntsen, T. O. Espelid, and A. Genz, "An adaptive algorithm
     for the approximate calculation of multiple integrals,"
     ACM Trans. Math. Soft. 17 (4), 437-451 (1991).

Note, however, that we do not use any of code from the above authors
(in part because their code is Fortran 77, but mostly because it is
under the restrictive ACM copyright license).  I did make use of some
GPL code from Rudolf Schuerer's HIntLib and from the GNU Scientific
Library as noted above, however.

I am also grateful to Dmitry Turbiner <dturbiner@alum.mit.edu>, who
implemented an initial prototype of the "vectorized" functionality for
evaluating multiple points in a single call (as opposed to multiple
functions in a single call).  (Although Dmitry implemented a working
version, I ended up re-implementing this feature from scratch as part
of a larger code-cleanup, and in order to have a single code path for
the vectorized and non-vectorized APIs.)

I subsequently extended the "vectorized" interface to use an algorithm
by Gladwell to increase the potential parallelism:

     I. Gladwell, "Vectorization of one dimensional quadrature codes,"
     pp. 230--238 in _Numerical Integration. Recent Developments,
     Software and Applications_, G. Fairweather and P. M. Keast, eds.,
     NATO ASI Series C203, Dordrecht (1987).

as described in:

     J. M. Bull and T. L. Freeman, "Parallel globally adaptive
     algorithms for multi-dimensional integration,"
     http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.6638
     (1994).

-----------------------------------------------------------------------
Usage:

You should compile cubature.c and link it with your program, and
#include the header file cubature.h.

The central subroutine you will be calling is:

int adapt_integrate(unsigned fdim, integrand f, void *fdata,
                    unsigned dim, const double *xmin, const double *xmax,
                    unsigned maxEval, double reqAbsError, double reqRelError,
                    double *val, double *err);

This integrates a function F(x), returning a vector of FDIM
integrands, where x is a DIM-dimensional vector ranging from XMIN to
XMAX (i.e. in a hypercube XMIN[i] <= x[i] <= XMAX[i]).

MAXEVAL specifies a maximum number of function evaluations (0 for no
limit).  Otherwise, the integration stops when the estimated |error|
is less than REQABSERROR (the absolute error requested), or when the
estimated |error| is less than REQRELERROR * |integral value| (the
relative error requested).

VAL and ERR are arrays of length FDIM, which upon return are the
computed integral values and estimated errors, respectively.  (The
estimated errors are based on an embedded cubature rule of lower
order; for smooth functions, this estimate is usually conservative.)

The return value of adapt_integrate is 0 on success and nonzero if
there was an error (currently, only out-of-memory situations).

The integrand function F should be a function of the form:

void F(unsigned ndim, const double *x, void *fdata,
       unsigned fdim, double *fval);

Here, the input is an array X of length NDIM (the point to be
evaluated), the output is an array FVAL of length FDIM (the vector of
function values at the point X).

The FDATA argument of F is equal to the FDATA argument passed to
adapt_integrate -- this can be used by the caller to pass any
additional information through to F as needed (rather than using
global variables, which are not re-entrant).  If F does not need any
additional data, you can just pass FDATA = NULL and ignore the FDATA
argument to F.

-----------------------------------------------------------------------
"Vectorized" interface:

For parallelization and other purposes, it is useful to call a single
integrand function with an array of points to evaluate, rather than a
single point at a time.  (e.g. you may wish to evaluate the different
points in parallel.)  This is accomplished by calling:

int adapt_integrate_v(unsigned fdim, integrand_v f, void *fdata,
                      unsigned dim, const double *xmin, const double *xmax,
                      unsigned maxEval, double reqAbsError, double reqRelError,
                      double *val, double *err);

where all the arguments are the same as above except that the
integrand function F should now be a function of the form:

void F(unsigned ndim, unsigned npts, const double *x, void *fdata,
       unsigned fdim, double *fval);

Now, x is not a single point, but an array of npts points, and upon
return the values of all fdim integrands at all npts points should be
stored in fval.  In particular, x[i*ndim + j] is the j-th coordinate
of the i-th point (i < npts and j < ndim), and the k-th function
evaluation (k < fdim) for the i-th point is returned in fval[k*npt + i].

The size of npts will vary with the dimensionality of the problem;
higher-dimensional problems will have (exponentially) larger npts,
allowing for the possibility of more parallelism.  Also npts will vary
between calls to your integrand---subsequent calls (if the first batch
of points did not achieve the desired accuracy) will have larger and
larger values of npts.

-----------------------------------------------------------------------
Test cases:

To compile a test case, just compile cubature.c while #defining
TEST_INTEGRATOR, e.g. (on Unix or GNU/Linux) via:

 cc -DTEST_INTEGRATOR -o cubature_test cubature.c -lm

The usage is then:

    ./cubature_test <dim> <tol> <integrand> <maxeval>

where <dim> = # dimensions, <tol> = relative tolerance, <integrand> is
0-7 for one of eight possible test integrands (see below) and
<maxeval> is the maximum # function evaluations (0 for none, the default).

The different test integrands are:

0: a product of cosine functions
1: a Gaussian integral of exp(-x^2), remapped to [0,infinity) limits
2: volume of a hypersphere (integrating a discontinuous function!)
3: a simple polynomial (product of coordinates)
4: a Gaussian centered in the middle of the integration volume
5: a sum of two Gaussians
6: an example function by Tsuda, a product of terms with near poles
7: a test integrand by Morokoff and Caflisch, a simple product of
   dim-th roots of the coordinates (weakly singular at the boundary)

For example:

    ./cubature_test 3 1e-5 4

integrates the Gaussian function (4) to a desired relative error
tolerance of 1e-5 in 3 dimensions.  The output is:

3-dim integral, tolerance = 1e-05
integrand 4: integral = 1, est err = 9.99952e-06, true err = 2.54397e-08
#evals = 82203

Notice that it finds the integral after 82203 function evaluations
with an estimated error of about 1e-5, but the true error (compared to
the exact result) is much smaller (2.5e-8): the error estimation is
typically conservative when applied to smooth functions like this.
