德国工程师 Simon Boehm 用六周时间，从零开始优化 CUDA 矩阵乘法内核，通过内存合并、共享内存分块、寄存器复用、向量化、自动调优和 Warp 分块等技术，将性能从 cuBLAS 的 1.3% 提升至 93.7%，揭示了 GPU 高性能计算的核心原理。
从零手撕 cuBLAS！一个程序员的 CUDA 矩阵乘法极限优化之旅
为什么大模型训练动不动就要几千张 A100？为什么英伟达的 cuBLAS 库被奉为“神级加速器”？其实，背后最核心的操作，就是矩阵乘法（GEMM）。
今天我们要聊的，不是什么高大上的 AI 前沿，而是一个硬核技术宅，用整整六周时间，从一行“天真的” CUDA 内核代码起步，硬生生把性能从 cuBLAS 的 1.3% 拉到 93.7% 的实战记录。
这位叫 Simon Boehm 的德国工程师，在他的个人网站 siboehm.com 上，用一篇“工作日志”风格的文章，毫无保留地展示了自己如何一步步优化 CUDA 矩阵乘法内核。这篇文章没有花里胡哨的术语包装，没有商业化的 KPI 压力，只有一张张手绘示意图、一行行带注释的代码、一次次性能测试数据的对比——它就像是一个老派黑客的工程笔记，既有学术深度，又有极客精神。而我们今天，就把他这段“硬核炼金术”，用最接地气的方式，掰开揉碎讲给你听。
初学者的“天真”写法，性能只有 cuBLAS 的 1.3%
故事的起点，是一个再简单不过的想法：每个 CUDA 线程负责计算结果矩阵 C 中的一个元素。比如，在 4092×4092 的矩阵相乘中，我们就启动 4092×4092 个线程，每个线程做一个点积。代码看起来清爽又直观，如下所示：
<strong>global</strong> void sgemm_naive(int M, int N, int K, float alpha, const float *A, const float *B, float beta, float *C) { const uint x = blockIdx.x * blockDim.x + threadIdx.x; const uint y = blockIdx.y * blockDim.y + threadIdx.y; if (x < M && y < N) { float tmp = 0.0; for (int i = 0; i < K; ++i) { tmp += A[x * K + i] * B[i * N + y]; } C[x * N + y] = alpha * tmp + beta * C[x * N + y]; } } |
但问题来了——这种写法的内存访问模式极其糟糕。GPU 的全局内存（GMEM）是以 128 字节为单位进行读取的，只有当同一线程束（warp，32 个线程）访问连续内存地址时，才能“合并”成一次高效读取，这叫“内存访问合并”（coalescing）。而在这个 naive 内核里，相邻线程读的是 A 的不同行、B 的不同列，地址完全不连续，导致大量带宽浪费。实测结果惨不忍睹：在 RTX A6000 上，性能只有 309 GFLOPs，而 cuBLAS 能跑 23249 GFLOPs，差距近 75 倍！这就像你开着法拉利去菜市场买菜，却只敢挂一档慢悠悠地开。
第一次优化：让内存访问“排好队”，性能暴涨 6 倍
Simon 没有急着换算法，而是先解决最基础的内存瓶颈。他意识到，问题的关键在于线程映射方式。原本用二维线程块（32×32），导致同一 warp 内的线程在内存中“东一榔头西一棒槌”。于是他改成一维线程块（1024 个线程排成一列），并通过整除和取余重新计算每个线程对应的 C 矩阵坐标：
const int x = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE); |
这一改，奇迹发生了：A 矩阵的读取变成连续的了！因为 threadIdx.x 连续变化，x 不变而 i 递增，A[x*K + i] 就是连续地址；B 矩阵虽然还是按列读，但只要 K 足够大，L2 缓存也能缓解一部分压力。更重要的是，全局内存带宽利用率从 15 GB/s 飙升到 110 GB/s！性能直接干到 1986 GFLOPs，是原来的 6.4 倍。这告诉我们：在 GPU 编程里，正确的数据布局比算法本身更重要。你写的不是代码，是内存访问的艺术。
引入共享内存：把热数据“搬上芯片”，性能再翻 1.5 倍
光靠全球内存还不够，因为 GMEM 延迟太高。于是 Simon 祭出了 GPU 的“高速缓存”——共享内存（SMEM）。SMEM 是每个流多处理器（SM）上的一小块片上存储，带宽高达 12 TB/s（是 GMEM 的 16 倍以上！）。他的思路是：每个线程块先从 GMEM 加载一小块 A 和 B 到 SMEM，然后在 SMEM 里做局部矩阵乘，重复滑动这个“块”直到算完整个结果。这就是经典的“分块算法”（tiling）。关键代码如下：
As[threadRow * BLOCKSIZE + threadCol] = A[threadRow * K + threadCol]; |
这样，每个元素被加载一次后，就能被块内多个线程复用，大幅减少 GMEM 访问次数。性能提升到 2980 GFLOPs，虽然增幅不如上一步大，但为后续优化铺平了道路。不过这时新的瓶颈出现了：SMEM 访问本身也开始成为瓶颈，因为每个 FMA（融合乘加）操作都要等两次 SMEM 读取，线程经常“饿着肚子等数据”。
让每个线程多干点活：从 1D 分块到 2D 分块，算力翻倍再翻倍
既然 SMEM 成了新瓶颈，那就减少对它的依赖——让每个线程一次计算多个 C 元素，把中间结果存在寄存器里。这就是“寄存器分块”（register tiling）。在 Kernel 4 中，每个线程算一列 8 个结果；到了 Kernel 5，直接升级成 8×8 的小方块！为什么是方块？因为这样 A 和 B 的元素能被最大程度复用。比如，一个 A 元素可以和 8 个 B 元素相乘，产出 8 个结果。代码结构也变得更复杂：
float threadResults[TM * TN] = {0.0}; |
这一招直接把算术强度（每字节数据能干多少次计算）拉高，让 GPU 从“内存受限”转向“计算受限”。性能从 8474 GFLOPs 飙到 15971 GFLOPs，几乎翻倍！Simon 在文中感慨：优化的本质，就是不断把数据从慢速存储搬到快速存储，最后塞进寄存器——那里才是算力的天堂。
向硬件极限冲锋：向量化 + 自动调优 + Warp 分块
到了 16 TFLOPs，离 cuBLAS 的 23.2 TFLOPs 还有不小差距。Simon 决定榨干最后一滴性能。
他做了三件事：
第一，向量化内存访问。通过 float4 类型一次性读 128 位数据，并在加载 A 时顺便转置，让 SMEM 访问也能向量化，省下 500 GFLOPs 开销。
第二，暴力自动调优（autotuning）。他写了个 Bash 脚本，把分块大小 BM/BN/BK、寄存器块大小 TM/TN 等参数组合全试一遍，找到最适合 A6000 的配置（BM=BN=128, BK=16, TM=TN=8），又榨出 5% 性能。
第三，也是最关键的——Warp 分块（warptiling）。Warp 是 GPU 调度的最小单位（32 线程），但之前的代码没显式利用它。Simon 在寄存器分块和线程块之间插入一层 Warp 级分块，让同一 Warp 内的线程协作计算一个更大的子块。这不仅提升寄存器局部性，还为未来使用 Tensor Core 铺路。
最终，性能达到 21779 GFLOPs，是 cuBLAS 的 93.7%！
为什么 cuBLAS 仍是不可逾越的高峰？
尽管 Simon 的内核在大矩阵上接近 cuBLAS，但在小矩阵（如 256×256）上差距巨大。
原因很简单：cuBLAS 不是一个内核，而是一个内核库！它内置了上百种针对不同尺寸、不同架构优化的 SGEMM 实现。比如在 256 尺寸，它会用“Split-K”策略：把 K 维拆给多个线程块并行计算，最后再用一个归约内核合并结果。
这种策略在小矩阵上能更好利用 GPU 资源。而 Simon 的内核只针对大矩阵优化，没做这么多特化。这也说明：工业级高性能库的背后，是海量的工程细节和硬件知识的沉淀。我们能手撕 93%，但最后那 7%，可能需要英伟达十几年的积累。
硬核优化的尽头，是理解硬件本身
Simon 在文末说：“优化 SGEMM 是理解 GPU 硬件性能特性的最佳方式。”这句话太对了。从内存合并到共享内存分块，从寄存器复用到 Warp 调度，每一步优化都迫使你深入思考：GPU 的内存层级是怎样的？SM 的资源限制是什么？Warp 调度器如何工作？甚至编译器如何优化循环？这种“第一性原理”式的实践，远比调用现成库更有价值。尤其在今天，当所有人都在追 AI 大模型时，反而更需要这种硬核底层能力——因为只有理解硬件，才能真正驾驭它。Simon 虽然没写出 cuBLAS，但他获得的，是对计算本质的深刻洞察。