第 4 讲 · GPU 加速：分块、内存并行、稀疏矩阵与 cuBLAS · GPU Acceleration

对应 Maps to：CMU 11-868 / 11-968 Large Language Model Systems（Spring 2026，Lei Li），Lecture 4 · Mon 1/26 · GPU Acceleration。 阅读材料：PMPP 第 4 版第 5 章（存储架构与数据局部性：分块）、第 6 章（性能考虑：合并访存、线程粗化）、第 14 章（稀疏矩阵计算）；NVIDIA CUDA Programming Guide；cuBLAS 文档。 依据说明：本讲以原课第 4 讲公开的 slides（llmsys-04-gpu-acceleration.pdf，50 页）为主线：计算–全局内存访问比与 A100 上 2.48% 的计算、存储层次的延迟、CUDA 的内存模型、分块矩阵乘法的步骤与 MatMulTiledKernel、内存的限制、DRAM 突发与合并访存、矩阵转置的三个版本与 bank 冲突、CSR 与 SpMV 内核、cuBLAS 的接口。模拟器上的全部计数、TILE 的扫描、漏掉栅栏 / 边界处理的反例、线程粗化、“一个线程束一行”的 SpMV、2:4 稀疏、行优先交换技巧和 roofline 阶梯都是本讲义自己设计的，不代表课堂上讲了什么。 配套代码 Code：code/lecture04/（5 个脚本 + accel.py；复用第 2 讲的 simt.py、gpu_model.py；每个脚本 6 秒以内，输出重复运行逐位相同） 预计阅读时间：4–5 小时（含运行代码）

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

0.1 学习目标 Learning goals

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

1. 算出一个内核的计算–全局内存访问比，并用它和 GPU 的平衡点解释”朴素矩阵乘法为什么只能用到 A100 峰值的 2.48%”。
2. 说出 CUDA 的内存模型（寄存器、本地、共享、全局、常量）各自的延迟、作用域、生命周期，以及变量声明怎样决定它放在哪里。
3. 写出分块矩阵乘法，解释它为什么把全局内存读取减少到 1/TILE、两个 __syncthreads() 各防什么、边界怎样处理、TILE 受哪些资源限制，以及线程粗化再省在哪里。
4. 用”一个线程束在同一条指令里访问了哪些地址”判断合并访存；理解矩阵转置的三个版本，算出共享内存的 bank 冲突，并解释”+1 填充”为什么有效。
5. 用 CSR 表示稀疏矩阵，写出 SpMV 内核，判断多稀疏才划算，并理解非结构化稀疏和 2:4 结构化稀疏的区别。
6. 正确调用 cuBLAS：列优先、leading dimension、transa / transb，以及让行优先张量直接调用 cuBLAS 的”交换技巧”；知道为什么几乎总是用库而不是自己写 GEMM。

0.2 先修与本讲范围 Prerequisites & scope

0.3 本讲地图 Map

§1  访存效率决定一切：计算–全局访存比、2.48%、延迟、内存模型、重复读取    ← 代码 01
§2  分块：共享内存里的复用；两个栅栏；边界；资源限制；线程粗化            ← 代码 02
§3  内存并行：DRAM 突发与合并访存；转置的三个版本；bank 冲突与 +1         ← 代码 03
§4  稀疏矩阵：CSR、SpMV、多稀疏才划算、2:4 结构化稀疏                     ← 代码 04
§5  cuBLAS：接口、列优先、leading dimension、行优先技巧、roofline 阶梯    ← 代码 05
§6  小结   §7  常见误区   §8  练习   §9  延伸阅读   §10 下一讲预告

slides 把本讲概括为四种 GPU 加速技术：分块（PMPP 第 5 章）、内存并行（第 5、6 章）、稀疏矩阵乘法、cuBLAS。前三种有一个共同点：计算本身一点没变，变的只是”哪个线程在什么时候从哪里取数”。

1. 访存效率决定一切 · Memory access efficiency is critical

1.1 计算–全局内存访问比

第 3 讲的朴素矩阵乘法内核，每个线程算 C 的一个元素：

__global__ void MatMulKernel(float *a, 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) return;
    float Pvalue = 0.0;
    for (int k = 0; k < N; k++)
        Pvalue += a[row * N + k] * b[k * N + col];     // 1 次乘 + 1 次加；读全局内存 2 次（各 4 字节）
    c[row * N + col] = Pvalue;
}

slides 定义了一个量：计算–全局内存访问比（compute-to-global-memory-access ratio）= 每从全局内存读一个字节，做多少次浮点运算。循环体做 2 次运算、读 8 字节，比值是 0.25 FLOP/字节。它就是第 1 讲的”算术强度”，只是专门按全局内存的访问来数。

1.2 A100 上的 2.48%

slides 用 A100 80GB PCIe 的规格算了一笔账（01 号脚本第 1 部分复现）：

显存每秒只能送来 1,935 GB 数据，每字节只能支撑 0.25 次运算，所以无论计算单元多快，这个内核最多跑 483.75 GFLOPS——访存受限（memory-bound）。要让 FP32 核心跑满，比值要提高到平衡点 10.1；要让 TF32 张量核心跑满，要 81。换句话说：每从显存读一个字节，要设法让它参与几十次运算。这一讲的前三节都是在想办法做到这一点。

1.3 读数据比计算慢得多

slides 用一条加法说明为什么”读”这么贵。C[i] = A[i] + B[i] 编译成 PTX（GPU 的中间汇编）是四条指令：

slides 给出的 CUDA 存储层次延迟（量级）：寄存器约 1 个周期；共享内存、L1 缓存、常量缓存约 5 个周期；全局内存约 500 个周期。共享内存比全局内存快约 100 倍，而且就在 SM 上。第 2 讲 §4.7 讲过 GPU 用大量线程束掩盖这 500 个周期；但掩盖得再好，带宽也是硬上限（§1.2）。真正的出路是少读全局内存。

1.4 CUDA 的内存模型

slides 第 11 页给出变量声明和内存的对应（01 号脚本第 3 部分）：

两个容易踩的坑：局部数组通常放不进寄存器，会被放到”本地内存”——名字叫本地，其实在显存里，和全局内存一样慢；共享内存是程序员手动管理的缓存——搬什么、什么时候搬、什么时候可以覆盖，都要自己决定（§2.5 的两个栅栏）。

1.5 机会在哪里：重复读取

slides 第 12–13 页指出了优化的机会：线程块里的线程会用到相同的数据。C 的同一行的所有元素都要 A 的同一行，同一列的都要 B 的同一列。01 号脚本第 4 部分在模拟器上把朴素矩阵乘法（N = 32）的字节数按三种口径数了一遍：

A 的每个元素被读了 N = 32 次。按 slides 的口径（每次读都去显存），比值是 0.25；如果每个数只从显存读一次，比值是 N/6 = 5.33，而且 N 越大越高。“要求读的”和”数据本身”差了 22 倍——这就是重复读取，也是优化的空间。真实 GPU 的 L1 / L2 缓存会自动接住一部分重复读取，但容量有限、命中与否不受控制；分块用共享内存把”重复使用”变成程序里明确的、有保证的事情。

2. 分块 · Tiling

2.1 直觉：把要反复用的书搬到书桌上

想象你要写一份报告，需要反复查阅图书馆书库里的几十本书。每查一个数就跑一趟书库（全局内存，500 个周期），太慢了。更好的办法：先把接下来一段时间要用的几本书一次搬到书桌上（共享内存，5 个周期），查完这几本再换下一批。书桌只有这么大（共享内存容量有限），所以要分批（分块）；几个人共用一张书桌（线程块共享），换书之前要确认大家都查完了（栅栏）。

2.2 分块矩阵乘法的步骤

slides 第 14–18 页一步步画出了这个过程。C 被分成 TILE × TILE 的小块，每个线程块负责一块，每个线程负责其中一个元素。要算这一块，需要 A 的一条”横带”（TILE 行）和 B 的一条”竖带”（TILE 列）；把它们沿 k 方向切成 TILE × TILE 的小块，一相（phase）处理一对：

1. （同时） 把 A、B 的第一对小块搬进共享内存：块内每个线程搬一个元素；__syncthreads() 等大家搬完。
2. 每个线程用共享内存里的小块算出部分和（As 的一行 · Bs 的一列，TILE 次乘加）；__syncthreads() 等大家算完。
3. （同时） 搬下一对小块……
4. 继续累加部分和……直到走完整条横带和竖带。

slides 第 36 页的内核（原课的现场编程环节）：

#define TILE_WIDTH 2
__global__ void MatMulTiledKernel(float* d_A, float* d_B, float* d_C, int N) {
    __shared__ float As[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Bs[TILE_WIDTH][TILE_WIDTH];
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    float Cvalue = 0;
    for (int ph = 0; ph < N / TILE_WIDTH; ++ph) {
        As[threadIdx.y][threadIdx.x] = d_A[row * N + ph * TILE_WIDTH + threadIdx.x];
        Bs[threadIdx.y][threadIdx.x] = d_B[(ph * TILE_WIDTH + threadIdx.y) * N + col];
        __syncthreads();
        for (int k = 0; k < TILE_WIDTH; ++k)
            Cvalue += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        __syncthreads();
    }
    d_C[row * N + col] = Cvalue;
}

为什么能省？块里同一行的 TILE 个线程都要 A 的同一个元素——现在由一个线程把它搬进共享内存，TILE 个线程从共享内存里读。每个从全局内存读进来的数被用了 TILE 次，全局内存读取减少到 1/TILE，计算–全局访存比提高到 0.25 × TILE。

2.3 实验：TILE 从 1 到 32

02 号脚本第 1 部分（N = 32，accel.matmul_tiled）：

• 线程读全局内存的次数正好是 2N³ / TILE。运算总量没变，读的地方变了：TILE = 16 时共享内存的读取 65,536 次（和朴素版的全局读取一样多），外加搬进来时写共享内存的 4,096 次。
• 结果和朴素版逐位相同：每个线程仍然按 k = 0, 1, …, N−1 的顺序累加。分块是纯粹的性能优化，不改变数值。
• 一个反常的数：TILE = 2 时线程读全局内存的次数减半，全局读事务却是朴素版的 4 倍。2 × 2 的块只有 4 个线程，凑不满一个线程束，每条读指令只取 2 个数，合并访存完全失效。块不能太小——一行至少要有 8 个线程（第 3 讲 §2.5）。
• 按 slides 的口径，TILE = 16 时能用到 A100 FP32 峰值的约 40%，TILE = 32 时约 80%。要超过平衡点还需要更多手段（§2.8、§5.5）。

2.4 两个栅栏各防什么

02 号脚本第 2 部分故意去掉其中一个栅栏（N = 32，TILE = 8）：

• 第一个栅栏（搬完之后）保证”读之前大家都写完了”（read-after-write）。漏掉它，线程会读到同伴还没搬进来的位置。模拟器的共享内存初始是 NaN，所以结果里出现了 NaN；真 GPU 上读到的是上一相留下的旧值，错得更隐蔽。
• 第二个栅栏（算完之后）保证”写之前大家都读完了”（write-after-read）。漏掉它，跑得快的线程会进入下一相、覆盖 As / Bs，而别的线程还在用旧的那一块。结果是错的数，而不是 NaN。
• 1024 个元素里有 16 个”碰巧对了”：模拟器按编号依次执行线程，每个块（共 16 个）里编号最大的那个线程运行时，它要的数都已经被前面的线程搬好了。真 GPU 上的执行顺序不确定，同样的 bug 可能时有时无——这正是并发 bug 难调的原因，也是模拟器另外做竞争检测的原因。

2.5 边界：N 不是 TILE 的倍数

slides 的内核用 ph < N / TILE_WIDTH（整数除法），而且不检查越界——slides 的例子 N 正好是 TILE 的倍数。02 号脚本第 3 部分取 N = 10、TILE = 4：

修正的要点：相数向上取整；搬运时如果行或列超出 N，就往共享内存写 0（0 乘任何数都是 0）；最后只让矩阵内的线程写 C。越界的线程不能提前 return——它们还要参与搬运和栅栏（第 2 讲的 BarrierError：__syncthreads() 必须被块内所有线程执行到）。

2.6 TILE 能有多大：资源限制

slides 第 20–21 页列出了限制：H100 每个 SM 有 256 KB 的 L1 / 共享内存、64K 个寄存器；A100 是 192 KB 和 64K（slides 的表里 A100 写的是 104 个 SM，NVIDIA 白皮书给出的是 108 个）。分块用的共享内存 As + Bs = 2 × TILE² × 4 字节：TILE = 32 时 As、Bs 各 32 × 32 × 4 = 4 KB（slides 的例子），合计 8 KB，离每块的上限（A100 163 KB、H100 227 KB）还很远。02 号脚本第 4 部分：

真正先撞上的是每块最多 1024 个线程：只要还是”一个线程一个输出”，TILE 就不能超过 32，比值就停在 8。要让每个块覆盖更大的输出块，必须让每个线程算好几个输出，用寄存器再做一层复用（§2.7、§5.5）。

关于 slides 第 21 页的寄存器算术。 slides 写”1024 个线程时，每个线程只能用 16384 / 1024 = 16 个寄存器”。16384 是一个分区的寄存器数（64 KB）；整个 SM 有 4 个分区、65,536 个寄存器，1024 个线程的块分布在 4 个分区上，所以每个线程最多能用 64 个；SM 驻留满 2048 个线程时每个 32 个（第 2 讲 §3.5）。结论的方向不变：线程越多，每个线程能用的寄存器越少。

2.7 线程粗化：一个块算好几个输出块

在分块版本里，横向相邻的两个块要的是同一条 A 的横带，却各自从全局内存读一遍。线程粗化（thread coarsening，PMPP 6.3 节）让一个块负责横向相邻的 COARSE 个输出块、每个线程算 COARSE 个输出：A 的小块每相只读一次，被 COARSE 个输出块共用。02 号脚本第 5 部分（N = 32，TILE = 8）：

A 的读取降到 1/COARSE，结果仍与朴素版逐位相同。代价是每个线程要 COARSE 个累加器（更多寄存器 → 占用率可能下降），块数减少（并行度降低；矩阵小时 GPU 可能喂不饱）。粗化是”用更少的线程做更多的事”，在并行度过剩时才划算。把它推到底——块在共享内存里分块、线程在寄存器里再分块（每个线程算 8 × 8 之类的小块）——就是 cuBLAS 这类库的做法（§5.5）。

3. 内存并行：合并访存与 bank 冲突 · Memory parallelism

3.1 DRAM 突发与合并访存

分块解决的是”读多少次”，这一节解决”怎么读”。slides 第 22 页：显存每次至少传 32 字节；一个线程束里连续的访问会被合并（coalesced）在一起；C 和 CUDA 的多维数组按行优先存放，这让 DRAM 的突发传输（burst）成为可能，比一个个单独访问快得多。

03 号脚本第 1 部分让一个线程束的 32 个线程按不同的步长和起点各读一个 float：

• 步长 1、起点对齐：32 个 float = 128 字节 = 4 个 sector，一个字节都不浪费（slides 第 23 页”一次读取服务一个线程束的多个线程”）。
• 起点偏移 1 个元素就跨了 5 个 sector；偏移 8 个元素（32 字节）又对齐了。cudaMalloc 返回的地址至少 256 字节对齐，所以”从数组开头按线程编号读”总是对齐的。
• 步长越大越浪费（slides 第 24 页”相隔 32 字节，两个线程要两次读取”）；步长 ≥ 8 后每个线程独占一个 sector。所有线程读同一个数是广播，1 个事务。

3.2 矩阵转置：读和写总有一边是按列走的

slides 用矩阵转置作为例子。朴素版本（第 26 页）：

#define INDX(row, col, ld) (((row) * (ld)) + (col))       // 行优先的二维下标
__global__ void naive_cuda_transpose(int m, float *a, float *c) {
    int myCol = blockDim.x * blockIdx.x + threadIdx.x;
    int myRow = blockDim.y * blockIdx.y + threadIdx.y;
    if (myRow < m && myCol < m)
        c[INDX(myCol, myRow, m)] = a[INDX(myRow, myCol, m)];
}

slides 第 27–28 页问：原数组里的位置是多少？转置后的位置是多少？访问合并吗？答案：读 a[myRow * m + myCol] 时，一个线程束（threadIdx.x = 0..31）读同一行的 32 个相邻元素——合并；写 c[myCol * m + myRow] 时，myCol 随 threadIdx.x 变化，相邻线程写同一列、相隔 m 个元素——不合并。转置的本质矛盾：读和写，总有一边是按列走的。

3.3 用共享内存把”按列走”挪进片上

slides 第 29–30 页的关键思路：先把一个小块按行（合并地）读进共享内存，__syncthreads()，再按行（合并地）写出它的转置——“按列走”的那一步发生在共享内存里，而共享内存不需要合并。

__global__ void smem_cuda_transpose(int m, float *a, float *c) {
    __shared__ float smemArray[THREADS_PER_BLOCK_X][THREADS_PER_BLOCK_Y];
    const int tileCol = blockDim.x * blockIdx.x;
    const int tileRow = blockDim.y * blockIdx.y;
    smemArray[threadIdx.x][threadIdx.y] = a[INDX(tileRow + threadIdx.y, tileCol + threadIdx.x, m)];
    __syncthreads();
    c[INDX(tileCol + threadIdx.y, tileRow + threadIdx.x, m)] = smemArray[threadIdx.y][threadIdx.x];
}

03 号脚本第 2 部分（64 × 64 矩阵，块 32 × 32）：

全局事务从 4,608 降到 1,024——读写都合并了。但第二行出现了新问题：共享内存 256 个请求却要 4,224 轮。

3.4 bank 冲突与 “+1”

slides 第 31 页：共享内存分成 32 个 bank，每个 bank 宽 4 字节；访问同一个 bank 的不同元素要依次进行。连续的 4 字节字依次落在 bank 0、1、……、31、0、1……，即第 k 个字在第 k mod 32 号 bank。一个线程束的 32 个线程如果访问 32 个不同 bank 的字，一轮完成；落在同一个 bank 的不同字要排队（读同一个字是广播，不排队）。03 号脚本第 3 部分把 bank 编号算了出来：

写 smemArray[threadIdx.x][threadIdx.y] 时，一个线程束写的是同一列的 32 个位置：线程 tx 访问第 tx × 32 号字，全都是 bank 0——32 路冲突，要排 32 轮。slides 第 33 页的修正只改了一个数：把数组声明成 [THREADS_PER_BLOCK_X][THREADS_PER_BLOCK_Y + 1]，多开一列不用。现在线程 tx 访问第 tx × 33 号字，bank = tx × 33 mod 32 = tx，32 个线程正好错开到 32 个 bank。多花 128 字节的共享内存，换掉了 32 路冲突（slides 第 34 页画出了加一列之后的新排布）。

规律和全局内存的合并访存是对应的：线程束里相邻的线程访问相邻的字，两种内存都最高兴。不同的是，全局内存怕”跨太多 sector”，共享内存怕”挤在同一个 bank”。

3.5 放到真实尺寸上

03 号脚本第 4 部分把模拟出的事务数等比放大到 8192 × 8192，按 A100 的 1,935 GB/s 估算（不计缓存、不计 bank 冲突）：

转置没有任何运算，数据量是读写各一遍 537 MB；naive 版本要搬 4.5 倍的字节，有效带宽只有峰值的 22%。两个共享内存版本在这个纸面模型里一样快——差别在 bank 冲突，模型没有计入；真实 GPU 上带冲突的版本会因为 32 路排队明显变慢。真实的 L2 缓存会接住 naive 版本的一部分写，所以实测差距通常没有 4.5 倍那么大，但结论不变：纯搬运的内核，性能完全取决于访存模式。

4. 稀疏矩阵 · Sparse matrix multiplication

4.1 为什么关心稀疏

很多矩阵大部分元素是 0：图的邻接矩阵、推荐系统的用户–物品矩阵、剪枝后的神经网络权重、MoE 里”每个 token 只用几个专家”（第 17 讲）。把 0 也存下来、也拿来乘，既浪费内存又浪费带宽。稀疏格式只存非零元。

4.2 CSR：压缩稀疏行

slides 第 41 页的 CSR（Compressed Sparse Row）用三个数组。04 号脚本第 1 部分：

M = [[3, 0, 1, 0],          data      = [3, 1, 2, 4, 1, 1, 1]   非零元的值，按行优先顺序
     [0, 0, 0, 0],          col_index = [0, 2, 1, 2, 3, 0, 3]   每个非零元在第几列
     [0, 2, 4, 1],          row_ptr   = [0, 2, 2, 5, 7]         第 r 行的非零元是 data[row_ptr[r] : row_ptr[r+1]]
     [1, 0, 0, 1]]

row_ptr 比行数多一个元素。第 1 行是空行，row_ptr[1] = row_ptr[2] = 2，内层循环一次也不执行——空行不需要特殊处理。

4.3 SpMV：串行与 CUDA

slides 第 42–43 页的稀疏矩阵-向量乘法 y = M · x：

// 串行
for (int row = 0; row < n; row++) {
    float dot = 0;
    for (int el = row_ptr[row]; el < row_ptr[row + 1]; el++)
        dot += x[col_index[el]] * data[el];
    y[row] += dot;
}
// CUDA：一个线程负责一行
__global__ void SpMVCSRKernel(float *data, int *col_index, int *row_ptr, float *x, float *y, int num_rows) {
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < num_rows) { /* 同上的内层循环 */ }
}

x = [1, 2, 3, 4] 时，串行循环、模拟器上的内核、稠密乘法都得到 [6, 0, 20, 5]。注意 x[col_index[el]]：x 是按 col_index 间接取的（gather）。这是稀疏计算和稠密计算最大的不同——访问哪里由数据决定，编译器和硬件都没法提前知道，合并访存和缓存都难以保证。

4.4 多稀疏才划算

04 号脚本第 2 部分（n = 4096，A100 纸面估算）：

CSR 每个非零元要存值和列号两样（各 4 字节），所以非零元比例超过 50% 时比稠密存储还大。SpMV 和稠密 matvec 都是访存受限的（每个矩阵元素只用一次），时间约等于搬运的字节 ÷ 带宽：x 能全部留在缓存里时，分界点在 50% 左右；gather 次次落空时，分界点降到 10% 左右。稀疏要足够稀疏才值得：图和推荐矩阵常常在 0.1% 以下，而剪枝到 50% 的神经网络权重用 CSR 反而更慢。

4.5 一个线程一行，还是一个线程束一行？

slides 的内核”一个线程一行”有两个问题：相邻线程同时读各自那一行的第 k 个非零元，位置相隔一整行——不合并；各行长度不同时，线程束要等它里面最长的那一行——分化（第 2 讲 §4.5）。PMPP 第 14 章的改进是”一个线程束一行”：32 个线程读同一行连续的 32 个非零元，最后在共享内存里树形归约。04 号脚本第 3 部分（256 行 × 1024 列）：

结果不是”一个线程束一行更好”那么简单：

• 访存：一个线程束一行读矩阵本身的事务少 8 倍（每行 64 个时 32,768 vs 4,096）；算上 x 的 gather（两种写法都是随机访问），全部读事务少 2.3 倍。
• 指令：一个线程束一行的发射次数多一个数量级——每一行都要一个线程束、都要做一次 5 轮的树形归约；行短时大部分车道闲着（每行 8 个只用 8 个车道读数）。
• 分化：行长不均匀时，一个线程一行的 SIMT 效率只有 10%。

SpMV 是访存受限的，真 GPU 上通常是访存那一项说了算，所以长行一般用一个线程束一行；行很短时，事务只省了 1.8 倍却要多付几十倍的指令，常用”几个线程一行”的折中。实际的库（cuSPARSE）会按行长分组、对不同的组用不同的内核，或者换别的格式（ELL、COO、混合格式）。没有一种分工在所有矩阵上都最好——这句话在第 3 讲的归约里已经出现过一次。

4.6 2:4 结构化稀疏：规格表上”含稀疏”的数字

第 2 讲读规格表时说，张量核心的峰值默认”含 2:4 结构化稀疏”。04 号脚本第 4 部分演示了它：每 4 个连续的权重里只保留绝对值最大的 2 个，其余置 0；压缩存储时只存这 2 个值，外加每个值在组内的位置（2 位）。用压缩形式算 W·x 和用剪枝后的稠密矩阵算，最大差 2.4e-07。FP16 权重压缩后每个原始权重占 9 位，是稠密的 56.2%。

2:4 稀疏的非零元比例固定是 50%——按 §4.4 的算账，CSR 在这个比例下毫无优势。但它的模式是固定的：位置只要 2 位就能记下，硬件（Ampere 起的张量核心）可以直接跳过一半的乘法，同样的时间做 2 倍的有效运算。代价是模型要按这个模式剪枝（通常还要微调恢复精度），而且只对权重适用。

稀疏的两条路：非结构化稀疏（CSR 等）适合极稀疏的数据，靠省存储和省访存取胜；结构化稀疏（2:4、块稀疏）适合神经网络，靠和硬件配合取胜。

5. cuBLAS · CUDA Basic Linear Algebra Subroutines

5.1 为什么用库

slides 第 45 页：cuBLAS 是 CUDA 的基础线性代数子程序库，一个专门做通用矩阵乘法（GEMM）等运算的轻量库。§1–§3 把朴素矩阵乘法从 2.5% 提到了约 80%（FP32）；要用上张量核心、在每一代 GPU、每一种数据类型、每一种矩阵形状上都接近峰值，需要的东西远不止分块（§5.5）。实践中，矩阵乘法几乎总是调用库；自己写 GEMM 是为了理解它。

5.2 接口

slides 第 46–48 页：

cublasStatus_t cublasCreate(cublasHandle_t *handle);         // 使用前必须调用：创建库的上下文
cublasStatus_t cublasDestroy(cublasHandle_t handle);         // 使用后必须调用
cublasStatus_t cublasSdot(handle, int n, const float *x, int incx, const float *y, int incy, float *result);
cublasStatus_t cublasSgemv(handle, cublasOperation_t trans, int m, int n, const float *alpha,
                           const float *A, int lda, const float *x, int incx, const float *beta, float *y, int incy);
                           // y = α · op(A) · x + β · y
cublasStatus_t cublasSgemm(handle, cublasOperation_t transa, cublasOperation_t transb, int m, int n, int k,
                           const float *alpha, const float *A, int lda, const float *B, int ldb,
                           const float *beta, float *C, int ldc);
                           // C = α · op(A) · op(B) + β · C

函数名里的 S 表示单精度（float），D 双精度，H 半精度；op(X) 是 X（CUBLAS_OP_N）或 Xᵀ（CUBLAS_OP_T）。handle 保存库的状态（使用的流、工作空间等），创建一次、反复使用。

5.3 列优先与 leading dimension

BLAS 诞生于 Fortran 的年代，Fortran 的数组是列优先的，所以 cuBLAS 也按列优先解释矩阵。05 号脚本第 1 部分：缓冲区 [1, 2, 3, 4, 5, 6] 按行优先读是 [[1,2,3],[4,5,6]]，按列优先读成 3 × 2 是 [[1,4],[2,5],[3,6]]——同一段内存，行优先看是 X，列优先看就是 Xᵀ。

leading dimension（lda、ldb、ldc）是”相邻两列的起点相隔几个元素”。整块矩阵时 ld 等于行数；对一个大矩阵的子块运算时，指针指向子块的起点，ld 给大矩阵的行数——不用拷贝。例如 5 × 4 列优先矩阵从 (1, 1) 开始的 3 × 2 子块：指针偏移 1 + 1 × 5 = 6，ld = 5。这和 P1 §5 的步长是同一个思想（列优先时行步长 1、列步长 ld）；incx、incy 则是向量的步长。BLAS 的接口从一开始就是”指针 + 步长”。

05 号脚本第 2 部分用 NumPy 按 BLAS 的语义实现了 sgemm / sgemv / sdot：四种 transa / transb 组合、子块（偏移 13、lda = 6）、incx = incy = 2 都与 NumPy 一致（最大差 ≤ 2.4e-07）；lda 写得比行数小时报参数错误。传错一个整数，库只会报”第几个参数非法”，或者更糟——静默地算错。

5.4 行优先张量怎样调用列优先的库

我们的张量（NumPy、PyTorch、miniTorch）都是行优先的。想算 C = A · B，又不想做任何转置拷贝，怎么办？05 号脚本第 3 部分：

• 行优先存储的 A（M × K）在列优先的库眼里是 Aᵀ（K × M），B 是 Bᵀ；我们想要行优先的 C，在库眼里就是 Cᵀ。
• 而 Cᵀ = (A·B)ᵀ = Bᵀ · Aᵀ。

所以只要把 A、B 的位置交换、m 和 n 交换：sgemm('N', 'N', m = N, n = M, k = K, B, ldb = N, A, lda = K, C, ldc = N)，结果缓冲区按行优先读出来就是 A · B（与 NumPy 差 4.6e-08）；不交换直接调用则错（差 2.98）。PyTorch 在 GPU 上调用 cuBLAS 时就是这样处理行优先张量的；张量是转置视图时，改用 ‘T’ 参数，同样不拷贝。

5.5 从朴素到库：roofline 阶梯

05 号脚本第 4 部分把本讲的各种写法放在同一张表里（A100 80GB PCIe，FP32 数据）：

块分块的比值：一个块算 BM × BN 的输出，每一步读 A 的 BM × BK 和 B 的 BK × BN，做 2·BM·BN·BK 次运算，比值 = 2·BM·BN / ((BM + BN) × 4 字节)，和 BK 无关——块越大、越接近正方形，比值越高。每个线程用寄存器算一个 8 × 8 之类的小块，块就能突破”每块 1024 个线程”的限制（§2.6）。到 128 × 128，FP32 已经计算受限。

但张量核心的平衡点是 81：同样的分块只能用到约 40%–53%。要把张量核心喂饱，cuBLAS（以及 CUTLASS、Triton 生成的内核）还要做很多事：用张量核心的矩阵指令；为每一代 GPU、每种数据类型、每类矩阵形状挑选分块大小（第 2 讲的波次量化也在考虑之内）；用双缓冲和异步拷贝（Hopper 上是 TMA，第 3 讲 §1.6）让搬运和计算重叠；让 L2 缓存在块之间复用数据；把偏置、激活函数融合进结果的写回（cuBLASLt 的 epilogue，第 10–11 讲的算子融合）。这就是 slides 总结里那句”现成可用的向量、矩阵-向量、矩阵-矩阵运算”背后的工作量。

6. 小结 · Summary

1. 计算–全局内存访问比决定访存受限的内核能跑多快。朴素矩阵乘法 0.25 FLOP/字节，在 A100 上只能用到 FP32 峰值的 2.48%。
2. 全局内存约 500 个周期，共享内存约 5 个，寄存器约 1 个。局部数组会落到”本地内存”（其实在显存里）。
3. 分块：线程块把 A、B 的小块搬进共享内存，块内反复使用，全局读取降到 1/TILE、比值升到 0.25 × TILE；结果与朴素版逐位相同。两个 __syncthreads() 分别防”读之前没写完”和”写之前没读完”；边界要向上取整、补 0、不能提前 return；块不能太小（合并失效）也不能超过 1024 个线程。
4. 线程粗化 / 寄存器分块：让每个线程算多个输出，再省一层读取，突破线程数的限制。
5. 合并访存：一个线程束在同一条指令里访问连续、对齐的地址。转置的读和写总有一边按列走，用共享内存把它挪到片上。
6. bank 冲突：共享内存 32 个 bank，第 k 个字在 k mod 32 号；同一 bank 的不同字要排队。按列访问 [32][32] 是 32 路冲突，声明成 [32][33] 即可消除。
7. 稀疏：CSR 用 data / col_index / row_ptr；比例超过约 50%（gather 落空时约 10%）就不划算。一个线程一行 vs 一个线程束一行是访存和指令的权衡。2:4 结构化稀疏靠固定模式和张量核心配合。
8. cuBLAS：列优先、leading dimension、op(X)；行优先张量用”交换 A、B”的技巧直接调用。生产中的 GEMM 几乎总是用库。

7. 常见误区 · Common pitfalls

• “分块减少了计算量。” 运算总量一次没少，减少的是全局内存的读取次数；共享内存的读取和原来的全局读取一样多。
• “TILE 越大越好。” 一个线程一个输出时 TILE 不能超过 32（每块 1024 个线程）；共享内存、寄存器用多了占用率会掉；TILE 太小（例如 2）反而让合并访存失效。
• “一个 __syncthreads() 就够了。” 搬完之后和算完之后各需要一个，分别防两种竞争；漏掉任何一个，结果都会错，而且在真 GPU 上可能时对时错。
• “越界的线程直接 return 就好。” 在用栅栏的内核里，越界的线程也要参与搬运（写 0）和栅栏。
• “共享内存很快，怎么访问都行。” 同一个 bank 的不同字要排队；按列访问 32 × 32 的数组是 32 路冲突。
• “合并访存只看一个线程依次访问的地址是否连续。” 看的是一个线程束在同一条指令里一起访问的地址（第 2 讲 §4.8）。
• “稀疏总比稠密快。” 非零元比例要低于约 50%（gather 不命中时约 10%）才划算；50% 稀疏的权重要靠 2:4 这样的结构化稀疏和专门的硬件。
• “cuBLAS 和我的张量一样是行优先的。” cuBLAS 是列优先的；直接把行优先缓冲区传进去会算错，要用交换技巧或 transa / transb。
• “自己写的分块 GEMM 能比 cuBLAS 快。” 在一般形状上几乎不可能；库用了张量核心、异步流水线、按形状调优等一整套技术。

8. 练习 · Exercises

1. 向量加法 C = A + B 和 AXPY（y = a·x + y）的计算–全局内存访问比各是多少（float32）？
2. 在 H100 SXM（3.35 TB/s、FP32 67 TFLOPS）上，朴素矩阵乘法最多能达到 FP32 峰值的多少？
3. N = 1024 的矩阵乘法，朴素版和 TILE = 16 的分块版各读全局内存多少次（按线程计）？
4. TILE = 32 时 As、Bs 合计多少字节共享内存？如果输入改成 FP16 呢？
5. 用自己的话解释分块矩阵乘法里两个 __syncthreads() 各防的是什么竞争。
6. N = 100、TILE = 16 时，正确的相数是多少？slides 的 ph < N / TILE_WIDTH 会做几相、漏掉几项？
7. 一个线程束的线程 tx 读 A[2 * tx]（float，A 从对齐的地址开始），产生几个 32 字节事务？读 A[tx + 3] 呢？
8. __shared__ float s[32][32]; 一个线程束（threadIdx.x = 0..31）读 s[threadIdx.x][0] 几路冲突？读 s[0][threadIdx.x] 呢？如果声明成 s[32][33]，读 s[threadIdx.x][0] 呢？
9. __shared__ float s[64]; 线程 tx（0..31）读 s[2 * tx]，有 bank 冲突吗？几路？
10. 写出矩阵 [[0, 5, 0], [7, 0, 0], [0, 0, 0], [1, 2, 3]] 的 CSR 三个数组。
11. 如果值用 8 字节的 double、列号用 4 字节的 int，CSR 在非零元比例多少以下才比稠密存储省空间？
12. A 是行优先的 M × K，B 是行优先的 N × K，想得到行优先的 C = A · Bᵀ（M × N）。sgemm 的参数应该怎样填？

参考答案

1. 向量加法：1 次运算，读 8 字节 + 写 4 字节 → 1/12 ≈ 0.083 FLOP/字节；AXPY：1 次乘 + 1 次加 = 2 次运算，读 x、y 写 y 共 12 字节 → 1/6 ≈ 0.17。两者都严重访存受限。
2. 3.35×10¹² × 0.25 = 837.5 GFLOPS，占 67 TFLOPS 的 1.25%（H100 的算力涨得比带宽快，比例比 A100 更低）。
3. 朴素：2N³ = 2,147,483,648 次；TILE = 16：2N³ / 16 = 134,217,728 次。
4. 2 × 32 × 32 × 4 = 8 KB；FP16（2 字节）时 4 KB。
5. 第一个（搬完之后）：保证块内所有线程都把自己负责的元素写进共享内存之后，才有人开始读——防”写后读”。第二个（算完之后）：保证所有线程都用完这一相的小块之后，才有人开始写下一相——防”读后写”。
6. 正确的相数 ⌈100/16⌉ = 7；slides 的写法做 100 // 16 = 6 相，漏掉最后 100 − 96 = 4 项（而且网格向上取整时会越界）。
7. A[2·tx]：32 个 float 跨 256 字节 → 8 个事务（有用 50%）。A[tx + 3]：起点偏移 12 字节，128 字节跨 5 个 sector。
8. s[tx][0]：字号 tx × 32，全在 bank 0 → 32 路冲突。s[0][tx]：连续的 32 个字 → 无冲突。s[32][33] 时 s[tx][0] 的字号 tx × 33，bank = tx → 无冲突。
9. 字号 0, 2, 4, …, 62 → bank 0, 2, …, 30, 0, 2, …, 30：每个偶数 bank 被两个不同的字占用 → 2 路冲突。
10. data = [5, 7, 1, 2, 3]，col_index = [1, 0, 0, 1, 2]，row_ptr = [0, 1, 2, 2, 5]（第 2 行是空行）。
11. 每个非零元 8 + 4 = 12 字节，稠密每个元素 8 字节：12·d < 8 → d < 2/3 ≈ 66.7%（忽略 row_ptr）。
12. 想要的 Cᵀ = (A·Bᵀ)ᵀ = B·Aᵀ。库眼里 A 的缓冲区是 Aᵀ（K × M，op = ‘N’），B 的缓冲区是 Bᵀ（K × N），而我们需要 B（N × K），所以 op = ‘T’：sgemm(‘T’, ‘N’, m = N, n = M, k = K, B, ldb = K, A, lda = K, C, ldc = N)。（用 05 号脚本的 sgemm 验证过，最大差 1.2e-07。）

9. 延伸阅读 · Further reading

• PMPP 第 4 版第 5 章（存储类型、分块、边界检查、存储用量对占用率的影响）、第 6 章（合并访存、隐藏延迟、线程粗化、优化清单）、第 14 章（COO、CSR、ELL、JDS 等稀疏格式与它们的 SpMV 内核）。
• NVIDIA CUDA C++ Best Practices Guide：“Memory Optimizations”（合并访存、共享内存与 bank 冲突、矩阵转置的例子）。
• Mark Harris, An Efficient Matrix Transpose in CUDA C/C++（NVIDIA 技术博客）：转置三个版本的实测带宽。
• Simon Boehm, How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance（博客）：从朴素到接近 cuBLAS 的十步，每一步都有实测数字，是 §5.5 阶梯的实践版。
• NVIDIA CUTLASS 文档：GEMM 的分层分块（线程块 → 线程束 → 指令）在模板库里的实现。
• Mishra et al. (2021) Accelerating Sparse Deep Neural Networks（NVIDIA）：2:4 结构化稀疏的硬件支持与”剪枝—微调”的流程。

10. 下一讲预告 · Next lecture

前四讲都在”底层”：GPU 怎样执行、内核怎样写快。第 5 讲往上走一层：深度学习框架与自动微分。loss.backward() 是怎样对任意网络算出每个参数的梯度的？slides 从通用的迭代学习算法（梯度下降）出发，讲计算图、拓扑排序、前向与反向模式的自动微分、向量–雅可比积，再到一个框架的组成（TensorFlow 式的占位符、变量、运算节点、会话）。我们会自己写一个迷你框架：标量和张量两层的自动微分（原课 HW2 的 topological_sort 和 backpropagate）、静态图和动态图两种风格，并把第 3 讲的 SIMT 内核接成它的一个”后端”，最后用它训练 P3 的影评分类网络。