c++ 由推力变换组成的for循环的并行化

5cnsuln7  于 2023-01-06  发布在  其他
关注(0)|答案(1)|浏览(144)

我实现了一个包含多个Thrust转换的for循环,目的是为i从0到N的每个值计算r[i],简单地说,r是一个列向量,它的每个元素都可以独立计算。
因此,我正在寻找一种并行化for循环的方法,如下所示:

for(int i=0; i < N; i++) {
   thrust::device_vector<float> P(N, 0.0);  
   thrust::device_vector<int> corr_col_indices_d(col_indices.begin() + row_begin[i], col_indices.begin() + row_begin[i+1]); // indices of the columns
   thrust::device_vector<float> corr_values_d(values_d.begin() + row_begin[i], values_d.begin() + row_begin[i+1]); // values of the columns

    // P[j] = corr_values_d[k] if j is in corr_col_indices_d, else 0  (increment k if j is in corr_col_indices_d)
    thrust::scatter(corr_values_d.begin(), corr_values_d.end(), corr_col_indices_d.begin(), P.begin());
    r2[i] = thrust::inner_product(P.begin(), P.end(), r1.begin(), 0.0f);
}
    • 1)**在大量的谷歌搜索,漫游Stackoverflow和NVIDIA之后,我试图用循环变量i将所有连续的转换放入一个更大的"转换"中。
auto counting_iter = thrust::make_counting_iterator(0);
thrust::transform(counting_iter, counting_iter + N, r2.begin(), [&](int i) {
    thrust::device_vector<float> P(N, 0.0);  
    thrust::device_vector<int> corr_col_indices_d(col_indices.begin() + row_begin[i], col_indices.begin() + row_begin[i+1]); /
    thrust::device_vector<float> corr_values_d(values_d.begin() + row_begin[i], values_d.begin() + row_begin[i+1]); 
    thrust::scatter(corr_values_d.begin(), corr_values_d.end(), corr_col_indices_d.begin(), P.begin());
    thrust::transform(P.begin(), P.end(), r1.begin(), P.begin(), thrust::multiplies<float>());
    return thrust::reduce(P.begin(), P.end());
});

不幸的是,它不起作用,要么是没有这样的转换,要么是我的语法错了。

    • 2)**然后我试着创建一个函子,把所有这些device_vectors作为输入并对它们进行操作。正如here所述,不可能从外部将device_vectors传递给一个函子-因此我试着将它们作为原始指针。
struct loop {
    // constructor that takes a vector as a parameter
    __host__ __device__
    loop(int *t_row_begin, int *t_col_indices, float*t_values, float *r1): 
        t_row_begin_(t_row_begin), t_col_indices_(t_col_indices), t_values_(t_values), r1_(r1) {}

    // member variable to store the vector
    int *t_row_begin_;
    int *t_col_indices_;
    float *t_values_;
    float *r1_;

    __host__ __device__
    float operator()(int i) const {
        thrust::device_vector<float> P(N, 0.0);  
        thrust::device_vector<int> corr_col_indices_d(t_col_indices_ + t_row_begin_[i], t_col_indices_ + t_row_begin_[i + 1]); // indices of the columns
        thrust::device_vector<float> corr_values_d(t_values_ + t_row_begin_[i], t_values_ + t_row_begin_[i+1]); // values of the columns

        thrust::scatter(corr_values_d.begin(), corr_values_d.end(), corr_col_indices_d.begin(), P.begin());
        return thrust::inner_product(P.begin(), P.end(), r1.begin(), 0.0f);
    }
};

以及循环本身:

loop lp(thrust::raw_pointer_cast(row_begin_d.data()), 
            thrust::raw_pointer_cast(col_indices_d.data()), 
            thrust::raw_pointer_cast(values_d.data()), 
            thrust::raw_pointer_cast(r1.data()));
auto iter = thrust::make_counting_iterator(0);
// perform the operations for each iteration of the loop using transform
thrust::transform(iter, iter + N, r2.begin(), lp);
    • 3)**我甚至尝试过将参数传递给operator而不是函子的构造函数:
struct loop {
    __host__ __device__
    float operator()(int i, thrust::device_vector<int>& col_indices, thrust::device_vector<float>& values_d, thrust::device_vector<int>& row_begin, thrust::device_vector<float>& r1) const {
        thrust::device_vector<float> P(N, 0.0);  
        thrust::device_vector<int> corr_col_indices_d(col_indices.begin() + row_begin[i], col_indices.begin() + row_begin[i+1]); // indices of the columns
        thrust::device_vector<float> corr_values_d(values_d.begin() + row_begin[i], values_d.begin() + row_begin[i+1]); // values of the columns
        thrust::scatter(corr_values_d.begin(), corr_values_d.end(), corr_col_indices_d.begin(), P.begin());
        return thrust::inner_product(P.begin(), P.end(), r1.begin(), 0.0f);
    }
};
auto iter = thrust::make_counting_iterator(0);
thrust::transform(iter, iter + N, r2.begin(), 
   thrust::make_transform_iterator(iter, loop()), 
   thrust::make_zip_iterator(thrust::make_tuple(col_indices, values_d, row_begin, r1)));

它们都不能编译,所有那些复杂的错误信息也没有什么帮助。所以,我在这一点上寻求一些帮助。
CUDA版本:11.2
推力版本:1.10.0
编辑:如果你想知道,这些向量对应于CSR矩阵表示的分量:

vector<int> row_begin;
vector<float> values;
vector<int> col_indices;

更新

  • 根据@paleonix的建议,将transformreduce融合为inner_product
n53p2ov0

n53p2ov01#

1.去掉循环内部的分配,不需要复制行,P可以重用:

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

void foo(int N,
         thrust::host_vector<int> const &row_begin,
         thrust::device_vector<int> const &col_indices,
         thrust::device_vector<float> const &values_d,
         thrust::device_vector<float> const &r1,
         thrust::host_vector<float> &r2) {

    thrust::device_vector<float> P(N);
    for(int i = 0; i < N; ++i) {
        thrust::fill(P.begin(), P.end(), 0.0f);

        // P[j] = corr_values_d[k] if j is in corr_col_indices_d, else 0  (increment k if j is in corr_col_indices_d)
        thrust::scatter(values_d.cbegin() + row_begin[i],
                        values_d.cbegin() + row_begin[i+1],
                        col_indices.cbegin() + row_begin[i],
                        P.begin());

        r2[i] = thrust::inner_product(P.cbegin(), P.cend(),
                                      r1.cbegin(),
                                      0.0f);
    }
}

1.在r1上使用置换迭代器,而不是将值分散到P中,这样效率更高。

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/inner_product.h>
#include <thrust/iterator/permutation_iterator.h>

void foo(int N,
         thrust::host_vector<int> const &row_begin,
         thrust::device_vector<int> const &col_indices,
         thrust::device_vector<float> const &values_d,
         thrust::device_vector<float> const &r1,
         thrust::host_vector<float> &r2) {

    auto const r1_iter =
        thrust::make_permutation_iterator(
            r1.cbegin(),
            col_indices.cbegin());

    for(int i = 0; i < N; ++i) {
        r2[i] =
            thrust::inner_product(
                values_d.cbegin() + row_begin[i],
                values_d.cbegin() + row_begin[i+1],
                r1_iter + row_begin[i],
                0.0f);
    }
}
  1. inner_product中没有多少并行性了,所以要按顺序并行化外部循环:
#include <thrust/device_vector.h>
#include <thrust/inner_product.h>
#include <thrust/transform.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/counting_iterator.h>

void foo(int N,
         thrust::device_vector<int> const &row_begin,
         thrust::device_vector<int> const &col_indices,
         thrust::device_vector<float> const &values_d,
         thrust::device_vector<float> const &r1,
         thrust::device_vector<float> &r2) {

    auto const row_begin_ptr = row_begin.data();
    auto const col_indices_ptr = col_indices.data();
    auto const values_d_ptr = values_d.data();
    auto const r1_iter =
        thrust::make_permutation_iterator(
            r1.cbegin(),
            col_indices.cbegin());

    thrust::transform(
        thrust::make_counting_iterator(0),
        thrust::make_counting_iterator(0) + N,
        r2.begin(),
        [=] __host__ __device__ (int i){
            return thrust::inner_product(thrust::seq,
                                         values_d_ptr + row_begin_ptr[i],
                                         values_d_ptr + row_begin_ptr[i+1],
                                         r1_iter + row_begin_ptr[i],
                                         0.0f);
        });
}

1.虽然上述解决方案对于例如行非常小且规则的带状矩阵应该是足够的,但是像单个长行这样的不规则性将使该解决方案再次变得非常低效。替代方案是使用由thrust::reduce_by_key实现的分段/批处理约简。在此要使用reduce_by_key,需要“解压缩”CSR矩阵矩阵(将行偏移转换为关键点),即使在引擎盖下推力可能返回到引擎盖下的行偏移(在后端使用CUB)。为了避免这种低效,我直接通过cub::DeviceSegmentedReduce::Sum使用CUB。为了仍然将变换/乘法部分融合到约简中,我们可以使用变换迭代器,为了代码质量,我也抛弃了置换迭代器,直接在变换迭代器中实现了收集:

#include <cub/cub.cuh>

#include <thrust/device_vector.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>

void foo(int N,
         thrust::device_vector<int> const &row_begin,
         thrust::device_vector<int> const &col_indices,
         thrust::device_vector<float> const &values_d,
         thrust::device_vector<float> const &r1,
         thrust::device_vector<float> &r2) {

    auto const col_indices_ptr = col_indices.data();
    auto const values_d_ptr = values_d.data();
    auto const r1_ptr = r1.data();

    auto const corr_iter =
        thrust::make_transform_iterator(
            thrust::make_counting_iterator(0),
            [=] __host__ __device__ (int j){
                return values_d_ptr[j] * r1_ptr[col_indices_ptr[j]];
            });
    
    // Determine temporary storage
    size_t temp_storage_bytes = 0;
    cub::DeviceSegmentedReduce::Sum(nullptr, temp_storage_bytes,
                                    corr_iter,
                                    r2.begin(),
                                    N,
                                    row_begin.cbegin(), row_begin.cbegin() + 1);
    // Allocate temporary storage
    thrust::device_vector<char> d_temp_storage(temp_storage_bytes);
    // Run sum-reduction
    cub::DeviceSegmentedReduce::Sum(thrust::raw_pointer_cast(d_temp_storage.data()),
                                    temp_storage_bytes,
                                    corr_iter,
                                    r2.begin(),
                                    N,
                                    row_begin.cbegin(), row_begin.cbegin() + 1);
}
避免临时缓冲区初始化(和分配)

最后一个解决方案中唯一缺少的“理想”性能是临时存储被不必要地初始化。这可以通过使用一个自定义分配器来避免,如Thrust示例uninitialized_vector.cu所示。为了避免膨胀,我没有在上面的代码中包含它。
一个更好的解决方案是来自RAPIDS Memory Managerrmm::device_buffer,但CUDA工具包中不包括这个解决方案。
在未来,libcudac++有望给我们提供一个类似的C++选项,因为他们目前正在处理内存资源。
如果重复执行此操作,也可以重用临时内存。即使不直接使用CUB,也可以使用池内存资源来实现此目的。请参阅cuda/custom_temporary_allocation.cumr_basic.cu

相关问题