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

问题描述

我的问题是关于获取某些相同长度数组的“总和”。例如，我总共有一个 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 = 254 毫秒| 422ms | 1365ms| 4361ms| 5399ms | 10635ms

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

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

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

Cublas 方式

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,d_in,d_out); //faster
getSum<<<1,M>>>(M,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）比较慢。我只是出于某些可能的原因在这里发帖，注册号？延迟？

解决方法

使用 CuBLAS 通常是一个非常好的主意，如果有您想要的专用函数，尤其是对于大型数据集，则应该首选。话虽如此，对于处理如此小的数据集的 GPU 内核来说，您的时间安排非常糟糕。让我们了解原因。

gridSize 越大，速度越快。我不知道原因。
getSum<<<M,1>>>(M,N,d_in,d_out);
getSum<<<1,M>>>(M,d_out);

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

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

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

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

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

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

更新：

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

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

arrays arrays arrays cublas cuda cuda kernel kernel sum sum