【发布时间】：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 时，我看到推力允许主机减少。