A numerical library which uses a paralleled algorithm to do one dimensional integration?

Question

Is there a numerical library which can use a paralleled algorithm to do one dimensional integration (global adaptive method)? The infrastructure of my code decides that I cannot do multiple numerical integrations in parallel, but I have to use a paralleled algorithm to speed up.

Thanks!

Answer 1

Nag C numerical library does have a parallel version of adaptive quadrature (link here). Their trick is to request the user the following function

void (*f)(const double x[], Integer nx, double fv[], Integer *iflag, Nag_Comm *comm)

Here the function "f" evaluates the integrand at nx abscise points given by the vector x[]. This is where parallelization comes along, because you can use parallel_for (implemented in openmp for example) to evaluate f at those points concurrently. The integrator itself is single threaded.

Nag is a very expensive library, but if you code the integrator yourself using, for example, numerical recipes, it is not difficult to modify serial implementations to create parallel adaptive integrators using NAG idea.

I can't reproduce numerical recipes book to show where modifications are necessary due to license restriction. So let's take the simplest example of trapezoidal rule, where the implementation is quite simple and well known. The simplest way to create adaptive method using trapezoidal rule is to calculate the integral at a grid of points, then double the number of abscise points and compare the results. If the result changes by less than the requested accuracy, then there is convergence.

At each step, the trapezoidal rule can be computed using the following generic implementation

double trapezoidal( double (*f)(double x), double a, double b, int n)
{
  double h = (b - a)/n;
  double s = 0.5 * h * (f(a) + f(b));
  for( int i = 1; i < n; ++i ) s += h * f(a + i*h);
  return s;
}

Now you can make the following changes to implement NAG idea

double trapezoidal( void (*f)( double x[], int nx, double fv[] ), double a, double b, int n)
{
  double h = (b - a)/n;
  double x[n+1];
  double fv[n+1];
  for( int i = 0; i < n; ++i ) x[i+1] = (a + i * h);
  x[n] = b;

  f(x, n, fv); // inside f, use parallel_for to evaluate the integrand at x[i], i=0..n

  double s = 0.5 * h * ( fv[0] + fv[n] );
  for( int i = 1; i < n; ++i ) s += h * fv[i];
  return s;
}

This procedure, however, will only speed-up your code if the integrand is very expensive to compute. Otherwise, you should parallelize your code at higher loops and not inside the integrator.