使用CUDA的Romberg方法的数值积分

Question

我尝试在谷歌上搜索CUDA上的数值积分库，但没有找到任何库。

1)我想问一下，有没有库可以用来在CUDA上执行集成(函数)？

2)如果我在CUDA上编写自己的代码，例如实现Romberg集成，我应该如何继续？假设我有一个函数，比如说f(x)；我是否需要计算这个函数在不同区间的积分，例如0.0 - 0.1，...，0.2 - 0.3，...，1.3 - 2.3？如何并行计算所有这些参数？

在我看来，策略是如果我必须执行，例如，1000集成，我生成1000线程，每个线程计算trapzoids以及误差估计。但是，如果我想要并行计算其中一个积分区间的陷波，以及其他积分，我不知道如何编程。

Answer 1

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

下面是一个简单的例子，实现了龙贝格积分方法的推力原语。它是相应的site代码的“直接”翻译，所以这个例子还展示了如何通过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;
}