企业官网建设流程全解析

我最初以为“矩阵乘法跑赢 cuBLAS”是那种只有专家才能碰的题目，直到自己花了两周时间，把一块消费级显卡上的 SGEMM 从 200 GFLOPS 一路优化到 cuBLAS 的 83%，我才意识到这件事离普通 CUDA 开发者并没有那么远。这个优化过程的核心不是玄学，而是 Warp Tiling、共享内存布局和异步拷贝这三件事。这篇文章就是把我的完整实录整理出来，从为什么矩阵乘法值得手写，到每一版 kernel 的改法，再到最后怎么用 profiler 找瓶颈，全部说清楚。适合正在学 CUDA、准备高性能计算方向，或者面试被问到“你怎么优化 GEMM”的朋友。

1. 怎么想到要和 cuBLAS 比

1.1 矩阵乘法是显卡的试金石

矩阵乘法，也就是 C = A × B，是所有基础线性代数里最有代表性的一道题。它的计算量一眼就能算出来：如果 A 是 M×K，B 是 K×N，那么 C 是 M×N，总计算量是 2·M·N·K 次浮点运算。以 4096×4096×4096 为例，一次全尺寸 SGEMM 就是 2×4096³ = 137.4 GFLOP 的运算量。这个数字很大，但 GPU 本身擅长做并行乘加，所以它成了检验一块显卡真实算力的标准项目。

很多面试官喜欢问矩阵乘法，并不是因为它简单，而是因为它能同时考察线程编排、存储层次、访存模式、指令级并行和性能分析。你把你自己的 kernel 和官方库放在同一块 GPU 上跑，差距就是你和 GPU 优化高手之间的距离。我第一版 naive kernel 大概只有 0.2 TFLOPS，那时候和 cuBLAS 的差距是十倍还多，没有任何讨论的空间。

所以我把矩阵乘法当作“硬碰硬”的试金石：不用张量核心，不用 FP16，就用最普通的 FP32 乘加指令，看纯手写代码到底能贴近官方库到什么程度。这个目标听起来简单，但做起来要拆的东西非常多。

1.2 cuBLAS 到底赢在哪

一开始我也天真地以为 cuBLAS 只是“官方库所以优化得多”。真正看下去才发现，它赢在几个层面。第一是分块策略极其丰富，同一个尺寸它会测试多种 block tile、多种 BK、多种调度方式，甚至对不同的 K 段采用不同的 kernel。第二是共享内存和寄存器都用得非常狠，通常每个线程负责一个 8×8 甚至更大的输出子矩阵，而不是每个线程只算一个输出点。第三是它有非常好的软件流水线，能通过预取隐藏访存延迟，让访存单元和计算单元尽量同时工作。

如果不看源码，你会觉得这些都很神秘。但它们本质上就是我们手写 kernel 时做的那些事，只不过官方库把这些事情做到了极致，并且用了大量 profiling 数据去调参。拿 FP32 SGEMM 来说，官方库在没有张量核心参与的路径上，也能把一块高带宽显卡的浮点利用率推到很高的水平。所以“到达 cuBLAS 的 83%”实际上意味着，我已经把手工 GEMM 做得像模像样了，剩下的差距更多来自参数枚举、汇编级微调这类体力活。

2. 设计一个 Warp Tiling 的 SGEMM

2.1 从线程、Warp 到 Block 的三层映射

矩阵乘法并行化的第一步，是把输出矩阵 C 切成很多小块，让不同的线程块去算不同的小块。我选的 block tile 是 128×64，也就是一个线程块负责 C 上的 128 行、64 列范围。整个线程块一共有 128 个线程，正好 4 个 warp。网格划分就是 (N/64, M/128)，如果有 4096×4096 的输出矩阵，就会生成 32×64 = 2048 个 block，每个 block 独立跑。

关键的一步是“Warp Tiling”：不要把线程块里的任务直接切成 128 个独立线程去算，而是先用 4 个 warp 去切分行方向，每个 warp 负责 32×64 的输出区域。再往下一层，每个 warp 里的 32 个 lane 又分成 4 组，每组 8 个 lane，每个 lane 负责一个 8×8 的小块。这样一层层套下来，最终每个线程算 64 个输出值。

这样做的好处非常明显。每个线程的 64 个结果都在寄存器里复用，同一个 A 值要和 8 个不同的 B 值做 8 次乘加，同一个 B 值也要和 8 个不同的 A 值做 8 次乘加。和“每个线程只算一个输出”的朴素做法相比，数据从寄存器读一遍就被复用了很多次，共享内存和全局内存的访存次数都大幅下降。GPU 最怕的不是计算，而是数据搬来搬去，Warp Tiling 本质上就是把“搬数据”这件事尽量挡在寄存器外面。

2.2 共享内存里的轨道布局

有了线程映射，还要解决数据供应的问题。全局内存虽然大，但延迟很高，所以矩阵 A 和 B 的块要先进共享内存。我选的内层分块大小是 BK=16，也就是循环里每次取 A 的 128×16 一小块，取 B 的 16×64 一小块。共享内存里安排两块缓冲区：

__shared__ float As[2][BK][TILE_M + 1]; __shared__ float Bs[2][BK][TILE_N + 1];

这里 TILE_M=128，TILE_N=64，加 1 是为了避免 bank conflict。共享存储器一共有 32 个 bank，每隔 4 字节换一个 bank。如果行没有任何 padding，那同一行相邻元素刚好落在相邻 bank 上，看起来没问题；但当你按行方向取不同数据时，可能多个线程访问同一个 bank，造成串行化。加一个 float 的 padding 后，行的步长从 128 变成 129，129 对 32 取余是 1，这样连续行的起始 bank 位置就不一样了，冲突概率会低很多。

共享内存大小也顺便算一下：As 是 16×129 个 float，Bs 是 16×65 个 float，一个 stage 加起来是 3104 个 float，约 12 KB。双缓冲就是约 24 KB，完全能塞进常见的 48 KB 静态共享内存范围。这个尺寸设计不是随便拍的，它让一个 SM 上还能放好几个 block，保证足够的并行度去掩盖延迟。

2.3 双缓冲与异步拷贝

矩阵乘法循环的长度是 K/BK 次，每次都要把新的 A、B 小块从全局内存搬进共享内存。最笨的办法是“先搬完，再算，再搬下一块”，但这会让计算单元在等数据时发呆。我的做法是双缓冲加异步拷贝：算第 i 块的同时，用 cp.async 把第 i+1 块从全局内存预取到另一块共享内存中。

cp.async 是 Ampere 之后架构提供的异步拷贝指令，它和普通 global load 最大的区别是：数据从全局内存到共享内存的过程不需要经过寄存器，也不会把当前线程的计算流卡住。你可以把它想象成在厨房做菜时，水龙头已经在另一个水池里放水，而不是你站在水池边等水满了才开始切菜。

这里有个容易搞错的点：双缓冲不只是把数组加一个[2]就完事了。同步逻辑必须处理好。我用的是两阶段 pipeline：每个迭代先把下一 stage 的拷贝发出去，再等待当前 stage 的数据就绪，然后用__syncthreads()保证整个 block 里的线程都看到这一 stage。如果同步写错，轻则计算结果错，重则死锁。后面实操部分我会把伪代码写出来。

3. 实操落地：内核代码与关键参数选择

3.1 第一版：先保证正确

我没有一上来就写花里胡哨的版本，而是先写了一个最简单的共享内存版本：每个线程只算 C 的一个点。线程块仍然是 128×64，每个线程加载 A 的一个值和 B 的一个值，乘加后写回 C。它的优点是正确性极好验证，一旦跑出来的结果和 cuBLAS 一致，就说明分块逻辑、索引换算、循环边界都没问题。

这个版本性能当然很惨。因为每个输出点只做一次乘加，却要从共享内存读两次数据，访存开销远大于计算。我实测同一块卡上大约只能跑到 0.8 TFLOPS 左右，naive 版本就更不用说了。不过这个版本是后面所有优化的“正确性底座”，后面每次改动都会和它做同一份随机矩阵的 diff。

很多人会跳过这一步直接上 register tiling，结果出了 bug 也不好定位。我的建议是，第一版哪怕再慢也要写完整，因为矩阵乘法调索引本身就很烦，如果一开始就混合了 tile 映射、双缓冲、cp.async，任何一个小错都会让你排查到崩溃。

3.2 第二版：把一个 Warp 当成计算单元

第二版开始引入真正的 Warp Tiling。我给每个线程分配 8×8 的输出矩阵，用 8 个寄存器存 A 的 8 个值，8 个寄存器存 B 的 8 个值，64 个寄存器存 C 的输出。这 64 个寄存器是性能的关键，因为它们在最内层循环里反复被读和累加。

内层循环的代码逻辑大概是：

for (int k = 0; k < BK; ++k) { float a[TM]; // 从 As 读 8 个 A 值 float b[TN]; // 从 Bs 读 8 个 B 值 for (int i = 0; i < TM; ++i) for (int j = 0; j < TN; ++j) c[i][j] += a[i] * b[j]; }

每读一组 a 和 b，就做 8×8=64 次乘加。在这 64 次乘加中，a 和 b 都在寄存器里，不再产生新的访存。这就是“寄存器级数据复用”的威力。第二版只是把朴素版本改成这个结构，性能就立刻翻了好几倍，实测能到 1.6 TFLOPS 左右。但这时候还没有加双缓冲，每次迭代都要等共享内存加载完才能继续算，卡在访存延迟上的时间很多。

3.3 第三版：铺平内存访问与控制 Bank Conflict

第三版要做的事情很具体：让共享内存的访问尽量不冲突，同时用向量化 load 把带宽打满。

先说 bank conflict。共享内存一次周期只能服务若干个 bank，如果同一 warp 里的 32 个线程恰好访问了同一个 bank 上的不同地址，硬件就会把它们拆成多次访问，性能直接按倍数下降。我在 2.2 里说的 padding 就是为这个服务的。具体到我的布局里，As 的行数是 128，加一个 padding 后变成 129，B 同理。这样在访问连续行时，不会所有线程都撞到同一组 bank 上。

还有一点容易被忽略：同一 warp 里如果有好几个线程在读同一个共享内存地址，硬件是可以做广播的，不产生冲突。所以在写索引的时候，我会刻意让某些线程组访问同一个地址，减少真正被访问的“地址数”。

另外，我还在加载阶段用 float4 一次性读 4 个 float。对齐好的 float4 可以把 load 指令数降到原来的四分之一，指令发射压力小很多。代价是 padding 和对齐要小心，如果行宽不是 4 的倍数，float4 就会错位。实际我选择让每一行都按 4 的倍数起步，宁可浪费一点点空间，也要保住向量化。

3.4 第四版：预取数据，把等待藏起来

第四版加入异步拷贝和双缓冲。这一步让性能第一次有了“跑起来很爽”的感觉。原先每个 K 循环迭代都是“等数据 → 计算 → 等数据”，现在变成“一边算第 k 块，一边加载第 k+1 块”。计算单元大多数时候都在做乘加，而不是空等。

代码层级的伪代码如下：

// 预取第一块 cp_async(As[0], A + offset0, ...); cp_async(Bs[0], B + offset0, ...); cp_async_commit(); for (int k0 = 0; k0 < K; k0 += BK) { int cur = (k0 / BK) & 1; int next = cur ^ 1; if (k0 + BK < K) { cp_async(As[next], A + offset_next, ...); cp_async(Bs[next], B + offset_next, ...); } cp_async_commit(); cp_async_wait<1>(); __syncthreads(); // 用 As[cur], Bs[cur] 计算这一个 128×64 小块 compute_warp_tile(As[cur], Bs[cur], ...); __syncthreads(); }

这里cp_async_wait<1>表示允许还在路上的最新一组异步拷贝未完成，但更早的那组必须完成。这正好保证当前要用的 cur 块已经就绪。第二次__syncthreads()是防止某个 warp 提前进入下一轮，把别人还在用的共享内存覆盖掉。

加了这版之后，我在 4096³ 的 FP32 用例上跑到了大约 2.6 TFLOPS，离这块卡上 cuBLAS 的 3.15 TFLOPS 已经不远了。83% 的成绩基本就是这个版本测出来的。

4. 性能实测：83% 是怎么来的

4.1 跑分环境与口径

先说跑分环境。我是在 WSL2 里配的 CUDA 环境，编译用 nvcc，运行时用同一个 CUDA 工具包，所以没有出现 Windows 驱动和 Linux 工具链打架的问题。看性能之前，建议先确认两件事：nvidia-smi里显示的驱动版本，以及nvcc --version或者cudart的版本，这两个版本不必完全一致，但你的代码编译和运行要使用同一套 runtime。

矩阵尺寸选的是常规的 M=N=K=4096，数据类型为 FP32。性能计算方式是：

TFLOPS = 2 * M * N * K / 时间(秒) / 1e12

我用同一尺寸跑 100 次取中位数，避免频率波动带来的干扰。对正确性，我还会和 cuBLAS 的结果做相对误差比较，绝对误差控制在 1e-3 以内，不然优化半天跑出错误结果就白忙了。

4.2 几个版本的数据对比

我把自己手上几个版本的实测数据整理了一下，一张表看下来就清楚每一步优化到底值多少钱：

注意，这些数字都是同一个运行环境下测的，不同显卡、不同尺寸、不同架构结果会不一样。比如矩阵尺寸换成 1024³，cuBLAS 的优势依然存在，但比例可能有变化；换成非对齐的 1000 之类尺寸，官方库还会表现出更大的调参优势。所以别把 83% 当成一个固定答案，它更多是告诉你可以做到“很接近”。

4.3 用分析器找瓶颈

跑完分之后，我并没有因为“看起来还行”就停手，而是用分析器把每个版本的瓶颈找出来。这里最常用的工具是 ncu。我会重点看几个指标：共享内存 bank conflict 数量、寄存器溢出情况、SM 上 warp 的平均占用率，以及访存等待占比。

第一版 naive 的问题很容易看，很多周期都花在 memory 等待上，SM 里的计算单元利用率很低。Warp Tiling 版本之后，计算占比明显上升。但我也遇到过很诡异的现象：代码看起来应该很快，实际帧数却突然掉一半，最后用 profiler 一查，发现寄存器溢出，大量局部变量被压到本地内存里，性能直接崩塌。

另一个值得关注的指标是stall类型。如果大量 warp 都是long_scoreboard这种等待全局访存返回的停顿，说明预取还不够；如果barrier停顿占比高，说明__syncthreads()太频繁或者调度粒度不对。这些指标读起来稍微有点门槛，但绝对是手写 kernel 的“仪表盘”，比瞎猜靠谱得多。

5. 常见问题与排查心得

5.1 共享内存 Bank Conflict 的“隐形降速”

共享内存 bank conflict 是最容易隐蔽的问题。它的可怕之处在于，程序不会报错，结果也正确，就是性能提不上去。我第一次从 1650 涨到 2050 的跨越，就是靠修 bank conflict 拿到的。

怎么排查？如果你不想装 profiler，可以先手工推索引：假设 warp 里有 32 个线程同时访问共享内存地址addr(tid)，看addr % 32的分布，重复越多冲突越重。最典型的错误是行宽恰好是 32 的整数倍，导致同一列上的线程全落在同一个 bank。padding 到“和 32 互质”是成本最低的修法。

不过也要注意，加了 padding 之后，共享内存容量会稍微变大。如果 BK 很大、双缓冲又开着，可能超过静态共享内存限制。这时可以用动态共享内存，或者在 launch 时设置更大共享内存配置。我一般先看 profiler 的l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum，那个值如果成百上千，肯定要先修这里。

5.2 寄存器溢出：性能直接从山峰落回盆地

寄存器溢出是我这次踩过最大的坑。每个线程 8×8 的输出意味着至少 64 个寄存器要用来保存 C 的值，再加上 A、B、地址、循环变量，很容易逼近 128 甚至 255 个寄存器上限。寄存器是有限的，如果被压到本地内存，所有看似漂亮的优化都会瞬间作废。

我当时的症状是：加上 cp.async 之后，性能不升反降。排查发现，为了写双缓冲逻辑，我引入了太多临时变量，编译器又使用了比较激进的展开策略，把寄存器用爆了。解决方法有三个方向：一是减少线程块的 tile，让每线程输出从 8×8 降到 4×8；二是把某些地址计算放到循环外复用；三是用__launch_bounds__(128, 4)这类约束限制编译器不要放飞自我，必要时-maxrregcount也要会用。

如果你想看寄存器到底用了几路，编译时加-Xptxas -v，它会直接打印每个 kernel 的寄存器数、shared memory 数和 spill 情况。我看到spill stores和spill loads不为 0 的时候，就会回去调整。

5.3 CUDA 多版本、WSL2 环境带来的坑

我在 WSL2 里做测试时，确实也遇到过 CUDA 版本不匹配的问题。最常见的是：机器上装了多个 CUDA 工具包，软链指到一个版本，但某些第三方库或者 Python 包是拿另一个版本编译的，跑起来就报版本不一致错误。遇到这种问题，不要急着重装驱动，先确认你的PATH和LD_LIBRARY_PATH指向哪个版本，再用nvcc --version和实际运行的cudart对一下。

另外，很多人装 CUDA 是为了编译 OpenCV 或者跑深度学习框架。这些库对 CUDA 版本和配套库版本有严格绑定，比如某个框架版本可能要求特定的 CUDA 和配套数学库版本。我的经验是：能不动环境就不动环境，项目尽量用容器或者独立工具链隔离版本。WSL2 的好处是可以用 Windows 侧驱动，Linux 侧只装 CUDA 工具包，省去了非常多驱动兼容问题，但坏处是如果你从网上复制安装命令，很容易装出一堆半吊子环境。

5.4 别把线程数和性能划等号

最后一类“问题”其实来自思维惯性。很多新手觉得，GPU 优化就是堆线程数，线程越多越快。我在实测里验证过：当共享内存和寄存器资源被占满后，再多开线程反而会因为 block 调度不过来、资源排队严重而变慢。

正确的思路是“让每个线程做更多有用功”，而不是“让更多线程各做一点点”。把 C 矩阵拆成 128×64 的块，每线程 8×8 输出，已经让总共 4096 个 block 在 2048 个线程的 GPU 上有了足够的调度余地。真正决定上限的，是每个 SM 上的寄存器、共享内存和调度带宽。我最终停留在这个配置，不是因为不能把 tile 调大，而是再调大就撞到寄存器墙，性价比不高了。

6. 关于 Warp Tiling 的延续方向

如果你也想把这套思路继续往前推，有几个方向是可以再走的。一是把每线程 tile 从 8×8 改成 16×8 甚至 16×16，进一步加大寄存器复用，但同时对共享内存 load 的向量化要求会更高。二是引入更复杂的 swizzle 模式，在共享内存索引上做 XOR 操作，让 float4 向量化 load 时也能规避 bank conflict，这比 padding 更省空间。三是切分 K 维度，把一个大矩阵按 K 方向拆成几段，交给不同线程块并行计算后再做累加，这在 K 特别大时很有用，也是 cuBLAS 的常用策略之一。

还有一个方向是尝试使用张量核心路径，比如 TF32 或者 FP16 累加。不过那就是另一个世界了，指令、布局和精度策略都不一样。这篇实录里讲的 83%，是在纯 FP32 乘加指令下得到的，如果你把张量核心打开，手写 kernel 和官方库的差距会重新拉开，优化目标也会变成如何把数据喂得更快更满。

我现在写矩阵乘法，已经不会因为它只到 cuBLAS 的 83% 而觉得失败。真正值钱的，是你清楚自己的内核慢在哪、快在哪，也知道下一个百分点要去哪个方向找。如果非要说还有一条最值得分享的经验，那就是：先把 kernel 跑对，再开 profiler，最后才去抠指令级优化。顺序反过来，你会在第一天就怀疑人生。