如何在cuda中获得"sum"个并行数组？
How to get "sum" of parallel arrays in cuda?
我的问题是关于获取一些相同长度数组的“总和”。比如我一共有一个M*N(100 * 2000)长度的float数组。我想获得每个 N(2000) 个浮点数的 M(100) 个总和值。我找到了两种方法来完成这项工作。一种是在 M 的 for 循环中使用 Cublas 函数，例如 cublasSasum。另一种是自写核函数，循环加数。我的问题是这两种方式的速度以及如何在它们之间进行选择。
对于Cublas方法，无论N(4000~2E6)有多大，耗时主要取决于循环次数M
自写的kennel函数，速度随N的变化很大。具体来说，如果N比较小，在5000以下，比Cublas的方式跑得快很多。然后时间消耗随着N的增加而增加。
N = 4000 |10000 | 40000 | 80000 | 1E6 | 2E6
t = 254 毫秒| 422 毫秒 | 1365毫秒| 4361毫秒| 5399 毫秒 | 10635 毫秒
如果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, N, d_in, d_out); //faster
getSum<<<1,M>>>(M, N, d_in, d_out);
这是blockSize时间参数扫描结果。 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。我只是想知道为什么 full occupancy(blockSize 1024) 比较慢。我只是 post 出于某些可能的原因，注册号？延迟？
使用 CuBLAS 通常是一个非常好的主意，如果有专门的功能可以直接满足您的需求，那么应该首选它，尤其是对于大型数据集。话虽如此，对于在如此小的数据集上工作的 GPU 内核来说，你的时间安排非常糟糕。让我们明白为什么。
Bigger gridSize is faster. I don't know the reason.
getSum<<<M,1>>>(M, N, d_in, d_out);
getSum<<<1,M>>>(M, N, d_in, d_out);
调用CUDA内核的语法是kernel<<<numBlocks, threadsPerBlock>>>。因此，第一行提交了一个包含 M 个线程块的内核。 不要那样做：这效率很低。事实上，CUDA programming manual 说：
The NVIDIA GPU architecture is built around a scalable array of multithreaded Streaming Multiprocessors (SMs). When a CUDA program on the host CPU invokes a kernel grid, the blocks of the grid are enumerated and distributed to multiprocessors with available execution capacity. The threads of a thread block execute concurrently on one multiprocessor, and multiple thread blocks can execute concurrently on one multiprocessor. As thread blocks terminate, new blocks are launched on the vacated multiprocessors. [...]
The multiprocessor creates, manages, schedules, and executes threads in groups of 32 parallel threads called warps. [...]
A warp executes one common instruction at a time, so full efficiency is realized when all 32 threads of a warp agree on their execution path. If threads of a warp diverge via a data-dependent conditional branch, the warp executes each branch path taken, disabling threads that are not on that path. Branch divergence occurs only within a warp; different warps execute independently regardless of whether they are executing common or disjoint code paths.
因此，第一次调用创建 M 个线程块，浪费了每个 warp 中 32 个可用的 31 个 CUDA 内核。这意味着您可能只会读取 GPU 峰值性能的 3%...
第二次调用创建一个 M 线程块。因为 M 不是 32 的倍数，所以很少有 CUDA 核心被浪费。此外，它只使用 1 个 SM，因为你只有一个块。现代 GPU 有几十个 SM（我的 GTX-1660S 有 22 个 SM）。这意味着您将仅使用一小部分 GPU 功能（几个 %）。更不用说内存访问模式不是连续的，这会进一步减慢计算速度...
如果您想更有效地使用 GPU，您需要提供更多的并行性并浪费更少的资源。您可以首先编写一个在 2D 网格上运行的内核，使用原子 执行 缩减。这并不完美，但比您的初始代码要好得多。您还应该注意连续读取内存（共享同一个warp 的线程应该read/write 一个连续的内存块）。
请在编写 CUDA 代码之前准确阅读 CUDA manual 或教程。它非常准确地描述了所有这些。
更新：
根据新信息，您正在试验 blockSize 的问题可能是由于内核中的 跨步内存访问 （更具体地说 N*i).跨步内存访问模式很慢，并且当跨步变大时通常更慢。在您的内核中，每个线程将访问内存中的不同块。如前所述，GPU（实际上是大多数硬件计算单元）针对访问 连续 数据块进行了优化。如果你想解决这个问题并获得更快的结果，你需要并行处理另一个维度（所以不是M而是N）。
此外，BLAS 调用效率低下，因为 CPU 上循环的每次迭代都会调用 GPU 上的内核。 调用内核会带来相当大的开销（通常从几微秒到大约 100 微秒）。这样循环调用几万次会很慢
我的问题是关于获取一些相同长度数组的“总和”。比如我一共有一个M*N(100 * 2000)长度的float数组。我想获得每个 N(2000) 个浮点数的 M(100) 个总和值。我找到了两种方法来完成这项工作。一种是在 M 的 for 循环中使用 Cublas 函数，例如 cublasSasum。另一种是自写核函数，循环加数。我的问题是这两种方式的速度以及如何在它们之间进行选择。
对于Cublas方法，无论N(4000~2E6)有多大，耗时主要取决于循环次数M
自写的kennel函数，速度随N的变化很大。具体来说，如果N比较小，在5000以下，比Cublas的方式跑得快很多。然后时间消耗随着N的增加而增加。
N = 4000 |10000 | 40000 | 80000 | 1E6 | 2E6
t = 254 毫秒| 422 毫秒 | 1365毫秒| 4361毫秒| 5399 毫秒 | 10635 毫秒
如果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, N, d_in, d_out); //faster
getSum<<<1,M>>>(M, N, d_in, d_out);
这是blockSize时间参数扫描结果。 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。我只是想知道为什么 full occupancy(blockSize 1024) 比较慢。我只是 post 出于某些可能的原因，注册号？延迟？
使用 CuBLAS 通常是一个非常好的主意，如果有专门的功能可以直接满足您的需求，那么应该首选它，尤其是对于大型数据集。话虽如此，对于在如此小的数据集上工作的 GPU 内核来说，你的时间安排非常糟糕。让我们明白为什么。
Bigger gridSize is faster. I don't know the reason.
getSum<<<M,1>>>(M, N, d_in, d_out);
getSum<<<1,M>>>(M, N, d_in, d_out);
调用CUDA内核的语法是kernel<<<numBlocks, threadsPerBlock>>>。因此，第一行提交了一个包含 M 个线程块的内核。 不要那样做：这效率很低。事实上，CUDA programming manual 说：
The NVIDIA GPU architecture is built around a scalable array of multithreaded Streaming Multiprocessors (SMs). When a CUDA program on the host CPU invokes a kernel grid, the blocks of the grid are enumerated and distributed to multiprocessors with available execution capacity. The threads of a thread block execute concurrently on one multiprocessor, and multiple thread blocks can execute concurrently on one multiprocessor. As thread blocks terminate, new blocks are launched on the vacated multiprocessors. [...]
The multiprocessor creates, manages, schedules, and executes threads in groups of 32 parallel threads called warps. [...]
A warp executes one common instruction at a time, so full efficiency is realized when all 32 threads of a warp agree on their execution path. If threads of a warp diverge via a data-dependent conditional branch, the warp executes each branch path taken, disabling threads that are not on that path. Branch divergence occurs only within a warp; different warps execute independently regardless of whether they are executing common or disjoint code paths.
因此，第一次调用创建 M 个线程块，浪费了每个 warp 中 32 个可用的 31 个 CUDA 内核。这意味着您可能只会读取 GPU 峰值性能的 3%...
第二次调用创建一个 M 线程块。因为 M 不是 32 的倍数，所以很少有 CUDA 核心被浪费。此外，它只使用 1 个 SM，因为你只有一个块。现代 GPU 有几十个 SM（我的 GTX-1660S 有 22 个 SM）。这意味着您将仅使用一小部分 GPU 功能（几个 %）。更不用说内存访问模式不是连续的，这会进一步减慢计算速度...
如果您想更有效地使用 GPU，您需要提供更多的并行性并浪费更少的资源。您可以首先编写一个在 2D 网格上运行的内核，使用原子 执行 缩减。这并不完美，但比您的初始代码要好得多。您还应该注意连续读取内存（共享同一个warp 的线程应该read/write 一个连续的内存块）。
请在编写 CUDA 代码之前准确阅读 CUDA manual 或教程。它非常准确地描述了所有这些。
更新：
根据新信息，您正在试验 blockSize 的问题可能是由于内核中的 跨步内存访问 （更具体地说 N*i).跨步内存访问模式很慢，并且当跨步变大时通常更慢。在您的内核中，每个线程将访问内存中的不同块。如前所述，GPU（实际上是大多数硬件计算单元）针对访问 连续 数据块进行了优化。如果你想解决这个问题并获得更快的结果，你需要并行处理另一个维度（所以不是M而是N）。
此外，BLAS 调用效率低下，因为 CPU 上循环的每次迭代都会调用 GPU 上的内核。 调用内核会带来相当大的开销（通常从几微秒到大约 100 微秒）。这样循环调用几万次会很慢