如何在 cuda 中获得并行数组的“总和”？

Posted 2023-03-23

【中文标题】如何在 cuda 中获得并行数组的“总和”？

【英文标题】：How to get "sum" of parallel arrays in cuda?

【发布时间】：2021-07-18 20:33:27

【问题描述】：

我的问题是关于获得一些相同长度数组的“总和”。例如，我总共有一个 M*N(100 * 2000) 长度的浮点数组。我想获得每个 N(2000) 浮点数的 M(100) 和值。我找到了两种方法来完成这项工作。一种是在 M 的 for 循环中使用 Cublas 函数，例如 cublasSasum。另一种是自写的核函数，循环加数。我的问题是这两种方式的速度以及如何在它们之间进行选择。

对于Cublas方法，无论N(4000~2E6)有多大，耗时主要取决于循环数M。

对于自写的犬舍功能，速度随N变化很大。具体来说，如果N很小，在5000以下，它比Cublas方式运行得快得多。则时间消耗随着N的增加而增加。

N = 4000 |10000 | 40000 | 80000 | 1E6 | 2E6

t = 254ms| 422ms | 1365毫秒| 4361ms| 5399 毫秒 | 10635毫秒

如果 N 足够大，它的运行速度会比 Cublas 方式慢得多。我的问题是我怎么能用 M 或 N 来决定我应该使用哪种方式？我的代码可能用于不同的 GPU 设备。我必须在扫描的参数中比较速度，然后“猜测”以在每个 GPU 设备中做出选择，还是可以从 GPU 设备信息中推断？

另外，对于核函数方式，我在决定blockSize 和gridSize 时也有问题。我必须在这里指出，我更关心的是速度而不是效率。因为内存有限。例如，如果我有 8G 内存。我的数据格式是 4 个字节的浮点数。 N是1E5。那么M最多为2E4，小于MaxGridSize。所以我有两种方法如下。我发现有一个更大的 gridSize 总是更好，我不知道原因。是关于每个线程的寄存器号的使用吗？但我认为在这个内核函数中每个线程不需要很多寄存器。

任何建议或信息将不胜感激。谢谢。

库布拉斯方式

for (int j = 0;j<M;j++)
    cublasStatus = cublasSasum(cublasHandle,N,d_in+N*j,1,d_out+j);

自写内核方式

__global__ void getSum(int M, int N, float* in, float * out)

    int i = threadIdx.x + blockIdx.x * blockDim.x;
    if(i<M)
        float tmp = 0;
        for(int ii = 0; ii<N; ii++)
            tmp += *(in+N*i+ii);
        
        out[i] = tmp;
    

更大的 gridSize 更快。不知道是什么原因。

getSum<<<M,1>>>(M, N, d_in, d_out); //faster
getSum<<<1,M>>>(M, N, d_in, d_out); 

这是一个blockSize-time参数扫描结果。 M = 1E4.N = 1E5。

cudaEventRecord(start, 0);
//blockSize = 1:1024;
int gridSize = (M + blockSize - 1) / blockSize;
getSum<<<gridSize1,blockSize1>>>...
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&time, start, stop);

看来我应该选择一个相对较小的blockSize，比如10~200。我只是想知道为什么完全占用（blockSize 1024）比较慢。我只是出于一些可能的原因在这里发帖，注册号码？延迟？

【问题讨论】：

如果要求和的值在内存中是连续的，您可以尝试 Thrusts reduce_by_key（请参阅 here）。它将在一次内核调用中完成操作，而不是 M，并且还会进行一些优化。对于理想的 Thrust 性能，您可能不会为键创建自己的数组，而是使用一些 Thrusts “花哨的迭代器”来创建您的键序列，这只是每个值的行号或列号。我想象transform_iterator 采用counting_iterator 并将其除以列/行大小。

此选项仅在您被允许使用 C++ 时可用。

实际上有一个官方的 Thrust 示例完全按照我的描述 here

而this 的答案展示了如何将 Thrust 与原始 CUDA 指针一起使用

@PaulG 很抱歉我无法使用thrust，因为我使用的是 Matlab mex。它现在不支持thrust。感谢您的意见。我可以用thrust 以exe 格式编写代码。

【参考方案1】：

使用 CuBLAS 通常是一个非常好的主意，如果有专门的功能可以直接满足您的需求，则应该首选，尤其是对于大型数据集。话虽如此，对于在如此小的数据集上工作的 GPU 内核而言，您的时间安排非常糟糕。让我们了解原因。

更大的 gridSize 更快。我不知道原因。 getSum<<<M,1>>>(M, N, d_in, d_out); getSum<<<1,M>>>(M, N, d_in, d_out);

调用 CUDA 内核的语法是kernel<<<numBlocks, threadsPerBlock>>>。因此第一行提交了一个带有 M 个 1 线程块的内核。 不要这样做：这效率很低。确实，CUDA programming manual 说：

NVIDIA GPU 架构围绕多线程流式多处理器 (SM) 的可扩展阵列构建。当主机 CPU 上的 CUDA 程序调用内核网格时，网格的块被枚举并分发到具有可用执行能力的多处理器。 一个线程块的线程在一个多处理器上并发执行，多个线程块可以在一个多处理器上并发执行。当线程块终止时，新块在空出的多处理器上启动。 [...] 多处理器以称为 warps 的 32 个并行线程组的形式创建、管理、调度和执行线程。 [...] 一个warp 一次执行一条公共指令，因此当一个warp 的所有32 个线程 就它们的执行路径达成一致时，就可以实现全部效率。如果一个warp的线程分歧通过一个依赖于数据的条件分支，那么warp会执行每个分支路径，禁用不在那个分支上的线程小路。分支分歧只发生在一个扭曲内；不同的 warp 独立执行，无论它们执行的是公共代码路径还是不相交的代码路径。

因此，第一次调用创建 M 1 个线程的块，浪费了每个 warp 中可用的 31 个 CUDA 内核，共 32 个。这意味着您可能只会读取 GPU 峰值性能的 3%...

第二次调用创建一个M 线程块。因为M 不是 32 的倍数，所以浪费了很少的 CUDA 内核。此外，它仅使用 1 个 SM，而不是 GPU 上的许多可用块，因为您只有一个块。现代 GPU 有几十个 SM（我的 GTX-1660S 有 22 个 SM）。这意味着您将只使用一小部分 GPU 功能（几 %）。更不用说内存访问模式不是连续减慢计算速度......

如果您想更有效地使用 GPU，则需要提供更多并行度并减少资源浪费。您可以首先编写一个在 2D 网格上工作的内核，该内核执行 使用原子的缩减。这并不完美，但比您的初始代码要好得多。您还应该注意连续读取内存（共享同一个 warp 的线程应该读/写一个连续的内存块）。

在编写 CUDA 代码之前，请仔细阅读 CUDA manual 或教程。它非常准确地描述了这一切。

---

更新：

根据新信息，您正在尝试使用blockSize 的问题可能是由于内核中的跨步内存访问（更具体地说是N*i）。跨步内存访问模式很慢，并且当跨度变大时通常会更慢。在您的内核中，每个线程将访问内存中的不同块。 GPU（实际上是大多数硬件计算单元）针对访问连续块数据进行了优化，如前所述。如果你想解决这个问题并获得更快的结果，你需要在另一个维度上并行工作（所以不是M，而是N）。

此外，BLAS 调用效率低下，因为 CPU 上循环的每次迭代都会调用 GPU 上的内核。 调用内核会带来相当大的开销（通常从几微秒到大约 100 微秒）。因此，在称为数万次的循环中执行此操作将非常慢。

【讨论】：

感谢您的建议。对于 Cublas，它有自己的内部优化方法，例如并行归约以获得稳定的性能。速度仅与循环时间有关。至于核函数，我关心的是速度，而不是效率。例如，如果我有一个 8Gb 的内存，所有的数据都是浮点格式。 N是1E5。那么M最多为2E4，小于允许的MaxGridSize。所以“调用创建 1 个线程的 M 个块”对我来说是可以的。这是我一次可以处理的最大数据量。在这种情况下，“更大的 gridSize 更快”对吗？

@ZhangWei “我关心的是速度，而不是效率”它们是密切相关的。如果您使用 3% 的 GPU 使用慢速标量操作，那么代码可能会慢 30 倍。只是为了让您了解在 CPU（而不是 GPU）上可以获得的速度，计算 8GB 数据应该在普通台式机上不到一秒的时间内完成（我的为 0.25 秒）。我的中端 GPU 应该能够在不到 0.05 秒的时间内完成计算，而我的理解显然需要几秒钟的时间。因此，在速度方面，效率非常重要。

@ZhangWei 至于你的问题，这实际上非常依赖于目标GPU。更大的网格大小通常会更快，因为您提供了更多的并行性，但并不总是更快，因为必须考虑其他参数。我认为你的情况应该没问题。

@ZhangWei 如果只是为了执行这个计算而将数据逐块发送到GPU，数据传输会大大降低计算速度。在这种情况下，CPU（无论是什么 GPU）上的计算总是会更快，因为 CPU 足够快，并且您使用具有良好算法的快速编译代码（现代台式计算机就足够了）。

再次感谢您。我同意你的观点，我可以在 CPU 中完成这项工作。数据量太小。困扰我的问题是这些8G数据是一个中端数据。它是在 GPU 中生成的。如果我将其全部从 GPU 转移到 CPU，则需要很长时间。所以我需要它来在 GPU 中求和。我做了一个blockSize参数扫过。看来我应该采取更小的 blockSize(