CUDA 阵列缩减优化

Posted 2023-03-23

【中文标题】CUDA 阵列缩减优化

【英文标题】：CUDA array reduction optimisation

【发布时间】：2021-11-17 15:07:52

【问题描述】：

我有两个数组 x（大小为 N ~1-1 亿）和 a（非常小 Na ~1000-10000），我想使用 x 将 a 定义为

for(int j = 0; j < N; j++) 
  float  i = floor( x[j] / da); // in principle i < size(a)

  a[(int)i] += 0.5;
  a[(int)i+1] += 0.5; // I simplify the problem

对于上下文，x 是粒子位置，a 是每个单元格的粒子数。

我想在 CUDA 中执行这个功能。主要问题是我可以同时对同一内存进行多次修改，因为 x 未排序。

我找到了以下解决方案，但我发现它很慢。 我定义了一个临时数组d_temp_a，大小为 Na * 使用的线程数。然后，我将其缩减为我的完整数组。

这里是代码（使用nvcc -std=c++11 example_reduce.cu -o example_reduce.out）

#include "stdio.h"
#include <cuda.h>
#include <random>
using namespace std;


__global__ void getA(float *d_x, float *d_a, float *d_temp_a, int N, int Na, float da)

// Get our global thread ID
  int index = blockIdx.x * blockDim.x + threadIdx.x;
  int stride = blockDim.x * gridDim.x;

  float ix ;

  // Compute a
  for(int x = index; x < N; x += stride) 
      ix = floor( d_x[x] / da );

      d_temp_a[((int)ix) + Na * index] += 0.5;
      d_temp_a[((int)ix + 1) + Na * index] += 0.5;
  
  __syncthreads();


  // Reduce
  for(int l = index; l < Na; l += stride) 
      for(int m = 0; m < stride; m += 1) 
          d_a[l] += d_temp_a[l + Na * m];
      
  
  __syncthreads();



int main(int argc, char **argv)


  int N = 1000000;   
  int Na = 4096;   

  float L = 50; // box size
  float dxMesh = L / Na; // cell size

  float *h_x, *h_a;  // host data

  h_x = (float *)malloc(N * sizeof(float));
  h_a = (float *)malloc(Na * sizeof(float));

  /* Initialize random seed: */
  std::default_random_engine generator;
  std::uniform_real_distribution<float> generate_unif_dist(0.0,1.0);

  // h_x random initialisation
  for(int x = 0; x < N; x++) 
        float random = generate_unif_dist(generator);
        h_x[x] = random * L;
    

 
  int blockSize = 512; // Number of threads in each thread block
  int gridSize = (int)ceil((float) N /blockSize); // Number of thread blocks in grid

  float *d_x, *d_a;  // device data

  cudaMalloc((void **) &d_x, N * sizeof(float));
  cudaMalloc((void **) &d_a, Na * sizeof(float));

  cudaMemcpy(d_x, h_x, N * sizeof(float), cudaMemcpyHostToDevice);

  // Create temp d_a array
  float *d_temp_a;
  cudaMalloc((void **) &d_temp_a, Na * blockSize * gridSize * sizeof(float));

  getA<<<gridSize,blockSize>>>(d_x, d_a, d_temp_a, N, Na, da);

  cudaMemcpy(h_a, d_a, Na * sizeof(float), cudaMemcpyDeviceToHost);

  free(h_x);
  free(h_a);

  cudaFree(d_x);
  cudaFree(d_a);
  cudaFree(d_temp_a);

  return 0;

这很慢，因为我只为数组的每个元素使用 1 个线程。 我的问题：有没有办法优化这种减少？我还发现拥有这个非常大的 Na * 线程数数组效率低下。有没有办法避免使用它？

请注意，我打算稍后编写一个 2D 版本，其中 x 和 y 定义 a[i][j]。

【问题讨论】：

建议您按照通常的方法进行共享内存扫描式并行缩减。学习教程here。不建议仅使用全局内存进行缩减。在cuda标签讨论共享内存并行减少中已经有很多问题，并且有一个CUDA示例代码与之前链接的教程材料一起使用。

请注意，CUB 可能会帮助您做到这一点。还要注意除法很昂贵，即使在 GPU 上也是如此，我认为您可以安全地将其替换为 da 乘以预先计算的 1 / da。如果浮点值始终为正，floor 也可以优化。最后，最后一个__syncthreads 没用。

感谢 Robert 和 Jérôme 的回答。我正在使用全局内存，因为我需要在 d_a 上执行 FFT（而 cuFFT 是一个主机 API）。在检查 CUB 时，我看到推力允许主机减少。

【参考方案1】：

我认为你的方法对于这个问题可能有点矫枉过正。

和 cmets 中的其他人一样，我也认为您可以实施您的推力还原想法。但是，我的方法涉及计算每个 idx 的出现次数，然后插入这些计数（参见 Counting occurrences of numbers in a CUDA array）

推力减小方法

所以我几乎完全使用推力方法（填充、变换、排序、reduce_by_key）实现了这一点，并在最终结果上运行了一个最终内核，以在两个相邻单元格之间拆分值。这很有效，并且比您的 CUDA 方法快得多，但它仍然比简单的 CPU 实现慢得多。最大的问题是 N 个值的排序和 reduce_by_key。

struct custom_functor
    float factor;
    custom_functor(float _factor)
      factor = _factor;
    
    __host__ __device__ int operator()(float &x) const 
        return (int) floor(x / factor);
    
;

__global__ void thrust_reduce_kernel(float *d_a, int* d_a_idxs, int* d_a_cnts, int N, int Na, int n_entries)

  int index = blockIdx.x * blockDim.x + threadIdx.x;

  if (index >= n_entries)
    return;

  int a_idx = d_a_idxs[index];
  int a_cnt = d_a_cnts[index];

  if ((a_idx + 1) >= Na || a_idx < 0 || a_idx >= Na || (a_idx + 1) < 0)
  
    printf("Should not happen according to you!\n");
    return;
  

  atomicAdd(&d_a[a_idx], a_cnt * 0.5f);
  atomicAdd(&d_a[a_idx+1], a_cnt * 0.5f);


void test_thrust_reduce(float *d_x, float *d_a, float *h_a, int N, int Na, float da)

  int *d_xi, *d_ones;
  int *d_a_cnt_keys, *d_a_cnt_vals;

  cudaMalloc((void**) &d_xi, N * sizeof(int));
  cudaMalloc((void**) &d_ones, N * sizeof(float));

  cudaMalloc((void**) &d_a_cnt_keys, Na * sizeof(int));
  cudaMalloc((void**) &d_a_cnt_vals, Na * sizeof(int));
  CUDA_CHECK;

  thrust::device_ptr<float> dt_x(d_x);
  thrust::device_ptr<float> dt_a(d_a);
  thrust::device_ptr<int> dt_xi(d_xi);
  thrust::device_ptr<int> dt_ones(d_ones);
  thrust::device_ptr<int> dt_a_cnt_keys(d_a_cnt_keys);
  thrust::device_ptr<int> dt_a_cnt_vals(d_a_cnt_vals);

  custom_functor f(da);
  thrust::fill(thrust::device, dt_a, dt_a + Na, 0.0f);
  thrust::fill(thrust::device, dt_ones, dt_ones + N, 1);
  thrust::fill(thrust::device, dt_a_cnt_keys, dt_a_cnt_keys + Na, -1);
  thrust::fill(thrust::device, dt_a_cnt_vals, dt_a_cnt_vals + Na, 0);

  thrust::transform(thrust::device, dt_x, dt_x + N, dt_xi, f);
  thrust::sort(thrust::device, dt_xi, dt_xi + N);

  thrust::pair<thrust::device_ptr<int>,thrust::device_ptr<int>> new_end;
  new_end = thrust::reduce_by_key(thrust::device, dt_xi, dt_xi + N, dt_ones, 
                                  dt_a_cnt_keys, dt_a_cnt_vals);

  int n_entries = new_end.first - dt_a_cnt_keys;
  int n_entries_2 = new_end.first - dt_a_cnt_keys;

  dim3 dimBlock(256);
  dim3 dimGrid((n_entries + dimBlock.x - 1) / dimBlock.x);
  thrust_reduce_kernel<<<dimGrid, dimBlock>>>(d_a, d_a_cnt_keys, d_a_cnt_vals, N, Na, n_entries);
  cudaMemcpy(h_a, d_a, Na * sizeof(float), cudaMemcpyDeviceToHost);

  cudaFree(d_xi);
  cudaFree(d_ones);
  cudaFree(d_a_cnt_keys);
  cudaFree(d_a_cnt_vals);

简单的 atomicAdd 方法

所以我很好奇你是否可以在 d_x 中为每个条目使用一个简单的 atomicAdd，这被证明是所有解决方案中最快的。

__global__ void simple_atomicAdd_kernel(const float *d_x, float *d_a, float da, int N, int Na)

  int index = blockIdx.x * blockDim.x + threadIdx.x;

  if (index >= N)
    return;

  int a_idx = floor(d_x[index] / da); // in principle i < size(a)

  atomicAdd(&d_a[a_idx], 0.5f);
  atomicAdd(&d_a[a_idx+1], 0.5f);
 

void test_simple_atomicAdd(float *d_x, float *d_a, float *h_a, int N, int Na, float da)

  cudaMemset(d_a, 0, Na * sizeof(float));

  dim3 dimBlock(256);
  dim3 dimGrid((N + dimBlock.x - 1) / dimBlock.x);
  simple_atomicAdd_kernel<<<dimGrid, dimBlock>>>(d_x, d_a, da, N, Na);
  cudaMemcpy(h_a, d_a, Na * sizeof(float), cudaMemcpyDeviceToHost);

结果

您可以在下面看到我的 N=100,000 和 da=0.1 的时间。您的初始值 N = 1,000,000 导致我出现 out_of_memory 异常。全部

Times: 
- CPU Reference:         912 us
- CUDA Custom reduce:    34275 us
- CUDA Thrust reduce:    2144 us
- CUDA Simple atomicAdd: 59 us

查看更高的 N 值，Thrust reduce 方法开始变得更好，因为我们在 atomicAdd 方法中有更多的冲突。这在很大程度上取决于您的 x 值和 da 的值：

Times (N=1,000,000, da=0.1): 
- CPU Reference:         9398 us
- CUDA Thrust reduce:    1287 us
- CUDA Simple atomicAdd: 409 us

Times (N=10,000,000, da=0.1): 
- CPU Reference:         92068 us
- CUDA Thrust reduce:    3879 us
- CUDA Simple atomicAdd: 3851 us

Times (N=100,000,000, da=0.1): 
- CPU Reference:         918950 us
- CUDA Thrust reduce:    21051 us
- CUDA Simple atomicAdd: 38583 us

免责声明：我远不是 CUDA 编程方面的专家，我可能遗漏了一些重要的东西。这些只是我的发现，我确信在你的情况下存在更好的方法。但是，简单的 atomicAdd 方法可能是解决您的问题的快速简便的方法。

您可以在此处查看完整代码：https://github.com/steimich96/cuda_reduction_experiments

我希望这会有所帮助。 干杯，迈克尔

【讨论】：

谢谢迈克尔。这是一个超级答案！

确实，在我的第一种方法中，我需要一个巨大的临时矩阵，因此不适用于大量粒子。对我来说，atomicAdd 显示了更好的结果Times(N=100,000,000, da=0.1): - CPU Reference: 3.38543e+06 ms - CUDA Thrust reduce: 44554 ms - CUDA Simple atomicAdd: 8887 ms。因此，我将使用 atomicAdd。我认为推力方法可以改进，特别是因为它也依赖于原子操作。