maths-cs-ai-compendium-zh/chapter 16: SIMD and GPU programming/04. GPU architecture and CUDA.md

# GPU架构与CUDA

*GPU通过提供数千个核心用于大规模并行计算，改变了AI。本文涵盖GPU与CPU的设计哲学对比、GPU存储层次、C++中的CUDA编程、SIMT执行模型、内存访问模式、同步、流、性能分析以及NVIDIA GPU代次——编写和理解GPU核函数所需的知识。*

- 有关带有完整工作示例的实践CUDA教程，请参见配套仓库：[github.com/HenryNdubuaku/cuda-tutorials](https://github.com/HenryNdubuaku/cuda-tutorials)。

- 现代NVIDIA GPU有超过10,000个CUDA核心。CPU有4-128个核心。这100-1000倍的核心优势是GPU主导ML的原因：训练一个Transformer需要数万亿次乘加操作，GPU以CPU无法匹敌的规模并行处理它们。

- 即使你从不自己编写CUDA核函数，理解GPU架构也能解释：为什么批次大小很重要（需要足够的工作来饱和GPU），为什么内存通常是瓶颈（而非计算），以及为什么某些操作（分散、条件分支）在GPU上很慢。

## GPU vs CPU：根本不同的设计

- CPU是为**延迟**设计的：最小化完成一个任务的时间。它将其晶体管预算的大部分用于缓存、分支预测器和乱序执行——所有让单一线程快速运行的技巧。

- GPU是为**吞吐量**设计的：最大化每秒完成的任务数量。它将大部分晶体管用于执行单元（ALU）。单个线程很慢，但有数千个。

| | CPU | GPU |
|--|-----|-----|
| 核心 | 4-128（复杂、快速） | 1,000-20,000（简单、慢速） |
| 时钟频率 | 3-5 GHz | 1-2.5 GHz |
| 缓存 | 大（32 MB+ L3） | 小（每SM共享内存） |
| 分支预测 | 精密 | 无（所有线程遵循相同路径） |
| 最适合 | 低延迟、复杂控制流 | 高吞吐量、数据并行工作 |
| 典型FLOPS（FP32） | 1-5 TFLOPS | 30-80 TFLOPS |
| 内存带宽 | 50-100 GB/s | 1-3 TB/s |

- GPU的内存带宽优势（10-30倍）通常比其计算优势更重要。许多ML操作是内存受限的（逐元素操作、归一化、注意力），GPU的带宽使其能够足够快地向核心输送数据。

## GPU存储层次

- 理解GPU内存至关重要，因为**内存访问是主要瓶颈**，而非计算。

| 内存 | 大小 | 延迟 | 带宽 | 作用域 |
|--------|------|---------|-----------|-------|
| 寄存器 | 每SM约256 KB | 0周期 | 最高 | 每线程 |
| 共享内存 | 每SM 48-228 KB | 约5周期 | 约20 TB/s | 每线程块 |
| L1缓存 | 每SM 128-256 KB | 约30周期 | | 每SM |
| L2缓存 | 4-96 MB | 约200周期 | 约6 TB/s | 全局 |
| 全局内存（HBM） | 24-192 GB | 约400周期 | 1-3.3 TB/s | 全局 |

- **寄存器**是最快但最有限的。每个线程有一组私有寄存器（通常最多255个）。每线程使用过多寄存器会降低**占用率**（可同时运行的线程更少）。

- **共享内存**是由程序员管理的缓存，由块中的所有线程共享。它是编写快速CUDA核函数的关键：将数据瓦片从慢速全局内存加载到快速共享内存，然后进行计算。这是主导GPU编程的**分块**模式。

- **全局内存（HBM）**：主GPU内存（VRAM）。大但慢（400周期延迟）。所有数据起始和结束于此。核函数优化的目标是尽量减少全局内存访问。

## CUDA编程模型

- CUDA（统一计算设备架构）是NVIDIA的GPU编程模型。你编写**核函数**：在GPU上运行的函数，由数千个线程同时执行。

### 层次结构：网格、块、线程

```
网格（整个启动）
├── 块 (0,0)
│   ├── 线程 (0,0)
│   ├── 线程 (1,0)
│   ├── 线程 (2,0)
│   └── ... （每块最多1024线程）
├── 块 (1,0)
│   ├── 线程 (0,0)
│   └── ...
└── ... （可能有数百万个块）
```

- **线程**：最小单位。每个线程在其块内有唯一ID（`threadIdx.x`）。
- **块**：一组可以共享内存和同步的线程。块ID：`blockIdx.x`。块大小：`blockDim.x`（最多1024线程）。
- **网格**：单个核函数启动的所有块。可以是1D、2D或3D。

- 每个线程计算其全局索引：`int idx = blockIdx.x * blockDim.x + threadIdx.x;`

### 你的第一个CUDA核函数

```cpp
// vector_add.cu — CUDA源文件（.cu扩展名）

#include <stdio.h>

// __global__ 标记此为GPU核函数（从CPU调用，在GPU上运行）
__global__ void vector_add(const float* a, const float* b, float* c, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {           // 边界检查（网格可能大于数据）
        c[idx] = a[idx] + b[idx];
    }
}

int main() {
    int n = 1 << 20;  // 约100万个元素
    size_t bytes = n * sizeof(float);

    // 分配主机（CPU）内存
    float *h_a = new float[n];
    float *h_b = new float[n];
    float *h_c = new float[n];

    // 初始化
    for (int i = 0; i < n; i++) {
        h_a[i] = 1.0f;
        h_b[i] = 2.0f;
    }

    // 分配设备（GPU）内存
    float *d_a, *d_b, *d_c;
    cudaMalloc(&d_a, bytes);
    cudaMalloc(&d_b, bytes);
    cudaMalloc(&d_c, bytes);

    // 将数据从CPU拷贝到GPU
    cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);

    // 启动核函数：每块256线程，足够的块覆盖n个元素
    int block_size = 256;
    int grid_size = (n + block_size - 1) / block_size;  // 上取整除法
    vector_add<<<grid_size, block_size>>>(d_a, d_b, d_c, n);

    // 将结果从GPU拷贝到CPU
    cudaMemcpy(h_c, d_a, bytes, cudaMemcpyDeviceToHost);

    // 验证
    printf("c[0] = %f（期望值 3.0）\n", h_c[0]);

    // 释放内存
    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
    delete[] h_a; delete[] h_b; delete[] h_c;

    return 0;
}
```

```bash
# 用NVIDIA编译器编译
nvcc -O3 -o vector_add vector_add.cu
./vector_add
```

- **CUDA中的关键C++概念**：
    - `__global__`：CUDA关键字，标记核函数。从CPU（主机）调用，在GPU（设备）上运行。
    - `<<<grid_size, block_size>>>`：核函数启动语法。指定使用多少块和线程。
    - `cudaMalloc` / `cudaFree`：分配/释放GPU内存（类似于`new`/`delete`，但针对GPU）。
    - `cudaMemcpy`：在CPU和GPU之间拷贝数据。这通常是最大的瓶颈（PCIe带宽约32 GB/s，而GPU内存带宽约3 TB/s）。

### 线程束与SIMT

- GPU以32个为一组称为**线程束**的组执行线程。一个线程束中的所有32个线程同时执行**相同指令**（单指令多线程——SIMT）。这是GPU的SIMD等效，但在线程级别。

- **线程束分歧**发生在同一线程束中的线程在`if`语句中走不同分支时。GPU不能在一个线程束中同时执行两条不同指令，因此它顺序执行两个分支，屏蔽掉不应参与的线程。这使性能减半（或更差）。

```cpp
// 糟糕：线程束分歧（同一线程束中的线程走不同路径）
if (threadIdx.x % 2 == 0) {
    c[idx] = a[idx] + b[idx];    // 偶数线程做这个
} else {
    c[idx] = a[idx] - b[idx];    // 奇数线程做这个（同一线程束，串行化）
}

// 更好：无分支（无分歧）
float sign = (threadIdx.x % 2 == 0) ? 1.0f : -1.0f;
c[idx] = a[idx] + sign * b[idx];  // 所有线程执行相同指令
```

### 内存合并

- **合并访问**：当连续的线程访问连续的内存地址时，GPU将它们组合成单个内存事务。这对性能至关重要。

```cpp
// 好：合并——线程0读a[0]，线程1读a[1]，...
c[idx] = a[idx] + b[idx];

// 坏：跨步——线程0读a[0]，线程1读a[步长]，...
c[idx] = a[idx * stride] + b[idx * stride];  // 步长 > 1 浪费带宽
```

- 对于一个32线程的线程束，合并访问在单次事务中加载128字节（32 × 4字节用于float32）。跨步访问需要多次事务，每次加载128字节但只使用一小部分。步长为32是最坏情况：每次事务加载128字节，但只有一个线程使用4字节（3%的利用率）。

### 共享内存与分块

- **分块模式**是最重要的GPU优化技术。其想法：将数据块从慢速全局内存加载到快速共享内存，进行计算，然后将结果写回。

```cpp
// 使用共享内存分块的矩阵乘法（简化版）
__global__ void matmul_tiled(const float* A, const float* B, float* C,
                              int M, int N, int K) {
    // A的一个瓦片和B的一个瓦片的共享内存
    __shared__ float tile_A[TILE_SIZE][TILE_SIZE];
    __shared__ float tile_B[TILE_SIZE][TILE_SIZE];

    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    float sum = 0.0f;

    // 遍历瓦片
    for (int t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; t++) {
        // 将A和B的一个瓦片加载到共享内存
        if (row < M && t * TILE_SIZE + threadIdx.x < K)
            tile_A[threadIdx.y][threadIdx.x] = A[row * K + t * TILE_SIZE + threadIdx.x];
        else
            tile_A[threadIdx.y][threadIdx.x] = 0.0f;

        if (col < N && t * TILE_SIZE + threadIdx.y < K)
            tile_B[threadIdx.y][threadIdx.x] = B[(t * TILE_SIZE + threadIdx.y) * N + col];
        else
            tile_B[threadIdx.y][threadIdx.x] = 0.0f;

        __syncthreads();  // 等待所有线程完成加载

        // 计算此瓦片的部分点积
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += tile_A[threadIdx.y][k] * tile_B[k][threadIdx.x];
        }

        __syncthreads();  // 在加载下一个瓦片前等待
    }

    if (row < M && col < N)
        C[row * N + col] = sum;
}
```

- **`__shared__`**：声明块内所有线程共享的内存（快速、片上）。
- **`__syncthreads()`**：一个屏障，等待块中所有线程到达此点。在写入共享内存和读取它之间必须使用（否则某些线程读取到过期数据）。
- **为什么分块有效**：没有它，每个线程每次乘法都从全局内存加载。有了分块，一个TILE_SIZE × TILE_SIZE的数据块被加载到共享内存一次，并被块中所有线程重用。重用因子为TILE_SIZE，将全局内存流量减少该因子。

## 流与并发

- 默认情况下，CUDA操作是顺序的：CPU启动一个核函数，等待它完成，然后启动下一个。**流**允许重叠：

```cpp
cudaStream_t stream1, stream2;
cudaStreamCreate(&stream1);
cudaStreamCreate(&stream2);

// 这些操作可以重叠：不同流并发执行
cudaMemcpyAsync(d_a, h_a, bytes, cudaMemcpyHostToDevice, stream1);
cudaMemcpyAsync(d_b, h_b, bytes, cudaMemcpyHostToDevice, stream2);

kernel1<<<grid, block, 0, stream1>>>(d_a, d_c);
kernel2<<<grid, block, 0, stream2>>>(d_b, d_d);
```

- 流将数据传输与计算重叠：当流1的核函数运行时，流2在拷贝数据。这隐藏了PCIe传输延迟，并保持GPU忙碌。

## 分析CUDA代码

```bash
# NVIDIA Nsight Compute：核函数级分析
ncu --set full ./my_program

# NVIDIA Nsight Systems：系统级时间线
nsys profile ./my_program

# 快速指标
ncu --metrics sm__throughput,dram__throughput ./my_program
```

- **需要关注什么**：
    - **占用率**：SM容量中被使用的比例。低占用率（< 50%）意味着线程太少，无法隐藏内存延迟。原因：每线程寄存器过多、每块共享内存过多。
    - **内存吞吐量**：与峰值带宽比较。如果你达到峰值带宽的50%以下，内存访问模式低效（非合并、存储体冲突）。
    - **计算吞吐量**：与峰值FLOPS比较。如果内存和计算吞吐量都低，核函数是延迟受限的（并行度不够）。

## 高级优化技术

- 除了合并和共享内存分块的基础知识外，高性能GPU（和CPU）代码还使用几种高级技术：

### 数据布局：AoS vs SoA

- **结构体数组（AoS）**：每个元素将所有字段存储在一起。`[{x,y,z}, {x,y,z}, {x,y,z}]`。
- **数组结构体（SoA）**：每个字段存储在自己的连续数组中。`{[x,x,x], [y,y,y], [z,z,z]}`。

```cpp
// AoS：对于SIMD/GPU不好（访问所有x值触及非连续内存）
struct Particle { float x, y, z, mass; };
Particle particles[N];
// particles[0].x, particles[1].x 相隔16字节

// SoA：对于SIMD/GPU好（所有x值连续）
struct Particles {
    float x[N], y[N], z[N], mass[N];
};
// x[0], x[1] 相隔4字节——非常适合合并访问和SIMD
```

- SoA对于数据并行工作负载（SIMD、GPU）几乎总是更快。AoS在你总是同时访问一个元素的所有字段时更好（在数值代码中很少见）。PyTorch张量本质上是SoA：每个特征是一个连续维度。

### 软件预取

- 可以告诉CPU在需要之前开始加载数据，隐藏内存延迟：

```cpp
#include <xmmintrin.h>  // for _mm_prefetch

for (int i = 0; i < n; i += 4) {
    _mm_prefetch((char*)(a + i + 64), _MM_HINT_T0);  // 预取之前64个元素
    // 用SIMD处理 a[i:i+4]
    __m128 va = _mm_load_ps(a + i);
    // ...
}
```

- 预取指令是一个提示：如果数据已在缓存中，它是空操作。如果不是，CPU在执行其他指令的同时开始在后台获取数据。预取距离（此示例中向前64个元素）应根据内存延迟和循环迭代时间进行调整。

### 核函数融合

- **核函数融合**将多个操作组合成一个核函数，以避免将中间结果写入内存。这是ML中最有影响力的单个GPU优化：

```
// 未融合：3次核函数启动，3次全局内存往返
y = matmul(x, W)     // 写y到全局内存
z = y + bias          // 读y，写z
out = relu(z)         // 读z，写out

// 融合：1次核函数启动，1次全局内存写入
out = fused_matmul_bias_relu(x, W, bias)  // y和z永不离开SRAM
```

- 对于内存受限操作（偏置加法、ReLU、层归一化），内存流量主导执行时间。融合完全消除了流量。PyTorch的`torch.compile`和Triton可以自动或通过最少努力实现融合。

### 混合精度核函数

- 使用较低精度（FP16、BF16、INT8）进行计算和较高精度（FP32）进行累加，达到两全其美：

```cpp
// 张量核心：乘FP16矩阵，在FP32中累加
// 每条张量核心指令：D（FP32）= A（FP16）× B（FP16）+ C（FP32）
nvcuda::wmma::mma_sync(c_frag, a_frag, b_frag, c_frag);
```

- FP16比FP32小2倍，因此它使内存带宽加倍（通常的瓶颈），并在缓存中容纳2倍的数据。张量核心以FP32 CUDA核心8-16倍的速度处理FP16。这就是为什么混合精度训练（第6章）提供2-3倍加速且精度损失最小。

### 内存池分配器

- `cudaMalloc` 很慢（每次调用约1毫秒），因为它与GPU同步。在每次迭代分配临时缓冲区的训练循环中，这会累积起来。

- **内存池**（PyTorch的缓存分配器、CUDA内存池）预先分配一大块GPU内存，并从其中子分配而无需系统调用：

```python
# PyTorch自动执行此操作——但理解原因很重要
# 每个 torch.empty() 从池中重用内存，无需cudaMalloc
temp = torch.empty(1024, 1024, device='cuda')  # 微秒，而非毫秒
```

- 这就是为什么PyTorch的 `torch.cuda.memory_allocated()` 和 `torch.cuda.max_memory_allocated()` 不同：allocated是当前使用的，max是峰值（池可能持有比当前使用更多的内存）。

### 分析指导的优化

- 不要盲目优化。**先分析**，识别瓶颈，优化那个，然后重新分析。屋顶线模型（文件01）告诉你瓶颈是内存还是计算：

    - **内存受限**（低算术强度）：优化数据布局（SoA）、融合核函数、使用较低精度、预取。
    - **计算受限**（高算术强度）：使用张量核心、增加并行度、使用更快指令（FMA）。
    - **延迟受限**（并行度不足）：增加占用率、减少寄存器使用、启动更多线程。

- 大多数ML工作负载是**内存受限的**。令人惊讶的推论：更快的GPU（更多FLOPS）通常没有帮助。更快的内存（HBM3 vs HBM2e）更有帮助。这就是为什么A100→H100升级不只是关于FLOPS——H100也有2倍的内存带宽。

## NVIDIA GPU代次

| 代次 | 年份 | 关键创新 | AI相关性 |
|---------|------|-------------|--------------|
| Pascal（P100） | 2016 | HBM2、NVLink | 第一代严肃的深度学习GPU |
| Volta（V100） | 2017 | **张量核心**（混合精度矩阵乘法） | 实现FP16训练，125 TFLOPS TF32 |
| Ampere（A100） | 2020 | TF32、稀疏性、第三代张量核心 | 312 TFLOPS TF32，结构性稀疏2:4 |
| Hopper（H100） | 2022 | **Transformer引擎**（FP8）、HBM3 | 989 TFLOPS FP8，动态精度切换 |
| Blackwell（B200） | 2024 | 第二代Transformer引擎、NVLink 5 | 2.5 PFLOPS FP4，多芯片设计 |

- **张量核心**是专用的矩阵乘法单元。单个张量核心指令在一个周期内计算4×4矩阵乘法（D = A×B + C）。常规CUDA核心需要64次FMA操作。张量核心就是为什么混合精度训练（float16计算，float32累加）如此快速。

- **Transformer引擎**（Hopper+）在单层内动态切换FP8和FP16精度，只在需要时选择更高精度。这最大化吞吐量而不牺牲模型质量。它专为Transformer架构（注意力+MLP）设计，后者主导现代AI。

## 编程任务（用nvcc编译）

1. 编写一个对数组应用ReLU的CUDA核函数。测量包括内存传输在内的时间。这教授核函数编写、cudaMalloc/cudaMemcpy以及主机↔设备传输瓶颈。
```cpp
// task1_relu.cu
// 编译：nvcc -O3 -o task1_relu task1_relu.cu

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

__global__ void relu_kernel(const float* input, float* output, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        output[idx] = input[idx] > 0.0f ? input[idx] : 0.0f;
    }
}

int main() {
    const int N = 1 << 24;  // 约1600万元素
    size_t bytes = N * sizeof(float);

    // 分配主机内存
    float* h_input  = (float*)malloc(bytes);
    float* h_output = (float*)malloc(bytes);
    for (int i = 0; i < N; i++) {
        h_input[i] = (float)(i % 100) - 50.0f;  // 正负混合
    }

    // 分配设备内存
    float *d_input, *d_output;
    cudaMalloc(&d_input, bytes);
    cudaMalloc(&d_output, bytes);

    // 计时完整流水线：拷贝到GPU、计算、拷贝回
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    cudaEventRecord(start);
    cudaMemcpy(d_input, h_input, bytes, cudaMemcpyHostToDevice);

    int block_size = 256;
    int grid_size = (N + block_size - 1) / block_size;
    relu_kernel<<<grid_size, block_size>>>(d_input, d_output, N);

    cudaMemcpy(h_output, d_output, bytes, cudaMemcpyDeviceToHost);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float ms = 0;
    cudaEventElapsedTime(&ms, start, stop);

    // 验证
    int errors = 0;
    for (int i = 0; i < N; i++) {
        float expected = h_input[i] > 0.0f ? h_input[i] : 0.0f;
        if (h_output[i] != expected) errors++;
    }

    printf("时间（含传输）: %.2f ms\n", ms);
    printf("带宽: %.1f GB/s\n", 2.0 * bytes / ms / 1e6);  // 读取+写入
    printf("错误: %d / %d\n", errors, N);

    cudaFree(d_input); cudaFree(d_output);
    free(h_input); free(h_output);
    return 0;
}
```

2. 在CUDA中使用共享内存编写分块矩阵乘法。将性能与朴素（非分块）版本进行比较。这教授共享内存、`__syncthreads`以及为什么分块重要。
```cpp
// task2_matmul.cu
// 编译：nvcc -O3 -o task2_matmul task2_matmul.cu

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

#define TILE 16

// 朴素矩阵乘法：每个线程计算C的一个元素
__global__ void matmul_naive(const float* A, const float* B, float* C, 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;
    }
}

// 分块矩阵乘法：使用共享内存减少全局内存访问
__global__ void matmul_tiled(const float* A, const float* B, float* C, int N) {
    __shared__ float sA[TILE][TILE];
    __shared__ float sB[TILE][TILE];

    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float sum = 0.0f;

    for (int t = 0; t < (N + TILE - 1) / TILE; t++) {
        sA[threadIdx.y][threadIdx.x] = (row < N && t*TILE+threadIdx.x < N)
            ? A[row * N + t*TILE + threadIdx.x] : 0.0f;
        sB[threadIdx.y][threadIdx.x] = (t*TILE+threadIdx.y < N && col < N)
            ? B[(t*TILE + threadIdx.y) * N + col] : 0.0f;

        __syncthreads();
        for (int k = 0; k < TILE; k++)
            sum += sA[threadIdx.y][k] * sB[k][threadIdx.x];
        __syncthreads();
    }

    if (row < N && col < N)
        C[row * N + col] = sum;
}

int main() {
    const int N = 1024;
    size_t bytes = N * N * sizeof(float);

    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A, bytes); cudaMalloc(&d_B, bytes); cudaMalloc(&d_C, bytes);

    // 初始化为1（容易验证：C应全为N）
    float* h_A = new float[N*N];
    for (int i = 0; i < N*N; i++) h_A[i] = 1.0f;
    cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_A, bytes, cudaMemcpyHostToDevice);

    dim3 block(TILE, TILE);
    dim3 grid((N+TILE-1)/TILE, (N+TILE-1)/TILE);

    // 基准测试朴素版
    cudaEvent_t start, stop;
    cudaEventCreate(&start); cudaEventCreate(&stop);

    cudaEventRecord(start);
    for (int i = 0; i < 10; i++)
        matmul_naive<<<grid, block>>>(d_A, d_B, d_C, N);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float naive_ms; cudaEventElapsedTime(&naive_ms, start, stop);

    // 基准测试分块版
    cudaEventRecord(start);
    for (int i = 0; i < 10; i++)
        matmul_tiled<<<grid, block>>>(d_A, d_B, d_C, N);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float tiled_ms; cudaEventElapsedTime(&tiled_ms, start, stop);

    double gflops_naive = 2.0 * N * N * N * 10 / naive_ms / 1e6;
    double gflops_tiled = 2.0 * N * N * N * 10 / tiled_ms / 1e6;

    printf("朴素版:  %.2f ms, %.1f GFLOPS\n", naive_ms/10, gflops_naive);
    printf("分块版:  %.2f ms, %.1f GFLOPS\n", tiled_ms/10, gflops_tiled);
    printf("加速比: %.1fx\n", naive_ms / tiled_ms);

    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    delete[] h_A;
    return 0;
}
```

3. 演示线程束分歧。编写一个核函数，其中同一线程束中的线程走不同分支，并与无分支版本比较。
```cpp
// task3_divergence.cu
// 编译：nvcc -O3 -o task3_diverge task3_divergence.cu

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

// 糟糕：线程束分歧——偶数/奇数线程走不同路径
__global__ void divergent_kernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        if (idx % 2 == 0) {
            data[idx] = data[idx] * 2.0f + 1.0f;
        } else {
            data[idx] = data[idx] * 0.5f - 1.0f;
        }
    }
}

// 好：无分支——所有线程执行相同指令
__global__ void branchless_kernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        float scale = (idx % 2 == 0) ? 2.0f : 0.5f;
        float offset = (idx % 2 == 0) ? 1.0f : -1.0f;
        data[idx] = data[idx] * scale + offset;
    }
}

int main() {
    const int N = 1 << 24;
    float* d_data;
    cudaMalloc(&d_data, N * sizeof(float));
    cudaMemset(d_data, 0, N * sizeof(float));

    int block = 256, grid = (N + block - 1) / block;

    cudaEvent_t start, stop;
    cudaEventCreate(&start); cudaEventCreate(&stop);

    // 分歧版
    cudaEventRecord(start);
    for (int i = 0; i < 100; i++)
        divergent_kernel<<<grid, block>>>(d_data, N);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float div_ms; cudaEventElapsedTime(&div_ms, start, stop);

    // 无分支版
    cudaEventRecord(start);
    for (int i = 0; i < 100; i++)
        branchless_kernel<<<grid, block>>>(d_data, N);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    float nodiv_ms; cudaEventElapsedTime(&nodiv_ms, start, stop);

    printf("分歧版:  %.2f ms\n", div_ms / 100);
    printf("无分支版: %.2f ms\n", nodiv_ms / 100);
    printf("加速比:    %.2fx\n", div_ms / nodiv_ms);

    cudaFree(d_data);
    return 0;
}
```