CUDA 数组缩减优化
CUDA array reduction optimisation
我有两个数组 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 未排序。
我找到了以下解决方案，但我发现它很慢。
我定义了一个大小为 Na * 使用的线程数的临时数组 d_temp_a。然后，我将它减少到我的完整数组。
这是代码（使用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]。
我认为您的方法对于这个问题可能有点矫枉过正。
和评论中的其他人一样，我也认为你可以实现你的推力减少的想法。但是，我的方法包括计算每个 idx 的 appea运行ces，然后插入这些计数（参见 Counting occurrences of numbers in a CUDA array）
推力减小方法
所以我几乎完全使用推力方法（填充、t运行sform、排序、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
希望对您有所帮助。
干杯，迈克尔
我有两个数组 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 未排序。
我找到了以下解决方案，但我发现它很慢。
我定义了一个大小为 Na * 使用的线程数的临时数组 d_temp_a。然后，我将它减少到我的完整数组。
这是代码（使用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]。
我认为您的方法对于这个问题可能有点矫枉过正。
和评论中的其他人一样，我也认为你可以实现你的推力减少的想法。但是，我的方法包括计算每个 idx 的 appea运行ces，然后插入这些计数（参见 Counting occurrences of numbers in a CUDA array）
推力减小方法
所以我几乎完全使用推力方法（填充、t运行sform、排序、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
希望对您有所帮助。 干杯，迈克尔