前阵子我把一个SGEMM kernel从八百 GFLOPS 调到了一万零六百 GFLOPS。终端里打印出10623.47 GFLOPS的时候我盯着屏幕停了几秒。这个数字在深度学习大模型遍地走的今天不算炸裂Tensor Core 随便一开就是百 T 级别但它是靠纯手工 SIMT 风格优化、用 float4 向量化、共享内存分块一步步怼上去的。对我来说这比直接调 cuBLAS 爽多了。先交代一下背景我用的是一张 FP32 理论峰值在三四十 TFLOPS 级别的消费级 GPU矩阵规模是 4096×4096数据类型是 FP32。全程没有用 Tensor Core没有用 cuBLAS所有内核都是自己写的 CUDA C一步步从千行 profile 数据里抠性能。这篇文章不是教科书是我把矩阵乘法从 800 GFLOPS 干到 10600 GFLOPS 的完整记录代码能抄坑也给你标出来了。如果你在学 CUDA、或者正在为什么 SGEMM 跑不满发愁看完这篇应该能少走不少弯路。1. 10000 GFLOPS 到底是个什么水平1.1 先算明白 GFLOPS 这个数是怎么来的FLOPS 全称是 Floating-point Operations Per SecondGFLOPS 就是每秒十亿次浮点运算。10000 GFLOPS 等于 10 TFLOPS也就是一秒内完成一万亿次浮点运算。这个单位在显卡参数页上特别常见但很多人没搞清楚它到底怎么测出来的。矩阵乘法的浮点运算量有固定公式一次规模为 M×N×K 的矩阵乘需要2 × M × N × K次浮点运算。为什么是 2因为每算一个输出元素要执行一次乘法和一次加法。拿我的测试规模来说MNK4096总运算量就是2 × 4096 × 4096 × 4096 137,438,953,472 ≈ 137.4 GFLOPS注意单位换算这个数是 GFLOP总运算量不是 GFLOPS每秒速率。速率需要用总运算量除以实际耗时。假如这 137.4 GFLOP 的任务在 1 秒内跑完速率就是 137.4 GFLOPS如果 13.7 毫秒跑完速率就是 137.4 / 0.0137 ≈ 10030 GFLOPS也就是我刚跨过的那条线。我见过不少人直接在 kernel 结束后用clock()测时间这个精度在毫秒级任务上不够看。正确做法是用 CUDA event也就是cudaEventRecord包住 kernel 启动前后再cudaEventElapsedTime拿到精确到微秒的耗时。另外还有个很多人忽略的点测之前必须确认计算结果是对的。NVCC 在开启优化后有可能把没用的计算整段删掉导致测出来的“性能”高得离谱。我自己的习惯是每版 kernel 跑完后做一次残差校验最大误差超过 1e-4 就直接判定优化无效。1.2 矩阵乘法为什么被当成性能标尺这不只是因为矩阵乘法在数值计算里出现频率高更重要的是它的算术强度Arithmetic Intensity极高。算术强度 总运算量 / 总访存量它决定了一个计算任务在计算机里到底是被算力卡住还是被带宽卡住。以 4096³ 的 FP32 矩阵乘为例总数据量大约是三个矩阵加一起约 201 MB总运算量 137.4 GFLOP算下来算术强度在 680 FLOP/Byte 左右。就算 GPU 的全局内存带宽是 1 TB/s要达到 10 TFLOPS 你也只需要大约 14.7 GB/s 的访存流量带宽余量相当充裕。换句话说这是个典型的“计算受限”问题——理论上只要把计算单元喂饱性能就能拉满。正因为它的计算密度极高矩阵乘法成了检验一个并行体系结构“底线效率”的试金石。深度学习里的全连接层、卷积的矩阵展开、Transformer 的 QKV 投影和注意力打分底层本质都是 GEMM。如果你能把 GEMM 写好你在 AI 推理和训练的数据流优化上就成功了一大半。2. 向量化不是玄学是一笔经济账2.1 FMA一次乘加两条流水线一起跑CUDA 内核里最常见的一个计算操作是这种模式acc a[k] * b[k];这个表达式在 CPU 上会被拆成一条乘法和一条加法但在 GPU 上编译器会尽可能把它融合成一条 FFMA 指令。FFMA 的全称是 Fused Multiply-Add它在一个指令周期内完成d a * b c的完整操作。为什么这很重要因为 GPU 的 FP32 吞吐量通常就是按 FMA 来标定的。一张声称“35 TFLOPS”的卡真实含义是它每秒能执行 17.5 万亿次 FFMA因为每次 FMA 包含两次浮点运算所以对外宣传按 2 倍来算。如果你的代码里乘法和加法没有成功融合实际吞吐会直接掉一半。怎么确认融合发生了把编译后的 SASS 拿出来看就行。编译时加-keep选项或者用nvdisasm打开生成的 cubin 文件搜FFMA字样。如果看到大量FFMA而不是单独的FMUL和FADD说明编译器干得不错。我早期一版 kernel 性能一直上不去查了半天发现某些分支路径里编译器没有做 FMA 融合改写代码结构之后立刻好了。2.2 float4 的账本指令少四倍带宽不白费做 CUDA 优化的人一提到向量化第一个想到的就是float4。这玩意本质是 CUDA 内置的一种向量类型内含四个 float内存中连续排列要求 16 字节对齐。为什么它能提速要点有三个。第一是减少指令数量。同一个线程读 4 个 float用标量方式要发 4 条 LDG 指令用float4只要 1 条 LDG.128。指令发射数量少了计算流水线就有更多slot留给 FFMA。第二是充分利用缓存行和访存带宽。GPU 的全局内存访问以 128 字节为一个事务一个 warp 有 32 个线程如果每个线程按 float4 读取正好是 32 × 16 512 字节也就是 4 个整齐的内存事务如果每个线程只读一个 float则需要 32 个线程共同拼凑才能填满 cache line。虽然连续标量访问在现代 GPU 上也能合并但 float4 让合并变得更整齐、更容易命中问题的本质。第三是流水线压力显著降低。在矩阵乘法的内层循环里如果每个线程每次只取一个 float那么为了凑够足够的计算量内存流水线会被频繁塞满请求。换成 float4 之后同样数量的数据只需要四分之一的内存请求次数memory pipeline 的排队延迟随之下降。我习惯把 float4 想象成快递一个地址同时到了四个包裹你非得分四次下楼去拿这不是折腾自己吗一次拿上来后面处理就有余力了。2.3 对齐是个隐形门槛不注意直接崩用float4最烦的一点是它要求 16 字节对齐。malloc和cudaMalloc分配的内存默认对齐到 256 字节一般没问题。但当你用指针偏移去访问一个二维数组的某一行时就需要小心了。比如你有一个float* A行数是 K行主序布局下访问A[row * K col]。如果你想把这个地址当成float4*用就必须保证(row * K col) * 4字节能被 16 整除也就是row * K col是 4 的倍数。如果 K 不是 4 的倍数那恭喜你某些行会天然错位。解决方式一般是两种一是在分配缓冲区时按 K 向上取整到 4 的倍数也就是 padding二是干脆把矩阵拷贝到一个对齐的新缓冲区里。千万别硬着头皮 reinterpret_cast未对齐的 float4 访问在 CUDA 里会直接导致非法地址错误而且这种错误还不容易定位。我后来在代码里养成了一个习惯凡是涉及 float4 的地方先写一个断言assert(reinterpret_castuintptr_t(ptr) % 16 0);别小看这一行断言它帮我抓到了无数潜在的边界 bug。3. 从 800 GFLOPS 到 10000 GFLOPS 的一条可复现路线3.1 基线版本列主序 B 矩阵让访存模式崩了我的出发点是一个“最朴素”的 kernel每个线程算一个输出元素三重循环直接怼没有任何分块。由于我的数据接口里 B 矩阵是列主序存放的也就是B[k col * K]这种索引方式问题一下子就来了。在同一 warp 里线程的col是连续的但访问 B 的时候地址步长是K × 4字节。这意味着 32 个线程访问的地址彼此相距巨大根本无法合并。GPU 的全局内存系统为每个线程单独发起了内存事务一个 warp 的数据请求被拆成了一堆碎片带宽利用率惨不忍睹。这版代码的性能惨到什么程度呢只有大约 800 GFLOPS。虽然这个数字对很多 CPU 来说已经是天文数字但是对 GPU 来说是在暴殄天物。它浪费的不是算力而是内存系统把所有时间都花在搬运没用的数据上了。这版失败的 kernel 给了我很深一个教训性能优化永远先看访存模式再看计算模式。一个错误的数据布局能让千万亿次算力原地趴窝。3.2 第一步优化先改布局别急着上高深技巧面对这个惨淡的基线我做的第一件事不是引入共享内存也不是上 float4而是先把 B 矩阵的存放方式改成行主序。改写之后最内层循环里线程访问 B 的方式变成了B[kk * N col]同一 warp 的 32 个线程地址连续全局内存访问立刻重新合并。就这一下性能从 800 GFLOPS 飙到了约 2600 GFLOPS三倍多的收益代码逻辑几乎没变。这一步的意义在于把问题分成两层数据布局问题和计算组织问题。数据布局不合理的时候后面所有优化都建立在沙子上。就像做饭之前先确认食材是新鲜的否则刀工再好也白搭。顺带说一句如果你没法改调用方传进来的 B 矩阵布局还有一个常见做法是在 kernel 外面用一个cudaMemcpy2D把列主序矩阵转成行主序或者专门写一个转置 kernel。转置本身也有开销但一次转置能换来几十倍的计算加速这笔账永远划算。3.3 第二步优化共享内存 Tiling float4 向量化数据布局理顺以后下一步要解决的是“数据复用”的问题。朴素 kernel 里每个线程直接访问全局内存它算完一个输出就不再复用任何数据。但实际上矩阵乘法的数据重用机会极大A 矩阵的一行可以被多个线程复用B 矩阵的一列也是。所以我引入了共享内存分块把 A 和 B 按小块先搬到共享内存里再由所有线程从共享内存读取计算。共享内存在 GPU 片内延迟比全局内存低一个数量级而且带宽高得多。所谓 tiling就是每次只处理一个 64×64 的输入块把它切成小片搬进共享内存算完再换下一片。同时我把搬数据的过程向量化了一次float4拷贝四个 float减少了加载指令数量。大致的 kernel 骨架如下#define BLOCK_M 64 #define BLOCK_N 64 #define BLOCK_K 16 __global__ void sgemm_tile_v2(const float* A, const float* B, float* C, int M, int N, int K, float alpha, float beta) { __shared__ float As[BLOCK_M][BLOCK_K]; __shared__ float Bs[BLOCK_K][BLOCK_N]; int tx threadIdx.x; // 0..63 int ty threadIdx.y; // 0..63 int rowBase blockIdx.y * BLOCK_M; int colBase blockIdx.x * BLOCK_N; float acc 0.0f; for (int k0 0; k0 K; k0 BLOCK_K) { __syncthreads(); // 向量化加载 A 和 B 的子块到共享内存 float4 aVec *reinterpret_castconst float4*(A[rowBase ty][k0 tx * 4]); float4 bVec *reinterpret_castconst float4*(B[k0 ty * 4][colBase tx]); // 这里需要保证 K 和索引按 4 对齐示意代码省略边界处理 // ... __syncthreads(); #pragma unroll for (int kk 0; kk BLOCK_K; kk) { acc As[ty][kk] * Bs[kk][tx]; } } // 写回 C省略 beta 处理 C[(rowBase ty) * N (colBase tx)] acc; }这个版本跑出来约 5300 GFLOPS。注意我全程用__syncthreads()控制共享内存读写顺序这个 barrier 放错位置会导致数据竞争计算结果直接花掉。我在调试这个版本的时候就吃过亏——一个__syncthreads()漏放结果矩阵右上角总有随机值排查了整整一下午。这一步的收益来自两方面一是全局内存访问次数大幅减少因为数据从全局内存搬进共享内存之后可以被 64×64 个线程反复复用二是 float4 让搬数据本身更高效了。3.4 第三步优化寄存器微块、循环展开与 FMA 指令密度到了 5300 GFLOPS计算瓶颈开始转移。用nsight compute一查issue active已经比较高了说明指令发射速率开始饱和。这时候再用“一个线程算一个元素”的结构就不行了——每个 FFMA 都对应着两次共享内存读取共享内存 load 指令占据了太多发射槽位。解法是寄存器微块。具体来说让每个线程算一个 4×4 或者 8×8 的微块也就是把多个输出元素放在寄存器里。这样每次从共享内存取一个 A 元素和一个 B 元素可以同时更新多个累加器大幅提高数据复用率。比如一个线程算 4×4 的微块它从共享内存读一个 A 向量和 B 向量到寄存器就能更新 16 个累加器。共享内存的 load 指令数量降到原来的八分之一FFMA 的密度也上来了。float acc[4][4] {}; float aReg[4], bReg[4]; #pragma unroll for (int kk 0; kk BLOCK_K; kk) { // 从共享内存按 float4 加载 *reinterpret_castfloat4*(aReg) *reinterpret_castconst float4*(As[ty][kk * 4]); *reinterpret_castfloat4*(bReg) *reinterpret_castconst float4*(Bs[kk * 4][tx]); #pragma unroll for (int i 0; i 4; i) { #pragma unroll for (int j 0; j 4; j) { acc[i][j] aReg[i] * bReg[j]; } } }寄存器微块还有一个额外好处多个独立累加器增加了指令级并行度。FFMA 虽然吞吐高但它有大约四五个周期的延迟。如果只有一个累加器每步计算都得等上一步结果回来有 16 个独立累加器的话流水线可以同时塞满多个 FFMA延迟被彻底隐藏。这版性能达到了约 8800 GFLOPS。此时计算和访存的平衡已经很接近上限了剩下的空间更像是在和流水线细节搏斗。3.5 第四步优化双缓冲与异步拷贝跨过万级门槛最后这 1800 GFLOPS 的提升是最费劲的。瓶颈从“计算密度”变成了“等待共享内存填充”。虽然 tiling 已经把全局内存流量压低了很多但是每次计算新的 K 分块时线程还是得停下来等__syncthreads()之后共享内存的数据就绪。解决办法是双缓冲。核心思想是在计算当前分块的同时让下一个分块的数据提前从全局内存搬到共享内存。这样计算和访存重叠把隐藏在同步点后的空闲时间全利用起来。现代 CUDA 架构提供了cp.async指令也就是异步拷贝。在 Ampere 及之后架构上我直接用__pipeline_memcpy_async把全局内存的数据异步搬到共享内存手动控制__pipeline_commit和__pipeline_wait_prior。__shared__ float As[2][BLOCK_M][BLOCK_K]; __shared__ float Bs[2][BLOCK_K][BLOCK_N]; // 预取第0块 __pipeline_memcpy_async(As[0][ty][tx * 4], A[...], 16); __pipeline_commit(); for (int k0 0; k0 K; k0 BLOCK_K) { int cur (k0 / BLOCK_K) 1; int next cur ^ 1; __pipeline_wait_prior(0); __syncthreads(); // 预取下一块 if (k0 BLOCK_K K) { __pipeline_memcpy_async(As[next][ty][tx * 4], A[...], 16); __pipeline_commit(); } // 使用当前块计算 ... __syncthreads(); }这一步跑完成绩实打实地站上了 10600 GFLOPS。相比最初的 800 GFLOPS整体提升了十三倍。而且快跑进 4096³ 单次 13 毫秒以内了这个规模下已经是相当可观的手工优化水准。3.6 关键版本的性能与瓶颈变化记录这是整个调优过程中的各阶段数据每版 kernel 都做了正确性验证误差控制均在 1e-4 以内。版本关键优化点GFLOPS相对基线提升主要瓶颈v0 基线无列主序 B 直接硬算8001x非合并访存带宽浪费v1 布局调整B 转行主序消除非合并访存26003.25x数据复用不足计算密度低v2 共享内存分块Tiling float4 向量化搬数据53006.6x共享内存 load 指令过多v3 寄存器微块每线程 4×4 累加器FMA 指令密度提升880011x同步等待访存未与计算重叠v4 双缓冲异步拷贝cp.async 预取下一分块1060013.3x逼近纯手工 SIMT 优化天花板NSight Compute 的数据也能看出瓶颈的变迁v0 是stall_long_scoreboard等待全局内存回包占绝对主导v2 变成了stall_short_scoreboard等待共享内存v3 是stall_barrier同步等待到 v4各种 stall 占比变得非常均匀说明已经很难再从简单手段里榨出性能了。4. CUDA 环境与向量化实践的高频坑4.1 换 CUDA 版本前先想清楚三件事做 CUDA 优化绕不开环境配置。我见过太多人一上来就装最新版 CUDA结果项目编不过、运行报错最后发现是版本匹配问题。决定能不能用某个 CUDA Toolkit 版本核心有三件事显卡的计算能力Compute Capability、驱动支持的 CUDA 版本、Toolkit 自身的版本号。先说驱动。nvidia-smi右上角会显示一行CUDA Version: 13.0之类的信息这代表当前驱动能支持的最高 CUDA 运行时版本。驱动是向后兼容的只要你的 Toolkit 版本不高于这个数字通常就能跑。反过来驱动太老、Toolkit 太新运行时就会报找不到 CUDA 入口的错误。再说计算能力。每个 NVIDIA GPU 都有一个算力代号比如 sm_863060Ti/3080、sm_894060/4070、sm_1205070。编译 kernel 时-archsm_89里的数字必须能被 Toolkit 识别。如果你的 GPU 是 5070sm_120而 CUDA Toolkit 还是 12.0那么即使驱动支持编译也会报Unsupported gpu architecture compute_120。这种情况不是你的代码有问题是 Toolkit 太老必须升级到能识别 sm_120 的版本。最后是 Toolkit 到底装哪个版本。我的建议是别追新先看你的深度学习框架、第三方库的要求。PyTorch 官方轮子自带 CUDA runtime对系统级的 Toolkit 其实依赖不强但如果你要自己编译 cuDNN 相关代码就需要配套版本。4.2 nvcc -V 和 nvidia-smi 显示不一致不是 Bug这是个极其常见的迷惑点。nvcc -V显示的是你安装的 CUDA Toolkit 版本比如 11.8而nvidia-smi右上角显示的是驱动支持的最高 CUDA 版本比如 13.4。两个数字不一样太正常了。如果你没在nvcc在 PATH 里配置好更会出现这种情况编译时用的是系统自带的旧版 Toolkit运行时驱动却能支持新版本。我见过有同学因此把驱动卸了重装折腾了半天结果根本不是问题。确认编译环境时用which nvcc和nvcc -V确认运行环境时用nvidia-smi两条链别混在一起。4.3 几个高频报错的排查思路我把平时遇到的典型报错整理成了一个速查表每一项都是自己或者身边的人踩过的坑报错特征常见原因处理办法nvcc fatal: Unsupported gpu architecture compute_120Toolkit 版本过老不认识新算力升级到支持对应 sm_xx 的 Toolkit 版本/usr/local/cuda/include/cuda_runtime.h: No such fileCUDA_HOME没设置或头文件路径不对检查CUDA_HOME是否指向实际安装目录比如/usr/local/cuda-12.4libcudart.so.12: cannot open shared object file运行时找不到 CUDA 动态库把/usr/local/cuda/lib64加进LD_LIBRARY_PATH并跑一次ldconfigPyTorch 启动时UserWarning: GPU is not compatible with the PyTorch versionPyTorch 自带的 CUDA runtime 不支持新显卡算力升级 PyTorch 到更新的 CUDA 编译轮子而不是去重装系统 CUDACUDA_MALLOC一直失败或者 allocation 为 0显存被占满或环境变量CUDA_CACHE_DISABLE被设置查nvidia-smi显存占用清掉僵尸进程确认缓存目录可写编译gpu_burn这类项目时报错指向/usr/local/cuda/include/cuda找不到编译脚本里的 CUDA 路径和你实际安装路径不一致设置CUDA_HOME或修改 Makefile 里的 CUDA 路径变量OpenCV 装了opencv-python但cv2.cuda模块不存在官方预编译 wheel 默认不启用 CUDA 支持需要自己从源码构建带 CUDA 的 OpenCV或使用第三方预编译包并核对版本搭配另外提醒一个反直觉的坑很多人以为装了 CUDA Toolkit 就等于装了显卡驱动其实 Toolkit 只是编译器加运行时库不包含驱动。驱动必须单独从 NVIDIA 官网装而且先装驱动再装 Toolkit 是更稳妥的顺序。5. 冲过万级 GFLOPS 之后的几点心得如果你也想冲 10000 GFLOPS 这条线我给的最实际建议是别一上来就抄 CUTLASS 或者 cuBLAS 的源码。那些库为了适配所有架构做了大量模板抽象直接读会一头雾水。先用手写 kernel 把一个简单的 SGEMM 优化链路完整走一遍再回头看那些工业级代码你会发现处处都能对上号。我在这轮优化里最大的体会是性能瓶颈会随着优化不断移动。最初是访存模式然后是数据复用接着是指令密度最后是流水线等待。每个阶段都要用 profiler 重新确认瓶颈在哪不能靠猜。NSight Compute 里的 stall 原因分布是非常可靠的诊断指标比肉眼猜性能强太多了。还有一个小技巧在做性能对比时每个优化只改一个变量。我 v1 到 v2 的版本只加了共享内存分块v2 到 v3 只加了寄存器微块v3 到 v4 只引入了异步拷贝。如果一次性把所有优化都堆上去出了问题你根本不知道是哪一步导致的。控制变量法在进行性能优化时同样成立。这次跑过 10600 GFLOPS 之后我又顺手试了一下把分块尺寸从 64×64 调成 128×128结果不升反降因为共享内存占用太大导致 occupancy 掉得厉害。想冲更高性能下一步还得从 Tensor Core 和更精细的 warp 调度入手。但那是另一条路的故事了。
