c++ - 使用Cuda进行并行尺寸缩减(3D到2D求和)

Question

在CUDA应用程序中，我有一个N x N x D矩阵，我想通过在整个第一(或第二)轴上求和来简化为N x D。如何最有效地做到这一点？

通常，N大于10000，D为2或3。

使用atomicAdd的快速而简单的解决方案如下:

namespace kernel {
    __global__ void sumNND(float* devPtrIn, float* devPtrOut, const int N, const int D) {
        int index = blockIdx.x * blockDim.x + threadIdx.x;
        int stride = blockDim.x * gridDim.x;

        for (int id = index; id < N * N * D; id += stride) {
            const unsigned int d = id % D;
            const unsigned int i = (id - d) / D;
            const unsigned int n = i / N;
            const unsigned int m = i % N;

            atomicAdd(&devPtrOut[d + D * n], devPtrIn[d + D * n + N * m]);
        }
    }
}

void sumNND(const int numBlocks, const int blockSize, float* devPtrIn, float* devPtrOut, const int N, const int D) {
    HANDLE_ERROR(cudaMemset(devPtrOut, 0, N * D * sizeof(float)));
    kernel::sumNND<<<numBlocks, blockSize>>>(devPtrIn, devPtrOut, N, D);
    HANDLE_ERROR(cudaDeviceSynchronize());
}

在其中调用 sumNND的地方
loopSize = N * N * D， blockSize = 768和 numBlocks = (loopSize + blockSize - 1) / blockSize。

这是我的时间轴上的瓶颈(不足为奇)，但是我不知道如何有效地并行化降维。有指针吗？

Answer 1

任何CUDA程序员的前两个优化优先级是:

使用很多线程

有效使用内存

对于您的问题，您不会遇到第一个问题-它很容易分解为一系列独立的问题，可以分配给许多并行线程。然后，第二优先级是您要关注的地方。关于全局内存，这意味着我们应该尽可能地争取合并访问。我们应该特别注意阅读内容。

我需要做一些假设。我假设您的维度组织为ROW，COLUMN，DEPTH，并且您的数据存储在通常的C样式(即行主存储)中。然后，利用这些假设，请求(在整个第一(或第二)轴上求和)实际上是在整个行上求和或在整个列上求和。如果您在此处在 cuda标记上进行了一些搜索，则会找到这两个示例的有效示例( here是这样的示例之一)。尽管它们不一定全部涵盖3D情况，但它们应该提供一个很好的路线图。您会发现这两种情况应该以不同的方式处理，着眼于合并的全局内存访问，即已经提到的优化优先级。行方向也是合并方向，因此，如果需要对行求和，则需要使用经典的并行约简技术，以便可以读取行并将元素求和在一起。如果我们需要对列求和，那么高效的内核更容易编写；每个线程可以负责一列，并且可以只将一个运行中的总和保持在for循环中。

就您而言，您似乎正在对列求和(但请参见下面的注释)。以下是一个有效的示例，将您的方法与运行速度更快的column-sum方法进行了比较，并结合了访问方式(相邻线程读取内存中的相邻元素):

$ cat t1263.cu
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

const int my_N = 10000;
const int my_D = 3;
const int my_blockSize = 768;
const int my_loopSize = my_N*my_N*my_D;
const int my_numBlocks = (my_loopSize + my_blockSize -1)/my_blockSize;
const int bsize = 512;
const float TOL = 0.1f;

#define HANDLE_ERROR(x) x

#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)

#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;
}

namespace kernel {
    __global__ void sumNND(float* devPtrIn, float* devPtrOut, const int N, const int D) {
        int index = blockIdx.x * blockDim.x + threadIdx.x;
        int stride = blockDim.x * gridDim.x;

        for (int id = index; id < N * N * D; id += stride) {
            const unsigned int d = id % D;
            const unsigned int i = (id - d) / D;
            const unsigned int n = i / N;
            const unsigned int m = i % N;

            atomicAdd(&devPtrOut[d + D * n], devPtrIn[d + D * n + N * m]);
        }
    }
}

void sumNND(const int numBlocks, const int blockSize, float* devPtrIn, float* devPtrOut, const int N, const int D) {
    HANDLE_ERROR(cudaMemset(devPtrOut, 0, N * D * sizeof(float)));
    kernel::sumNND<<<numBlocks, blockSize>>>(devPtrIn, devPtrOut, N, D);
    HANDLE_ERROR(cudaDeviceSynchronize());
}

// kernel assumes 1 block assigned per row, use block-striding methodology
// assumes block size is a power of 2
__global__ void sum_rows_NND(const float * __restrict__  devPtrIn, float * __restrict__  devPtrOut, const int N, const int D) {
  __shared__ float sdata[bsize];
  sdata[threadIdx.x] = 0;
  for (int i = threadIdx.x; i < N; i += blockDim.x) // block-stride
    sdata[threadIdx.x] += devPtrIn[(blockIdx.x * N) + i];
  __syncthreads();
  for (int i = blockDim.x>>1; i > 0; i>>=1){
    if (threadIdx.x < i) sdata[threadIdx.x] += sdata[threadIdx.x+i];
    __syncthreads();}
  if (!threadIdx.x) devPtrOut[blockIdx.x] = sdata[0];
}



// kernel assumes one thread assigned per column sum
// launch N threads
 __global__ void sum_cols_NND(const float * __restrict__  devPtrIn, float * __restrict__  devPtrOut, const int N, const int D) {
  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  int ido = idx;
  if (idx < N){
    for (int j = 0; j < D; j++){
      float temp = 0;
      for (int i = 0; i < N; i++) temp += devPtrIn[idx + (i*N)];
      devPtrOut[ido] = temp;
      ido += N;
      idx += N*N;}}
}

int main(){

  float *h_data, *d_data, *h_res1, *h_res2, *d_res;

  h_data = new float[my_loopSize];
  cudaMalloc(&d_data, my_loopSize*sizeof(d_data[0]));
  h_res1 = new float[my_N*my_D];
  h_res2 = new float[my_N*my_D];
  cudaMalloc(&d_res, my_N*my_D*sizeof(d_res[0]));
  for (int i = 0; i < my_loopSize; i++) h_data[i] = rand()/(float)RAND_MAX;
  cudaCheckErrors("CUDA failure");
  cudaMemcpy(d_data, h_data, my_loopSize*sizeof(d_data[0]), cudaMemcpyHostToDevice);
  // test original approach
  cudaMemset(d_res, 0, my_N*my_D*sizeof(d_res[0]));
  unsigned long long dt1 = dtime_usec(0);
  kernel::sumNND<<<my_numBlocks, my_blockSize>>>(d_data, d_res, my_N, my_D);
  cudaDeviceSynchronize();
  dt1 = dtime_usec(dt1);
  cudaMemcpy(h_res1, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);

  //test columnwise reduction
  unsigned long long dt2 = dtime_usec(0);
  //sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
  sum_cols_NND<<<(my_N + bsize -1)/bsize, bsize>>>(d_data, d_res, my_N, my_D);
  cudaDeviceSynchronize();
  dt2 = dtime_usec(dt2);
  cudaMemcpy(h_res2, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);

  // validate results
  for (int i = 0; i < my_N; i++)
    if (fabsf(h_res1[i] - h_res2[i]) > TOL) {printf("mismatch at %d, was %f, should be %f\n", i, h_res2[i], h_res1[i]); return -1;}
  cudaCheckErrors("program error");

  printf("results match,  kernel 1 time: %fs, kernel 2 time: %fs\n", dt1/(float)USECPSEC, dt2/(float)USECPSEC);
  // time row reduction kernel
  unsigned long long dt3 = dtime_usec(0);
  sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
  cudaDeviceSynchronize();
  dt3 = dtime_usec(dt3);
  printf("row reduction kernel time: %fs\n", dt3/(float)USECPSEC);
  cudaCheckErrors("program error");
}
$ nvcc -arch=sm_52 -o t1263 t1263.cu
$ ./t1263
results match,  kernel 1 time: 0.459971s, kernel 2 time: 0.013678s
row reduction kernel time: 0.013724s
$

笔记:

经过优化的内核比您的朴素原子内核快30倍左右。我怀疑其中很大一部分实际上不是原子的使用，而是未分批访问。新型GPU上的全局原子可能很快。

我的内核和您的内核之间元素列总和的第一个“页面”(NxN)匹配(即，前N个结果匹配)。第一页之后(前N个结果)，我们的结果有所不同。我很确定我的索引编制是正确的，但是花了一段时间试图弄清楚您的索引编制之后，我放弃了。如果您尝试对列求和，我怀疑您的内核索引中有一个错误，并且所有上述假设都是正确的。

我还包括了行求和内核的时序测量，它看起来有很大不同，但是产生的时序几乎相同。这是可以预料的，因为针对这些类型问题的最佳内核将受到内存带宽的限制，这在两种情况下都是相同的。最佳内核将以合并的方式一次加载所有数据。之后，行和与列和机制对内核时间的影响相对较小。

对数据的初始化进行少量修改，我认为很容易证明您的内核没有创建正确的索引，因此没有在第一个“页面”之后(即在第一个N结果之后)产生正确的行总和。 )。在对索引进行了更多研究之后，我对出了什么问题有了一些想法。一个示例问题是，对于无法被N整除的D，您的内核d变量在第一个“页面”之后不会重置为零，但这不是唯一的问题。

根据项目4，这是修改了数据初始化的代码版本，并对所有 N * D结果进行了全面测试。数据初始化为，第一页的第一列将全部为零，下一列的全部为1，下一列的全部为2，依此类推。在第二页上，我们将所有内容加1，因此第一列将全部为1，第二列全为2，依此类推。因此，应该很容易就列的总和达成一致。对于第一页，列的总和应为0、10000、20000等。对于第二页，它们的应为10000、20000、30000等。在第二页的第一列上，我的代码生成10000，您的代码生成1.在注释中更改索引后，第一页的第一列将产生0，而您的代码将产生9999。根据我描述的数据初始化，1和9999可能不是有效的列总和:

$ cat t1263.cu
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

const int my_N = 10000;
const int my_D = 3;
const int my_blockSize = 768;
const int my_loopSize = my_N*my_N*my_D;
const int my_numBlocks = (my_loopSize + my_blockSize -1)/my_blockSize;
const int bsize = 512;
const float TOL = 0.1f;

#define HANDLE_ERROR(x) x

#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)

#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;
}

namespace kernel {
    __global__ void sumNND(float* devPtrIn, float* devPtrOut, const int N, const int D) {
        int index = blockIdx.x * blockDim.x + threadIdx.x;
        int stride = blockDim.x * gridDim.x;

        for (int id = index; id < N * N * D; id += stride) {
            const unsigned int d = id % D;       // 0 1 2 0 1 2 0 1 2
            const unsigned int i = (id - d) / D; // 0 0 0 1 1 1 2 2 2
            const unsigned int n = i / N;        // 0 0 0 0 0 0 0 0 0
            const unsigned int m = i % N;        // 0 0 0 1 1 1 2 2 2

            atomicAdd(&devPtrOut[d + D * n],    //  0 1 2 0 1 2 0 1 2
              devPtrIn[d + D * n + N * m]);     //  0 1 2 0+N 1+N 2+N 0+2N 1+2N 2+2N
        }
    }
}

void sumNND(const int numBlocks, const int blockSize, float* devPtrIn, float* devPtrOut, const int N, const int D) {
    HANDLE_ERROR(cudaMemset(devPtrOut, 0, N * D * sizeof(float)));
    kernel::sumNND<<<numBlocks, blockSize>>>(devPtrIn, devPtrOut, N, D);
    HANDLE_ERROR(cudaDeviceSynchronize());
}

// kernel assumes 1 block assigned per row, use block-striding methodology
// assumes block size is a power of 2
__global__ void sum_rows_NND(const float * __restrict__  devPtrIn, float * __restrict__  devPtrOut, const int N, const int D) {
  __shared__ float sdata[bsize];
  sdata[threadIdx.x] = 0;
  for (int i = threadIdx.x; i < N; i += blockDim.x) // block-stride
    sdata[threadIdx.x] += devPtrIn[(blockIdx.x * N) + i];
  __syncthreads();
  for (int i = blockDim.x>>1; i > 0; i>>=1){
    if (threadIdx.x < i) sdata[threadIdx.x] += sdata[threadIdx.x+i];
    __syncthreads();}
  if (!threadIdx.x) devPtrOut[blockIdx.x] = sdata[0];
}



// kernel assumes one thread assigned per column sum
// launch N threads
 __global__ void sum_cols_NND(const float * __restrict__  devPtrIn, float * __restrict__  devPtrOut, const int N, const int D) {
  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  int ido = idx;
  if (idx < N){
    for (int j = 0; j < D; j++){
      float temp = 0;
      for (int i = 0; i < N; i++) temp += devPtrIn[idx + (i*N)];
      devPtrOut[ido] = temp;
      ido += N;
      idx += N*N;}}
}

int main(){

  float *h_data, *d_data, *h_res1, *h_res2, *d_res;

  h_data = new float[my_loopSize];
  cudaMalloc(&d_data, my_loopSize*sizeof(d_data[0]));
  h_res1 = new float[my_N*my_D];
  h_res2 = new float[my_N*my_D];
  cudaMalloc(&d_res, my_N*my_D*sizeof(d_res[0]));
  for (int i = 0; i < my_loopSize; i++) h_data[i] = i%my_N + i/(my_N*my_N); //rand()/(float)RAND_MAX;
  cudaCheckErrors("CUDA failure");
  cudaMemcpy(d_data, h_data, my_loopSize*sizeof(d_data[0]), cudaMemcpyHostToDevice);
  // test original approach
  cudaMemset(d_res, 0, my_N*my_D*sizeof(d_res[0]));
  unsigned long long dt1 = dtime_usec(0);
  kernel::sumNND<<<my_numBlocks, my_blockSize>>>(d_data, d_res, my_N, my_D);
  cudaDeviceSynchronize();
  dt1 = dtime_usec(dt1);
  cudaMemcpy(h_res1, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);

  //test columnwise reduction
  unsigned long long dt2 = dtime_usec(0);
  //sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
  sum_cols_NND<<<(my_N + bsize -1)/bsize, bsize>>>(d_data, d_res, my_N, my_D);
  cudaDeviceSynchronize();
  dt2 = dtime_usec(dt2);
  cudaMemcpy(h_res2, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);

  // validate results
  for (int i = 0; i < my_N*my_D; i++)
    if (fabsf(h_res1[i] - h_res2[i]) > TOL) {printf("mismatch at %d, was %f, should be %f\n", i, h_res2[i], h_res1[i]); return -1;}
  cudaCheckErrors("program error");

  printf("results match,  kernel 1 time: %fs, kernel 2 time: %fs\n", dt1/(float)USECPSEC, dt2/(float)USECPSEC);
  // time row reduction kernel
  unsigned long long dt3 = dtime_usec(0);
  sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
  cudaDeviceSynchronize();
  dt3 = dtime_usec(dt3);
  printf("row reduction kernel time: %fs\n", dt3/(float)USECPSEC);
  cudaCheckErrors("program error");
}
$ nvcc -arch=sm_52 -o t1263 t1263.cu
$ ./t1263
mismatch at 10000, was 10000.000000, should be 1.000000
$