继续学CUDA编程（三）——CUDA执行模型（上）

这一部分是 CUDA 的核心部分，涉及到了硬件和程序的执行模型。

CUDA 的执行层级是 Grid → Block → Warp → Thread，而真正被硬件调度的基本单位其实是 Warp，而不是 Thread。

然后必须先记住 CUDA 的约定：

threadIdx.x → 列（col）
threadIdx.y → 行（row）

SM

流式多处理器（Streaming Multiprocessor，SM）是构建整个 GPU 的核心模块。GPU 的硬件并行，是通过复制了多个 SM 来实现的。

一个 Block 只能在一个 SM 上被调度，Block 不会被拆分到多个 SM。

下图包含了 SM 的关键组件：

CUDA 核心（Core）
共享内存/一级缓存（Shared Memory/L1 Cache）
寄存器文件（Register File）
加载/存储单元（LD/ST）
特殊功能单元（SFU）
线程束调度器（Warp Scheduler）

CUDA 采用单指令多线程（SIMT）架构，每 32 个线程组成一个 Warp。在理想情况下，一个 Warp 内所有活跃线程在同一时刻执行同一条指令，但各自处理不同的数据。当程序出现条件分支时，Warp 调度器会依次执行各个分支路径，并临时屏蔽未参与当前路径的线程，这种现象称为 Warp 分化（Warp Divergence）。

Warp

一个 Warp 由 32 个连续的 Thread 组成，在一个 Warp 中，所有 Thread 都按照 SIMT 的方式执行。虽然 Block 可以是一维、二维、三维的，但是在硬件的角度看，所有 Thread 都是一维的。

在一个二维 Block 中，每个 Thread 的唯一索引都可以计算：

threadIdx.y * blockDim.x + threadIdx.x

注意在 Grid 中，每个 Thread 的全局坐标分量与 Block 内索引的计算方式不同：

blockIdx.x * blockDim.x + threadIdx.x
blockIdx.y * blockDim.y + threadIdx.y

因此我们重新计算一个 Thread 分别在 Grid 和 Block 中的唯一索引，在二维 Grid 和二维 Block 的情况下：

每个 Thread 在 Block 中的位置为：threadIdx.y 行、threadIdx.x 列。
每个 Block 在 Grid 中的位置为：blockIdx.y 行和 blockIdx.x 列。
任一个 Thread 在 Block 中的唯一索引为：threadIdx.y * blockDim.x + threadIdx.x^{[1]}
任一个 Block 在 Grid 中的唯一索引为：blockIdx.y * gridDim.x + blockIdx.x^{[2]}
任一个 Thread 在 Grid 中的位置为：blockIdx.y * blockDim.y + threadIdx.y 行，blockIdx.x * blockDim.x+ threadIdx.x 列。
随后计算 Thread 在 Grid 中的唯一索引：[blockDim.x * blockDim.y] * [blockIdx.x + gridDim.x * blockIdx.y] + blockDim.x * threadIdx.y + threadIdx.x

上式可以理解为，当某 Block 位于 Grid 中的某个位置时，其前面有 [blockDim.x * blockDim.y] * bIdx 个 Thread，这个 bIdx 就是上面提到的 Block 在 Grid 中的唯一索引^{[2]}，随后再加上 Thread 在 Block 中的唯一索引，也就是公式^{[1]}中的索引。

线程束分化

Warp 分化的含义就是：在 同一个 Warp 中的 Thread 执行不同的指令，分化会导致性能明显下降。当不得不在算法中加入其他分支的时候，确定一个合理的分支粒度可以有效避免 Warp 分化。

比如下面两个函数：

__global__ void mathKernel1(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    float a = 0.0;
    float b = 0.0;
    if (tid % 2 == 0){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

__global__ void mathKernel2(float *c){
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    float a = 0.0;
    float b = 0.0;
    if ((tid / warpSize) % 2 == 0){
        a = 100.0f;
    }
    else{
        b = 200.0f;
    }
    c[tid] = a + b;
}

在 mathKernel1 中，tid 为奇数的线程执行 else，tid 为偶数的线程执行 if。由于相邻线程的 tid 奇偶交替，同一 Warp 内约一半线程走 if、另一半走 else，Warp 存在严重的分化。

然而在 mathKernel2 中，分支粒度是 Warp 大小的倍数。当有两个 Warp 时，第一个 Warp 内的线程编号 tid 为从 0 到 31，因此 tid / warpSize 都等于 0，执行 if。第二个 Warp 内的线程编号 tid 是从 32 到 63，tid / warpSize 都不等于 0，执行 else。当一个线程束中所有的线程都执行 if 或者都执行 else 时，不会导致性能下降。

有一个指标，分支效率，定义为未分化的分支和全部分支之比：

\textstyle \text{BranchEfficiency} = 100 \times \frac{\text{Branches} - \text{DivergentBranches}}{\text{Branches}}

当分支效率低于 100%，并不一定会导致程序效率降低，也就是必要条件，而非充分条件。CUDA 编译器是有优化功能的，很短的分支并不会对程序效率产生明显影响，但是很长的代码路径必定会导致 Warp 分化和明显的效率降低。

资源分配

Warp 的本地执行上下文包含以下资源：

程序计数器
寄存器（Register File）
共享内存（Shared Memory）

执行上下文（Context）指设备与特定进程相关联的所有状态，是管理 CUDA 程序中所有对象生命周期的容器。

SM 处理的每个 Warp 的执行上下文（程序计数器、寄存器等）在 Warp 的整个生命周期内都在芯片上。因此，从一个执行上下文切换到另一个执行上下文是没有成本的，并且在每个指令发出时，Warp 调度器都会选择一个线程准备好执行其下一条指令（Warp 的活动线程）并将指令发布给这些线程。

每个 SM 都有一组 32 位宽的寄存器（寄存器文件），可在各 Thread 之间分配。

同时有固定数量的共享内存，可以在 Block 中进行分配。

因此对一个 Kernel，同时存在于一个 SM 中的 Block 和 Warp 数量取决于 SM 中 可用且所需 的寄存器和共享内存数量。

每个 Thread 需要的寄存器越多，那么 SM 中的 Warp 就越少。即减少 Thread 所需寄存器数量，即可增加 SM 中的 Warp 数。

每个 Block 需要的共享内存越多，那么 SM 中可以被同时处理的 Block 就会变少。即减少每个 Block 所需的共享内存，即可同时处理更多 Block。SM 内的资源没办法处理一个完整 Block，那么 Kernel 将无法启动。

上面提到的计算资源会限制 SM 中常驻 Block 的数量。当资源被分配给 Block 时，这个 Block 就变成活跃的，其中的 Warp 也是活跃的。

根据执行情况，活跃的 Warp 分为三类：

选定的线程束：活跃执行的 Warp
阻塞的线程束：没有做好执行准备
符合条件的线程束：准备执行，但未执行

当 32 个 CUDA Core 可用且当前指令中所有参数已就绪时，Warp 满足执行条件。

延迟

CPU 核心是为了最小化少数几个线程的延迟而设计，GPU 则是为了处理大量并发且轻量级的线程以最大化吞吐量。指令延迟被定义为：指令发出到指令完成的时钟周期。

指令延迟可以分为两种：

算术指令延迟：一个算术操作开始到产生输出之间的时间
内存指令延迟：发送出加载或存储操作和数据到达目的地之间的时间

上图描述了 Warp0 阻塞执行流水线的情况，Warp 调度器选取其他 Warp 执行，当 Warp0 符合条件时再执行。

当每个时钟周期中的所有线程调度器，都有一个符合条件的 Warp，可以达到计算资源的完全利用，通过在其他常驻 Warp 中发布其他指令，可以隐藏指令的延迟。如果想估算隐藏延迟所需的活跃 Warp 数量，Little's Law 可以估算一个近似值，即延迟和吞吐量的乘积：

\textstyle \text{Warp Num} = \text{Delay} \times \text{Throughput}

吞吐量是已经达到的值，描述单位时间内任何形式的信息和操作的执行速度；带宽指理论峰值，描述单位时间内最大可能达到的数据传输量。

图中绿色的箭头是 Warp，可以理解为，只要 Warp 足够多，那么吞吐量就不会下降。

占用率

占用率指每个 SM 中活跃的 Warp 占最大 Warp 数量的比值：

\textstyle \text{占用率} = \frac{\text{活跃 Warp 数}}{\text{最大 Warp 数}}

使用 cudaGetDeviceProperties() 函数可以获取设备中每个 SM 的最大 Warp 数。

通过 \textstyle \cfrac{\text{maxThreadsPerMultiProcessor}}{32} 获得最大 Warp 数量。

#include <stdio.h>
#include <cuda_runtime.h>

int main(int argc, char* argv[]){
    int iDev = 0;
    cudaDeviceProp iProp;
    cudaGetDeviceProperties(&iProp, iDev);
    printf("----------------------------------------------------------\n");
    printf("Number of multiprocessors:                      %d\n", iProp.multiProcessorCount);
    printf("Total amount of constant memory:                %4.2f KB\n",
           iProp.totalConstMem / 1024.0);
    printf("Total amount of shared memory per block:        %4.2f KB\n",
           iProp.sharedMemPerBlock / 1024.0);
    printf("Total number of registers available per block:  %d\n",
           iProp.regsPerBlock);
    printf("Warp size                                       %d\n", iProp.warpSize);
    printf("Maximum number of threads per block:            %d\n", iProp.maxThreadsPerBlock);
    printf("Maximum number of threads per multiprocessor:  %d\n",
           iProp.maxThreadsPerMultiProcessor);
    printf("Maximum number of warps per multiprocessor:     %d\n",
           iProp.maxThreadsPerMultiProcessor / 32);
    return EXIT_SUCCESS;
}

返回结果为：

----------------------------------------------------------
Number of multiprocessors:                      8
Total amount of constant memory:                64.00 KB
Total amount of shared memory per block:        48.00 KB
Total number of registers available per block:  65536
Warp size                                       32
Maximum number of threads per block:            1024
Maximum number of threads per multiprocessor:  1536
Maximum number of warps per multiprocessor:     48

记住一些准则：

每个 Block 中的 Thread 数量是 32 的整数倍
每个 Block 要有 128 或 256 个 Thread，即不要太小
根据 Kernel 资源调整 Block Size
Block 的数量要远远多于 SM 数量，保证足够并行，减少指令延迟

CUDA Toolkit 中包含了一个电子表格，名为《CUDA GPU Occupancy Calculator》，但是这个东西已经不能用了，目前推荐使用 Nsight。VS Code 有 Nsight 插件，后面会花些时间专门研究这个调试工具怎么用。（GDB 你还不会啊喂！）

栅栏同步

共享内存可以被 Block 中的多个 Thread 访问，CUDA 假设设备是一个弱序（Weakly-ordered）的内存模型，即一个 CUDA 线程将数据写入共享内存的顺序，与另一个 CUDA 或主机线程观察到的该数据被写入内存的顺序不一定相同。那么，两个线程在没有同步的情况下对同一个内存位置进行读写将出现未定义的行为。

CUDA 提供障碍（Barrier）和内存栅栏（Memory Fences）来实现块内同步。在障碍中，所有调用的线程等待其余调用的线程到达障碍点。在内存栅栏中，所有调用的线程必须等到全部内存修改对其余调用线程可见时才能继续执行。

如果让多个线程互相合作完成一项任务，这要求线程间可以进行协调。栅栏相当于程序中的一个集合点，当结果需要在中间进行整合的时候经常需要使用；当一个线程需要等待其他线程时候，可以让线程运行到栅栏处，一旦所有线程到达这个栅栏，栅栏就撤销。

同步在两个级别进行：

System 级：等待 Host 和 Device 完成所有工作。对 Host 来说，许多 CUDA API 调用和所有 Kernel 启动不是同步的，需要使用 cudaDeviceSynchronize() 函数来阻塞 Host 程序，直到所有 CUDA 操作完成。
Block 级：等待一个 Block 中的所有 Thread 达到同一个点。由于一个 Block 中的 Warp 会以未定义的顺序执行，使用 CUDA 的 Block 局部栅栏可以同步，__device__ void __syncthreads(void) 可以在 Kernel 中标记同步点。该函数被调用时，同一个 Block 中的 Thread 必须等待，直至 Block 中所有 Thread 都达到这个同步点。不过由于它强制 Thread 空闲，可能导致性能下降。

需要注意的是，不同 Block 间，无线程同步。因此唯一的办法是在每个 Kernel 执行结束时使用全局同步点。

参考文章

[1] GPU 编程 9：共享内存 3 → 共享内存线程同步

[2] 【CUDA 基础】3.2 理解线程束执行的本质（Part II）

[3] 极智开发 | CUDA 线程模型与全局索引计算方式

继续学CUDA编程（三）——CUDA执行模型（上）

SM

Warp

线程束分化

资源分配

延迟

占用率

栅栏同步

参考文章

推荐阅读

可变参数模板和折叠表达式的工程示例

C++模板

评论