使用Cuda平行降维（3D到2D，总和）
Posted
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了使用Cuda平行降维（3D到2D，总和）相关的知识，希望对你有一定的参考价值。
在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)
",
msg, cudaGetErrorString(__err),
__FILE__, __LINE__);
fprintf(stderr, "*** FAILED - ABORTING
");
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
", i, h_res2[i], h_res1[i]); return -1;}
cudaCheckErrors("program error");
printf("results match, kernel 1 time: %fs, kernel 2 time: %fs
", 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
", 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个结果）之后，我们的结果不同。我很确定我的索引是正确的，但在花了一段时间试图解开你的索引后，我放弃了。我怀疑你的内核索引中有一个错误，如果你试图对列进行求和，并且所有上述假设都是正确的。
- 我还包括行求和内核的时序测量，它看起来完全不同，但产生几乎相同的时序。这是可以预期的，因为这些类型问题的最佳内核将受到内存带宽的限制，这在两种情况下都是相同的。最佳内核将以合并的方式加载所有数据一次。之后，行和与列总和机制对内核时间的影响相对较小。
- 通过对数据初始化的一个小修改，我认为很容易证明你的内核没有创建正确的索引，因此在第一个“页面”之后（即在第一个
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)
",
msg, cudaGetErrorString(__err),
__FILE__, __LINE__);
fprintf(stderr, "*** FAILED - ABORTING
");
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;
以上是关于使用Cuda平行降维（3D到2D，总和）的主要内容，如果未能解决你的问题，请参考以下文章
是否可以通过 2d-3d 点关系获得图像或相机位置的正面平行视图？
使用Python，PCA对iris数据集降维2维/3维并进行2D，3D散点图绘制（包括有图例&无图例，有标题Label&无标题Label）