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

Question

我的问题是关于获得一些相同长度数组的“总和”。例如，我总共有一个 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）比较慢。我只是出于一些可能的原因在这里发帖，注册号码？延迟？

Answer 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 微秒）。因此，在称为数万次的循环中执行此操作将非常慢。