在 CUDA/Thrust 中,如何在 for-each 操作期间访问向量元素的邻居?

In CUDA / Thrust, how can I access a vector element's neighbor during a for-each operation?

我正在尝试使用 CUDA 中的 Thrust 库进行一些科学模拟,但我陷入了以下操作,这基本上是一个 for-each 循环:

device_vector<float> In(N);

for-each In(x) in In
      Out(x) = some_calculation(In(x-1),In(x),In(x+1));
end

我已经查过 whosebug.com 并找到了一些类似的问题: Similar questions 1

但似乎只有当some_calculation函数在2个参数之间完成时才可能使用变换迭代器,因为变换迭代器最多传递两个参数。

那么,对于问题2: Similar questions 2

讨论刚刚结束,没有结论。

我认为这是一个简单的问题,因为它是并行计算的自然要求。谁能告诉我该怎么做?

奇特的迭代器是这种操作的关键,它在推力上并不是那么直观。您可以使用 zip_iterator 创建可以迭代的值元组,因此对于典型的 f(x[i-1], x[i], x[i+1]) 类型函数,您会得到如下内容:

#include <iostream>
#include <cmath>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/tuple.h>
#include <thrust/transform.h>

struct divided_diff {
    float dx;
    divided_diff(float _dx) : dx(_dx) {};

    float operator()(const thrust::tuple<float, float, float> &in) const {
        float y0 = in.get<0>();
        float y1 = in.get<1>();
        float y2 = in.get<2>();

        return (y0 - 2.f * y1 + y2) / (dx * dx);
    }
};

int main() {
    const int N = 10;
    const float dx = 0.1f;
    float x[N], y[N], dydx[N];

    for (int i = 0; i < N; ++i) {
        x[i] = dx * float(i);
        y[i] = std::sin(x[i]);
        dydx[i] = 0.f;
    }

    auto begin = thrust::make_zip_iterator(thrust::make_tuple(&y[0], &y[1], &y[2]));
    auto end = thrust::make_zip_iterator(thrust::make_tuple(&y[N-2], &y[N-1], &y[N]));

    divided_diff f(dx);
    thrust::transform(begin, end, &dydx[1], f);

    for (int i = 0; i < N; ++i) {
        std::cout << i << " " << dydx[i] << std::endl;
    }

    return 0;
}

此处仿函数一次处理一个元组,其中该元组包含来自同一数组或迭代序列中三个不同起点的三个输入。


编辑:显然,将此代码的主机版本转换为使用设备构造对最初的发布者来说具有挑战性,因此这里是一个使用 thrust::device_vector 作为基本容器在设备上执行所有内容的版本:

#include <iostream>
#include <cmath>
#include <thrust/tuple.h>
#include <thrust/transform.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/device_vector.h>
#include <thrust/sequence.h>

struct divided_diff {
    float dx;
    divided_diff(float _dx) : dx(_dx) {};

    __device__
    float operator()(const thrust::tuple<float, float, float> &in) {
        float y0 = in.get<0>();
        float y1 = in.get<1>();
        float y2 = in.get<2>();

        return (y0 - 2.f*y1 + y2) / (dx * dx);
    }
};

struct mysinf {
    __device__
    float operator()(const float &x) { 
        return __sinf(x); 
    }
};

int main()
{

    const int N = 10;
    const float dx = 0.1f;
    thrust::device_vector<float> x(N), y(N), dydx(N-2);

    thrust::sequence(x.begin(), x.end(), 0.f, dx); 
    thrust::transform(x.begin(), x.end(), y.begin(), mysinf());

    auto start  = thrust::make_zip_iterator(thrust::make_tuple(y.begin(), y.begin()+1, y.begin()+2));
    auto finish = thrust::make_zip_iterator(thrust::make_tuple(y.end()-2, y.end()-1, y.end()));

    divided_diff f(dx);
    thrust::transform( start, finish, dydx.begin(), f);

    thrust::device_vector<float>::iterator it = dydx.begin();
    for(; it != dydx.end(); ++it) {
        float val = *it;
        std::cout << val << std::endl;
    }

    return 0;
}