第3章线程与执行模型 - SIMT 架构详解（Threads and Execution Model）

一、导读

1.1 本章学习目标

通过本章的学习，你将能够：

深入理解 SIMT（单指令多线程）执行模型的工作原理

掌握 warp 的概念及其在 GPU 执行中的作用

理解线程发散（Thread Divergence）的成因和性能影响

掌握多维线程层次结构的设计方法

学会选择合适的线程块大小以优化性能

理解 GPU 调度器如何管理大量并发线程

掌握 occupancy（占用率）的概念及其对性能的影响

能够分析和优化内核函数的并行效率

1.2 为什么本章如此重要

线程是 CUDA 编程的核心抽象。理解线程如何被组织、调度和执行，是编写高效 CUDA 程序的关键。许多性能问题——如线程发散、低占用率、负载不均衡——都源于对线程模型理解不足。

SIMT 架构是 GPU 高性能的秘密所在。它允许数千个线程同时执行，但又保持了编程模型的简洁性。然而，这种架构也带来了一些独特的挑战，特别是当线程的执行路径出现分歧时。

本章将帮助你建立对 GPU 执行模型的直觉，使你能够预测代码的性能特征，并做出明智的设计决策。无论你是要优化现有代码还是设计新的算法，这些知识都是不可或缺的。

1.3 核心问题引导

在深入学习之前，请思考以下问题：

什么是 warp？为什么 GPU 以 warp 为单位调度线程？

线程发散是什么？它为什么会降低性能？

如何选择最优的线程块大小？

什么是 occupancy？高 occupancy 一定意味着高性能吗？

GPU 如何在数千个线程之间快速切换？

线程块之间的执行顺序有保证吗？

为什么线程块的大小通常是 2 的幂次方？

如何避免或减轻线程发散的影响？

二、核心概念详解

3.1 SIMT 执行模型

SIMT（Single Instruction, Multiple Threads）是 NVIDIA GPU 的核心执行模型。它结合了 SIMD 的高吞吐量和多线程的灵活性。

SIMT 的基本原理

在 SIMT 模型中：

一组线程（称为 warp）同时执行相同的指令
每个线程有自己的程序计数器和寄存器状态
每个线程可以独立地处理不同的数据

这与传统的 SIMD 不同：在 SIMD 中，所有处理单元必须严格同步执行；而在 SIMT 中，线程可以有独立的控制流，尽管这会带来性能开销。

Warp：执行的基本单位

Warp 是 GPU 调度和执行的基本单位，包含 32 个线程。关键特性包括：

固定大小：每个 warp 始终包含 32 个线程，不可配置
同步执行：warp 内的所有线程同时执行同一条指令
共享程序计数器：warp 内所有线程共享一个程序计数器
硬件调度：warp 是 GPU 调度器分配给执行单元的最小单位

Warp 结构示意：
┌─────────────────────────────────────────┐
│              Warp (32 threads)           │
├──────┬──────┬──────┬──────┬─────┬──────┤
│Thread│Thread│Thread│Thread│ ... │Thread│
│  0   │  1   │  2   │  3   │     │  31  │
└──────┴──────┴──────┴──────┴─────┴──────┘
         All execute the same instruction

从线程到 Warp 的映射

线程块中的线程按顺序分组为 warps：

线程 0-31 组成 warp 0
线程 32-63 组成 warp 1
依此类推

如果线程块的线程数不是 32 的倍数，最后一个 warp 会有未使用的线程槽位，造成资源浪费。

3.2 Warp 调度与执行

Warp 调度器

现代 GPU 包含多个 warp 调度器，每个流多处理器（SM）通常有 2-4 个调度器。调度器的职责是：

选择就绪的 warp 进行执行
将指令分派给执行单元
处理内存访问请求
管理 warp 的状态转换

指令流水线

GPU 采用深度流水线来隐藏指令延迟：

指令执行阶段：
Fetch → Decode → Issue → Execute → Memory → Writeback

当一个 warp 等待内存访问时，调度器可以切换到另一个就绪的 warp，从而保持执行单元的利用率。

延迟隐藏

GPU 通过上下文切换来隐藏延迟：

内存访问可能需要数百个周期
在此期间，调度器执行其他 warp 的指令
当内存访问完成时，原 warp 恢复执行

这种机制要求有足够的活跃 warp 来填充空闲时间。这就是为什么高 occupancy 有助于性能的原因。

3.3 线程发散（Thread Divergence）

什么是线程发散？

当 warp 内的线程由于条件分支而执行不同的代码路径时，就发生了线程发散。例如：

cpp

__global__ void divergentKernel(float *data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n / 2) {
        data[idx] = data[idx] * 2.0f;  // 路径 A
    } else {
        data[idx] = data[idx] + 1.0f;  // 路径 B
    }
}

发散的性能代价

由于 warp 内的线程必须同步执行，当发生发散时：

GPU 串行执行每个分支路径

执行路径 A 时，路径 B 的线程处于空闲状态

执行路径 B 时，路径 A 的线程处于空闲状态

总执行时间是各分支时间之和

在最坏情况下（一半线程走路径 A，一半走路径 B），性能可能下降近 50%。

硬件如何处理发散

GPU 使用掩码（mask）来管理发散：

每个 warp 有一个活动线程掩码
执行某个分支时，只有该分支的线程被标记为活动
非活动线程不产生结果，但消耗时钟周期

避免发散的策略

重新组织数据：使条件相关的元素聚集在一起

使用算术代替分支：用三元运算符或数学技巧消除 if-else

平衡分支：确保各分支的计算量相近

接受必要发散：某些算法不可避免地需要分支，关键是最小化其影响

3.4 线程层次结构设计

一维 vs 多维线程配置

CUDA 支持最多三维的线程块和网格：

cpp

// 一维
dim3 block(256);
dim3 grid(n / 256);

// 二维
dim3 block(16, 16);
dim3 grid(width/16, height/16);

// 三维
dim3 block(8, 8, 8);
dim3 grid(nx/8, ny/8, nz/8);

多维配置的主要优势是：

更自然地映射到多维数据（图像、体积数据）
简化索引计算
提高代码可读性

线程索引计算

对于多维配置，全局索引的计算方式：

cpp

// 二维
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
int idx = row * width + col;

// 三维
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int z = blockIdx.z * blockDim.z + threadIdx.z;
int idx = z * width * height + y * width + x;

选择线程块大小的考虑因素

硬件限制：

- 最大线程数每块（通常 1024）

- 最大维度大小（通常 1024×1024×64）

资源约束：

- 每个线程的寄存器用量

- 共享内存需求

- SM 的资源上限

性能优化：

- 线程块大小应为 warp 大小（32）的倍数

- 考虑内存访问模式（合并访问）

- 平衡 occupancy 和资源使用

3.5 Occupancy（占用率）

定义

Occupancy 是指 SM 上活跃的 warp 数量与硬件支持的最大 warp 数量的比值：

Occupancy = 活跃 Warp 数 / 最大 Warp 数

影响 Occupancy 的因素

线程块大小：更大的线程块可以填满更多 warp 槽位

寄存器用量：每个线程使用的寄存器越多，能同时存在的线程越少

共享内存用量：每个块使用的共享内存越多，能同时驻留的块越少

Occupancy 与性能的关系

高 occupancy 有助于：

更好地隐藏内存延迟
提高指令吞吐量
容忍分支发散

但高 occupancy 不一定意味着最佳性能：

有时减少 occupancy 可以增加每个线程可用的资源
更好的内存访问模式可能比高 occupancy 更重要
某些算法不需要高 occupancy 就能达到良好性能

查询 Occupancy

cpp

int minGridSize;
int blockSize;
cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, 
                                    myKernel, 0, 0);

这个 API 可以帮助确定最优的线程块大小。

3.6 GPU 硬件架构概览

流多处理器（SM）

SM 是 GPU 的核心计算单元，包含：

CUDA 核心（流处理器）
Warp 调度器
寄存器文件
共享内存/L1 缓存
特殊功能单元（SFU）

内存层次结构

GPU 内存层次：
┌─────────────────────────────────────┐
│         Global Memory (GDDR/HBM)    │  ← 大容量，高带宽，高延迟
├─────────────────────────────────────┤
│      L2 Cache (shared across SMs)   │
├──────────┬──────────┬───────────────┤
│  SM 0    │  SM 1    │   ... SM N    │
│ ┌──────┐ │ ┌──────┐ │               │
│ │L1/   │ │ │L1/   │ │               │
│ │Shared│ │ │Shared│ │               │
│ └──────┘ │ └──────┘ │               │
│ ┌──────┐ │ ┌──────┐ │               │
│ │Regs  │ │ │Regs  │ │               │
│ └──────┘ │ └──────┘ │               │
└──────────┴──────────┴───────────────┘

执行资源分配

每个 SM 可以同时管理多个线程块，但受限于：

最大线程数（通常 2048）
最大线程块数（通常 32）
寄存器总数
共享内存总量

三、程序员的视角

3.1 在实际项目中的应用

图像处理

二维线程配置天然适合图像处理：

cpp

dim3 block(16, 16);
dim3 grid((width+15)/16, (height+15)/16);

每个线程处理一个像素，线程索引直接对应像素坐标。

矩阵运算

矩阵乘法可以使用二维或三维线程配置：

cpp

dim3 block(16, 16);
dim3 grid((N+15)/16, (N+15)/16);

粒子系统

一维线程配置适合粒子更新：

cpp

dim3 block(256);
dim3 grid((numParticles+255)/256);

3.2 常见陷阱与误区

误区一：线程块越大越好

事实：过大的线程块可能导致资源耗尽，反而降低 occupancy。需要根据寄存器和共享内存用量平衡。

误区二：必须填满所有 warp 槽位

事实：如果数据规模不是 32 的倍数，最后一个 warp 会有空闲线程。这是正常的，只需确保边界检查正确。

陷阱：忽视线程发散

复杂的条件分支可能导致严重的性能下降。尽量重构代码以减少发散。

陷阱：假设线程块执行顺序

线程块可以以任何顺序执行，甚至可能并行执行。不要依赖线程块之间的执行顺序。

3.3 与其他层级的联系

本章的线程模型概念与后续章节密切相关：

第4章：内存访问模式与线程索引的关系
第5章：线程同步的范围和机制
第6章：流与并发执行
第8章：基于线程模型的性能优化

四、代码实战

4.1 演示线程发散

cpp

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

// 存在线程发散的内核
__global__ void divergentKernel(float *data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        if (data[idx] > 0.5f) {
            data[idx] = data[idx] * 2.0f;
        } else {
            data[idx] = data[idx] + 1.0f;
        }
    }
}

// 重构后减少发散的内核
__global__ void optimizedKernel(float *data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        float val = data[idx];
        float result = (val > 0.5f) ? (val * 2.0f) : (val + 1.0f);
        data[idx] = result;
    }
}

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

int main() {
    int n = 1 << 20;
    size_t size = n * sizeof(float);
    
    float *h_data = (float *)malloc(size);
    for (int i = 0; i < n; i++) {
        h_data[i] = (float)rand() / RAND_MAX;
    }
    
    float *d_data;
    CUDA_CHECK(cudaMalloc((void **)&d_data, size));
    CUDA_CHECK(cudaMemcpy(d_data, h_data, size, cudaMemcpyHostToDevice));
    
    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
    
    // 测试发散版本
    cudaEvent_t start, stop;
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));
    
    CUDA_CHECK(cudaEventRecord(start));
    divergentKernel<<<blocksPerGrid, threadsPerBlock>>>(d_data, n);
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    
    float divergentTime = 0;
    CUDA_CHECK(cudaEventElapsedTime(&divergentTime, start, stop));
    printf("Divergent kernel time: %.3f ms\n", divergentTime);
    
    // 重置数据
    CUDA_CHECK(cudaMemcpy(d_data, h_data, size, cudaMemcpyHostToDevice));
    
    // 测试优化版本
    CUDA_CHECK(cudaEventRecord(start));
    optimizedKernel<<<blocksPerGrid, threadsPerBlock>>>(d_data, n);
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    
    float optimizedTime = 0;
    CUDA_CHECK(cudaEventElapsedTime(&optimizedTime, start, stop));
    printf("Optimized kernel time: %.3f ms\n", optimizedTime);
    
    printf("Speedup: %.2fx\n", divergentTime / optimizedTime);
    
    CUDA_CHECK(cudaFree(d_data));
    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaEventDestroy(stop));
    free(h_data);
    
    return 0;
}

4.2 探索不同线程块大小的性能

cpp

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

__global__ void matrixMulKernel(float *C, const float *A, const float *B, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            sum += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

void testBlockSize(int blockSize, int N) {
    size_t size = N * N * sizeof(float);
    
    float *h_A = (float *)malloc(size);
    float *h_B = (float *)malloc(size);
    float *h_C = (float *)malloc(size);
    
    for (int i = 0; i < N * N; i++) {
        h_A[i] = 1.0f;
        h_B[i] = 1.0f;
    }
    
    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc((void **)&d_A, size));
    CUDA_CHECK(cudaMalloc((void **)&d_B, size));
    CUDA_CHECK(cudaMalloc((void **)&d_C, size));
    
    CUDA_CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice));
    
    dim3 block(blockSize, blockSize);
    dim3 grid((N + blockSize - 1) / blockSize, 
              (N + blockSize - 1) / blockSize);
    
    cudaEvent_t start, stop;
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));
    
    CUDA_CHECK(cudaEventRecord(start));
    matrixMulKernel<<<grid, block>>>(d_C, d_A, d_B, N);
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    
    float elapsed = 0;
    CUDA_CHECK(cudaEventElapsedTime(&elapsed, start, stop));
    
    printf("Block size %dx%d (%d threads): %.3f ms\n", 
           blockSize, blockSize, blockSize * blockSize, elapsed);
    
    CUDA_CHECK(cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost));
    
    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));
    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaEventDestroy(stop));
    free(h_A);
    free(h_B);
    free(h_C);
}

int main() {
    int N = 512; // 512x512 矩阵
    
    printf("Testing different block sizes for %dx%d matrix multiplication:\n", N, N);
    
    int blockSizes[] = {4, 8, 16, 32};
    for (int i = 0; i < 4; i++) {
        testBlockSize(blockSizes[i], N);
    }
    
    return 0;
}

4.3 查询设备 Occupancy

cpp

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

__global__ void sampleKernel(float *data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        data[idx] = data[idx] * 2.0f;
    }
}

int main() {
    int device;
    cudaGetDevice(&device);
    
    cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, device);
    
    printf("Device: %s\n", prop.name);
    printf("Compute Capability: %d.%d\n", prop.major, prop.minor);
    printf("Max threads per block: %d\n", prop.maxThreadsPerBlock);
    printf("Max threads per SM: %d\n", prop.maxThreadsPerMultiProcessor);
    printf("Max blocks per SM: %d\n", prop.maxBlocksPerMultiProcessor);
    printf("Number of SMs: %d\n\n", prop.multiProcessorCount);
    
    // 查询不同线程块大小的 occupancy
    printf("Occupancy analysis:\n");
    printf("%-12s %-15s %-15s\n", "Block Size", "Blocks/SM", "Occupancy");
    printf("--------------------------------------------\n");
    
    int blockSizes[] = {32, 64, 128, 256, 512, 1024};
    
    for (int i = 0; i < 6; i++) {
        int blockSize = blockSizes[i];
        
        int minGridSize, actualBlockSize;
        cudaOccupancyMaxPotentialBlockSize(&minGridSize, &actualBlockSize, 
                                            sampleKernel, 0, 0);
        
        int maxBlocksPerSM;
        cudaOccupancyMaxActiveBlocksPerMultiprocessor(&maxBlocksPerSM, 
                                                       sampleKernel, 
                                                       blockSize, 0);
        
        int maxWarpsPerSM = prop.maxThreadsPerMultiProcessor / 32;
        int activeWarps = maxBlocksPerSM * blockSize / 32;
        float occupancy = (float)activeWarps / maxWarpsPerSM * 100.0f;
        
        printf("%-12d %-15d %.1f%%\n", blockSize, maxBlocksPerSM, occupancy);
    }
    
    return 0;
}

4.4 三维线程配置示例

cpp

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

// 处理三维体积数据的内核
__global__ void volumeProcessKernel(float *volume, int nx, int ny, int nz) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    int z = blockIdx.z * blockDim.z + threadIdx.z;
    
    if (x < nx && y < ny && z < nz) {
        int idx = z * ny * nx + y * nx + x;
        // 简单的平滑操作
        volume[idx] = volume[idx] * 0.5f + 0.25f;
    }
}

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

int main() {
    int nx = 64, ny = 64, nz = 64;
    size_t size = nx * ny * nz * sizeof(float);
    
    float *h_volume = (float *)malloc(size);
    for (int i = 0; i < nx * ny * nz; i++) {
        h_volume[i] = (float)i / (nx * ny * nz);
    }
    
    float *d_volume;
    CUDA_CHECK(cudaMalloc((void **)&d_volume, size));
    CUDA_CHECK(cudaMemcpy(d_volume, h_volume, size, cudaMemcpyHostToDevice));
    
    // 三维线程配置
    dim3 blockSize(8, 8, 4);
    dim3 gridSize((nx + 7) / 8, (ny + 7) / 8, (nz + 3) / 4);
    
    printf("Grid size: %d x %d x %d\n", gridSize.x, gridSize.y, gridSize.z);
    printf("Block size: %d x %d x %d\n", blockSize.x, blockSize.y, blockSize.z);
    printf("Total threads: %d\n", 
           gridSize.x * gridSize.y * gridSize.z * 
           blockSize.x * blockSize.y * blockSize.z);
    
    volumeProcessKernel<<<gridSize, blockSize>>>(d_volume, nx, ny, nz);
    CUDA_CHECK(cudaDeviceSynchronize());
    
    CUDA_CHECK(cudaMemcpy(h_volume, d_volume, size, cudaMemcpyDeviceToHost));
    
    printf("Volume processing completed.\n");
    
    CUDA_CHECK(cudaFree(d_volume));
    free(h_volume);
    
    return 0;
}

五、深入思考

5.1 历史背景

SIMT 架构的设计受到了多种早期并行计算模型的影响。向量超级计算机（如 Cray）展示了 SIMD 的威力，但也暴露了其缺乏灵活性的问题。图形着色器的发展证明了可编程并行处理的可行性。

NVIDIA 的创新在于将线程级的灵活性与 SIMD 的效率结合起来。warp 作为执行单位的设计是一个巧妙的折中：足够大以实现高效率，足够小以保持灵活性。

5.2 前沿方向

当前的研究正在探索更细粒度的并行模型：

Independent Thread Scheduling：Volta 架构引入的新特性，允许 warp 内的线程更独立地执行
Cooperative Groups：提供比线程块更灵活的线程分组机制
Dynamic Parallelism：允许内核函数启动子内核

5.3 延伸阅读

NVIDIA PTX ISA Documentation - 了解底层执行模型
"Programming Massively Parallel Processors" Chapter 4
NVIDIA Developer Blog 上的 warp 相关文章

六、实验笔记

实验 1：测量线程发散的影响

目标：量化线程发散对性能的影响。

步骤：

创建一个内核，其中包含根据数据值的条件分支

准备两组数据：一组导致严重发散，一组最小化发散

比较两种情况下的执行时间

使用 nsight compute 分析 warp 执行效率

预期结果：发散严重的情况执行时间明显更长。

实验 2：探索 Occupancy 与性能的关系

目标：验证高 occupancy 是否总是带来高性能。

步骤：

使用不同的线程块大小运行相同的内核

记录每种配置的 occupancy 和执行时间

绘制 occupancy 与执行时间的关系图

分析是否存在最优 occupancy

思考题：为什么有时较低的 occupancy 也能获得良好性能？

七、本章小结

本章深入探讨了 CUDA 的线程与执行模型。我们学习了 SIMT 架构的核心概念，包括 warp 的组织方式、调度机制和线程发散的影响。

关键要点包括：

Warp 是 GPU 执行的基本单位，包含 32 个线程
线程发散会导致性能下降，因为 warp 必须串行执行不同分支
Occupancy 衡量 SM 上活跃 warp 的比例，高 occupancy 有助于隐藏延迟
线程块大小的选择需要平衡资源使用和并行度
多维线程配置可以更自然地映射到多维数据

理解这些概念对于编写高效的 CUDA 程序至关重要。在下一章中，我们将探讨 GPU 内存层次结构，这是另一个影响性能的关键因素。

八、习题

解释 SIMT 与 SIMD 的区别。

为什么 warp 的大小固定为 32 个线程？

描述线程发散的发生过程及其性能影响。

如果一个线程块有 256 个线程，它包含多少个完整的 warp？

什么是 occupancy？列举影响 occupancy 的三个因素。

编写一个内核，展示如何通过重新组织数据来减少线程发散。

解释为什么线程块的执行顺序是不确定的。

对于一个 1024×1024 的图像处理任务，设计合适的线程配置。

查询你当前使用的 GPU 的最大 warp 数和 SM 数量。

分析以下代码是否存在线程发散，如果有，如何优化？

cpp

if (threadIdx.x % 2 == 0) {
    data[idx] = a[idx] + b[idx];
} else {
    data[idx] = a[idx] - b[idx];
}

第3章 线程与执行模型 - SIMT 架构详解（Threads and Execution Model）