极客工场

Appearance

共 4,979 字预计阅读约 17 分钟

矩阵乘法 tiling：把分块搬进 shared memory，让 GFLOPS 翻百倍 ​

更新时间：2026-09-25。矩阵乘法（GEMM）是 GPU 上被研究得最透的算子，几乎所有 CUDA 性能优化教材都从它开始——因为它把一个道理演示得最赤裸：算力从来不是白给的，喂不饱数据的 ALU 只会闲着。

写过 CUDA 的人多半有过同一段经历：照着教程敲出矩阵乘法，几千个线程跑得热火朝天，nsight 一看只有几个 GFLOPS，显卡标称的几十 TFLOPS 像是假的。问题不在线程开得不够多，而在每个线程都在反复穿越整条显存层级去取同一批数。本文把这段优化路径完整走一遍：朴素三重循环为什么慢、tiling（分块）利用了矩阵乘法什么样的数据复用结构、线程怎么映射到输出子块、内层 k 循环为什么需要两道 __syncthreads() 屏障，最后用算术强度把性能数字一笔笔算出来。

阅读本文前，建议先具备 GPU 内存层级的基本概念（全局显存、shared memory、寄存器各自的位置和带宽差异），可对照 GPU 内存层级与带宽；SIMT、warp 与合并访存的硬件背景见 GPU 架构深入：SIMT 与 warp 调度。文中代码用 CUDA C++ 书写，但同样的分块结构可以一一对应翻译到 ROCm、oneAPI 乃至 CPU 多线程加缓存分块，因为背后的矛盾是通用的：多层存储之间，数据搬运远比数据计算昂贵。

本文要回答的问题 ​

  • 朴素矩阵乘法核函数里，每个 C 元素都独立读一整行和一整列，全局显存到底被重复读了多少遍？
  • 为什么说矩阵乘法天然藏着一个 BLOCK 倍的复用机会？这个复用在数学上来自哪？
  • tiling 之后线程块（thread block）、线程和 C 的子块之间怎么对应？
  • 内层 k 循环里为什么要协作搬两块、又为什么必须放两道 __syncthreads()，放错位置会怎样？
  • BLOCK 取 16 还是 32？算术强度怎么算，性能为什么能从几个 GFLOPS 跳到几百个 GFLOPS？
  • shared memory 版之后，Tensor Core 版又凭什么再赢一个数量级？

一、朴素实现：O(N³) 的运算配 O(N³) 的访存 ​

先看最直白的 CUDA 写法。设 C = A × B，三个矩阵都是 N×N 的 float，每个线程负责 C 的一个元素：

cuda
__global__ void gemm_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) return;
    float acc = 0.0f;
    for (int k = 0; k < N; ++k)
        acc += A[row * N + k] * B[k * N + col];   // 每次迭代两次全局读
    C[row * N + col] = acc;
}

循环体做一次乘法和一次加法，共 2 次浮点运算，代价是两次全局显存读取。整个核的浮点运算量是 2N³，而全局内存加载次数也是 2N³ 量级——运算和访存一一对应，算术强度（arithmetic intensity，每字节访存换来多少次浮点运算）只有约 0.25 FLOP/B：8 次浮点运算对应 8 字节读取。

这个数字意味着什么？拿一块典型显卡算笔账：HBM/GDDR 显存带宽约 1 TB/s，峰值算力约 10 TFLOPS。带宽上限给这个核的性能天花板是 1e12 × 0.25 = 250 GFLOPS，而要让计算单元吃饱需要的喂数速率是 10e12 / 0.25 = 40 TB/s，四十倍于物理带宽。更糟的是实际访问模式更差：B 按列读取，同一 warp 内相邻线程访问的地址相隔 N 个元素，无法合并（coalesce），有效带宽再打折。实测朴素核通常只有几个 GFLOPS，不到峰值的千分之一。

病根是重复。A 的第 k 列要被计算第 k 个内积项的所有行使用，B 的第 k 行要被所有列使用；在朴素实现里，同一个 A 元素被不同线程分别从显存读 N 次，谁也不沾谁的光。N=4096 时，A、B 两个矩阵总共只有 128 MB，却支撑了约 137 GFLOP 的运算量——平均每字节数据被从显存搬出上千次。

有人会问：GPU 不是有 L1/L2 缓存吗，让缓存自动挡住重复读行不行？对这个模式基本指望不上。L1 容量只有几十到一两百 KB，而任意时刻在飞的工作集是整行整列、随 N 增长；不同线程进度错位，缓存里的数据还没等到被复用就被别的行挤出去，表现为持续的 capacity miss。缓存是被动的、对所有访问一视同仁；而矩阵乘法的复用是结构化的、程序员完全知道下一个要用谁——与其赌缓存的替换策略，不如显式控制一块程序员自己说了算的片上存储。这正是 tiling 和纯靠缓存的本质区别。

还要区分延迟和带宽这两种「慢」。单次显存读的延迟大约是几百个时钟周期，GPU 靠 warp 切换把它藏起来，所以朴素版的慢不是线程在干等一条读指令返回，而是切换完全部就绪 warp 之后，整体的取数吞吐仍然撑不起乘法器的胃口：每个周期从显存流进来的字节，只够配出峰值算力几百分之一的运算。nsight 里看到的显存带宽利用率可能已经不低，但计算单元（FMA pipe）利用率长期贴着地面——这是典型的带宽受限画像。

二、复用从哪来：一个 C 子块只依赖两条窄带 ​

矩阵乘法在数学上有个关键性质：输出 C 的任意一个子块，可以脱离矩阵其余部分独立算出来。把 C 切成 BLOCK×BLOCK 的小块，考虑行在 [m, m+BLOCK)、列在 [n, n+BLOCK) 的那个子块：

$$C_{tile} = \sum_{k} A_{m:m+B,,k} \cdot B_{k,,n:n+B}$$

随着 k 推进，需要的输入只是 A 的一条 BLOCK 行宽的横向条带和 B 的一条 BLOCK 列宽的纵向条带。把 k 也按 BLOCK 切段，每一段只用到 A 的一个 BLOCK×BLOCK 小块和 B 的一个 BLOCK×BLOCK 小块——两个小块一共 2·BLOCK² 个元素，却能产出子块内 BLOCK² 个输出元素各自的 BLOCK 次乘加，即 2·BLOCK³ 次浮点运算。

换句话说，每搬一个小块进来，块内每个元素平均被 BLOCK 个输出元素的计算各用一次：复用率正好是 BLOCK。朴素版让这份复用发生在显存里（每个线程自己去取），tiling 的思路就是把它挪到片上——搬进 shared memory，让一个线程块内的线程共享这批数据。

chip-design/intermediate/matrix-multiply-tiling.md D2 示意图

上图就是「分块数据复用」的核心视角：显存和计算单元之间插进一块由全块线程共享的高速存储，粗箭头从显存只走一次，细粒度复用全部在片上完成。shared memory 的带宽与延迟接近寄存器一两级，是这层 staging 的天然位置。

三、访存量怎么从 O(N³) 降到 O(N³/BLOCK) ​

把全矩阵的账算清楚。C 共有 (N/BLOCK)² 个子块；每个子块沿 k 方向要走 N/BLOCK 段，每段从显存搬 2·BLOCK² 个元素。于是全局访存元素总数为：

$$\left(\frac{N}{B}\right)^2 \cdot \frac{N}{B} \cdot 2B^2 = \frac{2N^3}{B}$$

运算量仍是 2N³ 不变，访存却除以一个 BLOCK。BLOCK=32 时，全局流量降为朴素版的 1/32；折合算术强度（float，每元素 4 字节）为：

$$I = \frac{2N^3}{(2N^3/B)\times 4\text{B} } = \frac{B}{4}\ \text{FLOP/B}$$

BLOCK=16 给出 4 FLOP/B，BLOCK=32 给出 8 FLOP/B。站在单个数据元素的视角更直观：一个 4 字节元素被搬进 shared memory 后，被 BLOCK 个输出计算各用一次、每次参与一次乘加，共服务 2·BLOCK 次浮点运算——BLOCK=32 时就是每个元素撑起 64 次 FLOP。这就是「32 块约 64 FLOP 复用」这个经验数字的来历：64 是每个元素对应的 FLOP 数，除以 4 字节即 16 FLOP/B 的单元素口径；两矩阵合计的整核口径为 8 FLOP/B。

8 FLOP/B 过没过山脊点（ridge point）取决于具体显卡，可参考 GPU roofline 性能模型 里带宽平台与计算平台的交汇分析。对主流数据中心卡，8 FLOP/B 已足以让算子从「纯带宽卡死」挪到接近计算受限区，这正是百倍加速的直接来源。

把抽象公式换成 N=4096、BLOCK=32 的实数：朴素版从显存搬运的元素数约 2N³，即约 1370 亿个、折合 549 GB 流量（这还没算非合并访问造成的额外事务）；tiling 版的全局流量是 2N³/B 个元素，约 42.9 亿个、17.2 GB，连前者的三十分之一。同一批 128 MB 的矩阵，被搬运的总字节从半张蓝光碟的量级缩到一部小安装包，而要做的 2740 亿次浮点运算一个不少。这组数字可以直接对照 profiler 里的 dram throughput 子项验证：tiling 版跑同样规模，显存子项读数应当显著下降，而 sm 计算活跃度反过来打满。

顺带一提，为什么复用率恰好是 BLOCK 而不是 BLOCK²：一个 A 小块里的单个元素，只参与子块内固定一行的 BLOCK 个输出，列方向的输出用不到它；同理 B 元素只喂固定一列。块级总计 2B² 个元素各被用 B 次，得到 2B³ FLOP。想进一步提高每字节的复用，就要让一个元素同时在更多输出里登场——这正是后面 Tensor Core 多级分块中「一个线程负责多个输出元素」「warp 级块再复用一次」的用武之地。

四、线程映射与内层 k 循环：两块一起搬，两道屏障 ​

映射规则很简单：一个线程块负责一个 C 子块，块内每个线程负责子块中的一个 C 元素（更进阶的实现里一个线程负责一小块，但入门版一对一最清楚）。二维线程网格 (BLOCK, BLOCK) 与二维子块天然对齐：

cuda
const int BLOCK = 32;
__global__ void gemm_tiled(const float* A, const float* B, float* C, int N) {
    __shared__ float As[BLOCK][BLOCK];
    __shared__ float Bs[BLOCK][BLOCK];

    int tx = threadIdx.x, ty = threadIdx.y;
    int row = blockIdx.y * BLOCK + ty;   // 本线程负责的 C 元素行
    int col = blockIdx.x * BLOCK + tx;   // 列
    float acc = 0.0f;

    for (int ph = 0; ph < (N + BLOCK - 1) / BLOCK; ++ph) {
        // ① 全块协作：把当前 k 段的 A、B 两个小块搬进 shared memory
        As[ty][tx] = A[row * N + ph * BLOCK + tx];
        Bs[ty][tx] = B[(ph * BLOCK + ty) * N + col];
        __syncthreads();                // 屏障一：搬完才算

        // ② 块内复用：每个线程用这两块数据累加自己的输出
        #pragma unroll
        for (int k = 0; k < BLOCK; ++k)
            acc += As[ty][k] * Bs[k][tx];
        __syncthreads();                // 屏障二：算完才允许下一轮覆盖
    }
    C[row * N + col] = acc;
}

两道屏障各自封死一种竞态，缺一不可：

  • 第一道（搬后算前）：线程之间进度不一，搬得快的线程开始做乘法时，邻居负责的元素可能还没写进 Bs。没有这道屏障，就是拿 shared memory 里的旧值或零值参与累加，结果随机出错。
  • 第二道（算后下一轮搬前）：进入下一个 ph 后所有线程会原地覆盖 As、Bs。若还有线程没读完当前段，它读到的就是下一段的数据。屏障保证「旧两块被所有人读完」才允许覆盖。

这是一个经典的双屏障 ping-pong 结构：共享缓冲区被反复复用，生产（搬运）和消费（计算）两个阶段各靠一道屏障隔开。搬运阶段还顺便把朴素版的 B 列访问改成了块内协作的合并读——同一 warp 的线程在 A[...] 和 B[...] 上访问连续地址，显存带宽利用率同时被修好。

chip-design/intermediate/matrix-multiply-tiling.md D2 示意图

第二张图给出的是线程层次的映射：网格上的块对应输出上的子块，块内线程对应子块内元素，所有线程沿 k 方向以锁步（lock-step）节奏推进各段。

值得花一行代码看懂搬运阶段为什么同时修好了合并访问。warp 沿 threadIdx.x 方向成组（即同一 ty、相邻 tx 的 32 个线程）：在 A[row*N + ph*B + tx] 上它们的地址只差 1，恰好合成一段连续的 128 字节事务；在 B[(ph*B + ty)*N + col] 上 col = blockIdx.x*B + tx，同样连续。也就是说，搬运时每个线程各取一个、拼成整块，取的过程在显存看来是整齐的长条形突发流；而到了块内计算阶段，访问模式可以随意（As[ty][k] 按行、Bs[k][tx] 按列），因为那已经是片上访问，与显存事务无关。对外合并、对内随意——这正是双缓冲区分工的精妙之处：合并性只需要在跨越内存层级边界的那一跳保证。

另外注意两个 shared 数组的下标方向：As[ty][tx] 与计算时的 As[ty][k] 都是行优先，同一 warp 固定 ty、遍历 k……不过内层 k 循环每个线程独立遍历，不存在跨线程合并问题，只需关注 bank 分布。这类下标方向一旦写反，程序结果仍可能正确（只是搬运映射错位到别的块），性能却悄悄回退，是 code review 时值得专门核对的一处。

五、性能数字：几 GFLOPS、几百 GFLOPS，再到 Tensor Core ​

把三个版本放在一起对比，能看清每一赢次赢在哪：

版本全局访存量算术强度（float）典型实测算力瓶颈
朴素三重循环O(N³)，B 列访问不合并~0.25 FLOP/B几个 GFLOPS（约 2–5）显存带宽 + 非合并访问
shared memory tiling（BLOCK=16/32）O(N³/BLOCK)4 / 8 FLOP/B数百 GFLOPS 起，调优后数 TFLOPS逐渐转为计算/延迟受限
Tensor Core（WMMA/CUTLASS，FP16/BF16）再配合多级分块与寄存器复用数十至数百 FLOP/B数十至上百 TFLOPS矩阵乘加阵列利用率

数字随显卡代际差别很大，量级关系是稳定的：tiling 这一步通常带来 50–200 倍提升，是整个优化链中回报率最高的一跳。一个朴素版跑到 3 GFLOPS 的核，换 16×16 分块后常见 150–300 GFLOPS，换 32×32 并处理好 bank conflict 后可再上一个台阶；同一张卡上 Tensor Core 版则直接进入几十 TFLOPS 的区间。这种「改数据放在哪、不改数学」的百倍跃迁，在系统优化里并不多见，也是 GEMM 成为教学经典的原因。

BLOCK 取 16 还是 32 不是越大越好——shared memory 容量、每块线程数上限（常见 1024）、寄存器占用决定的 occupancy 互相牵制；32×32 配一个线程一个输出元素在老卡上可能寄存器吃紧，实战中更常见 16×16 的块让每线程算多个元素（thread coarsening）。直觉上「BLOCK 越大、复用越多」，但块一大，同时驻留一个 SM 的块数就少，可切换的 warp 数下降，延迟隐藏能力变弱；块太小又让复利好不容易攒下的算术强度不够越过山脊点。实战选点是在这张三维约束网里实测出来的，不是公式推出来的。

此外 shared memory 被分成 32 个 bank，As[ty][k] 和 Bs[k][tx] 这两类访问在 warp 内可能落在同一 bank 产生冲突，访问被串行化，具体成因和 padding 解法见 shared memory 与 bank conflict。

走到这一步后，CUDA Core 版 GEMM 离峰值仍有距离：每条指令只做一对标量乘加。Tensor Core 把一条指令变成一个小矩阵（如 16×16×16）的乘加，配合 WMMA 接口或 CUTLASS 的多级流水分块，才能摸到标称的几十上百 TFLOPS，这条路线见 Tensor Core 矩阵加速原理。注意 tiling 的思想在那里原封不动，只是「块」的层次从「shared memory 块」扩展为「寄存器块、warp 块、Tensor Core 指令块」的金字塔。

六、动手时最容易踩的坑 ​

  • 漏第二道屏障：单测小矩阵、单段 k 时结果正确，矩阵一大、段数一多就出现随机误差——因为竞态只在覆盖共享缓冲区时才发作。
  • 边界处理偷懒：N 不是 BLOCK 整数倍时，越界搬运会踩到非法地址或读到垃圾值。入门做法是给 kernel 加行、列判断并在越界处向 shared memory 填零。
  • 只 tiling 不查合并：线程块内部是二维的，搬运用的下标要保证同一 warp（通常沿 x）访问连续地址，否则带宽红利吃不到。
  • 把 BLOCK 当性能旋钮硬调：块大小与 occupancy、bank conflict、寄存器压力联动，改动后应以 nsight 实测为准，而不是凭直觉。
  • 累加器不用寄存器：acc 必须留在局部标量（编译器分配寄存器），若写回 shared memory 或全局内存做累加，tiling 省下的访存又被加回来。

验证优化是否生效，不必靠感觉，对着 profiler 的三个指标即可：dram 读字节数应近似 2N³/B 个元素（而不是 2N³）；warp 级的非合并访问计数应趋近于零；FMA 吞吐占比应从几乎为零升到与 roofline 预测相符的区间。三者同时到位，才说明「搬得少、搬得齐、算得满」三件事都成立。如果字节数降了但 FMA 没涨，多半是 bank conflict 或 occupancy 不足把瓶颈从显存挪到了 shared memory 与调度上，接着按对应文章排查即可。

这一篇在知识体系里的位置 ​

本文处在 GPU 编程优化这条线的中段：往上，GPU 内存层级与带宽 给出了全局显存与 shared memory 的数量级差异，GPU 架构深入：SIMT 与 warp 调度 解释了线程与合并访问的硬件含义；往旁边，shared memory 与 bank conflict 解决块内复用的二次折损，GPU roofline 性能模型 提供算术强度的系统分析框架；往下，Tensor Core 矩阵加速原理 是同一套分块思想在矩阵乘加硬件上的延续。CPU 侧的缓存分块（loop tiling）与此同源，可以对照理解：片上存储换数据复用，是所有多层存储层级通用的优化母题。

总结 ​

矩阵乘法 tiling 的本质是把「同一批数据被读 BLOCK 次」从显存挪到片上：一个线程块认领一个 C 子块，沿 k 方向每段把 A、B 两个小块协作搬进 shared memory，靠两道 __syncthreads() 隔开搬运与计算，使全局访存从 O(N³) 降到 O(N³/BLOCK)、算术强度升到 BLOCK/4 FLOP/B，性能因此从几个 GFLOPS 跳到数百 GFLOPS——而这正是通向 Tensor Core 多级分块金字塔的第一级台阶。