使用 CUDA 的 Romberg 方法进行数值积分答案

【问题标题】：Numerical integration by the Romberg method using CUDA使用 CUDA 的 Romberg 方法进行数值积分
【发布时间】：2013-05-08 13:41:35
【问题描述】：

我尝试在 Google 上搜索库以在 CUDA 上进行数值集成，但找不到任何库。

1) 我想问一下，是否有任何库可用于在 CUDA 上执行（函数的）集成？

2) 如果我在 CUDA 上编写自己的代码，例如实施 Romberg 集成，我该如何进行？假设我有函数，比如f(x)；我需要为不同的时间间隔计算这个函数的积分吗？ 0.0 - 0.1, ..., 0.2 - 0.3, ..., 1.3 - 2.3?如何并行计算所有这些？

在我看来，策略是如果我必须执行，例如，1000 集成，我生成1000 线程，每个线程计算梯形以及错误估计。但是，如果我想计算其中一个积分区间的梯形以及其他积分，我不知道如何以编程方式处理这个问题。

【问题讨论】：

嗯，你知道 CUDA（和一般的 GPGPU 编程）是如何工作的吗（不是在 “并行运行一堆线程”的级别上，而是在“线程块”、“块内同步”、“共享内存与全局内存”、...）的级别？如果是，那么您的任务听起来相当容易并行化，就像您说的那样，启动一堆线程，每个线程集成一个间隔（如果您希望间隔在之后求和，请执行加法减少）。如果没有，请从适当的 CUDA 学习资源开始。
是的，我了解 CUDA 的工作原理。我知道我可以启动一个线程块，其中每个线程都可以在函数 f(x) 的间隔上执行积分。但是这种方式我使积分平行化。但是每个积分本身都是按顺序执行的。就像 CPU 上的任何串行程序都可以工作一样。是否有任何方法可以使每个积分也可以并行化。例如我可以为 X 次重复的给定间隔并行计算梯形。同时我可以并行执行其他积分。
是否有可用于执行集成的库？
嗯，如果我没记错的话，Romberg 是分层的，所以我认为各个术语不能很好地并行化。但是您可以只使用足够简单的东西，例如分段 Trapzoid 甚至 Rectangle 并使用数千个小间隔，然后进行缩减。你真正想要并行化的是什么，你想并行计算单个函数的积分，单个函数的许多不相关区间的积分，还是不同函数的多个积分？当然，单个线程是按顺序工作的，这就是为什么你有成千上万个线程。
积分只不过是对值求和（可能带有权重，尽管对于等距网格，该选项通常被高估了）。所以看看减少代码。对于 Romberg 集成，您可以在多个细化级别并行运行求和，在 CPU 上计算最终结果和误差估计。

标签： cuda numerical-integration

【解决方案1】：

正如上面 Tera 在他的评论中指出的那样，从并行编程的角度来看，集成基本上是一种简化，因此在 CUDA 中实现集成的一种非常简单的方法是利用 Thrust 库的原语（另请参阅我的回复Simpson's method to integrate real valued functions with CUDA)。

下面是一个通过 Thrust 原语实现 Romberg 集成方法的简单示例。它是site 提供的相应 Matlab 代码的“直接”翻译，因此该示例还展示了 Thurst 如何“简单”地将一些 Matlab 代码移植到 CUDA。

#include <thrust/sequence.h>

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>

#define pi_f  3.14159265358979f                 // Greek pi in single precision

struct sin_functor
{
    __host__ __device__
    float operator()(float x) const
    {
        return sin(2.f*pi_f*x);
    }
};

int main(void)
{
    int M = 5;                          // --- Maximum number of Romberg iterations

    float a     = 0.f;                  // --- Lower integration limit
    float b     = .5f;                  // --- Upper integration limit

    float hmin  = (b-a)/pow(2.f,M-1);   // --- Minimum integration step size 

    // --- Define the matrix for Romberg approximations and initialize to 1.f 
    thrust::host_vector<float> R(M*M,1.f);

    for (int k=0; k<M; k++) {

        float h = pow(2.f,k-1)*hmin;    // --- Step size for the k-th row of the Romberg matrix

        // --- Define integration nodes
        int N = (int)((b - a)/h) + 1;
        thrust::device_vector<float> d_x(N);
        thrust::sequence(d_x.begin(), d_x.end(), a, h);

        // --- Calculate function values
        thrust::device_vector<float> d_y(N);
        thrust::transform(d_x.begin(), d_x.end(), d_y.begin(), sin_functor());

        // --- Calculate integral
        R[k*M] = (.5f*h) * (d_y[0] + 2.f*thrust::reduce(d_y.begin() + 1, d_y.begin() + N - 1, 0.0f) + d_y[N-1]);

    }

    // --- Compute the k-th column of the Romberg matrix
    for (int k=1; k<M; k++) { 

        // --- The matrix of Romberg approximations is triangular!
        for (int kk=0; kk<(M-k+1); kk++) { 

            // --- See the Romberg integration algorithm
            R[kk*M+k] = R[kk*M+k-1] + (R[kk*M+k-1] - R[(kk+1)*M+k-1])/(pow(4.f,k)-1.f); 

        } 

    }

    // --- Define the vector Rnum for numerical approximations
    thrust::host_vector<float> Rnum(M); 
    thrust::copy(R.begin(), R.begin() + M, Rnum.begin());

    for (int i=0; i<M; i++) printf("%i %f\n",i,Rnum[i]);

    getchar();

    return 0;
}

【讨论】：