第 3 讲 · GPU 编程基础（二）：CUDA 操作与矩阵 / 张量计算 · GPU Programming Basics 2

对应 Maps to：CMU 11-868 / 11-968 Large Language Model Systems（Spring 2026，Lei Li），Lecture 3 · Wed 1/21 · GPU Programming Basics 2。 阅读材料：PMPP 第 4 版第 2 章（CUDA 程序结构）、第 3 章（多维网格与数据）；原课 HW1 说明（CUDA Programming：map、zip、reduce、矩阵乘法内核及其接入）；Mark Harris, Optimizing Parallel Reduction in CUDA（NVIDIA 教程）。 依据说明：本讲以原课第 3 讲公开的 slides（llmsys-03-gpu-programming2.pdf，32 页）为主线：CUDA 的基本操作（cudaMalloc、cudaMemcpy、发射内核、cudaFree）、GPU 的三种内存、H100 的 TMA 与异步执行、__global__ / __device__ / __host__、dim3 与内置变量、VecAdd / MatAdd / fullKernel 三个例子。§3 的四类内核按原课 HW1 的说明组织（slides 只给了 Colab 笔记本的链接）。模拟器上的全部实验、时间账、§4 的并行归约都是本讲义自己设计的，不代表课堂上讲了什么。 配套代码 Code：code/lecture03/（4 个脚本 + kernels.py；复用第 2 讲的 simt.py、gpu_model.py 和 P1 的 tensor_layout.py；每个脚本 4 秒以内，输出重复运行逐位相同） 预计阅读时间：4 小时（含运行代码）

0. 本讲导航 · How to read this lecture

0.1 学习目标 Learning goals

学完这一讲，你应该能够：

1. 写出一个完整 CUDA 程序的五步（分配显存、拷入、发射内核、拷出、释放），解释每个 API 的参数（特别是 cudaMalloc(&dA, size) 为什么要传 &dA），并能算出拷贝 vs 内核的时间账，判断一个运算值不值得搬到 GPU 上。
2. 说出 GPU 的三种内存（寄存器、共享内存、全局内存）各自的作用域和生命周期，以及 __global__、__device__、__host__ 的含义。
3. 用 dim3 配置一维 / 二维 / 三维的网格和线程块，写出任意维度下的全局线程编号，处理越界，并让线程束的访问合并（threadIdx.x 对应内存里连续的那一维）。
4. 用”存储 + 形状 + 步长”实现原课 HW1 的四类内核：map、zip（带广播）、reduce、矩阵乘法，并把它们接到主机端的张量上。
5. 理解并行归约：关键路径、树形归约的三种写法（分化与 bank 冲突的差别）、原子操作与数据竞争，以及为什么 GPU 上的浮点归约”每次可能不一样”。

0.2 先修与本讲范围 Prerequisites & scope

0.3 本讲地图 Map

§1  一个 CUDA 程序的五步：分配、拷入、发射、拷出、释放；三种内存；时间账；异步   ← 代码 01
§2  多维网格与线程编号：dim3、fullKernel、矩阵加法的 x / y、越界、块的形状       ← 代码 02
§3  张量计算：HW1 的 map、zip、reduce、矩阵乘法（步长与广播）；主机端封装      ← 代码 03
§4  并行归约：关键路径、树形归约的三种写法、原子操作与不确定性                   ← 代码 04
§5  小结   §6  常见误区   §7  练习   §8  延伸阅读   §9  下一讲预告

1. 一个 CUDA 程序的五步 · Basic CUDA operations

1.1 回顾：内核在哪里跑，数据在哪里放

第 2 讲讲了内核是怎样执行的：网格 → 线程块 → 线程束 → 线程，每个线程用 blockIdx、threadIdx 算出自己负责的数据。但还有一个问题没回答：数据是怎么到 GPU 上的？ 第 2 讲 §4.1 说过，主机（CPU）和设备（GPU）有各自的内存，内核只能读写显存。slides 第 11 页的图再一次提醒：CPU 的系统内存和 GPU 的显存之间隔着一条 PCIe（32 GB/s），而 GPU 访问自己显存的带宽是它的几十倍。

所以每个 CUDA 程序都要明确地安排数据的搬运。slides 把它总结成五个操作：

1. CPU 在 GPU 上分配显存：cudaMalloc
2. CPU 把数据拷到显存（host to device）：cudaMemcpy
3. CPU 发射内核
4. CPU 把结果拷回主机（device to host）：cudaMemcpy
5. 释放显存：cudaFree

1.2 五个 API

分配：cudaError_t cudaMalloc(void** devPtr, size_t size)。

int *dA;
cudaMalloc(&dA, n * sizeof(int));

为什么传 &dA？因为 cudaMalloc 要修改 dA 本身——把新分配的显存地址写进去；而它的返回值已经被用来返回错误码了。C 只有按值传参，要修改调用者的变量，就得传它的地址（P1 §3.4–§3.5）。dA 的类型是 int*，它的地址是 int**，统一转成 void** 传进去。slides 第 13 页画的正是这件事：CPU 内存里的变量 dA（地址 306）里存着一个显存地址（12）。

注意 dA 本身是主机上的一个变量，它的值是一个显存地址。主机代码不能解引用它（*dA 会读主机内存的地址 12，通常段错误）；只能把它传给内核，或传给 cudaMemcpy。

释放：cudaFree(void* devPtr)。slides 写着”不要忘了！“——显存不会因为函数返回而自动释放。

拷贝：cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind)。参数顺序是先目的、后来源（和 memcpy 一样），count 是字节数，kind 说明方向（cudaMemcpyHostToDevice、cudaMemcpyDeviceToHost）。

发射：VecAddKernel<<<num_blocks, threads_per_block>>>(dA, dB, dC, n);，其中 num_blocks = (n + threads_per_block - 1) / threads_per_block（向上取整，第 2 讲 §4.8）。

01 号脚本第 1 部分把 slides 第 29 页的 VecAddCUDA 一行行搬到模拟器上，每一步之后打印显存的状态：

结果和 NumPy 的 A + B 逐位相同，没有泄漏。

1.3 GPU 的三种内存

slides 第 15 页列出了内核能用的三种内存：

最后一行的”跨内核发射一直存在”很关键：一个内核写进全局内存的结果，下一个内核可以直接读。这让”数据留在显存里，被一串内核接力处理”成为可能（§1.5）。共享内存的用法是第 4 讲的主题。

1.4 时间账：什么样的运算值得搬到 GPU 上

01 号脚本第 2 部分按公开规格算了一笔账（每次拷贝的固定开销取 10 μs、每次发射内核取 5 μs，都是量级）。向量加法：拷入两个数组、内核读 8 字节写 4 字节、拷出一个数组：

拷贝占了九成以上。原因很简单：内核和拷贝搬的是同样多的字节，而显存带宽是 PCIe 的 24 倍（A6000）到 52 倍（H100）。只用一次的数据，专门搬到 GPU 上做一次加法，不如在 CPU 上直接做。

矩阵乘法呢？拷贝 $3n^2$ 个 float，计算 $2n^3$ 次——拷贝的每个字节对应约 $n/6$ 次运算。在 A6000 上按 FP32 峰值算：

即使是算术强度很高的矩阵乘法，n = 4096 时计算也只有拷贝的一半多；到 n = 16384 才是两倍多。用张量核心（BF16 快约 4 倍），拷贝的占比还会更大。PCIe 太慢了。

结论不是”GPU 没用”，而是：数据要留在显存里，被很多个内核反复使用。训练时权重、梯度、优化器状态全部常驻显存，每步从主机送进来的只有一小批 token 编号（第 2 讲 §2.2：66 KB）——正是这个道理。

1.5 数据留在显存里

01 号脚本第 3 部分做了一个对照：连续做 8 个逐元素运算 $y = a_j x + b_j$。

结果完全相同，时间差 6.2 倍。PyTorch 的设计就是第二种：张量一旦 .to('cuda')，之后所有运算的输入输出都在显存里；只有 .item()、.cpu()、print(tensor) 才会拷回主机——而且是同步的（见 §1.6），所以训练循环里随手写的 loss.item() 会让 CPU 停下来等 GPU。

再往前一步：8 个内核之间，中间结果还是要写回显存再读出来。如果把 8 个运算合成一个内核，每个元素只读一次、写一次——这就是算子融合（第 10–11 讲）。

1.6 异步执行与 TMA

上面的时间账假设每一步做完才开始下一步。真实的 CUDA 并不完全是这样：

• 内核发射对主机是异步的：kernel<<<...>>>(...) 只是把内核放进一个队列（流，stream）就立即返回，CPU 可以接着做别的事（比如准备下一批数据）。
• cudaMemcpy 是同步的：它会等之前队列里的内核全部完成、拷贝结束才返回。异步版本 cudaMemcpyAsync 配合多个流，可以让”拷贝下一批”和”计算这一批”重叠（P4 §3 的重叠，第 15 讲会用到）。

slides 第 17–18 页介绍了 Hopper（H100）的新硬件：张量内存加速器（Tensor Memory Accelerator，TMA）。在 H100 之前，把一块数据从全局内存搬到共享内存，要由线程们一条条地发读指令、再一条条地写共享内存；TMA 让一个线程发一条指令（cuda::memcpy_async(dst, src, shape, barrier)），由专门的硬件把一整块多维数据搬过去，搬完时通知一个屏障（barrier.wait()）。线程们在等待期间可以继续算上一块——搬运和计算在硬件层面异步重叠。FlashAttention-3（第 21 讲）就重度依赖 TMA 和这种异步流水线。本讲只需要记住它的动机：数据搬运太贵，要让专门的硬件去搬，算的单元专心算。

1.7 显存管理的错误，以及 __global__ / __device__ / __host__

01 号脚本第 4 部分在模拟器上犯了五个错误：

真 CUDA 里，这些错误大多不会立刻报出来：重复释放、释放后使用的后果往往是之后某个不相干的地方崩溃；方向写反时两个指针在 C 里都是 float*，编译器分辨不出。所以真实的 CUDA 代码会检查每一个 API 的返回值（cudaError_t），发射内核后调用 cudaGetLastError()。

slides 第 19 页的三个函数修饰符：

主机代码和设备代码写在同一个 .cu 文件里，nvcc 按修饰符把它们分开编译（第 2 讲 §4.9）。在模拟器里，__device__ 函数就是被内核调用的普通 Python 函数（脚本里的 device_helper）；HW1 的下标函数 to_index、index_to_position、broadcast_index 在 CUDA 版本里就是 __device__ 函数。

2. 多维网格与线程编号 · Multi-dimensional grids

2.1 dim3 与内置变量

slides 第 21–22 页：主机在发射时用 dim3 指定网格和线程块的大小，每一维都可以是 1–3 维：

dim3 Dg(4, 2, 1);           // 网格：4 × 2 × 1 = 8 个块
dim3 Db(8, 8, 1);           // 线程块：8 × 8 × 1 = 64 个线程（≤ 1024）
kernelFuncName<<<Dg, Db>>>(args);

内核里有四个内置变量，都有 .x .y .z 三个字段：gridDim（网格大小）、blockIdx（块在网格里的编号）、blockDim（块大小）、threadIdx（线程在块里的编号）。限制：每块最多 1024 个线程，blockDim.x、.y 最多 1024、.z 最多 64；网格的 x 方向最多 $2^{31}-1$，y、z 方向最多 65535。

2.2 fullKernel：三维编号的公式

slides 第 27–28 页的例子：blocks_per_grid(2, 3, 4)、threads_per_block(2, 4, 8)，一共 24 个块 × 64 个线程 = 1536 个线程。每个线程算出全局编号：

int block_id = blockIdx.x + blockIdx.y * gridDim.x + blockIdx.z * gridDim.x * gridDim.y;
int block_offset = block_id * blockDim.x * blockDim.y * blockDim.z;
int thread_offset = threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y;
int tid = block_offset + thread_offset;

02 号脚本第 1 部分让每个线程对 COUNT[tid] 做一次原子加，并打开了全局竞争检测：1536 个 tid 每个恰好被写了一次，0 次竞争——公式是一一对应的。它就是把 (blockIdx.z, .y, .x, threadIdx.z, .y, .x) 这六个数当作一个六维下标、按行优先展平（P1 §5）：x 变化最快，z 最慢。例如 blockIdx = (1, 2, 3)、threadIdx = (1, 3, 7) 的线程：block_id = 1 + 2·2 + 3·6 = 23，thread_offset = 1 + 3·2 + 7·8 = 63，tid = 23 × 64 + 63 = 1535，正是最后一个。

线程束按块内编号（thread_offset）每 32 个切一刀：块 (0,0,0) 的第 0 个线程束是 thread_offset 0–31，它们的 threadIdx 是 (0..1, 0..3, 0..3)——x 方向只有 2，所以一个线程束跨了 y 和 z 方向。块的哪一维放在 x 上，决定了哪些线程在同一个线程束里一起访存。

2.3 矩阵加法：x 对应行还是对应列？

slides 第 26 页的矩阵加法内核：

__global__ void MatAddKernel(float* A, float* B, float* C, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    C[i * N + j] = A[i * N + j] + B[i * N + j];
}

它让 x 方向的编号 i 当行号。02 号脚本第 2 部分把它和”行号来自 y、列号来自 x”的写法放在一起（N = 32，一个 32 × 32 的块）：

两种写法结果都对，但 slides 的写法搬运的字节是 8 倍。原因：一个线程束是 threadIdx.x = 0..31、threadIdx.y 相同的 32 个线程。slides 的写法里 x 是行号，这 32 个线程访问 32 个不同的行，地址相隔 N × 4 = 128 字节，每个线程独占一个 sector；改写后它们访问同一行的 32 个相邻元素，4 个 sector 就够了。

规则：行优先存储的矩阵，让 threadIdx.x（线程束里变化最快的那一维）对应列号（内存里连续的那一维）。这条规则在 HW1 的矩阵乘法提示和第 4 讲里都会反复出现。

2.4 越界：N 不是块大小的整数倍

slides 第 30 页的 MatAddCUDA 用 THREADS = 32、BLOCKS = (n + THREADS - 1) / THREADS，网格是 BLOCKS × BLOCKS。但 MatAddKernel 没有越界检查——slides 的例子 N = 32 正好整除。02 号脚本第 3 部分取 N = 40：网格 2 × 2 个块、4096 个线程，矩阵只有 1600 个元素。slides 的内核立刻越界（blockIdx=(1,0,0)、threadIdx=(8,0,0) 的线程访问 A[1600]）；加上 if (row < N && col < N) 之后结果正确，4096 个线程里有 2496 个（61%）什么也不做。

二维时要对行、列分别检查。多出来的线程是”向上取整”的代价，N 远大于块大小时可以忽略；矩阵很小或形状很”瘦”时（例如 1 × 4096 的矩阵用 32 × 32 的块）浪费会很大，这时应该换用一维的块，或者换一种分工。

2.5 二维线程块的形状

同样每块 256 个线程，块的形状不同，合并程度也不同。02 号脚本第 4 部分（64 × 64 的矩阵加法，行 = y、列 = x）：

在”按 32 字节 sector 计数”的模型里，x 方向至少 8 个 float（32 字节）时每个 sector 都被用满：blockDim.x = 16 时一个线程束跨两行、每行读 64 字节（2 个 sector），总数和 32 × 8 一样。x 方向只有 4 时每行只读 16 字节（半个 sector），事务翻倍；只有 1 时每个线程独占一个 sector。实践中 x 方向通常取 16 或 32（16 × 16、32 × 8、32 × 32）：x = 32 时一个线程束只碰一行、地址完全连续，对 128 字节的缓存行和显存的突发传输也最友好。另外，dim3(64, 32) = 2048 个线程会直接发射失败。

3. 张量计算：原课 HW1 的四类内核 · Tensor computation on GPU

3.1 张量 = 存储 + 形状 + 步长

第 2 讲 §1 说过，神经网络最后只剩几种底层运算。原课 HW1 要求在 combine.cu 里写其中四种的 CUDA 内核，并接到 miniTorch 的张量上：

难点不在运算本身，而在下标：miniTorch 的张量和 NumPy、PyTorch 一样，是”一维存储 + 形状 + 步长”（HW1 说明的 “Data layout and strides” 一节，P1 §5）：

转置、切片、广播都只改形状和步长，不搬数据（视图）。所以内核拿到的输入可能是任意步长的，同一个内核必须对连续的张量、转置之后的视图、广播之后的张量都算对。

kernels.py 用模拟器把这四个内核写了一遍（对应 HW1 的 combine.cu），并写了主机端的封装（对应 HW1 的 cuda_kernel_ops.py）。下标换算直接用 P1 的 tensor_layout.py：to_index（第 i 个元素 → 多维下标）、index_to_position（多维下标 → 一维位置）、broadcast_index（输出的下标 → 被广播的输入的下标）。

3.2 map：每个线程一个输出元素

def map_kernel(t, out, out_shape, out_strides, out_size, a, a_shape, a_strides, fn):
    i = t.blockIdx.x * t.blockDim.x + t.threadIdx.x
    if i < out_size:
        out_index = to_index(i, out_shape)                          # 我负责的输出的多维下标
        a_index = broadcast_index(out_index, out_shape, a_shape)    # 对应的输入下标（广播）
        out[index_to_position(out_index, out_strides)] = fn(a[index_to_position(a_index, a_strides)])

思路是 HW1 的提示：“每个线程处理输出的一个元素；用线程和块的编号算出全局编号；做越界检查；考虑基于步长的下标”。03 号脚本第 1 部分用同一个内核处理三种输入（32 × 32）：

内核只看 shape 和 strides，三种都算对了，没有一次额外的拷贝。但视图有代价：输出按行优先连续写，转置视图的输入却是按列读的，相邻线程读的地址相隔 128 字节，每个请求 32 个事务。PyTorch 里对转置后的张量调用 .contiguous()，就是先付一次拷贝（把数据按新的顺序重排），换后续运算的合并访存。什么时候值得？取决于这个张量之后会被读几次。

3.3 zip：两个输入各自广播

zip 的输出形状由广播规则决定（shape_broadcast：从最后一维对齐，每一维要么相等、要么有一个是 1），每个线程把自己的输出下标分别”缩”回 a 和 b 的形状。03 号脚本第 2 部分：

逐元素运算没有求和顺序的问题，所以和 NumPy 逐位相同。广播的那个输入里，同一个元素被很多线程读到，但内存里只有一份——这就是”加偏置”不需要先把偏置复制成整个矩阵的原因。

3.4 reduce：一个线程一个输出，还是一个块一个输出？

归约要把一整条”纤维”（沿 dim 的所有元素）合成一个数。有两种分工：

• 一个线程算一个输出：线程 i 负责第 i 个输出，自己沿 dim 从头加到尾。简单，不需要同步。
• 一个线程块算一个输出（HW1 的提示）：块号 = 第几个输出；块内线程沿 dim 跨步各加一部分，写进共享内存，再做树形归约（§4.2），最后线程 0 写结果。

03 号脚本第 3 部分对形状 (8, 16, 64) 的张量分别沿第 0 维和第 2 维求和：

两点观察：

1. 数值：所有写法都和”按同样顺序计算”的参考实现逐位相同，和 float64 的精确和相差 1e-6 量级。不同的分工就是不同的求和顺序，float32 的舍入就不同（P2 §3）。所以测试归约内核要用容差比较——原课 HW1 的测试也是这样。
2. 合并访存取决于”沿哪一维归约”和”怎样分工”的搭配。沿不连续的第 0 维时，“一个线程 / 输出”让相邻线程算相邻的输出、读相邻的地址，合并得很好；“一个块 / 输出”让块内线程沿第 0 维读，地址相隔 4 KB，每个线程独占一个 sector（这一维只有 8 个元素，32 个线程里还只有 8 个有活干）。沿连续的第 2 维时正好反过来。没有一种分工在所有情况下都最好——PyTorch 的归约内核会根据归约的维度和形状选择不同的策略。

同一个内核换一个 fn 和初始值就是别的归约：求最大值用 $-\infty$（与 NumPy 逐位相同），求积用 1。

3.5 矩阵乘法：一个线程一个输出元素

HW1 的”简单并行”提示：每个线程算输出矩阵的一个元素，读 A 的一行和 B 的一列做内积。kernels.matmul_kernel 用三维网格：blockIdx.z 是第几个批，x 方向对应输出的列、y 方向对应行（§2.3 的规则），批维长度为 1 时步长当作 0（广播）。03 号脚本第 4 部分：

• 和 NumPy 不是逐位相同，差在最后几位：NumPy 调用 BLAS，BLAS 会分块、换求和顺序、用 FMA（一次舍入的乘加）。我们的内核按 k 从 0 到 m−1 依次”先乘、再加”（两次舍入），和同顺序的参考逐位相同。两者都是正确的 float32 结果。（真 GPU 上 nvcc 默认会把 a*b + c 编译成 FMA，所以真 GPU 的结果又会差最后一位。）
• 线程读次数 = 2·n·p·m：每个输出元素读 m 个 A 的元素和 m 个 B 的元素。每做一次乘加（2 FLOPs）就要读两个 float（8 字节）——算术强度 0.25 FLOP/字节。按第 2 讲的平衡点（A100 FP32 约 10），这个内核只能用到 FP32 峰值的约 2.5%。为什么？怎么办？——这正是第 4 讲的开场。

3.6 主机端：把内核接到张量上

HW1 每道题的 Part B 是”集成”：在 cuda_kernel_ops.py 里写 CudaKernelOps.map(fn)、zip(fn)、reduce(fn, start)、matrix_multiply(a, b)，让 miniTorch 的张量运算调用你的内核。kernels.py 里的 tensor_map 等函数做的是同样的事：

tensor_map(dev, fn, a):
    建一个连续的输出张量（形状由输入或广播决定）
    把 a 的存储 cudaMemcpy 到显存，为输出 cudaMalloc
    块数 = ⌈输出元素数 / 每块线程数⌉
    发射 map_kernel，参数是两个设备指针 + 两组 (shape, strides) + fn
    把输出 cudaMemcpy 回来，cudaFree 两块显存

注意它和 §1.4–§1.5 的矛盾：每次调用都拷入拷出，正是”每次都拷来拷去”的写法。原课的 miniTorch 为了教学清晰也是这样做的（所以作业里的 CUDA 后端在小张量上不一定比 CPU 快）；PyTorch 则让张量常驻显存，只在需要时拷回。第 5 讲会把这几个内核接成迷你框架的一个”后端”，看清”前端—自动微分—后端”的分层。

4. 并行归约 · Parallel reduction

4.1 关键路径：并行算法和串行算法的根本区别

把 16384 个数加成一个数，04 号脚本第 1 部分比较了三种分工：

三种写法读的数一样多（16,384 次），总的加法次数也差不多，差别全在关键路径：每个线程先依次累加自己那份（n / 线程数步），树形归约每轮一步（log₂256 = 8 轮），打在同一个地址上的原子操作只能排队（16 步）。只要线程足够多，GPU 的时间由关键路径决定——这是写并行算法时首先要问的问题。

树形归约的思想来自结合律：$(a+b)+(c+d)$ 和 $((a+b)+c)+d$ 数学上相等，但前者的两次加法可以同时做。n 个数的树有 $\log_2 n$ 层，每层内的加法互不依赖。顺便，它在数值上也更好：和 P2 的成对求和一样，树形求和的误差（5.16e-06）比从头加到尾（3.76e-04）小两个数量级。

4.2 树形归约的三种写法

块内树形归约有几种经典写法，原课 HW1 的可选提示给出的是第一种（这几种写法出自 Mark Harris 的教程 Optimizing Parallel Reduction in CUDA）：

// interleaved（HW1 的 reduce0）：第 s 轮，tid 是 2s 的倍数的线程干活
for (s = 1; s < blockDim.x; s *= 2) { if (tid % (2*s) == 0) sdata[tid] += sdata[tid + s]; __syncthreads(); }
// strided：第 s 轮，前 blockDim/(2s) 个线程干活，各算 index = 2·s·tid
for (s = 1; s < blockDim.x; s *= 2) { int index = 2*s*tid; if (index < blockDim.x) sdata[index] += sdata[index + s]; __syncthreads(); }
// sequential：步长从 blockDim/2 开始减半，前 s 个线程干活
for (s = blockDim.x/2; s > 0; s >>= 1) { if (tid < s) sdata[tid] += sdata[tid + s]; __syncthreads(); }

interleaved 和 strided 合并的是完全相同的配对（第 s 轮把位置 2s·k 和 2s·k + s 合并），只是换了”由哪个线程来做”；sequential 换了一种配对（位置 k 和 k + s）。04 号脚本第 2 部分（8 行 × 1024 列，每行一个 256 线程的块，开了指令计数）：

• interleaved：每一轮都让每个线程束里的一部分线程干活（s = 1 时每个线程束里偶数号的 16 个线程）。这是线程束分化（第 2 讲 §4.5），SIMT 效率 87%；而且每个线程束每一轮都要发共享内存请求（1200 个，另两种只有 360 个）。
• strided：让编号连续的线程干活，分化少了，请求也降到 360 个。但线程 tid 访问 sdata[2s·tid]：s = 1 时一个线程束访问第 0, 2, 4, …, 62 号字，两两落在同一个 bank 上（共享内存有 32 个 bank，第 k 个字在第 k mod 32 号 bank；同一 bank 的不同字要排队），s 越大冲突越严重。平均每请求 3.33 轮，总轮数还是 1200——省下来的全被 bank 冲突吃回去了；它的发射次数也没比 interleaved 少（多算了一行 index）。
• sequential：干活的线程编号连续（分化只在最后几轮 s < 32 时出现），访问的地址也连续（没有 bank 冲突）。发射次数最少，共享内存总轮数只有另两种的 30%。

（这里的计数包含了每个线程先累加 4 个数的阶段，三种写法完全相同，所以整体的发射次数差别看起来不大。）这个例子把第 2 讲的分化和第 4 讲要正式讲的 bank 冲突放到了一起：同样的计算，“让哪个线程做哪件事”决定了效率。

4.3 原子操作、竞争与不确定性

多级归约的最后一步，是把每个块的部分和加到同一个位置上。04 号脚本第 3 部分（2048 个数）：

第一行最有教育意义：不用原子操作的写法在模拟器上”算对了”——因为模拟器让线程一个接一个地执行。真 GPU 上，两个线程可能同时读到旧值、各自加完再写回，其中一个的加法就丢了（P4 §3 的”读-加-写”竞争）。这正是”测试通过不代表没有竞争”的典型例子，模拟器的竞争检测器把它抓了出来。

atomicAdd 把”读-加-写”变成硬件保证不可分割的一步，结果正确；代价是打在同一个地址上的原子操作只能排队。先在块内归约、每块只做一次原子操作，原子操作从 2048 次降到 8 次。

但原子操作只保证”不丢”，不保证顺序：真 GPU 上每次运行，各个块完成的先后不同，加法的顺序就不同。float32 不满足结合律，所以结果每次可能差最后几位。脚本对照：同样 2048 个 float32 按 20 种随机顺序依次相加，得到 15 个不同的结果（范围 −95.087830 到 −95.087662）。这是 GPU 训练”不能逐位复现”的常见来源之一。需要确定性时，要用固定顺序的归约（PyTorch 的 torch.use_deterministic_algorithms(True) 会换用确定性的内核，通常更慢）。

5. 小结 · Summary

1. CUDA 程序的五步：cudaMalloc（传 &dA，因为要修改指针本身）→ cudaMemcpy 拷入（先目的、后来源，字节数）→ 发射内核 → cudaMemcpy 拷出 → cudaFree。
2. PCIe 比显存慢几十倍：只用一次的数据不值得搬到 GPU 上；正确的用法是数据常驻显存、被一串内核反复使用（再进一步是算子融合）。
3. 三种内存：寄存器（线程私有）、共享内存（块内共享、块结束即消失）、全局内存（所有线程可见、跨内核发射一直存在）。__global__ 由主机调用在设备执行，__device__ 由设备调用在设备执行。
4. 多维编号就是把 (blockIdx, threadIdx) 按行优先展平，x 变化最快；线程束按块内编号切分。行优先矩阵要让 threadIdx.x 对应列号，否则事务数成倍增加；二维也要对行、列分别做越界检查。
5. HW1 的四类内核：每个线程（或每个块）负责一个输出位置，用 to_index → broadcast_index → index_to_position 找到输入。同一个内核对连续张量、转置视图、广播都正确，但视图会改变访存模式。
6. 归约：一个线程 / 输出 vs 一个块 / 输出，哪个合并取决于沿哪一维归约。浮点归约的结果依赖顺序，测试要用容差或同顺序参考。
7. 朴素矩阵乘法的算术强度只有 0.25 FLOP/字节——第 4 讲的问题。
8. 并行归约看关键路径：树形归约把 n 步变成 log₂n 层；sequential 写法同时避开了分化和 bank 冲突；跨块汇总用原子操作，它保证正确但不保证顺序。

6. 常见误区 · Common pitfalls

• “cudaMalloc(dA, size) 就行。” 必须传 &dA：函数要把地址写回你的变量。传 dA 只是把一个未初始化的值复制进去。
• “cudaMemcpy(src, dst, …)。” 顺序是先目的、后来源，和 memcpy 一样；count 是字节数，不是元素个数。
• “把计算搬到 GPU 上总会更快。” 对只用一次的数据，PCIe 拷贝比计算贵得多（向量加法拷贝占 96%）。
• “i = blockIdx.x * blockDim.x + threadIdx.x 当行号很自然。” 行优先矩阵里 x 方向应该对应列号，否则一个线程束访问 32 个不同的行。
• “例子里 N = 32 没有越界检查，我的代码也不用。” slides 的例子正好整除；换一个 N 就越界。
• “转置视图是免费的。” 不拷贝数据是免费的，但之后的读取变成了按列读，事务可能多 8 倍。
• “归约结果和 NumPy 不完全一样，所以我的内核错了。” 浮点加法不满足结合律，换一种分工就换了一种顺序；要和同顺序参考比，或用容差。
• “树形归约里写 if (tid < s) { …; __syncthreads(); }。” 栅栏不能放在只有部分线程执行的分支里（第 2 讲 BarrierError）；if 只包住”加”那一行。
• “不用原子操作也测对了。” 在模拟器上或碰巧的时序下”算对”不代表没有竞争。多线程往同一个位置写，要么用原子操作，要么先归约。
• “用了原子操作，结果就可复现。” 原子操作不保证顺序，浮点结果每次可能差最后几位。

7. 练习 · Exercises

1. 写出 cudaMalloc 的原型，解释为什么它的第一个参数是 void**，以及调用时为什么写 &dA。
2. 在 A6000 上（PCIe 32 GB/s、显存 768 GB/s，每次拷贝 + 10 μs、每次发射 + 5 μs）对 10⁷ 个 float32 做向量加法，拷贝和内核各要多久？把答案和 01 号脚本的表对照。
3. dim3 Dg(4, 2, 1); dim3 Db(8, 8, 1); 一共多少个块、多少个线程？每个块有几个线程束？
4. 在 slides 的 fullKernel（网格 (2,3,4)、块 (2,4,8)）里，blockIdx = (1, 2, 3)、threadIdx = (1, 3, 7) 的线程算出的 tid 是多少？
5. 1000 × 1000 的矩阵加法用 16 × 16 的块，网格多大？一共发射多少个线程？多少个线程什么也不做？
6. slides 的 MatAddKernel（x 方向当行号）用在 1024 × 1024 的矩阵上、块 32 × 32。一个线程束的一条读指令产生几个 32 字节事务？改成 x 方向当列号呢？
7. 一个形状 (2, 3, 4) 的连续张量做 permute(2, 0, 1)，新的形状和步长是什么？新张量的元素 [3, 1, 2] 在存储里的位置是多少？
8. a 的形状 (3, 1)、b 的形状 (1, 4)，zip(add) 的输出形状是什么？输出下标 (2, 3) 对应 a 和 b 的哪个下标？
9. 对形状 (8, 16, 64) 的张量沿第 2 维求和，用”一个块 / 输出、每块 32 个线程、sequential 树形归约”。一共几个块？每个线程先累加几个元素？树形归约几轮？每个块几次栅栏？全部块一共几次？和 03 号脚本的表对照。
10. 朴素矩阵乘法 C (n×p) = A (n×m) · B (m×p)，每个线程算一个输出。一共读多少次全局内存？做多少 FLOPs？算术强度（按线程的读取字节算）是多少？
11. interleaved 树形归约，块大小 256。第 s = 1 轮每个线程束有几个线程干活？s = 16 轮呢？s = 32 轮有几个线程束里有线程干活？
12. 2048 个线程各对同一个地址做一次 atomicAdd，和”8 个块各自树形归约、再各做一次 atomicAdd”相比，原子操作那一段的关键路径各是多少步？为什么后者在真 GPU 上快得多？

参考答案

1. cudaError_t cudaMalloc(void** devPtr, size_t size)。返回值被错误码占用了，新分配的显存地址只能通过参数”带回来”；C 是按值传参，要修改调用者的指针变量 dA，就要传它的地址 &dA（类型 float**，转成通用的 void**）。
2. 拷贝：3 × (10 μs + 4×10⁷ 字节 / 32×10⁹) = 3 × (0.01 + 1.25) ms = 3.78 ms；内核：5 μs + 12×10⁷ / 768×10⁹ ≈ 0.005 + 0.156 = 0.161 ms。拷贝占 95.9%，与脚本一致。
3. 块：4 × 2 × 1 = 8；每块 8 × 8 = 64 个线程，共 512 个线程；每块 64 / 32 = 2 个线程束。
4. block_id = 1 + 2 × 2 + 3 × 2 × 3 = 23；thread_offset = 1 + 3 × 2 + 7 × 2 × 4 = 63；tid = 23 × 64 + 63 = 1535（最后一个线程）。
5. ⌈1000/16⌉ = 63，网格 63 × 63；线程 63² × 256 = 1,016,064；空转 1,016,064 − 1,000,000 = 16,064。
6. 一个线程束是 threadIdx.x = 0..31（同一个 threadIdx.y）：x 当行号时 32 个线程访问 32 个不同的行，地址相隔 4096 字节，32 个事务；x 当列号时访问同一行的 32 个相邻 float（128 字节），4 个事务。
7. 原步长 (12, 4, 1)；新形状 (4, 2, 3)、新步长 (1, 12, 4)；位置 = 3 × 1 + 1 × 12 + 2 × 4 = 23。
8. 输出 (3, 4)；a 的下标 (2, 0)（第 1 维大小为 1，取 0），b 的下标 (0, 3)。
9. 输出 8 × 16 × 1 = 128 个块；每个线程累加 64 / 32 = 2 个元素；log₂32 = 5 轮；每块栅栏 1（写共享内存之后）+ 5 = 6 次；共 128 × 6 = 768 次，与脚本表中的 768 一致。
10. 读 2·n·p·m 次（每个输出 m 次 A、m 次 B）；FLOPs 2·n·p·m；算术强度 = 2npm / (8npm) = 0.25 FLOP/字节。
11. s = 1：tid 是 2 的倍数的线程，每个线程束 16 个；s = 16：tid 是 32 的倍数，每个线程束 1 个（SIMT 效率 1/32）；s = 32：tid 是 64 的倍数，8 个线程束里只有 4 个有线程干活（线程束 0、2、4、6），每个只有 1 个线程。
12. 前者 2048 步（同一地址的原子操作串行）；后者树形归约 8 轮 + 8 步原子操作。真 GPU 上同一地址的原子操作要在 L2 里排队处理，2048 次的排队远比 8 次长；而且块内归约用的是快得多的共享内存。

8. 延伸阅读 · Further reading

• PMPP 第 4 版第 2 章（向量加法、cudaMalloc / cudaMemcpy、错误检查）、第 3 章（多维网格、彩色转灰度与图像模糊、矩阵乘法）、第 10 章（归约：树形归约、分化、合并、线程粗化、原子操作）。
• Mark Harris, Optimizing Parallel Reduction in CUDA（NVIDIA 教程幻灯片）：从 interleaved 到 sequential，再到循环展开、线程粗化，一步步把归约做到接近带宽上限。§4.2 的三种写法出自这里。
• 原课 HW1 说明（https://llmsystem.github.io/llmsystemhomework/ ）与 CUDA_Code_Examples.ipynb（原课 llmsys_code_examples 仓库的 simple_cuda_demo，slides 推荐在 Colab 上运行）：本讲的四类内核在真 GPU 上的版本。
• NVIDIA CUDA C++ Programming Guide：“Programming Interface → CUDA Runtime”（设备内存、异步并发执行、流）；“Asynchronous Data Copies”与 Hopper 调优指南里的 TMA。
• PyTorch 文档 Reproducibility：哪些 CUDA 运算是非确定性的，以及 use_deterministic_algorithms 的代价。

9. 下一讲预告 · Next lecture

第 3 讲留下了一个数字：朴素矩阵乘法的算术强度只有 0.25 FLOP/字节。原课第 4 讲从这里出发：在 A100 上，这个内核受显存带宽限制，只能用到 FP32 峰值的 2.48%。怎样把它提高几十倍？答案是分块（tiling）：让一个线程块把 A 和 B 的一小块搬进共享内存，块内的线程反复使用它们。之后是内存并行：合并访存、转置的三个版本、共享内存的 bank 冲突和”+1 填充”的技巧；再之后是稀疏矩阵（CSR 格式与 SpMV 内核）和 cuBLAS（为什么通常不该自己写 GEMM，以及它”列优先”的接口怎么用）。模拟器会把每一步省下的访存都数出来。