c++ - 使用Cuda进行并行尺寸缩减(3D到2D求和)-6ren

c++ - 使用Cuda进行并行尺寸缩减(3D到2D求和)

转载作者：行者123 更新时间：2023-11-30 01:38:23

在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循环中。

就您而言，您似乎正在对列求和(但请参见下面的注释)。以下是一个有效的示例，将您的方法与运行速度更快的column-sum方法进行了比较，并结合了访问方式(相邻线程读取内存中的相邻元素):

$ 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
$

笔记:

经过优化的内核比您的朴素原子内核快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)\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
$

关于c++ - 使用Cuda进行并行尺寸缩减(3D到2D求和)，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/47993833/

文章推荐： c++ - 使用 QAudioOutput 进行实时流式传输

文章推荐： c++ - 使用 OpenCL，如何获取 GPU 内存使用情况？

文章推荐： c++ - 类型别名模板参数包

文章推荐： c++ - 从多字节字符集转换旧的 Visual Studio C++ 项目

c - 为什么这个 malloc 包装器对请求的内存大小执行此操作？尺寸 = (尺寸 + 3) & ~3;
我正在查看 DOOM 源代码，我找到了 this行。 void * Z_Malloc (int size, int tag, void *user)
window 尺寸
关闭。这个问题需要details or clarity .它目前不接受答案。想改进这个问题吗？通过 editing this post 添加细节并澄清问题. 关闭 9 年前。 Improve t
c++ - 如何获得特征矩阵的形状(尺寸)？
我正从 Python 和 Numpy 转向 C++ 和 Eigen。在 Python 中，我可以使用 .shape 属性获取 Numpy 数组/矩阵的形状(维度)，如下所示: import nump
youtube - 嵌入YouTube:尺寸
固定嵌入式YouTube视频的宽度并自己照顾自己的高度是否安全？我有一个应用程序，用户可以将通知发布到公告板上。这些通知主要是文本(带有有限的html标签)和嵌入式图像。我现在要添加对嵌入式YouT
javascript - 创建后如何更改盒子的尺寸/尺寸？
可以轻松创建一个 THREE.BoxGeometry，在创建宽度、高度和深度的三个独立参数时，您必须在其中传递参数。我想创建任何和所有不带参数的THREE[types]()，并在之后设置值。有没有
JavaScript Canvas 尺寸
我在 HTML 页面上有一个 Canvas : 属性width和height将 Canvas 拉伸(stretch)到某个字段，但不调整其大小。所以 var canvasElement = docu
CSS 定位/尺寸
我在我的 css 中使用 @media all 和 (max-width: 600px) {} 作为响应式菜单，问题是它没有正确显示。我想让橙色填充绿色空间……当然，还要将绿色空间变成透明的。基本上
c++ - 尺寸、长度换算？
(我知道我问了很多关于这个的问题!) 基本上，我正在尝试将一些代码从 Matlab 转换为 C++，我遇到了这个: n = sum(size(blocks)) - len; 现在我计算了 vector
javascript - 选择类别然后显示类别.尺寸
您好，我有一个用于创建产品的表单。用户应该能够选择类别(例如 T 恤)，然后 T 恤的所有尺码(例如 S、M、L)都会下拉。用户可以输入每种尺寸的数量。 Javascript 对此不起作用。用户可以选
html - 尺寸、对齐方式和列表样式等社交图标命令对页脚没有任何影响
我正在尝试在页脚中定位和调整我的社交图标链接的大小，但是，这些命令似乎都没有效果，尤其是当我尝试调整它们的大小时。我试过将宽度和高度标记为“!重要”，但这也没有效果。这是代码的 JSFiddle:h
canvas - 给定特定尺寸的最大潜在 PNG 尺寸
我目前正在创建一个 HTML5 canvas基于绘图程序。用户可以绘制一张图像或几张图像“页面”，并将其保存到云端以供日后快速检索。这是用于交互式白板的；老师不能总是确定他们计划类(class)使用的
image - 如何为网页保存图像？ (格式/尺寸)
为网站存储图像的最佳方式是什么？我不应该超过什么尺寸？现在，我将所有界面文件保存在 png(主要是 Sprite )中，并将常用图像保存在 jpg 中。一些图像大约为 100-150Kb。保存图
jquery - fancybox iframe 尺寸
在 fancybox 主页 ( http://fancybox.net/home ) 中，有一个打开尺寸为屏幕 75% 的 iFrame 的示例。我无法按照网站上的说明通过修改 .js 文件的宽度和
iPhone 4 Nib 尺寸
我想做一个仅适用于 iPhone 4 的应用程序，该应用程序使用 iAd AdBannerView。当我添加它时，它的固定大小为 320x50。在更高分辨率下这如何工作？有人可以解释一下 iPhon
cocoa - 旋转的 NSString 尺寸
我们有一个 NSString，我们使用 - (NSSize)sizeWithAttributes:(NSDictionary *)attributes 来测量边界框。一切都好。现在我们使用标准 NS
delphi - 如何仅通过知道其句柄来确定 Canvas 尺寸？
我想知道 Canvas 的宽度和高度，但我只知道它的 HDC。我尝试过这段代码: procedure TForm92.Button1Click(Sender: TObject); var hBi
javascript - 旋转路径后计算 SVG 尺寸
问题是如何使用数学从 START SVG 维度(不带旋转)和 END SVG 维度(带旋转)获取 >开始 SVG 信息。基本上，要从 START SVG 到 END SVG，我需要执行 -115.60
oracle - 鱿鱼进口 Parquet 尺寸
我的问题是，我有一个包含50万行的Oracle表。我设置了sqoop以将其作为 Parquet 文件导入到HDFS。我将--num-partition参数设置为32，得到了32个 Parquet 文件
wordpress - YouTube 尺寸 WordPress
是否可以更改 WordPress 中当前主题的 YouTube(或其他视频)的默认嵌入尺寸？我搜索了一个插件和一些代码，但似乎找不到。我的意思是当您将 YouTube 网址粘贴到帖子或页面中时使用的
JavaFX 缩小 ImageView 尺寸
我有一个组，其中包含一个矩形和顶部的图像。我希望矩形可以调整大小，并且图像应该具有固定大小，除非矩形小于图像的情况。然后图像应该随着矩形缩小。图像还应该始终居中并有一些填充。除了图像的缩小尺寸部分

行者123

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

c++ - 使用Cuda进行并行尺寸缩减(3D到2D求和)