c++ CUDA/推力：如何对交错数组的列求和？

使用Thrust可以直接对一个扁平的（即由向量支持）矩阵，如示例here所示。
我想做的是对数组的列求和。
我尝试使用类似的结构，即：

// convert a linear index to a column index
template <typename T>
struct linear_index_to_col_index : public thrust::unary_function<T,T>
{
  T C; // number of columns

  __host__ __device__
  linear_index_to_col_index(T C) : C(C) {}

  __host__ __device__
  T operator()(T i)
  {
    return i % C;
  }
};

// allocate storage for column sums and indices
thrust::device_vector<int> col_sums(C);
thrust::device_vector<int> col_indices(C);

// compute row sums by summing values with equal row indices
thrust::reduce_by_key
  (thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_col_index<int>(C)),
   thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_col_index<int>(C)) + (R*C),
   array.begin(),
   col_indices.begin(),
   col_sums.begin(),
   thrust::equal_to<int>(),
   thrust::plus<int>());

然而，这导致仅对第一列求和，其余列被忽略。我猜测为什么会发生这种情况，正如reduce_by_key文档中所指出的：
对于范围[ keys_first，keys_last）that are equal, reduce_by_key copies the first element of the group to the keys_output中的每组**连续**键。[* 强调我的 *] 如果我的理解是正确的，因为行迭代器中的键是连续的（即索引[0 - (C-1)]将给予0，然后[C - (2C-1)]将给出1，依此类推），它们最终被求和在一起。但是列迭代器将把索引[0 - (C-1)]Map到[0 - (C-1)]，然后再次开始，索引[C - (2C-1)]将Map到[0 - (C-1)]`，依此类推。使得所产生的值不连续。
这种行为对我来说是不直观的，我希望分配给同一个键的所有数据点都被分组在一起，但这是另一个讨论。
无论如何，我的问题是：如何使用Thrust对交错数组的列求和？

这些操作（对行求和、对列求和等）通常在GPU上受存储器带宽限制。因此，我们可能要考虑如何构建一个算法，使GPU内存带宽的最佳利用。特别是，如果可能的话，我们希望从thrust代码生成的底层内存访问被合并。简而言之，这意味着相邻的GPU线程将从内存中的相邻位置读取。
原始row-summing example显示此属性：由thrust产生的相邻线程将读取存储器中的相邻元素。例如，如果我们有R行，那么我们可以看到，在reduce_by_key操作期间，由thrust创建的第一个R线程将全部阅读矩阵的第一“行”。由于与第一行相关联的内存位置都分组在一起，因此我们得到合并访问。
解决这个问题（如何对列求和）的一种方法是使用与行求和示例类似的策略，但使用permutation_iterator使同一个键序列的所有线程读取 * 列 * 数据，而不是 * 行 * 数据。这个置换迭代器将接受底层数组和一个Map序列。这个Map序列是由transform_iterator使用一个特殊的函子创建的，该函子应用于counting_iterator，将线性（行为主）索引转换为列为主索引，因此第一个C线程将读取矩阵的第一列 * 的元素，而不是第一行。由于第一个C线程将属于相同的键序列，因此它们将在reduce_by_key操作中相加。这就是我在下面的代码中所称的方法1。
然而，这种方法的缺点是相邻线程不再阅读内存中相邻的值-我们破坏了合并，正如我们将看到的，性能影响是明显的。
对于以行优先顺序存储在内存中的大型矩阵（我们在这个问题中已经讨论过的顺序），对 * 列 * 求和的一个相当理想的方法是让每个线程使用for循环对单个列求和。这在CUDA C中实现起来相当简单，我们可以在Thrust中使用适当定义的函子来执行此操作。
我在下面的代码中将其称为方法2。此方法只会启动与矩阵中的列一样多的线程。对于具有足够大数量的列（例如10，000或更多）的矩阵，该方法将使GPU饱和并有效地使用可用的存储器带宽。如果你检查函子，你会发现这是一个有点“不寻常”的推力适应，但完全法律的。
下面是比较这两种方法的代码：

$ cat t994.cu
#include <thrust/device_vector.h>
#include <thrust/reduce.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/functional.h>
#include <thrust/sequence.h>
#include <thrust/transform.h>

#include <iostream>

#define NUMR 1000
#define NUMC 20000
#define TEST_VAL 1

#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}

typedef int mytype;

// from a linear (row-major) index, return column-major index
struct rm2cm_idx_functor : public thrust::unary_function<int, int>
{
  int r;
  int c;

  rm2cm_idx_functor(int _r, int _c) : r(_r), c(_c) {};

  __host__ __device__
  int operator() (int idx)  {
    unsigned my_r = idx/c;
    unsigned my_c = idx%c;
    return (my_c * r) + my_r;
  }
};

// convert a linear index to a column index
template <typename T>
struct linear_index_to_col_index : public thrust::unary_function<T,T>
{
  T R; // number of rows

  __host__ __device__
  linear_index_to_col_index(T R) : R(R) {}

  __host__ __device__
  T operator()(T i)
  {
    return i / R;
  }
};

struct sum_functor
{
  int R;
  int C;
  mytype *arr;

  sum_functor(int _R, int _C, mytype *_arr) : R(_R), C(_C), arr(_arr) {};

  __host__ __device__
  mytype operator()(int myC){
    mytype sum = 0;
      for (int i = 0; i < R; i++) sum += arr[i*C+myC];
    return sum;
    }
};


int main(){
  int C = NUMC;
  int R = NUMR;
  thrust::device_vector<mytype> array(R*C, TEST_VAL);

// method 1: permutation iterator

// allocate storage for column sums and indices
  thrust::device_vector<mytype> col_sums(C);
  thrust::device_vector<int> col_indices(C);

// compute column sums by summing values with equal column indices
  unsigned long long m1t = dtime_usec(0);
  thrust::reduce_by_key(thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_col_index<int>(R)),
   thrust::make_transform_iterator(thrust::counting_iterator<int>(R*C), linear_index_to_col_index<int>(R)),
   thrust::make_permutation_iterator(array.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator<int>(0), rm2cm_idx_functor(R, C))),
   col_indices.begin(),
   col_sums.begin(),
   thrust::equal_to<int>(),
   thrust::plus<int>());
  cudaDeviceSynchronize();
  m1t = dtime_usec(m1t);
  for (int i = 0; i < C; i++)
    if (col_sums[i] != R*TEST_VAL) {std::cout << "method 1 mismatch at: " << i << " was: " << col_sums[i] << " should be: " << R*TEST_VAL << std::endl; return 1;}
  std::cout << "Method1 time: " << m1t/(float)USECPSEC << "s" << std::endl;

// method 2: column-summing functor

  thrust::device_vector<mytype> fcol_sums(C);
  thrust::sequence(fcol_sums.begin(), fcol_sums.end());  // start with column index
  unsigned long long m2t = dtime_usec(0);
  thrust::transform(fcol_sums.begin(), fcol_sums.end(), fcol_sums.begin(), sum_functor(R, C, thrust::raw_pointer_cast(array.data())));
  cudaDeviceSynchronize();
  m2t = dtime_usec(m2t);
  for (int i = 0; i < C; i++)
    if (fcol_sums[i] != R*TEST_VAL) {std::cout << "method 2 mismatch at: " << i << " was: " << fcol_sums[i] << " should be: " << R*TEST_VAL << std::endl; return 1;}
  std::cout << "Method2 time: " << m2t/(float)USECPSEC << "s" << std::endl;
  return 0;
}
$ nvcc -O3 -o t994 t994.cu
$ ./t994
Method1 time: 0.034817s
Method2 time: 0.00082s
$

很明显，对于足够大的矩阵，方法2比方法1快得多。
如果你不熟悉置换迭代器，可以看看thrust quick start guide。

c++ CUDA/推力：如何对交错数组的列求和？

1条答案

相关问题

热门标签

最新问答