GEMM 优化是 CUDA 高性能计算里绕不开的一道坎。网上讲 CUDA GEMM 分块Tiling的教程不少但大多停在“能跑”的阶段真正把算力利用率从 1% 一路推到 95% 的路径讲清楚的并不多。这一篇把整个优化过程拆成九重天从最朴素的 Naive Kernel 开始一层层把全局内存复用、Bank Conflict、寄存器累加、Tensor Core 逐层打通并给出可以直接抄的代码和 Nsight Compute 排查思路。适合谁看刚写完 CUDA 入门程序、Profile 一看利用率不到 10% 的开发者以及在推理引擎里想把 GEMM、卷积、FlashAttention 这类算子性能再压一压的工程同学。整个过程偏实战我会把为什么这样做也一并讲清楚。1. 先搞清楚算力利用率为什么只有 1%1.1 算力利用率到底用什么衡量算力利用率的定义很简单内核实际达到的 FLOPs 除以硬件的理论峰值 FLOPs。GEMM 的浮点运算量是 2×M×N×K比如一个 4096×4096×4096 的矩阵乘浮点运算量大约 137 GFLOP。用 CUDA Event 测出内核耗时一除就能得到实测性能再和硬件理论峰值比就是算力利用率。困难的地方在于“理论峰值”在不同硬件、不同精度下差别很大。同样是 A100FP32 CUDA Core 峰值和 FP16 Tensor Core 峰值可以差一个数量级后者在 300T 这个量级。你拿 FP32 的 kernel 去和 Tensor Core 的理论峰值比那利用率当然惨不忍睹。所以做性能分析前先确认自己在拿什么精度、什么硬件单元和谁对比。1% 是什么概念一个 4096³ 的 FP32 GEMM如果硬件理论峰值为 20T那么 1% 利用率意味着实际只跑到 0.2T耗时接近 0.7 秒。这个数字听起来很离谱但当你写一个完全不考虑数据复用的 Naive 内核时结果真的会这么差原因不是计算单元慢而是数据根本喂不进去。1.2 Memory Bound 才是第一道墙GPU 的算力峰值可以想象成一条高速流水线计算单元是传送带内存系统是送料口。料供不上传送带再快也是空转。这种状态叫 Memory Bound。衡量一个计算任务是否受限于访存常用算术强度Arithmetic Intensity这个概念单位是 FLOP/Byte也就是每读入一个字节的数据能做多少次浮点运算。一个理想的 4096³ GEMMA、B、C 三个矩阵合计约 200MB对应 137 GFLOP算术强度大约 680 FLOP/Byte远高于 GPU 的平衡点理论上应该是一个纯计算受限的任务。但 Naive 实现的实际算术强度远没有这么高。每个输出元素都要重新读一整行 A 和一整列 B一个乘加操作对应两次全局内存读取真实访存量比“所有数据只读一遍”放大了上千倍直接跌进 Memory Bound 区域。这就是 1% 利用率的根本原因不是 GPU 不行是你的内核在拿带宽换计算而带宽是 GPU 上比算力更稀缺的资源。2. 九重天跃迁路线图从 Naive 到 Tensor Core在动手写代码前先把九重天的路径放在一张表里后面每一节对应一层跃迁层级关键动作核心收益预期利用率参考趋势第一重Naive 三重循环基线建立正确性参考1% 以下第二重向量化加载与布局调整提升访存带宽利用率2%-5%第三重共享内存 Tiling数据复用降低全局访存10%-20%第四重Bank Conflict 消除解除共享内存瓶颈20%-35%第五重寄存器 Tiling 双缓冲隐藏延迟提高指令级并行35%-50%第六重引入 Tensor Core启用硬件矩阵单元50%-70%第七重MMA PTX ldmatrix指令级精细控制70%-80%第八重Profiler 驱动调参找到当前架构最优配置80%-90%第九重特化、Padding、Split-K稳定冲击峰值90%-95%预期利用率是基于 V100/A100 这类常见数据中心的经验趋势不同架构和精度会有浮动但优化思路完全一致。2.1 第一重Naive GEMM先明确基线先把最简单的版本写出来做正确性参照__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; } }这个版本的问题一眼就能看出来每个线程独立计算一个输出元素线程之间没有任何数据复用。A 的一行被多个线程重复读取B 的一列也被多个线程重复读取当矩阵规模超过 L2 容量后这些重复读取都会变成真实的全局内存流量。内层循环每做一次乘加都要等两次全局内存访问返回计算管线基本处于停摆状态。这段代码的价值在于它是最容易和 CPU 参考实现对齐的基线。后续所有优化版本都要和它做正确性比对如果一开始基线就是错的后面每一步都白搭。2.2 第二重向量化加载与布局调整第一重里每个线程一次读一个 float4 字节。现代 GPU 的全局内存访问以 32 字节为最小事务单位如果线程束里 32 个线程连续读同一个行硬件能把访问合并成少数几次事务效率尚可但 B 的列访问是跨步的合并效果极差。一个很简单的改进是让数据访问尽量连续并借助向量化指令。float4 一次可以加载 16 字节相当于把四次标量加载合并成一次指令数也减少到四分之一。A 的行方向天然连续可以用 float4B 的列方向不连续最简单的办法是先对 B 做一次转置让计算变成 C^T B^T × A^T 的形式这样两个输入矩阵都能顺序访问。实际工程里CUTLASS 等库会直接要求输入按 K-Major 布局本质就是为这个目的服务的。这一层优化通常能把利用率从 1% 拉到几个百分点但还不够。因为数据复用问题没有解决全局内存流量依然巨大。2.3 第三重共享内存 Tiling复用才是核心分块Tiling的核心思想是把大矩阵切成小块先搬运到共享内存让同一个线程块内的线程反复读取这些小块。共享内存的带宽比全局内存高一个数量级延迟也低很多。以 16×16 的块为例线程块内 256 个线程协作先把 A 的一个 16×16 小块和 B 的一个 16×16 小块从全局内存搬运到共享内存然后每个线程只从共享内存里取数计算。内层循环每处理一个 K 维度上的块全局内存只需要读 2×16×16×42KB 数据但计算需要执行 16×16×16 次乘加数据复用倍数等于块大小。#define BLOCK_SIZE 16 __global__ void sgemm_tiled(const float* A, const float* B, float* C, int M, int N, int K) { __shared__ float As[BLOCK_SIZE][BLOCK_SIZE]; __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE]; int row blockIdx.y * BLOCK_SIZE threadIdx.y; int col blockIdx.x * BLOCK_SIZE threadIdx.x; float acc 0.0f; for (int kk 0; kk K; kk BLOCK_SIZE) { As[threadIdx.y][threadIdx.x] A[row * K kk threadIdx.x]; Bs[threadIdx.y][threadIdx.x] B[(kk threadIdx.y) * N col]; __syncthreads(); for (int k 0; k BLOCK_SIZE; k) { acc As[threadIdx.y][k] * Bs[k][threadIdx.x]; } __syncthreads(); } C[row * N col] acc; }这段代码有两个同步点缺一不可。第一个__syncthreads()保证 A、B 的小块都写入共享内存后才能开始计算第二个__syncthreads()保证所有线程都读完了当前块的数据才允许下一次迭代覆盖共享内存。如果漏掉第二个同步会出现某个线程还在读 As、Bs另一个线程已经往里面写新数据的情况结果随机出错。这层优化是从 Naive 到可用的关键一步通常能把利用率拉到 10%-20%。但如果块大小为 16你会发现共享内存访问还存在一个隐藏的性能杀手就是下一层要处理的 Bank Conflict。2.4 第四重Bank Conflict 与共享内存映射共享内存的硬件结构是 32 个 Bank每个 Bank 宽度 4 字节一次可以同时服务 32 个地址。理想情况下线程束内 32 个线程访问的地址应该均匀分布在 32 个不同 Bank 上。如果两个以上的线程命中同一个 Bank硬件就会把这些访问串行化这就是 Bank Conflict。在As[threadIdx.y][k]这样的内层循环里注意访问的是同一行、不同列。线程束内 32 个线程的 threadIdx.y 和 threadIdx.x 是交错的当它们同时访问 As 的不同行时如果共享内存数组按[16][16]存储每行 16 个 float 占 16 个 Bank第 0 行占 Bank 0-15第 1 行又从 Bank 0 开始两个不同行的同列数据会落在同一个 Bank 上产生冲突。最经典的解法是 Padding把As[BLOCK_SIZE][BLOCK_SIZE]改成As[BLOCK_SIZE][BLOCK_SIZE 1]也就是每行多出 1 个 float 的偏移。这样第 0 行占 Bank 0-16第 1 行从 Bank 1 开始线程束内不同线程访问的行就分散到了不同 Bank。__shared__ float As[BLOCK_SIZE][BLOCK_SIZE 1]; __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE 1];一行代码的变化可能带来接近翻倍的效果。Tensor Core 版本里为了同时兼顾 Bank Conflict 和向量化加载还会用到更复杂的 Swizzle 映射但核心逻辑都是让同一时刻线程束访问的地址散列到不同 Bank 上。2.5 第五重寄存器 Tiling、双缓冲与异步拷贝到这里访存模型已经基本健康但距离 95% 还很远。接下来要处理的是指令级并行和延迟隐藏。寄存器 Tiling 是指每个线程计算多个输出元素比如每个线程算 4×4 或 8×8 的小块。这样做有三个好处第一从共享内存读入一个数据后可以与多个累加结果复用第二每个线程的循环体中有多个独立的乘加链指令级并行度提高可以掩盖 FMA 指令的延迟第三所需的线程数变少同样的 SM 资源能容纳更多线程块提高占用率。双缓冲则是让“加载下一块数据”和“计算当前块”重叠。GEMM 内层循环里只有当前块计算完才能加载下一块否则共享内存会被覆盖。如果准备两块共享内存一个用于计算一个用于预加载下一轮数据就能把内存延迟隐藏掉。在 Ampere 及之后的架构上cp.async指令可以直接从全局内存异步拷贝到共享内存不需要过寄存器也不要求立即同步结合双缓冲效果更好。// 伪代码示意cp.async 发起异步拷贝 asm volatile(cp.async.ca.shared.global [%0], [%1], 16; :: r(smem_addr), l(gmem_addr));这一层是纯 CUDA Core 优化里最接近天花板的一步。到了这里一个 FP32 GEMM 通常能跑到理论峰值的一半左右接下来想再往上走就必须换一种思维方式用 Tensor Core。2.6 第六重Tensor Core 登场用 WMMA APITensor Core 是 GPU 上专门做小矩阵乘累加的硬件单元最早从 Volta 架构引入Ampere 上已经非常成熟。它的核心思想是让一条指令完成一次 4×4×4Volta或 16×16×16Ampere 及之后级别的矩阵乘加而不是像普通 CUDA Core 那样一条指令只做一个标量乘加。对开发者最友好的是 WMMAWarp Matrix Multiply-AccumulateAPI它把硬件细节封装起来编程模型是在 warp 粒度上操作 fragment。一个 16×16×16 的计算A 和 B 分别被切分成多个 fragment每个线程持有一部分数据mma_sync会在 warp 内完成整个矩阵乘累加。#include mma.h using namespace nvcuda; wmma::fragmentwmma::matrix_a, 16, 16, 16, __half, wmma::row_major a_frag; wmma::fragmentwmma::matrix_b, 16, 16, 16, __half, wmma::col_major b_frag; wmma::fragmentwmma::accumulator, 16, 16, 16, float c_frag; wmma::load_matrix_sync(a_frag, A_ptr, lda); wmma::load_matrix_sync(b_frag, B_ptr, ldb); wmma::fill_fragment(c_frag, 0.0f); wmma::mma_sync(c_frag, a_frag, b_frag, c_frag); wmma::store_matrix_sync(C_ptr, c_frag, ldc, wmma::mem_row_major);这段代码的底层效果是一次mma_sync把原本需要 16×16×16 次标量 FMA 的工作量变成一条硬件指令指令发射压力骤降数据搬运成为唯一制约因素。WMMA 对新手来说不用关心 fragment 内部到底怎么分布但使用时有几个硬性要求必须以整个 warp 为单位执行fragment 不能随意通过数组下标访问输入必须是半精度或 TF32 格式。这里面最容易踩的坑是忘记把 FP32 数据转成__half以及矩阵地址没有 16 字节对齐。2.7 第七重MMA PTX 与 ldmatrix 的精细控制WMMA API 已经很方便但封装也让一部分性能留在桌上。真正追求极致性能的库比如 CUTLASS会直接使用 PTX 级的指令。核心指令是mma.sync.aligned.m16n8k16。它的语义是一个 warp 内A 矩阵取 16×16、B 矩阵取 16×8K 维度 16累加到 8×8 的 C 矩阵上。A、B 数据不是从共享内存直接读的而是要求先分布在线程寄存器里。所以还需要另一个关键指令ldmatrix它可以把共享内存中的 16-bit 矩阵数据高效加载到线程寄存器的特定布局中。// 示意ldmatrix 加载 4 个 8x8 的 16bit 矩阵到寄存器组 asm volatile( ldmatrix.sync.aligned.m8n8.x4.b16 {%0,%1,%2,%3}, [%4]; : r(r0), r(r1), r(r2), r(r3) : l(smem_addr));为什么要绕开 WMMA 自己控制两方面的原因。第一WMMA 对 fragment 内部布局有固定选择而针对特定矩阵形状和 Swizzle 策略你可能想要不同的寄存器排列第二PTX 层面可以把ldmatrix、mma、cp.async的调度完全重叠起来减少同步等待。这个阶段的优化收益不是来自某一条指令本身而是来自整个数据流水线的精细编排。如果你不想手写 PTX直接读 CUTLASS 的 Cute 库也是一个很好的进阶路径它把数据布局和指令调度抽象成了更通用的描述但理解线下的 PTX 原理会让你更容易看懂它为什么那样设计。2.8 第八重调参与 Profiler 驱动优化当内核已经用上 Tensor Core剩下的事情就是调参。GEMM 的调参空间其实不大核心变量是这几个线程块维度BLOCK_M、BLOCK_N、K 维度块大小BLOCK_K、流水线级数NUM_STAGES、线程数、Swizzle 策略。组合起来可能上千种但真正有效的配置没几个。Nsight Compute 的 Speed of Light 是这一步最重要的工具。它会把内核运行时的计算管线和内存管线利用率同时显示出来。如果 Compute 利用率高而 Memory 利用率低说明计算正在吃满方向正确如果 Memory 利用率接近 100% 而 Compute 只有 50%说明数据搬运跟不上需要加大 Tiling 尺寸或增加流水线级数如果两者都不高可能是延迟没隐藏好需要提高占用率。以 A100 上的 TF32 GEMM 为例900 多个候选中能稳定跑到理论峰值 90% 以上的也就十几种配置而这十几种配置的共同特点是共享内存足够容纳 4-8 个 Stage 的流水线寄存器数量没有溢出Swizzle 没有引入额外的 Bank Conflict。调参不是随机试每一步都要用 Profiler 数据说明白到底在消除哪个瓶颈。2.9 第九重特化与边界处理最后一步是从“通用内核”走向“场景特化”。当 M、N、K 不是 Tile 的整数倍时要么在加载和存储时加边界判断要么在 Host 端 Padding 到对齐尺寸。边界判断会引入分支增加指令开销Padding 则会占用额外的显存和带宽。对大 GEMM 来说Padding 通常是更好的选择因为浪费很小而分支带来的性能损失是实时发生的。另一个维度的特化是 K 维切割。当 K 非常大时可以把 K 分成多段让不同线程块分别计算结果再用原子加或二阶归约汇总。这个方法叫 Split-K能显著提高小矩阵 GEMM 的并行度。更进一步的 Stream-K 则能精确控制每个线程块的工作量减少负载不均。到了这一层一个针对特定 shape 精心特化的 Tensor Core GEMM在 A100/H100 这类硬件上做到 95% 以上是很现实的。但也要清醒地认识到没有任何一个内核能同时吃遍所有 shape。95% 不是你写完一个通用 kernel 就有的而是针对你的实际生产 shape 反复裁剪出来的。3. 核心代码实现Tiled GEMM 怎么落地3.1 一个可直接编译的共享内存 Tiling 例程把第二重到第四重的思路合并可以得到一个可直接运行的教学版本。这段代码假设 M、N、K 都能被 16 整除方便你观察核心逻辑#include cuda_runtime.h #include cstdio #define BLOCK_SIZE 16 __global__ void sgemm_tiled_padded(const float* A, const float* B, float* C, int M, int N, int K) { __shared__ float As[BLOCK_SIZE][BLOCK_SIZE 1]; __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE 1]; int row blockIdx.y * BLOCK_SIZE threadIdx.y; int col blockIdx.x * BLOCK_SIZE threadIdx.x; float acc 0.0f; for (int kk 0; kk K; kk BLOCK_SIZE) { As[threadIdx.y][threadIdx.x] A[row * K kk threadIdx.x]; Bs[threadIdx.y][threadIdx.x] B[(kk threadIdx.y) * N col]; __syncthreads(); #pragma unroll for (int k 0; k BLOCK_SIZE; k) { acc As[threadIdx.y][k] * Bs[k][threadIdx.x]; } __syncthreads(); } C[row * N col] acc; } int main() { const int M 1024, N 1024, K 1024; size_t bytes M * K * sizeof(float); float* h_A (float*)malloc(bytes); float* h_B (float*)malloc(bytes); float* h_C (float*)malloc(M * N * sizeof(float)); for (int i 0; i M * K; i) h_A[i] 1.0f * (rand() % 10) / 10.0f; for (int i 0; i K * N; i) h_B[i] 1.0f * (rand() % 10) / 10.0f; float *d_A, *d_B, *d_C; cudaMalloc(d_A, bytes); cudaMalloc(d_B, bytes); cudaMalloc(d_C, M * N * sizeof(float)); cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice); cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice); dim3 block(BLOCK_SIZE, BLOCK_SIZE); dim3 grid(N / BLOCK_SIZE, M / BLOCK_SIZE); sgemm_tiled_paddedgrid, block(d_A, d_B, d_C, M, N, K); cudaDeviceSynchronize(); cudaMemcpy(h_C, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost); // 简单的正确性抽查 bool ok true; for (int i 0; i M; i) { for (int j 0; j N; j) { float ref 0.0f; for (int k 0; k K; k) ref h_A[i * K k] * h_B[k * N j]; if (fabs(h_C[i * N j] - ref) 1e-2f) { ok false; break; } } } printf(%s\n, ok ? PASS : FAIL); cudaFree(d_A); cudaFree(d_B); cudaFree(d_C); free(h_A); free(h_B); free(h_C); return 0; }注意这段代码里的共享内存行宽是BLOCK_SIZE 1很多人会奇怪为什么要浪费这 1 个 float其实就是为了消除 Bank Conflict上一节已经解释过。如果你把它改回BLOCK_SIZE功能完全一样但 Profiler 里会看到明显变高的 shared 冲突计数。3.2 共享内存容量与占用率手算选 Tile 大小不能只看算力利用率还要看硬件资源够不够。以 16×16 的块为例每个线程块需要两块 16×16×4B1KB 的共享内存总共 2KB。如果启用双缓冲就是 4KB四缓冲就是 8KB。A100 的一个 SM 有 228KB 共享内存但每 SM 最大线程数 2048。一个 16×16 块有 256 线程8 个块就能打满线程上限共享内存只用了 16KB很多空闲。如果改用 32×32 的块每个块 1024 线程2 个块就打满线程上限但共享内存需要 2×32×32×4B×216KB占用率也到 100%。你可能会想既然线程数先打满共享内存还有很多剩余那为什么不把双缓冲改成四缓冲完全正确多 Stage 流水线的代价就是共享内存线性增长8 个 Stage 的 32×32 块需要 8×8KB64KB这时共享内存先成为瓶颈并发块数量下降延迟隐藏能力反而可能变差。调参就是在“数据复用率”和“并发度”之间找平衡点。3.3 正确性验证不要用等号比较浮点数优化 GEMM 过程中最容易忽略的是正确性验证。线程块大小、累加顺序、向量化方式都会改变浮点运算的中间舍入顺序所以结果和 CPU 三重循环参考实现不可能逐位相等。验证时用相对误差更合理对于 FP32 计算通常允许 1e-4 到 1e-5 级别的差异如果使用了 TF32绝对值误差会放大到 1e-2 级别FP16 输入累加 FP32 大约在 1e-2 到 1e-3 之间。最好把验证阈值和精度绑定在一起不要拿同一个阈值测所有版本。性能测试时也有一个常见坑如果内核之后跟着cudaMemcpy而你没有做cudaDeviceSynchronize计时器的结束时间可能落在拷贝完成的前面得到的时间完全是乱的。正确做法是用 CUDA Event 包住内核同时保证两次 Event 之间的所有操作都完成。4. Tensor Core 思维从数据布局到指令发射4.1 为什么 Tensor Core 对数据布局这么苛刻Tensor Core 和普通 CUDA Core 最大的区别是它做乘法时数据不是从寄存器里随意拿两个数而是按照固定的 16×16 小块从一组线程的寄存器中批量取数。这就好比一个拼板机原材料必须在送进来之前就被裁剪成它规定的尺寸和形状。用 WMMA API 时load_matrix_sync会自动帮你完成从全局内存到 fragment 的布局转换但这个转换本身要消耗指令和寄存器。用 PTX 指令时你需要自己通过ldmatrix把共享内存中的数据按特定布局装进寄存器所以共享内存里的存储方式也要预先设计好这就是 Swizzle 存在的原因。一个成熟的 Tensor Core GEMM 内核整体数据流是全局内存 → 共享内存通过cp.async或 TMA→ 寄存器通过ldmatrix→ Tensor Core 执行 MMA。你会发现真正的 GEMM 优化早已不是“怎么算乘法”而是“怎么搬运数据”。4.2 一个完整的 WMMA 分块循环模式用 WMMA 写 GEMM 时建议的循环结构和共享内存 Tiling 类似只是内层计算从标量乘加变成了mma_sync。下面是一个片段式的模式#include mma.h using namespace nvcuda; // 每个线程块负责 64x64 的 C 块 // 每个 warp 负责 32x32 的 C 子块分 4 个 16x16 累加器 __global__ void sgemm_wmma(const __half* A, const __half* B, float* C, int M, int N, int K) { __shared__ __half As[64][64]; __shared__ __half Bs[64][64]; wmma::fragmentwmma::matrix_a, 16, 16, 16, __half, wmma::row_major a_frag; wmma::fragmentwmma::matrix_b, 16, 16, 16, __half, wmma::col_major b_frag; wmma::fragmentwmma::accumulator, 16, 16, 16, float c_frag[4]; int base_row blockIdx.y * 64; int base_col blockIdx.x * 64; for (int i 0; i 4; i) wmma::fill_fragment(c_frag[i], 0.0f); for (int kk 0; kk K; kk 64) { // 搬数据到 As/Bs... __syncthreads(); for (int kk_inner 0; kk_inner 64; kk_inner 16) { // 每个 warp 从 As/Bs 加载 16x16 的 fragment // 执行 4 次 mma分别对应 4 个输出 16x16 子块 __syncthreads(); } } // store c_frag[0..3] 到 C 的对应位置 }这段代码略去了坐标计算和边界判断但对理解 WMMA 分块模式足够了。实现时最容易出错的地方是 A、B 的矩阵布局。WMMA 的matrix_a和matrix_b有行主序和列主序之分一旦声明错了布局load_matrix_sync虽然不会报错但结果会是完全错误的值而且这种错误非常难排查。4.3 CUTLASS 给我们的启示提到 Tensor Core GEMM绕不开 CUTLASS。它的核心抽象是把 GEMM 的复制Copy、矩阵累加Mma、尾处理Epilogue彻底解耦每种操作都能独立选择和调优。Hopper 架构上CUTLASS 甚至可以借助 TMATensor Memory Accelerator和wgmma指令让数据从全局内存直接流进 Tensor Core中间不需要经过寄存器搬运。对普通开发者来说不一定要全部消化 CUTLASS但你可以观察它内部用的 Tile 形状、Swizzle 函数、Stage 数量再回到自己的内核里实验这些参数。很多时候把 CUTLASS 某个 SASS 级配置抄下来适配到自己的场景比从零瞎猜高效得多。5. 实测与调参怎么确认自己到了第几重5.1 Nsight Compute 的基本使用不要靠感觉优化用 Profiler 数据说话。最常用的是 Nsight Compute命令行就可以跑ncu --set full ./your_gemm或者只看关键摘要ncu --section SpeedOfLight ./your_gemmSpeed of Light 页面会给出两个核心百分比Compute Throughput 和 Memory Throughput。如果 Compute 高、Memory 低说明内核正在做正事如果 Memory 高、Compute 低说明它是 Memory Bound两者都不高通常意味着延迟隐藏不够需要提高并行度。其他关键指标还包括dram__bytes_read.sum看全局内存读取量shared_load/store看共享内存吞吐smsp__sass_inst_executed_op_shared_ld附近的计数器看是否有 Bank Conflict。在 Nsight Compute 的 Source 视图里甚至可以直接看到是源码的哪一行贡献了最多的 stall。5.2 从 1% 到 95% 的实测决策表我把每一次 Profiler 诊断对应的优化动作整理成表实战时直接对照诊断结果主要瓶颈下一层动作DRAM Throughput 高Compute 低全局内存访问太多做共享内存 TilingL1/TEX 吞吐高DRAM 不高数据复用不足或调度开销大加大 Tile / 寄存器 Tilingshared 指令存在 Bank Conflict共享内存访问模式冲突Padding / Swizzle计算管线活跃但指令槽浪费标量 FMA 指令过多改用 Tensor CoreTensor Core 活跃但等待数据搬运和计算未重叠双缓冲 / 多 Stage / cp.async整体利用率接近但差几个点边界分支 / 负载不均特化 shape / Padding / Split-K这套决策表基本就是九重天跃迁的操作手册。每一步都对应一个明确的诊断信号不是玄学调参。5.3 环境准备CUDA Toolkit、驱动与架构兼容很多人内核写好了性能也优化了结果一换机器就报 Illegal Address 或者内核启动失败最后发现是 Toolkit 版本和 GPU 架构不匹配。比如新一代显卡的算力标识是 SM_120旧版 CUDA Toolkit 编译出来的 cubin 根本没法加载运行时就会提示“no kernel image is available”。遇到这种情况先做三个检查nvcc --version查编译工具链nvidia-smi查驱动支持的 CUDA 版本ncu --version查 Profiler 版本。驱动和 Toolkit 不一定必须完全一致驱动版本要大于等于 Toolkit 要求的最低版本而内核编译时指定的-archsm_xx必须是目标 GPU 能支持的架构。写 PTX 或内联汇编时还要确认目标架构是否支持相关指令比如cp.async需要 Ampere 及以上wgmma需要 Hopper 及以上。6. 常见问题与排查技巧实录6.1 矩阵维度不是 Tile 整数倍怎么办最省事的方法是 Host 端 Padding把 M、N、K 向上取整到 16 或 64 的倍数多余区域填 0。这样做不会影响结果因为填充区域的乘加不影响有效输出。如果不方便 Padding可以给加载和存储加 Predicate 判断。比如共享内存加载时如果全局坐标超出 M 或 K 范围就写入 0存储到 C 时如果坐标超出范围就跳过。这种做法灵活但会在 warp 内引入分支分支收敛状态不同步会导致部分线程空转性能会打折扣。优先考虑 Padding。6.2 结果不对但看起来每个值都“差不多”如果输出矩阵的每个值都和参考实现对不上差的幅度不大不小先怀疑布局问题。WMMA 里把 row_major 和 col_major 写反、PTX 里 A/B 操作数顺序颠倒、ldmatrix 的地址不是 16 字节对齐都可能产生这种情况。这类错误用断点很难查最快的方法是写一个小矩阵测试比如 16×16×16把 fragment 里的值逐个 dump 出来和 CPU 参考对照哪个 EM 阶段出了问题一下就能看出来。另一个技巧是用compute-sanitizer跑一遍它能捕捉到未初始化值加载、越界访问这类隐藏错误。6.3 明白共享内存 Bank Conflict 但不知道哪些指令触发的Nsight Compute 的--set full里可以直接看l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum这类计数器。不为 0 就说明有冲突。再结合 Source 视图能精确定位到某条LDS指令。修复时先试 Padding如果 Padding 影响了 ldmatrix 需要 16 字节对齐那就要做一个完整的 Swizzle 映射把共享内存地址重新编码让线程束内访问的地址在 Bank 上错开同时保证每个 16 字节向量仍在连续的 4 个 Bank 内。这是 Tensor Core GEMM 里最考验功力的部分。6.4 内核工作一会儿就崩溃或挂起如果内核在大矩阵上崩溃先怀疑索引越界。Naive 和 Tiling 版本里最常见的 bug 是row或col超出矩阵范围读取了非法地址。解决办法是先把计算结果清零程序还能跑也方便对比。挂起则要优先检查__syncthreads()是否在条件分支里。如果线程块内某些线程进入了同步点另一些没有就会造成死锁。注意__syncthreads()必须让块内所有线程都能到达不能被if (tid 32)之类的分支隔离。6.5 性能始终上不去但 Profiler 一切正常这种情况我遇到过不止一次最后的原因往往很基础编译时没加足够高的优化选项。没有-O3的话编译器可能不会自动展开循环、向量化访问也无法生成最高效的 SASS。另一个原因是 host 端没有把输入布局转成内核期望的 K-Major 布局DataSource 不对GPU 端再优化也只是把错的数据搬来搬去。排查的时候先确认三件事编译选项带没带-O3输入矩阵是不是行主序且字节对齐cudaMalloc的地址是否满足 16 字节对齐要求。很多时候性能问题不是算法问题而是工程细节问题。这一路从 1% 到 95%本质上是不断追问“数据从哪来、到哪去、在哪个环节等了多久”。Naive 版本慢是因为数据从全局内存反复搬运Tiling 版本快是因为把数据复用在了共享内存Tensor Core 快是因为把计算从标量指令换成了矩阵指令。实际操作中我习惯每完成一层优化就把对应的 Profiler 截图和参数记录保存下来这样不仅知道当前卡在哪还能在后续改动中快速 revert 到某个稳定版本。如果你的内核已经能稳定跑到 90% 以上下一步可以试试针对具体 shape 做特化或者深入研究 CUTLASS 里面的 cute 抽象那又是一片新天地。 CUDA GEMM优化九重天：从1%到95%算力利用率实战指南
GEMM 优化是 CUDA 高性能计算里绕不开的一道坎。网上讲 CUDA GEMM 分块Tiling的教程不少但大多停在“能跑”的阶段真正把算力利用率从 1% 一路推到 95% 的路径讲清楚的并不多。这一篇把整个优化过程拆成九重天从最朴素的 Naive Kernel 开始一层层把全局内存复用、Bank Conflict、寄存器累加、Tensor Core 逐层打通并给出可以直接抄的代码和 Nsight Compute 排查思路。适合谁看刚写完 CUDA 入门程序、Profile 一看利用率不到 10% 的开发者以及在推理引擎里想把 GEMM、卷积、FlashAttention 这类算子性能再压一压的工程同学。整个过程偏实战我会把为什么这样做也一并讲清楚。1. 先搞清楚算力利用率为什么只有 1%1.1 算力利用率到底用什么衡量算力利用率的定义很简单内核实际达到的 FLOPs 除以硬件的理论峰值 FLOPs。GEMM 的浮点运算量是 2×M×N×K比如一个 4096×4096×4096 的矩阵乘浮点运算量大约 137 GFLOP。用 CUDA Event 测出内核耗时一除就能得到实测性能再和硬件理论峰值比就是算力利用率。困难的地方在于“理论峰值”在不同硬件、不同精度下差别很大。同样是 A100FP32 CUDA Core 峰值和 FP16 Tensor Core 峰值可以差一个数量级后者在 300T 这个量级。你拿 FP32 的 kernel 去和 Tensor Core 的理论峰值比那利用率当然惨不忍睹。所以做性能分析前先确认自己在拿什么精度、什么硬件单元和谁对比。1% 是什么概念一个 4096³ 的 FP32 GEMM如果硬件理论峰值为 20T那么 1% 利用率意味着实际只跑到 0.2T耗时接近 0.7 秒。这个数字听起来很离谱但当你写一个完全不考虑数据复用的 Naive 内核时结果真的会这么差原因不是计算单元慢而是数据根本喂不进去。1.2 Memory Bound 才是第一道墙GPU 的算力峰值可以想象成一条高速流水线计算单元是传送带内存系统是送料口。料供不上传送带再快也是空转。这种状态叫 Memory Bound。衡量一个计算任务是否受限于访存常用算术强度Arithmetic Intensity这个概念单位是 FLOP/Byte也就是每读入一个字节的数据能做多少次浮点运算。一个理想的 4096³ GEMMA、B、C 三个矩阵合计约 200MB对应 137 GFLOP算术强度大约 680 FLOP/Byte远高于 GPU 的平衡点理论上应该是一个纯计算受限的任务。但 Naive 实现的实际算术强度远没有这么高。每个输出元素都要重新读一整行 A 和一整列 B一个乘加操作对应两次全局内存读取真实访存量比“所有数据只读一遍”放大了上千倍直接跌进 Memory Bound 区域。这就是 1% 利用率的根本原因不是 GPU 不行是你的内核在拿带宽换计算而带宽是 GPU 上比算力更稀缺的资源。2. 九重天跃迁路线图从 Naive 到 Tensor Core在动手写代码前先把九重天的路径放在一张表里后面每一节对应一层跃迁层级关键动作核心收益预期利用率参考趋势第一重Naive 三重循环基线建立正确性参考1% 以下第二重向量化加载与布局调整提升访存带宽利用率2%-5%第三重共享内存 Tiling数据复用降低全局访存10%-20%第四重Bank Conflict 消除解除共享内存瓶颈20%-35%第五重寄存器 Tiling 双缓冲隐藏延迟提高指令级并行35%-50%第六重引入 Tensor Core启用硬件矩阵单元50%-70%第七重MMA PTX ldmatrix指令级精细控制70%-80%第八重Profiler 驱动调参找到当前架构最优配置80%-90%第九重特化、Padding、Split-K稳定冲击峰值90%-95%预期利用率是基于 V100/A100 这类常见数据中心的经验趋势不同架构和精度会有浮动但优化思路完全一致。2.1 第一重Naive GEMM先明确基线先把最简单的版本写出来做正确性参照__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; } }这个版本的问题一眼就能看出来每个线程独立计算一个输出元素线程之间没有任何数据复用。A 的一行被多个线程重复读取B 的一列也被多个线程重复读取当矩阵规模超过 L2 容量后这些重复读取都会变成真实的全局内存流量。内层循环每做一次乘加都要等两次全局内存访问返回计算管线基本处于停摆状态。这段代码的价值在于它是最容易和 CPU 参考实现对齐的基线。后续所有优化版本都要和它做正确性比对如果一开始基线就是错的后面每一步都白搭。2.2 第二重向量化加载与布局调整第一重里每个线程一次读一个 float4 字节。现代 GPU 的全局内存访问以 32 字节为最小事务单位如果线程束里 32 个线程连续读同一个行硬件能把访问合并成少数几次事务效率尚可但 B 的列访问是跨步的合并效果极差。一个很简单的改进是让数据访问尽量连续并借助向量化指令。float4 一次可以加载 16 字节相当于把四次标量加载合并成一次指令数也减少到四分之一。A 的行方向天然连续可以用 float4B 的列方向不连续最简单的办法是先对 B 做一次转置让计算变成 C^T B^T × A^T 的形式这样两个输入矩阵都能顺序访问。实际工程里CUTLASS 等库会直接要求输入按 K-Major 布局本质就是为这个目的服务的。这一层优化通常能把利用率从 1% 拉到几个百分点但还不够。因为数据复用问题没有解决全局内存流量依然巨大。2.3 第三重共享内存 Tiling复用才是核心分块Tiling的核心思想是把大矩阵切成小块先搬运到共享内存让同一个线程块内的线程反复读取这些小块。共享内存的带宽比全局内存高一个数量级延迟也低很多。以 16×16 的块为例线程块内 256 个线程协作先把 A 的一个 16×16 小块和 B 的一个 16×16 小块从全局内存搬运到共享内存然后每个线程只从共享内存里取数计算。内层循环每处理一个 K 维度上的块全局内存只需要读 2×16×16×42KB 数据但计算需要执行 16×16×16 次乘加数据复用倍数等于块大小。#define BLOCK_SIZE 16 __global__ void sgemm_tiled(const float* A, const float* B, float* C, int M, int N, int K) { __shared__ float As[BLOCK_SIZE][BLOCK_SIZE]; __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE]; int row blockIdx.y * BLOCK_SIZE threadIdx.y; int col blockIdx.x * BLOCK_SIZE threadIdx.x; float acc 0.0f; for (int kk 0; kk K; kk BLOCK_SIZE) { As[threadIdx.y][threadIdx.x] A[row * K kk threadIdx.x]; Bs[threadIdx.y][threadIdx.x] B[(kk threadIdx.y) * N col]; __syncthreads(); for (int k 0; k BLOCK_SIZE; k) { acc As[threadIdx.y][k] * Bs[k][threadIdx.x]; } __syncthreads(); } C[row * N col] acc; }这段代码有两个同步点缺一不可。第一个__syncthreads()保证 A、B 的小块都写入共享内存后才能开始计算第二个__syncthreads()保证所有线程都读完了当前块的数据才允许下一次迭代覆盖共享内存。如果漏掉第二个同步会出现某个线程还在读 As、Bs另一个线程已经往里面写新数据的情况结果随机出错。这层优化是从 Naive 到可用的关键一步通常能把利用率拉到 10%-20%。但如果块大小为 16你会发现共享内存访问还存在一个隐藏的性能杀手就是下一层要处理的 Bank Conflict。2.4 第四重Bank Conflict 与共享内存映射共享内存的硬件结构是 32 个 Bank每个 Bank 宽度 4 字节一次可以同时服务 32 个地址。理想情况下线程束内 32 个线程访问的地址应该均匀分布在 32 个不同 Bank 上。如果两个以上的线程命中同一个 Bank硬件就会把这些访问串行化这就是 Bank Conflict。在As[threadIdx.y][k]这样的内层循环里注意访问的是同一行、不同列。线程束内 32 个线程的 threadIdx.y 和 threadIdx.x 是交错的当它们同时访问 As 的不同行时如果共享内存数组按[16][16]存储每行 16 个 float 占 16 个 Bank第 0 行占 Bank 0-15第 1 行又从 Bank 0 开始两个不同行的同列数据会落在同一个 Bank 上产生冲突。最经典的解法是 Padding把As[BLOCK_SIZE][BLOCK_SIZE]改成As[BLOCK_SIZE][BLOCK_SIZE 1]也就是每行多出 1 个 float 的偏移。这样第 0 行占 Bank 0-16第 1 行从 Bank 1 开始线程束内不同线程访问的行就分散到了不同 Bank。__shared__ float As[BLOCK_SIZE][BLOCK_SIZE 1]; __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE 1];一行代码的变化可能带来接近翻倍的效果。Tensor Core 版本里为了同时兼顾 Bank Conflict 和向量化加载还会用到更复杂的 Swizzle 映射但核心逻辑都是让同一时刻线程束访问的地址散列到不同 Bank 上。2.5 第五重寄存器 Tiling、双缓冲与异步拷贝到这里访存模型已经基本健康但距离 95% 还很远。接下来要处理的是指令级并行和延迟隐藏。寄存器 Tiling 是指每个线程计算多个输出元素比如每个线程算 4×4 或 8×8 的小块。这样做有三个好处第一从共享内存读入一个数据后可以与多个累加结果复用第二每个线程的循环体中有多个独立的乘加链指令级并行度提高可以掩盖 FMA 指令的延迟第三所需的线程数变少同样的 SM 资源能容纳更多线程块提高占用率。双缓冲则是让“加载下一块数据”和“计算当前块”重叠。GEMM 内层循环里只有当前块计算完才能加载下一块否则共享内存会被覆盖。如果准备两块共享内存一个用于计算一个用于预加载下一轮数据就能把内存延迟隐藏掉。在 Ampere 及之后的架构上cp.async指令可以直接从全局内存异步拷贝到共享内存不需要过寄存器也不要求立即同步结合双缓冲效果更好。// 伪代码示意cp.async 发起异步拷贝 asm volatile(cp.async.ca.shared.global [%0], [%1], 16; :: r(smem_addr), l(gmem_addr));这一层是纯 CUDA Core 优化里最接近天花板的一步。到了这里一个 FP32 GEMM 通常能跑到理论峰值的一半左右接下来想再往上走就必须换一种思维方式用 Tensor Core。2.6 第六重Tensor Core 登场用 WMMA APITensor Core 是 GPU 上专门做小矩阵乘累加的硬件单元最早从 Volta 架构引入Ampere 上已经非常成熟。它的核心思想是让一条指令完成一次 4×4×4Volta或 16×16×16Ampere 及之后级别的矩阵乘加而不是像普通 CUDA Core 那样一条指令只做一个标量乘加。对开发者最友好的是 WMMAWarp Matrix Multiply-AccumulateAPI它把硬件细节封装起来编程模型是在 warp 粒度上操作 fragment。一个 16×16×16 的计算A 和 B 分别被切分成多个 fragment每个线程持有一部分数据mma_sync会在 warp 内完成整个矩阵乘累加。#include mma.h using namespace nvcuda; wmma::fragmentwmma::matrix_a, 16, 16, 16, __half, wmma::row_major a_frag; wmma::fragmentwmma::matrix_b, 16, 16, 16, __half, wmma::col_major b_frag; wmma::fragmentwmma::accumulator, 16, 16, 16, float c_frag; wmma::load_matrix_sync(a_frag, A_ptr, lda); wmma::load_matrix_sync(b_frag, B_ptr, ldb); wmma::fill_fragment(c_frag, 0.0f); wmma::mma_sync(c_frag, a_frag, b_frag, c_frag); wmma::store_matrix_sync(C_ptr, c_frag, ldc, wmma::mem_row_major);这段代码的底层效果是一次mma_sync把原本需要 16×16×16 次标量 FMA 的工作量变成一条硬件指令指令发射压力骤降数据搬运成为唯一制约因素。WMMA 对新手来说不用关心 fragment 内部到底怎么分布但使用时有几个硬性要求必须以整个 warp 为单位执行fragment 不能随意通过数组下标访问输入必须是半精度或 TF32 格式。这里面最容易踩的坑是忘记把 FP32 数据转成__half以及矩阵地址没有 16 字节对齐。2.7 第七重MMA PTX 与 ldmatrix 的精细控制WMMA API 已经很方便但封装也让一部分性能留在桌上。真正追求极致性能的库比如 CUTLASS会直接使用 PTX 级的指令。核心指令是mma.sync.aligned.m16n8k16。它的语义是一个 warp 内A 矩阵取 16×16、B 矩阵取 16×8K 维度 16累加到 8×8 的 C 矩阵上。A、B 数据不是从共享内存直接读的而是要求先分布在线程寄存器里。所以还需要另一个关键指令ldmatrix它可以把共享内存中的 16-bit 矩阵数据高效加载到线程寄存器的特定布局中。// 示意ldmatrix 加载 4 个 8x8 的 16bit 矩阵到寄存器组 asm volatile( ldmatrix.sync.aligned.m8n8.x4.b16 {%0,%1,%2,%3}, [%4]; : r(r0), r(r1), r(r2), r(r3) : l(smem_addr));为什么要绕开 WMMA 自己控制两方面的原因。第一WMMA 对 fragment 内部布局有固定选择而针对特定矩阵形状和 Swizzle 策略你可能想要不同的寄存器排列第二PTX 层面可以把ldmatrix、mma、cp.async的调度完全重叠起来减少同步等待。这个阶段的优化收益不是来自某一条指令本身而是来自整个数据流水线的精细编排。如果你不想手写 PTX直接读 CUTLASS 的 Cute 库也是一个很好的进阶路径它把数据布局和指令调度抽象成了更通用的描述但理解线下的 PTX 原理会让你更容易看懂它为什么那样设计。2.8 第八重调参与 Profiler 驱动优化当内核已经用上 Tensor Core剩下的事情就是调参。GEMM 的调参空间其实不大核心变量是这几个线程块维度BLOCK_M、BLOCK_N、K 维度块大小BLOCK_K、流水线级数NUM_STAGES、线程数、Swizzle 策略。组合起来可能上千种但真正有效的配置没几个。Nsight Compute 的 Speed of Light 是这一步最重要的工具。它会把内核运行时的计算管线和内存管线利用率同时显示出来。如果 Compute 利用率高而 Memory 利用率低说明计算正在吃满方向正确如果 Memory 利用率接近 100% 而 Compute 只有 50%说明数据搬运跟不上需要加大 Tiling 尺寸或增加流水线级数如果两者都不高可能是延迟没隐藏好需要提高占用率。以 A100 上的 TF32 GEMM 为例900 多个候选中能稳定跑到理论峰值 90% 以上的也就十几种配置而这十几种配置的共同特点是共享内存足够容纳 4-8 个 Stage 的流水线寄存器数量没有溢出Swizzle 没有引入额外的 Bank Conflict。调参不是随机试每一步都要用 Profiler 数据说明白到底在消除哪个瓶颈。2.9 第九重特化与边界处理最后一步是从“通用内核”走向“场景特化”。当 M、N、K 不是 Tile 的整数倍时要么在加载和存储时加边界判断要么在 Host 端 Padding 到对齐尺寸。边界判断会引入分支增加指令开销Padding 则会占用额外的显存和带宽。对大 GEMM 来说Padding 通常是更好的选择因为浪费很小而分支带来的性能损失是实时发生的。另一个维度的特化是 K 维切割。当 K 非常大时可以把 K 分成多段让不同线程块分别计算结果再用原子加或二阶归约汇总。这个方法叫 Split-K能显著提高小矩阵 GEMM 的并行度。更进一步的 Stream-K 则能精确控制每个线程块的工作量减少负载不均。到了这一层一个针对特定 shape 精心特化的 Tensor Core GEMM在 A100/H100 这类硬件上做到 95% 以上是很现实的。但也要清醒地认识到没有任何一个内核能同时吃遍所有 shape。95% 不是你写完一个通用 kernel 就有的而是针对你的实际生产 shape 反复裁剪出来的。3. 核心代码实现Tiled GEMM 怎么落地3.1 一个可直接编译的共享内存 Tiling 例程把第二重到第四重的思路合并可以得到一个可直接运行的教学版本。这段代码假设 M、N、K 都能被 16 整除方便你观察核心逻辑#include cuda_runtime.h #include cstdio #define BLOCK_SIZE 16 __global__ void sgemm_tiled_padded(const float* A, const float* B, float* C, int M, int N, int K) { __shared__ float As[BLOCK_SIZE][BLOCK_SIZE 1]; __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE 1]; int row blockIdx.y * BLOCK_SIZE threadIdx.y; int col blockIdx.x * BLOCK_SIZE threadIdx.x; float acc 0.0f; for (int kk 0; kk K; kk BLOCK_SIZE) { As[threadIdx.y][threadIdx.x] A[row * K kk threadIdx.x]; Bs[threadIdx.y][threadIdx.x] B[(kk threadIdx.y) * N col]; __syncthreads(); #pragma unroll for (int k 0; k BLOCK_SIZE; k) { acc As[threadIdx.y][k] * Bs[k][threadIdx.x]; } __syncthreads(); } C[row * N col] acc; } int main() { const int M 1024, N 1024, K 1024; size_t bytes M * K * sizeof(float); float* h_A (float*)malloc(bytes); float* h_B (float*)malloc(bytes); float* h_C (float*)malloc(M * N * sizeof(float)); for (int i 0; i M * K; i) h_A[i] 1.0f * (rand() % 10) / 10.0f; for (int i 0; i K * N; i) h_B[i] 1.0f * (rand() % 10) / 10.0f; float *d_A, *d_B, *d_C; cudaMalloc(d_A, bytes); cudaMalloc(d_B, bytes); cudaMalloc(d_C, M * N * sizeof(float)); cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice); cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice); dim3 block(BLOCK_SIZE, BLOCK_SIZE); dim3 grid(N / BLOCK_SIZE, M / BLOCK_SIZE); sgemm_tiled_paddedgrid, block(d_A, d_B, d_C, M, N, K); cudaDeviceSynchronize(); cudaMemcpy(h_C, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost); // 简单的正确性抽查 bool ok true; for (int i 0; i M; i) { for (int j 0; j N; j) { float ref 0.0f; for (int k 0; k K; k) ref h_A[i * K k] * h_B[k * N j]; if (fabs(h_C[i * N j] - ref) 1e-2f) { ok false; break; } } } printf(%s\n, ok ? PASS : FAIL); cudaFree(d_A); cudaFree(d_B); cudaFree(d_C); free(h_A); free(h_B); free(h_C); return 0; }注意这段代码里的共享内存行宽是BLOCK_SIZE 1很多人会奇怪为什么要浪费这 1 个 float其实就是为了消除 Bank Conflict上一节已经解释过。如果你把它改回BLOCK_SIZE功能完全一样但 Profiler 里会看到明显变高的 shared 冲突计数。3.2 共享内存容量与占用率手算选 Tile 大小不能只看算力利用率还要看硬件资源够不够。以 16×16 的块为例每个线程块需要两块 16×16×4B1KB 的共享内存总共 2KB。如果启用双缓冲就是 4KB四缓冲就是 8KB。A100 的一个 SM 有 228KB 共享内存但每 SM 最大线程数 2048。一个 16×16 块有 256 线程8 个块就能打满线程上限共享内存只用了 16KB很多空闲。如果改用 32×32 的块每个块 1024 线程2 个块就打满线程上限但共享内存需要 2×32×32×4B×216KB占用率也到 100%。你可能会想既然线程数先打满共享内存还有很多剩余那为什么不把双缓冲改成四缓冲完全正确多 Stage 流水线的代价就是共享内存线性增长8 个 Stage 的 32×32 块需要 8×8KB64KB这时共享内存先成为瓶颈并发块数量下降延迟隐藏能力反而可能变差。调参就是在“数据复用率”和“并发度”之间找平衡点。3.3 正确性验证不要用等号比较浮点数优化 GEMM 过程中最容易忽略的是正确性验证。线程块大小、累加顺序、向量化方式都会改变浮点运算的中间舍入顺序所以结果和 CPU 三重循环参考实现不可能逐位相等。验证时用相对误差更合理对于 FP32 计算通常允许 1e-4 到 1e-5 级别的差异如果使用了 TF32绝对值误差会放大到 1e-2 级别FP16 输入累加 FP32 大约在 1e-2 到 1e-3 之间。最好把验证阈值和精度绑定在一起不要拿同一个阈值测所有版本。性能测试时也有一个常见坑如果内核之后跟着cudaMemcpy而你没有做cudaDeviceSynchronize计时器的结束时间可能落在拷贝完成的前面得到的时间完全是乱的。正确做法是用 CUDA Event 包住内核同时保证两次 Event 之间的所有操作都完成。4. Tensor Core 思维从数据布局到指令发射4.1 为什么 Tensor Core 对数据布局这么苛刻Tensor Core 和普通 CUDA Core 最大的区别是它做乘法时数据不是从寄存器里随意拿两个数而是按照固定的 16×16 小块从一组线程的寄存器中批量取数。这就好比一个拼板机原材料必须在送进来之前就被裁剪成它规定的尺寸和形状。用 WMMA API 时load_matrix_sync会自动帮你完成从全局内存到 fragment 的布局转换但这个转换本身要消耗指令和寄存器。用 PTX 指令时你需要自己通过ldmatrix把共享内存中的数据按特定布局装进寄存器所以共享内存里的存储方式也要预先设计好这就是 Swizzle 存在的原因。一个成熟的 Tensor Core GEMM 内核整体数据流是全局内存 → 共享内存通过cp.async或 TMA→ 寄存器通过ldmatrix→ Tensor Core 执行 MMA。你会发现真正的 GEMM 优化早已不是“怎么算乘法”而是“怎么搬运数据”。4.2 一个完整的 WMMA 分块循环模式用 WMMA 写 GEMM 时建议的循环结构和共享内存 Tiling 类似只是内层计算从标量乘加变成了mma_sync。下面是一个片段式的模式#include mma.h using namespace nvcuda; // 每个线程块负责 64x64 的 C 块 // 每个 warp 负责 32x32 的 C 子块分 4 个 16x16 累加器 __global__ void sgemm_wmma(const __half* A, const __half* B, float* C, int M, int N, int K) { __shared__ __half As[64][64]; __shared__ __half Bs[64][64]; wmma::fragmentwmma::matrix_a, 16, 16, 16, __half, wmma::row_major a_frag; wmma::fragmentwmma::matrix_b, 16, 16, 16, __half, wmma::col_major b_frag; wmma::fragmentwmma::accumulator, 16, 16, 16, float c_frag[4]; int base_row blockIdx.y * 64; int base_col blockIdx.x * 64; for (int i 0; i 4; i) wmma::fill_fragment(c_frag[i], 0.0f); for (int kk 0; kk K; kk 64) { // 搬数据到 As/Bs... __syncthreads(); for (int kk_inner 0; kk_inner 64; kk_inner 16) { // 每个 warp 从 As/Bs 加载 16x16 的 fragment // 执行 4 次 mma分别对应 4 个输出 16x16 子块 __syncthreads(); } } // store c_frag[0..3] 到 C 的对应位置 }这段代码略去了坐标计算和边界判断但对理解 WMMA 分块模式足够了。实现时最容易出错的地方是 A、B 的矩阵布局。WMMA 的matrix_a和matrix_b有行主序和列主序之分一旦声明错了布局load_matrix_sync虽然不会报错但结果会是完全错误的值而且这种错误非常难排查。4.3 CUTLASS 给我们的启示提到 Tensor Core GEMM绕不开 CUTLASS。它的核心抽象是把 GEMM 的复制Copy、矩阵累加Mma、尾处理Epilogue彻底解耦每种操作都能独立选择和调优。Hopper 架构上CUTLASS 甚至可以借助 TMATensor Memory Accelerator和wgmma指令让数据从全局内存直接流进 Tensor Core中间不需要经过寄存器搬运。对普通开发者来说不一定要全部消化 CUTLASS但你可以观察它内部用的 Tile 形状、Swizzle 函数、Stage 数量再回到自己的内核里实验这些参数。很多时候把 CUTLASS 某个 SASS 级配置抄下来适配到自己的场景比从零瞎猜高效得多。5. 实测与调参怎么确认自己到了第几重5.1 Nsight Compute 的基本使用不要靠感觉优化用 Profiler 数据说话。最常用的是 Nsight Compute命令行就可以跑ncu --set full ./your_gemm或者只看关键摘要ncu --section SpeedOfLight ./your_gemmSpeed of Light 页面会给出两个核心百分比Compute Throughput 和 Memory Throughput。如果 Compute 高、Memory 低说明内核正在做正事如果 Memory 高、Compute 低说明它是 Memory Bound两者都不高通常意味着延迟隐藏不够需要提高并行度。其他关键指标还包括dram__bytes_read.sum看全局内存读取量shared_load/store看共享内存吞吐smsp__sass_inst_executed_op_shared_ld附近的计数器看是否有 Bank Conflict。在 Nsight Compute 的 Source 视图里甚至可以直接看到是源码的哪一行贡献了最多的 stall。5.2 从 1% 到 95% 的实测决策表我把每一次 Profiler 诊断对应的优化动作整理成表实战时直接对照诊断结果主要瓶颈下一层动作DRAM Throughput 高Compute 低全局内存访问太多做共享内存 TilingL1/TEX 吞吐高DRAM 不高数据复用不足或调度开销大加大 Tile / 寄存器 Tilingshared 指令存在 Bank Conflict共享内存访问模式冲突Padding / Swizzle计算管线活跃但指令槽浪费标量 FMA 指令过多改用 Tensor CoreTensor Core 活跃但等待数据搬运和计算未重叠双缓冲 / 多 Stage / cp.async整体利用率接近但差几个点边界分支 / 负载不均特化 shape / Padding / Split-K这套决策表基本就是九重天跃迁的操作手册。每一步都对应一个明确的诊断信号不是玄学调参。5.3 环境准备CUDA Toolkit、驱动与架构兼容很多人内核写好了性能也优化了结果一换机器就报 Illegal Address 或者内核启动失败最后发现是 Toolkit 版本和 GPU 架构不匹配。比如新一代显卡的算力标识是 SM_120旧版 CUDA Toolkit 编译出来的 cubin 根本没法加载运行时就会提示“no kernel image is available”。遇到这种情况先做三个检查nvcc --version查编译工具链nvidia-smi查驱动支持的 CUDA 版本ncu --version查 Profiler 版本。驱动和 Toolkit 不一定必须完全一致驱动版本要大于等于 Toolkit 要求的最低版本而内核编译时指定的-archsm_xx必须是目标 GPU 能支持的架构。写 PTX 或内联汇编时还要确认目标架构是否支持相关指令比如cp.async需要 Ampere 及以上wgmma需要 Hopper 及以上。6. 常见问题与排查技巧实录6.1 矩阵维度不是 Tile 整数倍怎么办最省事的方法是 Host 端 Padding把 M、N、K 向上取整到 16 或 64 的倍数多余区域填 0。这样做不会影响结果因为填充区域的乘加不影响有效输出。如果不方便 Padding可以给加载和存储加 Predicate 判断。比如共享内存加载时如果全局坐标超出 M 或 K 范围就写入 0存储到 C 时如果坐标超出范围就跳过。这种做法灵活但会在 warp 内引入分支分支收敛状态不同步会导致部分线程空转性能会打折扣。优先考虑 Padding。6.2 结果不对但看起来每个值都“差不多”如果输出矩阵的每个值都和参考实现对不上差的幅度不大不小先怀疑布局问题。WMMA 里把 row_major 和 col_major 写反、PTX 里 A/B 操作数顺序颠倒、ldmatrix 的地址不是 16 字节对齐都可能产生这种情况。这类错误用断点很难查最快的方法是写一个小矩阵测试比如 16×16×16把 fragment 里的值逐个 dump 出来和 CPU 参考对照哪个 EM 阶段出了问题一下就能看出来。另一个技巧是用compute-sanitizer跑一遍它能捕捉到未初始化值加载、越界访问这类隐藏错误。6.3 明白共享内存 Bank Conflict 但不知道哪些指令触发的Nsight Compute 的--set full里可以直接看l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum这类计数器。不为 0 就说明有冲突。再结合 Source 视图能精确定位到某条LDS指令。修复时先试 Padding如果 Padding 影响了 ldmatrix 需要 16 字节对齐那就要做一个完整的 Swizzle 映射把共享内存地址重新编码让线程束内访问的地址在 Bank 上错开同时保证每个 16 字节向量仍在连续的 4 个 Bank 内。这是 Tensor Core GEMM 里最考验功力的部分。6.4 内核工作一会儿就崩溃或挂起如果内核在大矩阵上崩溃先怀疑索引越界。Naive 和 Tiling 版本里最常见的 bug 是row或col超出矩阵范围读取了非法地址。解决办法是先把计算结果清零程序还能跑也方便对比。挂起则要优先检查__syncthreads()是否在条件分支里。如果线程块内某些线程进入了同步点另一些没有就会造成死锁。注意__syncthreads()必须让块内所有线程都能到达不能被if (tid 32)之类的分支隔离。6.5 性能始终上不去但 Profiler 一切正常这种情况我遇到过不止一次最后的原因往往很基础编译时没加足够高的优化选项。没有-O3的话编译器可能不会自动展开循环、向量化访问也无法生成最高效的 SASS。另一个原因是 host 端没有把输入布局转成内核期望的 K-Major 布局DataSource 不对GPU 端再优化也只是把错的数据搬来搬去。排查的时候先确认三件事编译选项带没带-O3输入矩阵是不是行主序且字节对齐cudaMalloc的地址是否满足 16 字节对齐要求。很多时候性能问题不是算法问题而是工程细节问题。这一路从 1% 到 95%本质上是不断追问“数据从哪来、到哪去、在哪个环节等了多久”。Naive 版本慢是因为数据从全局内存反复搬运Tiling 版本快是因为把数据复用在了共享内存Tensor Core 快是因为把计算从标量指令换成了矩阵指令。实际操作中我习惯每完成一层优化就把对应的 Profiler 截图和参数记录保存下来这样不仅知道当前卡在哪还能在后续改动中快速 revert 到某个稳定版本。如果你的内核已经能稳定跑到 90% 以上下一步可以试试针对具体 shape 做特化或者深入研究 CUTLASS 里面的 cute 抽象那又是一片新天地。