c++ - CUDA:减少算法

Question

我是 C++/CUDA 的新手。我尝试通过递归内核输出(在内核包装器).

例如Implementing Max Reduce in Cuda这个问题的最佳答案，当线程大小足够小时，他/她的实现基本上是顺序的。

但是，当我编译和运行它时，我总是收到“Segmentation fault”......？

>> nvcc -o mycode mycode.cu
>> ./mycode
Segmentail fault.

使用 cuda 6.5 在 K40 上编译

这是内核，与我链接“越界”检查器的 SO 帖子基本相同，不同之处在于:

#include <stdio.h>

/* -------- KERNEL -------- */
__global__ void reduce_kernel(float * d_out, float * d_in, const int size)
{
  // position and threadId
  int pos = blockIdx.x * blockDim.x + threadIdx.x;
  int tid = threadIdx.x;

  // do reduction in global memory
  for (unsigned int s = blockDim.x / 2; s>0; s>>=1)
  {
    if (tid < s)
    {
      if (pos+s < size) // Handling out of bounds
      {
        d_in[pos] = d_in[pos] + d_in[pos+s];
      }
    }
  }

  // only thread 0 writes result, as thread
  if (tid==0)
  {
    d_out[blockIdx.x] = d_in[pos];
  }
}

当 1 个 block 不包含所有数据时，我提到要处理的内核包装器。

/* -------- KERNEL WRAPPER -------- */
void reduce(float * d_out, float * d_in, const int size, int num_threads)
{
  // setting up blocks and intermediate result holder
  int num_blocks = ((size) / num_threads) + 1;
  float * d_intermediate;
  cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);

  // recursively solving, will run approximately log base num_threads times.
  do
  {
    reduce_kernel<<<num_blocks, num_threads>>>(d_intermediate, d_in, size);

    // updating input to intermediate
    cudaMemcpy(d_in, d_intermediate, sizeof(float)*num_blocks, cudaMemcpyDeviceToDevice);

    // Updating num_blocks to reflect how many blocks we now want to compute on
      num_blocks = num_blocks / num_threads + 1;

    // updating intermediate
    cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
  }
  while(num_blocks > num_threads); // if it is too small, compute rest.

  // computing rest
  reduce_kernel<<<1, num_blocks>>>(d_out, d_in, size);

}

初始化输入/输出并创建用于测试的虚假数据的主程序。

/* -------- MAIN -------- */
int main(int argc, char **argv)
{
  // Setting num_threads
  int num_threads = 512;
  // Making bogus data and setting it on the GPU
  const int size = 1024;
  const int size_out = 1;
  float * d_in;
  float * d_out;
  cudaMalloc(&d_in, sizeof(float)*size);
  cudaMalloc((void**)&d_out, sizeof(float)*size_out);
  const int value = 5;
  cudaMemset(d_in, value, sizeof(float)*size);

  // Running kernel wrapper
  reduce(d_out, d_in, size, num_threads);

  printf("sum is element is: %.f", d_out[0]);
}

Answer 1

我会在您的代码中指出一些事情。

1. 作为一般规则/样板，我总是建议使用 proper cuda error checking并使用 cuda-memcheck 运行您的代码，任何时候您在使用 cuda 代码时遇到问题。然而，这些方法对段错误没有多大帮助，尽管它们可能在以后有所帮助(见下文)。

2. 实际段错误发生在这条线上:

   printf("sum is element is: %.f", d_out[0]);
   

   您违反了一条基本的 CUDA 编程规则:主机指针不得在设备代码中取消引用，并且设备指针不得在主机代码中取消引用。后一种情况适用于此。 d_out 是一个设备指针(通过 cudaMalloc 分配)。如果您尝试在主机代码中取消引用这些指针，那么这些指针没有意义，并且这样做会导致段错误。

   解决方案是在打印出来之前将数据复制回主机:

   float result;
   cudaMemcpy(&result, d_out, sizeof(float), cudaMemcpyDeviceToHost);
   printf("sum is element is: %.f", result);
   

3. 在循环中对同一个变量使用 cudaMalloc，而不执行任何 cudaFree 操作，这不是好的做法，并且可能会导致错误长时间运行的循环中的内存错误，如果在较大的程序中使用这样的构造，也可能导致程序内存泄漏:

   do
   {
     ...
   
     cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
   }
   while...
   

   在这种情况下，我认为更好的方法和简单的修复方法是在重新分配之前使用 cudaFree d_intermediate:

   do
   {
     ...
     cudaFree(d_intermediate);
     cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
   }
   while...
   

4. 这可能不是您认为的那样:

   const int value = 5;
   cudaMemset(d_in, value, sizeof(float)*size);
   

   您可能已经意识到这一点，但是 cudaMemset 与 memset 一样，对字节数量进行操作。所以你用对应于 0x05050505 的值填充 d_in 数组(我不知道当解释为 float 数量)。由于您提到的是虚假值，您可能已经意识到这一点。但这是一个常见错误(例如，如果您实际上试图在每个 float 位置用值 5 初始化数组)，所以我想我会指出来。

您的代码还有其他问题(如果您进行上述修复，然后使用 cuda-memcheck 运行您的代码，您会发现这些问题)。要了解如何进行良好的并行缩减，我建议学习 CUDA 并行缩减 sample code和 presentation .出于性能原因，不建议并行减少全局内存。

为了完整起见，以下是我发现的一些其他问题:

1. 您的内核代码需要一个适当的__syncthreads() 语句来确保 block 中所有线程的工作在任何线程进入 for 循环的下一次迭代之前完成。

2. 您对内核中全局内存的最终写入也需要以读取位置在边界内为条件。否则，您始终启动额外 block 的策略将允许从此行读取越界(cuda-memcheck 将显示这一点)。

3. reduce 函数中循环中的缩减逻辑通常是困惑的，需要通过多种方式重新处理。

我并不是说这段代码没有缺陷，但它似乎适用于给定的测试用例并产生正确答案 (1024):

#include <stdio.h>

/* -------- KERNEL -------- */
__global__ void reduce_kernel(float * d_out, float * d_in, const int size)
{
  // position and threadId
  int pos = blockIdx.x * blockDim.x + threadIdx.x;
  int tid = threadIdx.x;

  // do reduction in global memory
  for (unsigned int s = blockDim.x / 2; s>0; s>>=1)
  {
    if (tid < s)
    {
      if (pos+s < size) // Handling out of bounds
      {
        d_in[pos] = d_in[pos] + d_in[pos+s];
      }
    }
    __syncthreads();
  }

  // only thread 0 writes result, as thread
  if ((tid==0) && (pos < size))
  {
    d_out[blockIdx.x] = d_in[pos];
  }
}

/* -------- KERNEL WRAPPER -------- */
void reduce(float * d_out, float * d_in, int size, int num_threads)
{
  // setting up blocks and intermediate result holder
  int num_blocks = ((size) / num_threads) + 1;
  float * d_intermediate;
  cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
  cudaMemset(d_intermediate, 0, sizeof(float)*num_blocks);
  int prev_num_blocks;
  // recursively solving, will run approximately log base num_threads times.
  do
  {
    reduce_kernel<<<num_blocks, num_threads>>>(d_intermediate, d_in, size);

    // updating input to intermediate
    cudaMemcpy(d_in, d_intermediate, sizeof(float)*num_blocks, cudaMemcpyDeviceToDevice);

    // Updating num_blocks to reflect how many blocks we now want to compute on
      prev_num_blocks = num_blocks;
      num_blocks = num_blocks / num_threads + 1;

    // updating intermediate
    cudaFree(d_intermediate);
    cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
    size = num_blocks*num_threads;
  }
  while(num_blocks > num_threads); // if it is too small, compute rest.

  // computing rest
  reduce_kernel<<<1, prev_num_blocks>>>(d_out, d_in, prev_num_blocks);

}

/* -------- MAIN -------- */
int main(int argc, char **argv)
{
  // Setting num_threads
  int num_threads = 512;
  // Making non-bogus data and setting it on the GPU
  const int size = 1024;
  const int size_out = 1;
  float * d_in;
  float * d_out;
  cudaMalloc(&d_in, sizeof(float)*size);
  cudaMalloc((void**)&d_out, sizeof(float)*size_out);
  //const int value = 5;
  //cudaMemset(d_in, value, sizeof(float)*size);
  float * h_in = (float *)malloc(size*sizeof(float));
  for (int i = 0; i <  size; i++) h_in[i] = 1.0f;
  cudaMemcpy(d_in, h_in, sizeof(float)*size, cudaMemcpyHostToDevice);

  // Running kernel wrapper
  reduce(d_out, d_in, size, num_threads);
  float result;
  cudaMemcpy(&result, d_out, sizeof(float), cudaMemcpyDeviceToHost);
  printf("sum is element is: %.f\n", result);
}