在CUDA应用程序中，我有一个N x N x D矩阵，我想通过对整个第一（或第二）轴求和来减少到N x D。我如何最有效地完成这项工作？
通常，N大于10000且D为2或3。
使用atomicAdd的快速而天真的解决方案如下：
namespace kernel {
__global__ void sumNND(float* devPtrIn, float* devPtrOut, const int N, const int D) {
int index = blockIdx.x * blockDim.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;
for (int id = index; id < N * N * D; id += stride) {
const unsigned int d = id % D;
const unsigned int i = (id - d) / D;
const unsigned int n = i / N;
const unsigned int m = i % N;
atomicAdd(&devPtrOut[d + D * n], devPtrIn[d + D * n + N * m]);
}
}
}
void sumNND(const int numBlocks, const int blockSize, float* devPtrIn, float* devPtrOut, const int N, const int D) {
HANDLE_ERROR(cudaMemset(devPtrOut, 0, N * D * sizeof(float)));
kernel::sumNND<<<numBlocks, blockSize>>>(devPtrIn, devPtrOut, N, D);
HANDLE_ERROR(cudaDeviceSynchronize());
}
在哪里召唤sumNND
loopSize = N * N * D，blockSize = 768和numBlocks = (loopSize + blockSize - 1) / blockSize。
这（毫不奇怪）是我的时间线中的瓶颈，但我无法弄清楚如何有效地并行化降维。有什么指针吗？
任何CUDA程序员的前两个优化优先级是：
对于你的问题，你可以毫不费力地使用第一个问题 - 它很容易分解成一组独立的问题，可以分配给很多并行线程。第二个优先级是你想要关注的地方。关于全局内存，这意味着我们应尽可能争取合并访问。我们应该特别注意阅读。
我需要做一些假设。我假设您的维度组织是ROW，COLUMN，DEPTH，并且您的数据存储在通常的C风格，即行主要存储中。然后，根据这些假设，请求（在整个第一（或第二）轴上求和）有效地对整行进行求和或对整列进行求和。如果你在cuda标签上做一些搜索，你会发现两者的工作实例（here就是这样一个例子）。虽然它们并不一定都涵盖3D案例，但它们应该提供相当不错的路线图。您将发现的是，这两种情况应该以不同的方式处理，着眼于合并的全局内存访问，即已经提到的优化优先级。行方向也是合并方向，所以如果我们需要对行进行求和，那么我们需要使用经典的并行缩减技术，这样我们就可以读取行，并将元素加在一起。如果我们需要对列进行求和，那么高效的内核就更容易编写;每个线程都可以负责一个列，并且可以在for循环中保持运行总和。
在您的情况下，您似乎是对列进行求和（但请参阅下面的注释）。接下来是一个有效的例子，比较你的方法与更快的运行列和方法，合并访问（相邻线程读取内存中的相邻元素）：
$ cat t1263.cu
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
const int my_N = 10000;
const int my_D = 3;
const int my_blockSize = 768;
const int my_loopSize = my_N*my_N*my_D;
const int my_numBlocks = (my_loopSize + my_blockSize -1)/my_blockSize;
const int bsize = 512;
const float TOL = 0.1f;
#define HANDLE_ERROR(x) x
#define cudaCheckErrors(msg) \
do { \
cudaError_t __err = cudaGetLastError(); \
if (__err != cudaSuccess) { \
fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
msg, cudaGetErrorString(__err), \
__FILE__, __LINE__); \
fprintf(stderr, "*** FAILED - ABORTING\n"); \
exit(1); \
} \
} while (0)
#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL
long long dtime_usec(unsigned long long start){
timeval tv;
gettimeofday(&tv, 0);
return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}
namespace kernel {
__global__ void sumNND(float* devPtrIn, float* devPtrOut, const int N, const int D) {
int index = blockIdx.x * blockDim.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;
for (int id = index; id < N * N * D; id += stride) {
const unsigned int d = id % D;
const unsigned int i = (id - d) / D;
const unsigned int n = i / N;
const unsigned int m = i % N;
atomicAdd(&devPtrOut[d + D * n], devPtrIn[d + D * n + N * m]);
}
}
}
void sumNND(const int numBlocks, const int blockSize, float* devPtrIn, float* devPtrOut, const int N, const int D) {
HANDLE_ERROR(cudaMemset(devPtrOut, 0, N * D * sizeof(float)));
kernel::sumNND<<<numBlocks, blockSize>>>(devPtrIn, devPtrOut, N, D);
HANDLE_ERROR(cudaDeviceSynchronize());
}
// kernel assumes 1 block assigned per row, use block-striding methodology
// assumes block size is a power of 2
__global__ void sum_rows_NND(const float * __restrict__ devPtrIn, float * __restrict__ devPtrOut, const int N, const int D) {
__shared__ float sdata[bsize];
sdata[threadIdx.x] = 0;
for (int i = threadIdx.x; i < N; i += blockDim.x) // block-stride
sdata[threadIdx.x] += devPtrIn[(blockIdx.x * N) + i];
__syncthreads();
for (int i = blockDim.x>>1; i > 0; i>>=1){
if (threadIdx.x < i) sdata[threadIdx.x] += sdata[threadIdx.x+i];
__syncthreads();}
if (!threadIdx.x) devPtrOut[blockIdx.x] = sdata[0];
}
// kernel assumes one thread assigned per column sum
// launch N threads
__global__ void sum_cols_NND(const float * __restrict__ devPtrIn, float * __restrict__ devPtrOut, const int N, const int D) {
int idx = threadIdx.x+blockDim.x*blockIdx.x;
int ido = idx;
if (idx < N){
for (int j = 0; j < D; j++){
float temp = 0;
for (int i = 0; i < N; i++) temp += devPtrIn[idx + (i*N)];
devPtrOut[ido] = temp;
ido += N;
idx += N*N;}}
}
int main(){
float *h_data, *d_data, *h_res1, *h_res2, *d_res;
h_data = new float[my_loopSize];
cudaMalloc(&d_data, my_loopSize*sizeof(d_data[0]));
h_res1 = new float[my_N*my_D];
h_res2 = new float[my_N*my_D];
cudaMalloc(&d_res, my_N*my_D*sizeof(d_res[0]));
for (int i = 0; i < my_loopSize; i++) h_data[i] = rand()/(float)RAND_MAX;
cudaCheckErrors("CUDA failure");
cudaMemcpy(d_data, h_data, my_loopSize*sizeof(d_data[0]), cudaMemcpyHostToDevice);
// test original approach
cudaMemset(d_res, 0, my_N*my_D*sizeof(d_res[0]));
unsigned long long dt1 = dtime_usec(0);
kernel::sumNND<<<my_numBlocks, my_blockSize>>>(d_data, d_res, my_N, my_D);
cudaDeviceSynchronize();
dt1 = dtime_usec(dt1);
cudaMemcpy(h_res1, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);
//test columnwise reduction
unsigned long long dt2 = dtime_usec(0);
//sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
sum_cols_NND<<<(my_N + bsize -1)/bsize, bsize>>>(d_data, d_res, my_N, my_D);
cudaDeviceSynchronize();
dt2 = dtime_usec(dt2);
cudaMemcpy(h_res2, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);
// validate results
for (int i = 0; i < my_N; i++)
if (fabsf(h_res1[i] - h_res2[i]) > TOL) {printf("mismatch at %d, was %f, should be %f\n", i, h_res2[i], h_res1[i]); return -1;}
cudaCheckErrors("program error");
printf("results match, kernel 1 time: %fs, kernel 2 time: %fs\n", dt1/(float)USECPSEC, dt2/(float)USECPSEC);
// time row reduction kernel
unsigned long long dt3 = dtime_usec(0);
sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
cudaDeviceSynchronize();
dt3 = dtime_usec(dt3);
printf("row reduction kernel time: %fs\n", dt3/(float)USECPSEC);
cudaCheckErrors("program error");
}
$ nvcc -arch=sm_52 -o t1263 t1263.cu
$ ./t1263
results match, kernel 1 time: 0.459971s, kernel 2 time: 0.013678s
row reduction kernel time: 0.013724s
$
笔记：
N结果之后）没有产生正确的行和。 。在对你的索引进行一点研究之后，我对出了什么问题有所了解。一个示例问题是，对于N不能被D整除，你的内核d变量在第一个“页面”之后不会重置为零，但这不是唯一的问题。根据第4项，这里是修改数据初始化的代码版本，以及所有N * D结果的完整测试。数据初始化使得第一页的第一列全部为零，下一列全部为1，下一列全部为2，等等。在第二页上，我们将所有内容增加1，因此第一列将全部1，第二列将全部为2，等等。因此，应该很容易就列总和应该是什么达成一致。对于第一页，列总和应为0,10000,20000等。对于第二页，它们应为10000,20000,30000等。在第二页的第一列，我的代码产生10000，您的代码产生1.通过注释中更改的索引，我为第一页的第一列生成0，并且您的代码生成9999.根据我描述的数据初始化，1和9999不可能是有效的列总和：
$ cat t1263.cu
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
const int my_N = 10000;
const int my_D = 3;
const int my_blockSize = 768;
const int my_loopSize = my_N*my_N*my_D;
const int my_numBlocks = (my_loopSize + my_blockSize -1)/my_blockSize;
const int bsize = 512;
const float TOL = 0.1f;
#define HANDLE_ERROR(x) x
#define cudaCheckErrors(msg) \
do { \
cudaError_t __err = cudaGetLastError(); \
if (__err != cudaSuccess) { \
fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
msg, cudaGetErrorString(__err), \
__FILE__, __LINE__); \
fprintf(stderr, "*** FAILED - ABORTING\n"); \
exit(1); \
} \
} while (0)
#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL
long long dtime_usec(unsigned long long start){
timeval tv;
gettimeofday(&tv, 0);
return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}
namespace kernel {
__global__ void sumNND(float* devPtrIn, float* devPtrOut, const int N, const int D) {
int index = blockIdx.x * blockDim.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;
for (int id = index; id < N * N * D; id += stride) {
const unsigned int d = id % D; // 0 1 2 0 1 2 0 1 2
const unsigned int i = (id - d) / D; // 0 0 0 1 1 1 2 2 2
const unsigned int n = i / N; // 0 0 0 0 0 0 0 0 0
const unsigned int m = i % N; // 0 0 0 1 1 1 2 2 2
atomicAdd(&devPtrOut[d + D * n], // 0 1 2 0 1 2 0 1 2
devPtrIn[d + D * n + N * m]); // 0 1 2 0+N 1+N 2+N 0+2N 1+2N 2+2N
}
}
}
void sumNND(const int numBlocks, const int blockSize, float* devPtrIn, float* devPtrOut, const int N, const int D) {
HANDLE_ERROR(cudaMemset(devPtrOut, 0, N * D * sizeof(float)));
kernel::sumNND<<<numBlocks, blockSize>>>(devPtrIn, devPtrOut, N, D);
HANDLE_ERROR(cudaDeviceSynchronize());
}
// kernel assumes 1 block assigned per row, use block-striding methodology
// assumes block size is a power of 2
__global__ void sum_rows_NND(const float * __restrict__ devPtrIn, float * __restrict__ devPtrOut, const int N, const int D) {
__shared__ float sdata[bsize];
sdata[threadIdx.x] = 0;
for (int i = threadIdx.x; i < N; i += blockDim.x) // block-stride
sdata[threadIdx.x] += devPtrIn[(blockIdx.x * N) + i];
__syncthreads();
for (int i = blockDim.x>>1; i > 0; i>>=1){
if (threadIdx.x < i) sdata[threadIdx.x] += sdata[threadIdx.x+i];
__syncthreads();}
if (!threadIdx.x) devPtrOut[blockIdx.x] = sdata[0];
}
// kernel assumes one thread assigned per column sum
// launch N threads
__global__ void sum_cols_NND(const float * __restrict__ devPtrIn, float * __restrict__ devPtrOut, const int N, const int D) {
int idx = threadIdx.x+blockDim.x*blockIdx.x;
int ido = idx;
if (idx < N){
for (int j = 0; j < D; j++){
float temp = 0;
for (int i = 0; i < N; i++) temp += devPtrIn[idx + (i*N)];
devPtrOut[ido] = temp;
ido += N;
idx += N*N;}}
}
int main(){
float *h_data, *d_data, *h_res1, *h_res2, *d_res;
h_data = new float[my_loopSize];
cudaMalloc(&d_data, my_loopSize*sizeof(d_data[0]));
h_res1 = new float[my_N*my_D];
h_res2 = new float[my_N*my_D];
cudaMalloc(&d_res, my_N*my_D*sizeof(d_res[0]));
for (int i = 0; i < my_loopSize; i++) h_data[i] = i%my_N + i/(my_N*my_N); //rand()/(float)RAND_MAX;
cudaCheckErrors("CUDA failure");
cudaMemcpy(d_data, h_data, my_loopSize*sizeof(d_data[0]), cudaMemcpyHostToDevice);
// test original approach
cudaMemset(d_res, 0, my_N*my_D*sizeof(d_res[0]));
unsigned long long dt1 = dtime_usec(0);
kernel::sumNND<<<my_numBlocks, my_blockSize>>>(d_data, d_res, my_N, my_D);
cudaDeviceSynchronize();
dt1 = dtime_usec(dt1);
cudaMemcpy(h_res1, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);
//test columnwise reduction
unsigned long long dt2 = dtime_usec(0);
//sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
sum_cols_NND<<<(my_N + bsize -1)/bsize, bsize>>>(d_data, d_res, my_N, my_D);
cudaDeviceSynchronize();
dt2 = dtime_usec(dt2);
cudaMemcpy(h_res2, d_res, my_N*my_D*sizeof(d_res[0]), cudaMemcpyDeviceToHost);
// validate results
for (int i = 0; i < my_N*my_D; i++)
if (fabsf(h_res1[i] - h_res2[i]) > TOL) {printf("mismatch at %d, was %f, should be %f\n", i, h_res2[i], h_res1[i]); return -1;}
cudaCheckErrors("program error");
printf("results match, kernel 1 time: %fs, kernel 2 time: %fs\n", dt1/(float)USECPSEC, dt2/(float)USECPSEC);
// time row reduction kernel
unsigned long long dt3 = dtime_usec(0);
sum_rows_NND<<<my_N*my_D, bsize>>>(d_data, d_res, my_N, my_D);
cudaDeviceSynchronize();
dt3 = dtime_usec(dt3);
printf("row reduction kernel time: %fs\n", dt3/(float)USECPSEC);
cudaCheckErrors("program error");
}
$ nvcc -arch=sm_52 -o t1263 t1263.cu
$ ./t1263
mismatch at 10000, was 10000.000000, should be 1.000000
$
这取决于您的矩阵存储在哪个顺序以及您希望减少哪个维度。
目前，我将忽略D维度，因为操作可以被认为是减少包含NxN条目的矩阵，其中每个条目包含多个浮点数。
如果您的矩阵以行主要顺序存储，并且您希望将每一行减少到其总和（或列主要和列减少），答案很简单：
const int row = blockIdx.x * blockDim.x + threadIdx.x;
if (row < N) { // necessary if N is not divisible by the thread block size
float sum = 0; // stores the partial sum in a register
for (int col = 0; col < N; ++col) {
sum += devPtrIn[col + N * row];
}
devPtrOut[row] = sum; // no atomic operation necessary
}
这样，每个线程都以聚合方式读取内存（有关全局内存访问模式的讨论，请参阅NVIDIA's parallel forall blog），除最终结果外，不需要共享或全局内存写入。
如果你想沿着一个较小的维度减少 - 让我们说行 - 主矩阵上的列减少 - 答案变得有点困难：由于大步幅，如果我们只使用一个，内存访问会或多或少地表现得像随机访问一次输入列。
因此，每个线程并行地沿矩阵减少少量列并将部分结果存储在共享内存中是有意义的：
constexpr int numCols = ...;
__shared__ float partial[numCols * blockDim.x];
const int threadId = blockIdx.x * blockDim.x + threadIdx.x;
const int begin_col = threadId * numCols;
const int end_col = min(N, (threadId + 1) * numCols);
// initialize partial to 0
...
for (int row = 0; row < N; ++row) {
for (int col = begin_col; col < end_col; ++col) {
partial[threadIdx.x * numCols + col] += devPtrIn[col + N * row];
}
}
// store partial to global memory
...
根据GPU所具有的每个线程的寄存器数量，也可以通过展开内部循环并使用局部变量而不是数组来将部分和存储在寄存器中，因为arrays are usually not stored in registers
这样，我们总是从内存中读取连续的numCols浮点块，这比带有大步幅的访问提供了更大的带宽。
您可能不得不尝试使用numCols的最佳值，但它应该足够大，至少GPU内存的内存宽度用于加载这样的块，同时足够小以至于所有共享内存单个线程块适合GPU（有关详细信息，请再次参阅parallel forall）