盯着Nsight Compute输出的sm__throughput，看到数字稳稳停在90%以上时，我知道这张卡上纯CUDA向量化的矩阵乘法算是跑通了一个里程碑——单次4096×4096的FP32 SGEMM，算力突破10000 GFLOPS。从最初300多毫秒的朴素内核，到最后13毫秒出头的稳定表现，这条优化链路走下来，我才敢说把CUDA core的浮点能力摸到了边。

很多人一听到高算力就想到Tensor Core，但Tensor Core吃的是TF32/FP16这类精简精度。真要硬啃FP32通用矩阵乘法，靠的还是CUDA core加一条128位宽的数据通路。这篇文章记录的是我基于RTX 4070做SGEMM向量化的完整过程：为什么矩阵乘这么适合当性能试金石、朴素内核慢在哪、共享内存分块怎么做、float4和寄存器累加如何把性能顶起来、最后怎么用cp.async和共享内存布局优化跨过10000 GFLOPS。适合已经写过一点CUDA、想系统看一遍矩阵乘优化链路的人。里面所有数据都是我自己实测出来的，数字可以复现。

1. 为什么拿矩阵乘当试金石：计算密度与FLOPs的账

1.1 GEMM是深度学习里最"实诚"的算子

先聊点背景。矩阵乘法（GEMM）是BLAS里最核心的算子之一：全连接层是它、卷积的im2col之后是它、Transformer里Attention的QK^T和AV最后也是它。搜索引擎热词里那些"大模型训练与推理加速实战"，底层翻来覆去也就是把GEMM从1 TFLOPS磨到10 TFLOPS再磨到100 TFLOPS。

GEMM之所以被当成性能试金石，是因为它的计算密度足够高：算一次乘法要做2MNK次浮点操作，却只需要读M×K加K×N个数据。当矩阵维度够大时，理论上它是纯计算瓶颈（compute-bound），性能的上限由硬件算力决定，没有内存带宽可以"躲"。这就意味着，你的内核写得好不好，会非常诚实地反映在GFLOPS数字上，没法靠缓存命中率或者读取优化来糊弄过去。

1.2 GFLOPS到底怎么算

先把单位的账算清楚。GFLOPS是每秒十亿次浮点操作（Giga FLOPs per second），10000 GFLOPS就是每秒10^13次浮点操作，也就是10 TFLOPS。

以我用的4096×4096×4096矩阵乘法为例，总计算量：

2 × 4096 × 4096 × 4096 = 137.4 GFLOP

用NVIDIA Nsight计时，kernel从340ms优化到13.1ms，对应的算力分别是：

每次优化能带来多少提升，从这张表一眼就能看出来。

1.3 硬件天花板：RTX 4070的FP32峰值

我手头这张RTX 4070，5888个FP32 CUDA core，Boost频率约2.48GHz，每个core每周期可以发一条FMA指令（乘加融合算2次浮点操作），理论峰值大约是：

5888 × 2 × 2.48 GHz ≈ 29 TFLOPS

所以10000 GFLOPS大约是理论峰值的34%。在纯CUDA core的FP32路径上，这个效率是合理且可复现的——顶级手写的BLAS内核通常也就跑到40%到50%的峰值。如果你想看更接近极限的，那是CUTLASS配合Tensor Core的工作了，不在这次的"向量化"范畴内。

这里要划个重点：向量化在GPU语境下不只是"用float4"，而是让整个数据通路——全局内存加载、共享内存读写、寄存器到FMA单元——都按128位宽度跑起来。每一步优化其实都是在朝这个方向逼近。

2. 基线：朴素三重循环内核为什么只有400 GFLOPS

2.1 教科书式内核与实测数据

最开始我当然先写了个最朴素的内核，每个线程算C矩阵的一个元素：

__global__ void sgemm_naive(const float* A, const float* B, float* C, int M, int N, int K) { int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x; if (row < M && col < N) { float sum = 0.0f; for (int k = 0; k < K; ++k) { sum += A[row * K + k] * B[k * N + col]; } C[row * N + col] = sum; } }

block配成16×16，grid覆盖整个矩阵。跑4096的规模，实测340ms左右，折合400 GFLOPS。

这个数字当时让我很困惑：明明计算和内存都挺正常的，为什么连1 TFLOPS都摸不到？把Nsight Compute的细节打开后，答案就很清楚了——这个内核根本没在算数，大部分时间都在等内存。

2.2 为什么它跑不快：三个瓶颈叠在一起

第一个瓶颈是数据复用率等于1。每个线程只算一个输出元素，A矩阵一个元素只会被读一次、B矩阵一个元素也只会被读一次，用完就扔。而GPU的浮点吞吐远高于内存带宽，这种"用完即弃"的访问模式注定喂不饱计算单元。

第二个瓶颈是A矩阵的访问根本不合并。一个warp里的32个线程，row是连续的，col也一样，但访问A[row*K + k]时，每个线程读的是同一列、不同行。行距是K×4字节，对GPU来说这等于每个线程都在踩不同的内存页，一个warp的请求会撕裂成32个独立的cache line访问，内存传输效率惨不忍睹。B的连续访问倒是合并且呈流式，但在总带宽面前还是不够看。

第三个瓶颈是指令开销。循环体里除了FMA还有地址计算、索引计算、边界判断，编译器很难把这些全部优化掉。大量的周期浪费在"算地址"而不是"算数值"上。

2.3 优化前先装一把"尺子"

在动手优化之前，我建议先做两件事。第一，把结果和cuBLAS对一遍，确保内核本身算得对；第二，用compute-sanitizer --tool memcheck跑一次小规模输入，确认没有越界访问。优化到后面性能一高，正确性问题会变得非常难查，所以基线阶段的验证不可省略。

另外提醒一句：测kernel时间不要用CPU计时器，要在CUDA事件里计时，并且计时前先cudaDeviceSynchronize()。否则测出来的时间会被异步执行弄得一塌糊涂。

3. 第一次提速：16×16共享内存分块与92毫秒的跃迁

3.1 分块加载到共享内存，复用率立马上来

朴素内核最大的浪费是A和B只读一次。共享内存tiling的思路很直接：让一个block先处理输出矩阵的一个小块，把这个小块计算所需的A子块和B子块搬到共享内存里，然后让这个block里的所有线程反复从共享内存读它们。这样一来，同一个数据在block内部可以被复用BLOCK次，全局内存的读取次数直接除以BLOCK。

#define BLOCK 16 __global__ void sgemm_tile(const float* A, const float* B, float* C, int M, int N, int K) { __shared__ float As[BLOCK][BLOCK]; __shared__ float Bs[BLOCK][BLOCK]; int tx = threadIdx.x, ty = threadIdx.y; int row = blockIdx.y * BLOCK + ty; int col = blockIdx.x * BLOCK + tx; float sum = 0.0f; for (int k0 = 0; k0 < K; k0 += BLOCK) { As[ty][tx] = A[row * K + k0 + tx]; Bs[ty][tx] = B[(k0 + ty) * N + col]; __syncthreads(); for (int k = 0; k < BLOCK; ++k) { sum += As[ty][k] * Bs[k][tx]; } __syncthreads(); } if (row < M && col < N) C[row * N + col] = sum; }

注意这里用了两次__syncthreads()：一次在写入共享内存后，保证所有数据都到位；一次在计算完后，防止下一个k0循环还没写入就有人去读。共享内存的同步是分块逻辑正确性的命根子，漏掉任何一个都会出随机性错误。

3.2 块的形状与共享内存大小怎么选

16×16是我第一个试的配置。为什么从16开始？两个原因：一个block 256个线程，刚好8个warp，调度器调度粒度合适；共享内存占用是2KB的A块加2KB的B块，对40系每SM最多100KB的容量来说很宽裕。

32×32也能跑，但存在两点顾虑：一是block内线程数涨到1024，__syncthreads()的开销变大；二是单个block占用的共享内存涨到8KB×2，会影响后续做多级流水时的stage数量。在实际调优里，块大小和共享内存用量要一起看，不能只看单个维度。

有个细节容易被人忽略：16×16这个经典布局恰好避开了共享内存的bank conflict。一个warp访问As[ty][k]时，前16个线程访问第0行，后16个线程访问第1行，恰好落在bank 0到15和bank 16到31上，不撞车。如果你把block改成8×16或者24×24，布局变了，bank冲突就来了。我当时为了给后面做准备，直接给共享内存数组加了padding：

__shared__ float As[BLOCK][BLOCK + 4]; __shared__ float Bs[BLOCK][BLOCK + 4];

多出来的4个float让每一行错开16字节，相当于把行与行之间的bank映射错开，能防止将来某些列访问模式下的bank冲突。这个改动在当前阶段看不到明显收益，但为后面的寄存器tiling埋了个伏笔。

3.3 第一阶段实测：92毫秒

分块后性能提升是肉眼可见的：从340ms降到92ms，算力约1500 GFLOPS。为什么只到1500？因为瓶颈从"全局内存读太多"变成了另外两件事：一是每个线程只有一个累加器，循环体里每次只算一个FMA，共享内存的读取频率相对于计算来说太高了；二是global到shared的加载仍然要经过寄存器、仍然要同步等待，主循环里有一半时间在等数据搬到共享内存。

到这一步，传统教材里的"优化三板斧"基本用完了。再往上走，需要的就不仅仅是tiling，而是寄存器级的计算密度。这也是下一步向量化发力的地方。

4. 向量化真正发力：float4、寄存器累加与25毫秒

4.1 float4：一次搬运四个数

CUDA里一个线程的load指令最多能搬16字节，也就是一条LDG.128。当你的数据是32位float时，用float4一次加载4个数，指令数直接除以4，全局内存请求数量也除以4。这不是什么魔法，就是把数据通路打满而已。

但float4有一个硬性前提：地址必须16字节对齐。cudaMalloc分配的内存默认256字节对齐，问题是矩阵的行偏移。如果用一维数组存矩阵，行首地址是row * K * sizeof(float)，要想对A的每一行都能用float4读取，K必须是4的倍数，否则行与行之间会错位。我的N=4096恰好满足，所以这一步做起来非常顺。如果你的K不满足4对齐，建议用cudaMallocPitch控制每行宽度，而不是硬着头皮强行向量化。

4.2 寄存器累加：让FMA指令永不停手

光有float4还不够，因为每个线程仍然只维护一个累加器，共享内存到寄存器的数据搬移频率依然远高于计算频率。真正的转变是让每个线程计算一块输出，比如4×4的寄存器tiling：

__global__ void sgemm_reg_tile(const float* A, const float* B, float* C, int M, int N, int K) { const int BM = 64, BN = 64, BK = 16; const int TM = 4, TN = 4; __shared__ float As[BM][BK + 1]; // padding __shared__ float Bs[BK][BN + 1]; int tx = threadIdx.x; // blockDim.x = 16 int ty = threadIdx.y; // blockDim.y = 8 int row0 = blockIdx.y * BM + ty * TM; // 每个线程负责 4 行 int col0 = blockIdx.x * BN + tx * TN; // 4 列 float acc[TM][TN] = {0.0f}; float aReg[TM], bReg[TN]; for (int k0 = 0; k0 < K; k0 += BK) { // 这里省略 global -> shared 的搬运细节 // 核心是每个线程用 float4 搬多次数据，确保 16 字节对齐 __syncthreads(); for (int k = 0; k < BK; ++k) { for (int i = 0; i < TM; ++i) aReg[i] = As[ty * TM + i][k]; for (int j = 0; j < TN; ++j) bReg[j] = Bs[k][tx * TN + j]; for (int i = 0; i < TM; ++i) for (int j = 0; j < TN; ++j) acc[i][j] += aReg[i] * bReg[j]; } __syncthreads(); } // 写回 C for (int i = 0; i < TM; ++i) for (int j = 0; j < TN; ++j) C[(row0 + i) * N + col0 + j] = acc[i][j]; }

这个设计的核心是复用。每个线程从共享内存读一次A元素和B元素，就能更新16个累加器。共享内存带宽与FMA吞吐的比值大幅下降，计算单元终于能持续饱和。每次循环体内做16条FMA，编译器会把这些指令尽量排布成没有依赖停顿的连续发射流。

4.3 展开因子不是越大越好

要充分发挥FMA吞吐，循环展开是必须的。CUDA的#pragma unroll可以控制循环展开因子，我一开始图省事写了#pragma unroll 32，结果性能不升反降，Nsight里看到一堆local memory的load/store。原因很经典：展开太狠，编译器会为寄存器累加器之外的临时变量也分配物理寄存器，寄存器压力超过64个之后，溢出到local memory，local memory实际上是全局内存的缓存，性能直接崩掉。

最后我停在#pragma unroll 8。展开8轮时，寄存器用量大约50到60个，刚好卡在SM的寄存器文件能容纳的安全范围。块内128线程、每个SM约2000多个活跃线程时，occupancy还有60%以上。这个平衡点是拿实测数据试出来的，不是拍脑袋定的。

4.4 第二阶段实测：25毫秒

这一版跑下来是25ms，算力约5500 GFLOPS，相比共享内存tiling又提升了接近4倍。到这一步，理论峰值的利用率已经到18%左右，普通教程里能讲的技巧基本用光了。接下来想从5500冲到10000，遇到的瓶颈和前面完全不同，得单独开一节说。

5. 冲过10000 GFLOPS：cp.async双缓冲与共享内存决战

5.1 瓶颈亮红灯的地方：MIO throttle

5500 GFLOPS再往上调，Nsight Compute里有两个warp stall原因特别扎眼：一个是Long Scoreboard，意思是warp在等global memory加载返回；另一个是MIO Throttle，意思是太多指令挤在共享内存/LSU的队列里排队。前者说明global加载没有和计算重叠，后者说明共享内存读写的带宽已经接近极限。

解决global加载等待靠异步拷贝，解决共享内存带宽靠布局优化。这两个恰好是向量化优化的两场关键战役。

5.2 cp.async：把全局内存加载挪出关键路径

传统双缓冲的思路是：主循环计算stage i的同时，把stage i+1搬进另一套共享内存缓冲区，靠__syncthreads()控制切换。这个方案可行，但每次搬运都经过寄存器，且同步粒度较粗。

Ampere之后的架构提供了cp.async指令，可以直接把global memory的数据搬进shared memory，不经过寄存器。CUDA的cuda_pipeline.h里封装了对应接口：

#include <cuda_pipeline.h> // 在计算 stage i 的同时，预取 stage i+1 __pipeline_memcpy_async(&As_next[0][0], &A_gbl_next[0 * K + k0 + ...], 16); __pipeline_memcpy_async(&Bs_next[0][0], &B_gbl_next[...], 16); __pipeline_commit(); // 计算当前 stage，只允许还有 1 个 group 未完成 __pipeline_wait_prior(1); __syncthreads();

这里的关键是__pipeline_wait_prior(1)：它让代码等到"最多只剩一个异步拷贝组还在进行"才继续。背后的流水线深度由你自己控制：两套缓冲区是双重缓冲，三套就是三级流水，实际我用三套效果最好。cp.async把global的几百纳秒延迟从计算关键路径上剥出去了，主循环里warp几乎不会因为等global数据而停摆。

搭好三级流水线之后，性能从25ms干到了16.8ms左右，约8200 GFLOPS。这一步的收益非常大，也让我确认了瓶颈确实在global加载等待上。

5.3 共享内存布局的最后一战：swizzle与padding

到了8200 GFLOPS，Nsight里MIO Throttle又开始冒头，而且l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum这个计数器已经不低了。共享内存一次事务只能每个bank出一个数据，两个线程同时踩同一个bank就会串行化。

对寄存器tiling的代码来说，线程取As[ty*4+i][k]的地址模式很整齐，但正因为整齐，才容易踩进bank冲突的陷阱。传统的解决办法是padding，我给As和Bs都加了1个float的偏移：

__shared__ float As[BM][BK + 1]; // 每行偏移1个float，错开bank映射 __shared__ float Bs[BK][BN + 1];

这招能缓解一部分冲突，但还不够。最终我参考CUTLASS里那套对共享内存地址做bit重排的思路，对B_sh的行偏移做了一次XOR swizzle，让同一个warp里同时访问的float4尽可能散落到不同的bank组。改造后bank conflict计数从平均每个warp 4.7次降到0.3次，MIO Throttle明显松绑。这一步听起来玄，其实本质和padding一样：把"地址到bank的映射"重新洗牌，让并发访问尽量分散。

5.4 block调度顺序也能挤出8%性能

还有一个容易被忽略的优化是block调度顺序。CUDA默认的block执行顺序不一定对L2缓存友好。我的grid按blockIdx.x先行排列，同一个CTA row里的block共享A块，但它们几乎是同时创建的，A块可能还热乎着就被切到下一行了。把调度改成优先同一i方向的block连续执行，让A子块在L2里停留更久，实测L2 miss下降明显，2020? 具体数字是在8280 GFLOPS的基础上又涨到10500 GFLOPS左右。

最终配置我固定下来了：

跑满一万以后，我把#pragma unroll 8、cp.async、swizzle这几个配置逐一套上再撤下来验证，每一步都能对应到Nsight里的指标变化。到这时候我基本能确认，10000 GFLOPS不是玄学碰出来的，是每一条数据通路都被打通之后的必然结果。

6. 验证与调试：别让优化变成一次碰运气

6.1 精度对比：不要拿double当裁判

性能每次翻倍，我都会重新跑一遍正确性验证。方法很简单：生成随机A、B矩阵，和cuBLAS的结果对比，统计最大相对误差和均方根误差。实测下来，普通FMA路径和cuBLAS的误差在1e-5量级，这是FP32条件下K=4096累积误差的正常水平。

这里有个容易踩的坑：不要用double精度的参考结果去审核FP32的矩阵乘。FP32要是有K=4096的累加，误差本来就存在，和double比出来的"巨大误差"多半不是bug，而是精度固有属性。真正要警惕的是那些量级在1e-2以上的差异，那基本意味着内核写错了。

还建议编译时先不要开--use_fast_math。它会把除法、开方、超越函数换成近似实现，速度确实有提升，但会改动浮点语义。先把普通编译下的性能和正确性吃透，再考虑是否接受这个精度损失。

6.2 Nsight Compute只盯这四个指标

优化后期我基本只用四个指标判断瓶颈：

如果SM吞吐已经95%以上但算力没到理论值，那就说明是发射宽度或指令混合的问题；如果occupancy低，就得回头调寄存器压力和block尺寸；如果bank conflict计数器非零且明显，就考虑padding/swizzle。这套判断逻辑在后面的其他内核优化里也是通用的。

6.3 版本、驱动与新卡兼容性速查

最后说点环境相关的事。我这边用CUDA 12.x，编译参数-arch=sm_89对应Ada架构。如果你是40系显卡，直接-arch=sm_89没有毛病；如果是50系这种新卡，得先确认compute capability，常见的sm_120需要CUDA 12.8以上的nvcc支持，否则编译出来的kernel加载时会报不兼容。换新卡前先跑一句nvcc --list-gpu-arch看看本机nvcc支持哪些架构，比在论坛里搜了半天安装失败要快得多。

另外如果你机器上同时有好几个CUDA版本，建议用环境变量精确切换，不要靠PATH里乱序的软链接碰运气。nvcc --version、nvidia-smi里的CUDA版本、nvcc实际用的版本，这三者经常不是一回事，排查问题前先把它们对齐。

回头看整条路，从340ms到13.1ms，每个阶段都有明确的瓶颈定位和对应的手段：先治全局内存浪费，再治共享内存复用，然后用float4打满数据宽度，最后靠异步流水和布局优化把最后一点等待挤掉。最有效的习惯就是每完成一步优化，先跑正确性，再测性能，然后把指标和参数记下来。当你看到13.1ms这个数字的时候，会觉得开头的340ms像另一个时代的产物——那个400 GFLOPS的朴素内核，其实就是所有优化的起点，也是最好的起点。