企业官网建设流程全解析

1. 别急着写kernel：先理解“取行”和“取列”到底在说什么

1.1 行主序矩阵的物理存储与访存方向

我做CUDA优化这些年，越来越觉得GEMM是所有并行编程课程里最适合用来练内功的题目。Stanford CS217里有一整套关于矩阵乘法的讨论，但我印象最深的不是某个trick，而是教授反复强调的一个基础问题：你的矩阵到底是怎么放进内存的。行主序下，A[i][j]的地址是i * ld + j，也就是说，你访问A[i][0]、A[i][1]、A[i][2]时，地址是连续递增的；但访问A[0][j]、A[1][j]、A[2][j]时，地址每跳一次就跨过一个完整的行宽。这个“取行”和“取列”的差别，看起来只是方向问题，实际上决定了GEMM后续所有优化手段的走向。

全局内存里，一个warp的32个线程去读32个连续float，只需要一个128字节事务；如果这32个线程读的是同一列的不同行，地址分散在32个不同的行区间里，可能就要产生几十个访存事务，带宽利用率直接掉一个数量级。共享内存更敏感，它有32个bank，同一指令里不同线程访问同一个bank的不同地址就会发生bank conflict，轻则多几个周期，重则把一次load变成32次串行操作。所以“取行还是取列”从来不是代码风格问题，而是算法复杂度级别的取舍。

用生活里的例子说，取行就像在书架上从左到右拿一排书，手顺着扫过去就行；取列就像每次从不同层拿第3本，每拿一本都要踮脚或弯腰。程序虽然不会抱怨“累”，但profiler会诚实地显示带宽利用率上不去，只是很多同学看不出来是因为访存方向错了。

1.2 向量内积：一个输出元素，牵涉两个方向

矩阵乘法的教科书定义就是向量内积：C[i][j]等于A[i][:]和B[:][j]这两个向量做点乘。按这个定义写kernel，每个线程负责一个输出元素，内层循环对k做累加。这个写法最直观，但有两个天生的问题。

第一个问题是，行主序下A[i][k]是连续取行，B[k][j]是跨步取列。虽然你可以通过调整线程映射，让一个warp里的线程分别处理相邻的j，这样同一时刻warp访问的是B的第k行里一段连续列，看起来合并了；但你的A访问变成了同一地址的广播。说得直白一点，内积公式里总有一个方向是逆着物理存储方向走的，你只能选择把惩罚放在哪个矩阵上，没法同时让两个矩阵都按照连续方向访问。

第二个问题是数据复用率低。A[i][k]会被所有需要第i行输出的线程重复读，B[k][j]会被所有需要第j列输出的线程重复读。如果不加共享内存tiling，每个输入元素从全局内存被读的次数会随着矩阵规模线性放大，访存流量直接到了2 * M * N * K的水平。矩阵稍微大一点，这个kernel就完全是带宽瓶颈，GPU的浮点算力根本用不上。

1.3 向量外积：把“点乘”改成“秩1更新”

换个角度看，矩阵乘法也可以写成C = Σ_{k} A[:,k] ⊗ B[k,:]。这就是向量外积视角：把A的第k列看成一个向量u，把B的第k行看成一个向量v，两者做外积得到一个M×N的秩1矩阵，把所有k的秩1更新累加得到C。

这个写法不会改变数学结果，但它把“取行”和“取列”分开了：B[k,:]天然是连续的一行，适合做合并访问和向量化加载；A[:,k]虽然在全局内存里是跨步的列，但我们可以先把它搬进共享内存或寄存器，再在计算阶段按需要的方向去读。也就是说，外积视角让数据搬运和计算可以分开设计，而不是像内积那样把跨步访问塞进最内层循环。

外积更核心的优势是复用。做秩1更新时，A[:,k]这一列可以被这一列对应行上的所有列复用，B[k,:]这一行可以被这一行对应列上的所有行复用。如果我们在一个block里同时加载了A的一小块和B的一小块，那么每个从全局内存读进来的元素都可能参与很多次乘加。这个复用程度，直接决定了你能把全局内存带宽需求压到多低，也是GEMM优化能做到多少TFLOPs的关键。

2. 内积与外积：GEMM并行化结构的分叉点

2.1 内积结构：线程数与输出元素一一对应

最朴素的CUDA kernel就是每个线程算一个C[i][j]，Thread与输出元素一一对应。这种结构的最大优点是线程之间完全独立，没有跨线程归约，不需要__syncthreads()，写出来的代码几乎不会有同步错误，非常适合做正确性baseline。

但它的缺点也非常明显。每个线程只有一个累加器，但内层循环每步要从A和B各读一个float，然后做一次乘加。你可以把它想象成一条流水线：每个线程像一个小工，每天只生产一个零件，却要跑遍整个仓库去提两种原料。线程数量虽然多，但单线程算力利用率很低，访存次数和乘加次数的比值高得惊人。

另一个问题是向量化困难。A行连续，理论上可以一次读4个float；但B列跨步，一次读4个B[k][j]就是4次独立访存，编译器很难自动优化。如果你写内积循环，建议不要期待编译器能帮你生成高效的float4指令，它通常会被跨步访问拖住，最终只能生成一堆标量load。

2.2 外积结构：粗粒度累加与寄存器复用

外积结构在CUDA里的典型形态是寄存器分块，也就是register tiling。每个线程不再只算一个C元素，而是负责C矩阵里一个TM×TN的子块。每个线程需要维护TM×TN个累加寄存器，内层循环每遍历一个k，就从共享内存或者寄存器里读入TM个A元素和TN个B元素，然后做TM×TN次乘加，更新整个子块。

这才是我理解的“外积路径”：一次读入一小段A列和一小段B行，做一次秩1更新。相比内积结构，A和B的每个元素会被更多输出重复使用，内存访问次数被摊薄了。假设TM=4, TN=4，每个线程一次内层迭代做16次乘加，但只需要分别读4个A元素和4个B元素，计算访存比立刻变得好看很多。

寄存器tiling还能让访存方向更可控。你可以把加载tile和计算tile拆开：加载阶段从全局内存按连续方向搬运数据，计算阶段再从共享内存按固定方向读取。想优化的取行就连续读，想优化的取列就通过padding、转置或者swizzle来把bank conflict压下去。所以在实际高性能GEMM里，几乎看不到一个线程只算一个输出元素的版本，大家都在做各种尺寸的register tiling。

2.3 什么时候该用谁：K维度、矩阵规模与硬件限制

我把内积和外积结构做一个对比，方便你看着选：

• 线程输出数：内积为1；外积为TM×TN
• 线程总数：内积需要M×N；外积只需要(M×N)/(TM×TN)
• 累加寄存器：内积1个；外积TM×TN个
• 访存复用：内积低，外积高
• 向量化友好度：内积的B列跨步很难向量化；外积连续读行+寄存器复用更友好
• bank conflict风险：内积更容易在取列方向踩坑；外积可以通过padding/swizzle控制
• 典型场景：内积适合小矩阵、正确性验证和动态形状；外积适合大矩阵、追求峰值性能的训练/推理GEMM

看起来外积全面占优，但实际工程里也不是无脑选。外积要占用大量寄存器，比如TM=4, TN=4就需要16个累加器，再加上中间变量，一个线程轻松用到40~80个寄存器，在部分GPU上会压低occupancy。矩阵规模特别小，或者K特别短的时候，tiling带来的搬运开销可能比计算收益还大，这时候简单内积反而更稳。真正的产品级GEMM通常会用TM/TN的各种组合适配不同的矩阵形状，再叠加double buffering、ldmatrix、swizzle等高级手段。

3. 一步步把kernel从内积改成外积

3.1 完整代码骨架与编译环境

我用CUDA 12.x，编译目标用-arch=sm_80（A100）或者sm_86/sm_89（Ampere/Ada），命令大概是：

nvcc -O3 -arch=sm_80 -o gemm gemm.cu

矩阵约定为行主序，先写一个计时函数，后面所有版本都用同一套计时逻辑：

cudaEvent_t start, stop; cudaEventCreate(&start); cudaEventCreate(&stop); cudaEventRecord(start); gemm_kernel<<<grid, block>>>(A, B, C, M, N, K); cudaEventRecord(stop); cudaEventSynchronize(stop); float ms = 0.0f; cudaEventElapsedTime(&ms, start, stop); double gflops = 2.0 * M * N * K / (ms * 1e6);

后面每个版本都建议先用小矩阵跟CPU计算做对比，确认没算错再谈性能。优化最忌讳一开始就上一堆tiling、double buffering，结果错了都不知道是哪一步引入的。

3.2 版本一：全局内存内积，先建立一个正确的baseline

第一版就是最简单直接的kernel：

__global__ void gemm_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; } }

启动配置一般用dim3 block(16, 16)，dim3 grid((N + 15) / 16, (M + 15) / 16)。这个版本性能很差，但它是后续所有版本的正确性参照。如果你连这个都写不对，后面加了共享内存更没法debug。我每次做新优化，都会先把naive版本跑一遍，确定输出矩阵没问题，再开始动结构。

如果你用Nsight Compute看这个kernel，会看到全局内存的访问事务特别多。原因就是前面说的：每个线程都要从A和B读大量元素，而且中间有跨步取列，访存流水线一直在等数据。这个阶段的瓶颈是内存带宽，不是计算。

3.3 版本二：共享内存tiling，观察取列带来的bank conflict

第二版加共享内存tiling。假设块大小是16×16，代码大概是：

#define BLOCK_SIZE 16 __global__ void gemm_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 bx = blockIdx.x, by = blockIdx.y; int tx = threadIdx.x, ty = threadIdx.y; int row = by * BLOCK_SIZE + ty; int col = bx * BLOCK_SIZE + tx; float sum = 0.0f; for (int kk = 0; kk < K; kk += BLOCK_SIZE) { As[ty][tx] = A[row * K + kk + tx]; Bs[ty][tx] = B[(kk + ty) * N + col]; __syncthreads(); for (int k = 0; k < BLOCK_SIZE; ++k) { sum += As[ty][k] * Bs[k][tx]; } __syncthreads(); } C[row * N + col] = sum; }

这个版本已经比naive快很多，因为A和B的tile都只从全局内存读一次，复用到了块内所有输出。但如果你仔细看计算阶段的访问，As[ty][k]是固定k列、不同ty访问不同行，在共享内存里属于“取列”访问；Bs[k][tx]是固定k行、不同tx访问不同列，属于“取行”访问。取行的方向通常没什么问题，取列的方向在不同线程组织下容易触发bank conflict。

这里我不给你一个“肯定冲突”的结论，因为bank conflict和线程布局强相关。最简单的自查方法是把共享内存的地址线性化，然后看同一条load指令里32个线程分别落在哪些bank。比如As[ty][k]，如果BLOCK_SIZE=32，ty从0到31，地址偏移是ty * 32 + k，对32取模之后全部等于k，这就会产生32路bank conflict。如果BLOCK_SIZE=16，一个warp里有两个不同的ty，冲突就变成2路。你只需要记住规则：不同线程在同一时刻访问同一个bank的不同地址就是冲突。加padding是解决这类问题的通用手段，把As[BLOCK_SIZE][BLOCK_SIZE+1]，让每一行多出一个float占位，列的方向就会被错开，很多冲突就消失了。

3.4 版本三：外积风格的寄存器tiling，让取行取列各归其位

第三版开始上register tiling，也就是外积结构。每个线程负责一个TM×TN子块，内层循环反复做秩1更新。我常用的一个演示配置是TILE_M=64, TILE_N=64, TILE_K=16, TM=4, TN=4，这样block里线程数是(64/4) * (64/4) = 16 * 16 = 256，比较均衡。

#define TILE_M 64 #define TILE_N 64 #define TILE_K 16 #define TM 4 #define TN 4 // blockDim.x = TILE_N / TN = 16 // blockDim.y = TILE_M / TM = 16 __global__ void gemm_outer(const float* A, const float* B, float* C, int M, int N, int K) { __shared__ float As[TILE_M][TILE_K + 1]; __shared__ float Bs[TILE_K][TILE_N + 1]; int tx = threadIdx.x; int ty = threadIdx.y; int rowBase = blockIdx.y * TILE_M + ty * TM; int colBase = blockIdx.x * TILE_N + tx * TN; float rC[TM][TN] = {0.0f}; for (int kk = 0; kk < K; kk += TILE_K) { for (int idx = tx + ty * blockDim.x; idx < TILE_M * TILE_K; idx += blockDim.x * blockDim.y) { int r = idx / TILE_K; int c = idx % TILE_K; As[r][c] = A[(blockIdx.y * TILE_M + r) * K + kk + c]; } for (int idx = tx + ty * blockDim.x; idx < TILE_K * TILE_N; idx += blockDim.x * blockDim.y) { int r = idx / TILE_N; int c = idx % TILE_N; Bs[r][c] = B[(kk + r) * N + blockIdx.x * TILE_N + c]; } __syncthreads(); for (int k = 0; k < TILE_K; ++k) { float rA[TM]; float rB[TN]; #pragma unroll for (int i = 0; i < TM; ++i) rA[i] = As[ty * TM + i][k]; #pragma unroll for (int j = 0; j < TN; ++j) rB[j] = Bs[k][tx * TN + j]; #pragma unroll for (int i = 0; i < TM; ++i) { #pragma unroll for (int j = 0; j < TN; ++j) { rC[i][j] += rA[i] * rB[j]; } } } __syncthreads(); } for (int i = 0; i < TM; ++i) { float* cptr = C + (rowBase + i) * N + colBase; for (int j = 0; j < TN; ++j) { cptr[j] = rC[i][j]; } } }

这个kernel的核心变化在计算循环里。每次k迭代，我先从共享内存把TM个A元素和TN个B元素读进寄存器，然后一次性做TM×TN次乘加。这其实就是外积视角里的秩1更新：用rA这一小段列和rB这一小段行，去更新rC这个子块。因为数据在寄存器里被反复复用，共享内存的访问次数大大减少，性能会明显高于前面版本。

代码里我给共享内存加了padding，As第二维是TILE_K+1，Bs第二维是TILE_N+1。这样做的目的是让取列方向的访问在多数布局下尽量分散到不同bank。需要说明的是，这个代码是教学级别的register tiling，不是CUTLASS那种极致优化版。比如Bs的rB[j]读取在TN=4的情况下可能还有少量bank冲突，工程上会继续用swizzle或者ldmatrix来彻底解决。但作为CS217课程作业级别的展示，它已经足够说明“外积+寄存器复用”带来的收益。

如果你用的GPU比较新，还可以把全局内存访问改成float4向量化，前提是矩阵leading dimension和地址都要对齐到16字节。这一步对带宽利用率提升非常明显，但会增加代码复杂度，尤其处理边界的时候很麻烦。

3.5 性能对比与调参要点

我把几个版本的性能趋势列一下，注意不要当成固定数值，因为不同GPU和矩阵尺寸差异很大：

调参的时候，我一般先固定TILE_K，从16开始，因为在大多数GPU上16×64的共享内存尺寸不会爆容量。然后调TM×TN，4×4是比较安全的起点；如果寄存器够用，可以试8×4或者8×8。TILE_M和TILE_N不要无脑放大，因为共享内存和寄存器的压力会同时上升，最终occupancy掉下去，性能反而下降。调参没有银弹，最好的方法是把矩阵规模固定，逐组参数跑一遍，再用Nsight Compute看瓶颈是内存还是计算。

4. 调试实录：踩过的坑与排查技巧

4.1 共享内存bank conflict的快速定位

我自己调试时，最常遇到的坑就是bank conflict。你可以用Nsight Compute跑一下，重点看这几个指标：shared memory的load/store bank conflict次数、全局内存的吞吐量、achieved occupancy。

如果怀疑某个load有bank conflict，最直接的办法是在纸上手算一次。取一个warp的32个线程，写出它们在同一指令下访问的共享内存地址，然后对32取模，看落到同一个bank的地址是否互不相同。注意：多个线程访问同一个地址属于broadcast，不算冲突；但访问同一个bank的不同地址就算冲突。很多教程把“同一bank”和“同一地址”混为一谈，排查时候容易判断失误。

解决bank conflict的常用手段有：给共享内存行加padding、把矩阵在共享内存里转置、使用swizzle模式重排地址。对初学者来说，padding最简单，通常一次就能解决大部分冲突。不过加padding会浪费一点共享内存容量，所以tile越大越需要仔细计算是否值得。

4.2 float4对齐、矩阵leading dimension和边界处理

用float4做向量化是GEMM提升带宽利用率的重要一步，但它要求地址16字节对齐。cudaMalloc返回的显存地址通常是256字节对齐，问题不大；坑往往出现在矩阵的leading dimension上。如果你的矩阵宽度不是4的倍数，那么每一行的起始地址可能会错开，导致float4访问越界。

解决方法有两个：一是用cudaMallocPitch对齐行宽；二是在分配矩阵时手动往每一行末尾多padding几个float，让宽度凑成16字节对齐。注意，这种padding会改变矩阵在内存里的布局，所以传给kernel的leading dimension也要跟着改。我在实际项目里更喜欢用ld参数单独传给kernel，而不是直接用N，这样代码更灵活，边界处理也更清晰。

矩阵的M、N、K不是tile整数倍时，主流做法是让tile稍微越界读，然后计算阶段用mask把越界位置过滤掉。这个操作很容易出错，尤其边界块的写回。我的习惯是先在host端把矩阵padding到tile整数倍，再做性能测试，等基本逻辑都稳定了再去处理复杂的边界分支，不要一开始就混合在一起。

4.3 寄存器压力、occupancy与launch bounds

register tiling的累加器数量直接决定寄存器占用。TM=4, TN=4有16个累加器，看起来不多；但加上指针计算、索引变量、rA/rB，一个线程用到40~64个寄存器很常见。如果GPU每SM的寄存器总数有限，线程块能跑多少个block就会变少，occupancy下降，反而可能掩盖掉计算效率的提升。

我的经验是，在GEMM这种计算密集型kernel里，适当地降低occupancy换取更多寄存器的指令级并行，通常是划算的。但你得用__launch_bounds__或者maxrregcount控制好上限，不能让它肆意占用导致kernel无法启动。比如一个block有256线程，每个线程用到64个寄存器，那block一共就需要16384个寄存器，在每SM 64K寄存器的GPU上刚好可以放4个block，这就是一个比较舒服的配置。

调寄存器最忌讳凭空猜测，直接用--ptxas-options=-v看编译报告，或者用cudaOccupancyMaxPotentialBlockSize这类API算一下。很多时候，性能瓶颈不是register pressure，而是你为了压低occupancy强行减少了寄存器数量，结果spill到本地内存，性能反而更烂。

4.4 我每次写GEMM前必做的3个检查

第一，先在纸上画出A、B、C的tile内存布局，把每个访问方向标清楚，哪些是取行，哪些是取列。这一步看起来麻烦，但能避免一半以上的方向性错误。第二，模拟一个warp的第一次内存访问，手算共享内存bank index，确认没有明显的bank conflict。第三，跑之前先写naive版本做正确性校验，跑完之后立刻用Nsight Compute看访存指标，不要只盯着总耗时。

说句实在话，CUDA GEMM优化这个题目，网上教程特别多，但很多人抄了tiling、抄了padding，性能还是上不去。问题往往出在最基本的方向判断上：你的算法结构到底是内积还是外积，每一步访问到底是取行还是取列。把这些底层逻辑想清楚，再去看那些高级kernel代码，你会突然明白为什么人家要这样排线程、这样存共享内存、这样写累加循环。这个“想清楚”的过程，比单纯背几个优化技巧值钱得多。