CUDA：减少算法

Question

我是 C++/CUDA 的新手。我尝试实现并行算法“reduce”，能够处理任何类型的输入大小和线程大小，而不会通过递归内核输出来增加渐近并行运行time（在内核包装器).

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

但是，我在编译时一直收到“分段错误”并且运行它..？

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

使用 cuda 6.5 在 K40 上编译

这是内核，与SO post基本相同post我为"out of bounds"链接的检查器不同：

#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 个块不包含所有数据的情况。

/* -------- 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);

}

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

/* -------- 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

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

作为一般 rule/boilerplate，我总是建议使用 proper cuda error checking 和运行你的代码 cuda-memcheck，任何时候你遇到问题一个cuda代码。然而，这些方法对段错误没有多大帮助，尽管它们稍后可能会有所帮助（见下文）。
实际段错误发生在这条线上：
```
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);
```
在循环中对同一个变量使用 cudaMalloc，而不进行任何 cudaFree 操作，这不是好的做法，可能会导致 out-of-memory 错误在 long-running 循环中，如果在更大的程序中使用这样的构造，也可能导致程序内存泄漏：
```
do
{
  ...

  cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
}
while...
```
在这种情况下，我认为更好的方法和简单的解决方法是 cudaFree d_intermediate 就在您 re-allocate:
之前
```
do
{
  ...
  cudaFree(d_intermediate);
  cudaMalloc(&d_intermediate, sizeof(float)*num_blocks);
}
while...
```
这可能与您认为的不同：
```
const int value = 5;
cudaMemset(d_in, value, sizeof(float)*size);
```
可能您已经意识到这一点，但是 cudaMemset 与 memset 一样，对字节数量进行操作。所以你用对应于 0x05050505 的值填充 d_in 数组（我不知道当解释为 float 数量时该位模式对应什么）。由于您提到的是虚假值，您可能已经意识到这一点。但这是一个常见的错误（例如，如果你实际上试图在每个 float 位置用值 5 初始化数组），所以我想我会指出它。

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

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

您的内核代码需要一个适当的 __syncthreads() 语句来确保块中所有线程的工作在任何线程进入 for-loop 的下一次迭代之前完成.
您对内核中全局内存的最终写入也需要以 read-location 为 in-bounds 为条件。否则，您始终启动额外块的策略将允许从该行读取 out-of-bounds （cuda-memcheck 将显示此）。
reduce 函数中循环中的归约逻辑通常是混乱的，需要以多种方式 re-worked。

我并不是说这段代码是 defect-free，但它似乎适用于给定的测试用例并产生正确答案 (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);
}

CUDA：减少算法

CUDA: Reduce algorithm

c++

algorithm

parallel-processing

reduce

cuda