使用cuda的平行尺寸降低(3d至2d）
Parallel dimension reduction (3D to 2D with sum) using Cuda
在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程序员的前两个优化优先级是：
- 使用很多线程
- 有效使用内存
出于您的问题，您对第一个问题没有任何麻烦 - 它很容易分解为一组独立的问题，可以分配给许多并行线程。第二个优先级是您要专注的地方。关于全球记忆，这意味着我们应该尽可能争取合并的访问。我们应该特别注意阅读。
我需要做一些假设。我将假设您的维度组织是行，列，深度，并且您的数据存储在通常的c风格中，即行 - 马约尔存储。然后，使用这些假设，请求（在整个第一个（或第二）轴 中总结）正在有效地求和整个行或求和整列。如果您在cuda标签上进行一些搜索，则会找到两者的工作示例（这就是一个这样的示例）。尽管它们不一定全部涵盖3D案例，但它们应该提供一个很好的路线图。您会发现的是，这两种情况应以不同的方式处理，并着眼于融合的全球内存访问，即已经提到的优化优先级。行向也是合并方向，因此，如果我们需要总和行，那么我们需要使用经典的并行还原技术，以便我们可以读取行并将元素汇总在一起。如果我们需要汇总列，则有效的内核更容易编写。每个线程都可以负责列，并且只能将运行总和保存在循环中。
在您的情况下，您似乎正在求和列（但请参见下面的注释）。接下来是一个有效的示例，将您的方法与更快的运行柱和方法进行比较，并结合访问（相邻的线程在内存中读取相邻元素）：
$ 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 - ABORTINGn");
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 %fn", i, h_res2[i], h_res1[i]); return -1;}
cudaCheckErrors("program error");
printf("results match, kernel 1 time: %fs, kernel 2 time: %fsn", 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: %fsn", 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
$
注意：
优化的内核比您的天真原子核快30倍。我怀疑其中的很大一部分实际上不是原子的使用，而是无法透露的访问。较新的GPU上的全球原子可以很快。
我的内核和内核之间的元素列总和的第一个"页面"（nxn）（即第一个n结果匹配）。在第一页之后（第一N结果）之后，我们的结果有所不同。我很确定我的索引是正确的，但是花了一段时间试图解开您的索引，我就放弃了。我怀疑您的内核索引中有一个错误，如果您要概括列，并且上述所有假设都是TUME。
我还包括了对行召唤内核的定时度量，该计算看起来完全不同，但产生了几乎相同的时机。这是可以预期的，因为这些类型的问题的最佳内核将受到内存带宽的限制，这两种情况下都是相同的。最佳内核将以合并的方式加载所有数据。之后，行 - 和列和列技术对内核时间的影响相对较小。
对数据的初始化进行了微小的修改，我认为很容易证明您的内核不会创建正确的索引，因此在第一个"页面"之后没有产生正确的行 - 和正确的索引（即。在第一个
N结果之后）。经过一些对您的索引的研究，我对出了什么问题有所了解。一个示例问题是，对于N，D不可分解的是，您的内核d变量不会在第一个"页面"之后重置为零，但这不是唯一的问题。
根据项目4，以下是已修改数据初始化的代码版本，以及所有N * D结果的完整测试。数据初始化是使第一页的第一列全部为零，下一个列全部为1，下一列全部2等。在第二页上，我们将所有内容递增1，因此第一列将全部1，第二列将是全部2等。因此，应该很容易就列总和应达成。对于第一页，列总和应为0、10000、20000等。对于第二页，它们应为10000、20000、30000等。在第二页的第一列中，我的代码生成100001.在评论中使用了您的索引，我为第一页的第一列生产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 - ABORTINGn");
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 %fn", i, h_res2[i], h_res1[i]); return -1;}
cudaCheckErrors("program error");
printf("results match, kernel 1 time: %fs, kernel 2 time: %fsn", 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: %fsn", 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条目的矩阵，其中每个条目都包含多个浮点。
如果您的矩阵存储在 row-major order中，并且要将每一行降低到其总和（或 column-major 和列降低），则答案是简单：
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的平行forall博客，以讨论全局内存访问模式），并且除最终结果外，不需要共享或全球记忆。
。如果要沿次要维度减少 - 假设列在行矩阵上的列减少 - 答案变得更加困难：由于大步很大，如果我们的记忆访问或多或少会像随机访问一样，如果我们一次仅使用了列的一个条目。
因此，每个线程并行减少矩阵的少量列并将部分结果存储在共享内存中是有意义的：
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每个线程的寄存器数量，也可以通过展开内部循环和使用局部变量而不是数组来存储部分总和，因为数组通常不存储在寄存器
以这种方式，我们总是读取numCols从内存中浮动的连续块，这比访问大步很大。
您可能必须对numCols的最佳值进行实验，但是它应该足够大，以至于至少使用GPU内存的内存宽度来加载这样的块，同时又足够小，以至于所有单个线程块的共享内存拟合到GPU上（有关详细信息，请参见Parallel Forall）
相关文章：
- 如何声明一个可以在整个程序中使用的全局 2d 3d 4d .. 数组(堆版本）变量?
- 如何在 3D OpenGL 场景上绘制 SDL 2D 矩形?
- 根据组件 1、2(和 3）对空间 (2D/3D） 向量进行排序
- 删除使用 new 创建的 2D 或 3D 指针
- 将 3D 变换转换为 2D
- CUDA，使用2D和3D阵列
- 用于将笛卡尔 (x，y，z） 转换为圆柱坐标 (ρ，θ，z） 坐标 2D/3D 的代码
- 从 C++ 中的 3D 矩阵中提取 2D 矩阵
- 使用 pcl 或 opencv 在 2D 场景中匹配 3D 模型
- 将 2D 图像坐标转换为 z = 0 的 3D 世界坐标
- 将每个对象渲染到空间中，并在3D场景上渲染2D场景
- 使用 2D 数据创建 3D 直方图(OpenCV?
- 使用 OpenCV 在 3D 矩阵中导入 2D 矩阵
- 将 2D 屏幕坐标转换为 3D 世界坐标
- OpenCV从2D像素获取3D坐标
- 使用cuda的平行尺寸降低(3d至2d）
- 将2D数组的实例分配给3D数组中的索引
- OPENCV(C ） - 从已知的3D对象和相机位置计算图像的2D坐标
- 如何使用3D char数组调整双指针2D char阵列的大小
- 在单目视觉里程计中使用 3D-2D 运动估计的相机位置
