CUDA 阵列缩减优化

问题描述 投票：0回答：1

我有两个数组

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]
。

arrays cuda reduction
1个回答
0
投票

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

和评论中的其他人一样，我也认为你可以用推力来实现你的减少想法。但是，我的方法涉及计算每个 idx 的出现次数，然后插入这些计数（请参阅计算 CUDA 数组中数字的出现次数）

减少推力的方法

因此，我几乎完全使用推力方法（填充、变换、排序、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 的较高值，推力减少方法开始变得更好，因为原子添加方法中存在更多冲突。这很大程度上取决于您的 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://pastebin.com/BZYJ2rYe

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

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.