使用并行算法进行一维积分的数值库？

Question

是否有可以使用并行算法进行一维积分（全局自适应方法）的数值库？我的代码的基础结构决定了我不能并行进行多个数值积分，但我必须使用并行算法来加快速度。

谢谢！

Answer 1

Nag C 数值库确实具有自适应正交的并行版本（链接here）。他们的技巧是向用户请求以下功能

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

这里的函数“f”计算向量x[]给出的nx横坐标处的被积函数。这就是并行化出现的地方，因为您可以使用parallel_for（例如在openmp 中实现）同时在这些点上评估f。集成器本身是单线程的。

Nag 是一个非常昂贵的库，但是如果您自己使用例如numerical recipes 编写积分器，则使用 NAG 的想法修改串行实现以创建并行自适应积分器并不困难。

由于许可限制，我无法复制数字食谱书来显示需要修改的地方。因此，让我们以trapezoidal rule 的最简单示例为例，该示例的实现非常简单且众所周知。使用梯形规则创建自适应方法的最简单方法是计算点网格处的积分，然后将横坐标点数加倍并比较结果。如果结果的变化小于要求的精度，则存在收敛。

在每一步，梯形规则都可以使用以下通用实现来计算

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;
}

现在您可以进行以下更改来实现 NAG 想法

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;
}

但是，如果被积函数的计算成本非常高，则此过程只会加速您的代码。否则，您应该在更高的循环中并行化您的代码，而不是在积分器内部。

Answer 2

为什么不简单地围绕一个单线程算法实现一个包装器，该算法将边界细分的积分分派到不同的线程，然后在最后将它们相加？例如

thread 0: i0 = integral(x0, (x0+x1)/2)
thread 1: i1 = integral((x0+x1)/2, x1)

i = i0 + i1