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

Question

我正在尝试使用 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

我已经查找了 stackoverflow.com 并发现了一些类似的问题： Similar questions 1

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

那么，对于问题 2： Similar questions 2

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

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

Answer 1

花哨的迭代器是这种操作的关键，但在推力上并不是那么直观。您可以使用 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;
}