cuBLAS GEMM/GEMV 基准测试：精度、参数与带宽避坑指南

简介：cublas_benchmarks 是一份面向 CUDA/C++ 开发者与高性能计算学习者的轻量级基准测试工程，用于测量并比较 NVIDIA CuBLAS 库执行 GEMM（通用矩阵乘法）和 GEMV（矩阵向量乘法）时的性能。工程围绕不同矩阵、向量规模组织测试，涉及数据初始化、CuBLAS 函数调用、GPU 同步与耗时记录等环节，适合用来观察特定硬件配置下的吞吐表现，并为代码优化或计算策略选择提供参考。

压缩包共 10 个文件，约 13KB，包含 4 个 .h 头文件、4 个 .cu CUDA 源文件，以及 1 个 README 和 1 个 Makefile；头文件承担接口与工具声明，CUDA 源文件承载矩阵工具及 GEMM/GEMV 测试逻辑，Makefile 便于编译链接。已有 1095 人学习下载，说明该主题在 CUDA 学习与性能调优中具有一定关注度。读者可从中了解基准测试程序的目录组织、CuBLAS 调用方式、计时与同步写法，并借鉴其测试不同规模矩阵性能、排查配置差异的思路。

1. 从一次 3 倍性能误判说起：cublas_benchmarks 到底在测什么

两块卡跑同一个 4096³ 的矩阵乘，报表上差了三倍，换台机器重跑又只剩 1.2 倍。这类玄学我碰过不止一次，根子往往不在硬件，而在基准本身：把列主序当行主序用、没预热就开表、把 TF32 的数字和 FP32 的数字放进同一张对比图。cublas_benchmarks 就是冲着这些坑写的一套 C++ 基准程序，围绕 cuBLAS 把 GEMM 和 GEMV 两条最常见的路径单独拆出来测，分别喂不同的 M/N/K、数据类型和转置标志，看吞吐怎么变。它不产出好看的跑分海报，产出的是能复现、能横向对齐的数字。适合正在做推理或训练性能评估、要拿自家卡和主流卡对齐，或者单纯想知道某个 shape 到底该走 gemm 还是 gemv 的人。

2. C++ 工程化跑通：CMake、cuBLAS 链接与第一次计时

拿到这类基准包，第一步不是急着看跑分，而是先确认自己能不能把它编译出来、跑出和作者一致量级的数字。cuBLAS 的基准程序对工具链很敏感：CUDA 版本、驱动版本、编译时是否开了 fast math、链接的是静态库还是动态库，任何一项不同都可能让结果偏移 10% 以上。所以这一章的动作顺序是固定的——先把工程搭起来，再把计时框架搭对，最后才谈参数扫描。

2.1 依赖确认与编译骨架

基准程序依赖三样东西：CUDA Toolkit（提供 nvcc 和 cudart）、cuBLAS（随 Toolkit 一起装，头文件是cublas_v2.h）、一个 C++17 的编译器。检查方式很直接：

nvcc --version # 确认 CUDA 版本，和驱动支持的上限对齐 ls /usr/local/cuda/include/cublas_v2.h ldconfig -p | grep cublas # 确认动态库在搜索路径里

如果 cuBLAS 头文件在，但链接时报undefined reference to cublasCreate_v2，八成是 CMake 里只写了-lcublas而没带库路径，或者 Toolkit 装在非默认位置。常见做法是用 CMake 的CUDAToolkit模块，让它自己去找，比手写-L稳得多：

cmake_minimum_required(VERSION 3.20) project(cublas_benchmarks LANGUAGES CXX CUDA) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CUDA_STANDARD 17) # 让 CMake 自己去定位 Toolkit，避免手写 -L/-I 在换机器后失效 find_package(CUDAToolkit REQUIRED) add_executable(gemm_bench src/gemm_bench.cu src/timer.cpp) target_link_libraries(gemm_bench PRIVATE CUDA::cublas CUDA::cudart) # 单独开一条 release 配置，避免 Debug 下 -G 把内核性能打回原形 target_compile_options(gemm_bench PRIVATE $<$<CONFIG:Release>:-O3> $<$<COMPILE_LANGUAGE:CUDA>:--generate-code=arch=compute_80,code=sm_80>)

这里有两个参数值得说明。--generate-code里的sm_80要换成自己卡的计算能力，写错了不会报错，只会让 cuBLAS 走通用内核，测出来的数字偏低且没有任何提示，属于典型的无声翻车。-O3只加在 Release 配置上，因为 Debug 配置下 nvcc 会插入大量边界检查，GEMM 小尺寸下甚至能慢一倍。编译指令就一句：

cmake -B build -DCMAKE_BUILD_TYPE=Release && cmake --build build -j

2.2 行主序与列主序：一次调用就把语义定死

cuBLAS 是列主序的，而绝大多数 C++ 代码里的矩阵是行主序的。基准程序如果在这一步搞混，结果矩阵会转置，性能数字也会因为访问模式改变而失真。标准做法是不动数据，只动参数：要算行主序的C(M×N) = A(M×K) * B(K×N)，就把它改写成列主序的C^T(N×M) = B^T(N×K) * A^T(K×M)，然后两个操作标志都传CUBLAS_OP_N。

cublasHandle_t handle; cublasCreate(&handle); const int M = 4096, N = 4096, K = 4096; const float alpha = 1.0f, beta = 0.0f; // 行主序 A(MxK) 按列主序看是 KxM，leading dim = K // 行主序 B(KxN) 按列主序看是 NxK，leading dim = N // 行主序 C(MxN) 按列主序看是 NxM，leading dim = N // 于是调用变成 C^T(NxM) = B^T * A^T，m=N, n=M, k=K cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, // 注意：第一个参数是列数，不是行数 &alpha, dB, N, // B 在前 dA, K, // A 在后 &beta, dC, N);

cublasSgemm的参数顺序是(transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)。最容易错的是m和n的位置——它对应的是列主序视角下的维度，正好和直觉相反。另一个高频错误是lda：它必须大于等于矩阵的行数（列主序下），写成M或K都能编译通过，但会读到越界内存，测出来的时间毫无意义。写完这一段，强烈建议先用 M=N=K=256 的小矩阵和 CPU 参考实现比一次数值，确认布局对了再上大尺寸。

2.3 计时框架：预热、迭代与事件同步

用cudaEvent计时本身不难，难的是把三个干扰项排掉：首次调用的初始化开销、内核启动的异步性、以及单次测量本身的抖动。

// 预热：首次调用会建 plan、选内核，个别版本还会触发 JIT for (int i = 0; i < 5; ++i) { cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, dA, K, &beta, dC, N); } cudaDeviceSynchronize(); cudaEvent_t start, stop; cudaEventCreate(&start); cudaEventCreate(&stop); const int iters = 50; cudaEventRecord(start); for (int i = 0; i < iters; ++i) { cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, dA, K, &beta, dC, N); } cudaEventRecord(stop); cudaEventSynchronize(stop); // 必须同步，否则 stop 记录的是入队时间 float ms = 0.0f; cudaEventElapsedTime(&ms, start, stop); // 2*M*N*K 是一次 GEMM 的浮点运算量（一乘一加算两次） double tflops = 2.0 * M * N * K * iters / (ms * 1e-3) / 1e12; printf("M=%d N=%d K=%d iters=%d ms=%.3f tflops=%.2f\n", M, N, K, iters, ms, tflops);

iters的取法有讲究：小尺寸（比如 256）单次可能只有几十微秒，事件精度不够，要取 500 甚至 1000 次；大尺寸（8192 以上）单次就要几十毫秒，取 20 次足够，再多会让卡撞温度墙。50 是个折中值，但扫描脚本里最好按FLOPs / 目标时长动态算——常见做法是让每组至少跑满 200ms。cudaEventSynchronize那一行不能省，省掉之后stop事件记的是内核入队时刻，测出来的时间会小得离谱。

3. GEMM 基准的参数面：M/N/K、tile 与数据类型怎么组合才测得准

GEMM 是 compute-bound 的典型，也是 cuBLAS 里优化最狠的一条路径。正因为优化狠，它对参数的敏感度也高：同一个 4096 的方阵，换成 4096×4096×64 的长条形，吞吐可能掉到十分之一。基准程序的价值就在于把这个变化曲线画出来，而不是只报一个峰值。

3.1 参数扫描：取点策略比取点数量重要

全量扫 M/N/K 三个维度是组合爆炸，8192 以内每个维度取 8 个点就是 512 组，每组 200ms 也要跑近两分钟，还不算编译和预热。合理的做法是分层：先扫方阵确定峰值区间，再在峰值附近扫长条形，最后补几个实际业务里的真实 shape。

import csv import re import subprocess # 第一层：方阵，找峰值区间 SQ = [256, 512, 1024, 2048, 4096, 8192] # 第二层：固定 K，放大 M/N，看 tile 效率 RECT = [(m, m, 64) for m in (1024, 4096, 8192)] + \ [(m, m, 256) for m in (1024, 4096, 8192)] # 第三层：业务真实 shape，比如 batch 推理里的 [B*S, H] * [H, H] REAL = [(4096, 4096, 11008), (2048, 2048, 11008)] rows = [] for m, n, k in SQ + RECT + REAL: out = subprocess.run( ["./build/gemm_bench", "-m", str(m), "-n", str(n), "-k", str(k), "--dtype", "fp32", "--iters", "50"], capture_output=True, text=True, check=True).stdout tf = float(re.search(r"tflops=([0-9.]+)", out).group(1)) rows.append((m, n, k, tf)) print(f"{m}x{n}x{k}: {tf:.2f} TFLOPS") with open("gemm_fp32.csv", "w", newline="") as f: csv.writer(f).writerows([("m", "n", "k", "tflops")] + rows)

这段脚本的逻辑是三次分层而不是一次全组合：SQ负责找峰值在哪个体量上，RECT负责暴露 K 很小的时候 tile 利用率不足的问题，REAL把结论拉回业务。check=True不能去掉，不然某组 OOM 或者参数非法时脚本会静默跳过，最后画出来的曲线缺一段自己都不知道。--iters 50对小尺寸偏少，可以在脚本里按m*n*k加个条件分支，小尺寸自动抬到 500。

3.2 数据类型矩阵：FP32、TF32、FP16 不能放进同一张表

cuBLAS 的数学模式会直接改变实际执行的指令。FP32 输入在默认模式下，某些版本会自动走 TF32 的张量核心路径，算出来的数字漂亮得不像话，但精度只有 10 位尾数。基准程序必须显式控制这一点：

// 显式关掉 TF32，保证测的是真 FP32 cublasSetMathMode(handle, CUBLAS_DEFAULT_MATH); // 需要 TF32 时才打开，并单独标记这组结果 cublasSetMathMode(handle, CUBLAS_TF32_TENSOR_OP_MATH);

对应到扫描脚本，--dtype参数要能区分fp32、tf32、fp16（FP16 输入 + FP32 累加）、fp16-accum（FP16 累加）。四组数各自成表，绝不合并。经验上，同一张卡同一尺寸，TF32 相对纯 FP32 常见有 2～4 倍差距，FP16 输入再翻一倍左右，混在一张图里得出的"这卡不行"结论基本是错的。

表里最后一列的相对误差是 max relative error 的量级参考，具体阈值要按自己的 K 大小调——K 越大累加误差越大，这个规律比绝对值本身更重要。

3.3 batch 与 strided 变体：别把 N 次启动当成一次 GEMM

很多业务场景是一堆小矩阵乘，直接用循环调cublasSgemm会把内核启动开销算进性能里，测出来的数字远低于硬件能力。cuBLAS 提供了cublasSgemmStridedBatched和cublasGemmBatchedEx两条路，基准程序里应该把它们和单次 GEMM 并列成一组：

const int batch = 32; long long strideA = (long long)M * K; long long strideB = (long long)K * N; long long strideC = (long long)M * N; // 一次提交 32 个小 GEMM，避免 32 次内核启动的固定开销 cublasSgemmStridedBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, strideB, dA, K, strideA, &beta, dC, N, strideC, batch);

strideA/B/C是相邻两个 batch 之间起始地址的间隔，单位是元素个数不是字节，写成字节会直接越界。batch从 1 加到 32 的过程中，如果单 batch 吞吐明显上升，说明前面循环版本的瓶颈在启动开销而不是算力——这个结论只有把两种写法并排放才能看出来。

4. GEMV 基准的带宽视角：为什么算力表在这里完全失效

GEMV 是y = A * x，M×K 的矩阵只和一个长度为 K 的向量相乘。它的算术强度只有 0.5 FLOP/Byte 左右，也就是说每做一次乘加要搬 4 个字节。这种形状下，SM 的计算单元绝大部分时间是闲着的，性能上限由显存带宽决定。用 FLOPS 去评价 GEMV，就像用马力评价一艘船的载重，指标本身就选错了。

4.1 用有效带宽替代 FLOPS

基准程序在 GEMV 路径上应该直接输出 GB/s，而不是 GFLOPS。计算方式是把实际搬运的字节数除以耗时：

// GEMV: y(M) = A(MxK) * x(K)，FP32 // 主要搬运量是矩阵本身，向量可以忽略但为了严谨还是算上 double bytes = 1.0 * M * K * sizeof(float) + 1.0 * (M + K) * sizeof(float); double gbps = bytes * iters / (ms * 1e-3) / 1e9; // 算术强度，用来判断到底该看带宽还是看算力 double intensity = (2.0 * M * K) / bytes; printf("M=%d K=%d gbps=%.1f intensity=%.2f FLOP/Byte\n", M, K, gbps, intensity);

intensity这个值很有用：它低于 10 就基本可以确定是带宽瓶颈，高于 100 才是算力瓶颈，中间地带两者都有影响。同一个矩阵规模，GEMV 测出来的 GB/s 应该接近卡的峰值带宽的 70%～90%，如果只有 30%，通常不是硬件问题，而是访问模式或转置标志选错了。

4.2 转置标志对 GEMV 的影响远大于 GEMM

GEMV 里A是不是转置，决定了读矩阵时是顺着走还是跳着走，对带宽利用率的影响在 GEMV 上被放大得特别明显。基准程序必须把四种组合都跑一遍：

// 1) y = A*x，A 列主序，读 A 时沿 K 方向连续 cublasSgemv(handle, CUBLAS_OP_N, M, K, &alpha, dA, M, dx, 1, &beta, dy, 1); // 2) y = A^T*x，A 的访问变成跨行，带宽利用率通常明显下降 cublasSgemv(handle, CUBLAS_OP_T, M, K, &alpha, dA, M, dx, 1, &beta, dy, 1);

注意lda在两个分支里都是M，因为它描述的是列主序下矩阵的行数，和是否转置无关。这一点经常被写反，写反之后OP_T那组的数字会异常地好，因为实际读的范围变小了。跑完四组，把OP_N和OP_T的 GB/s 并排看，差距超过 40% 就要怀疑是不是访问模式出了问题，而不是急着下"这卡转置性能差"的结论。

4.3 用规模扫描找到 L2 溢出的拐点

GEMV 的性能曲线在矩阵超出 L2 缓存时会出现明显拐点。基准程序按 M/K 从小到大扫一遍，把 GB/s 画出来，拐点位置基本就等于 L2 容量除以 4 字节再开方。

import subprocess, re for k in (512, 1024, 2048, 4096, 8192, 16384): out = subprocess.run( ["./build/gemv_bench", "-m", "4096", "-k", str(k), "--dtype", "fp32"], capture_output=True, text=True, check=True).stdout gbps = float(re.search(r"gbps=([0-9.]+)", out).group(1)) print(f"K={k:6d} {gbps:7.1f} GB/s")

跑完之后，如果 K 从 2048 涨到 4096 时 GB/s 掉了一大截，那就是缓存装不下了。这个拐点的位置本身就是一个有用的硬件特征：它应该和卡标称的 L2 容量对得上。如果拐点跑到很靠前的位置（比如 K=512 就掉了），说明基准程序里有其他内存流量在抢带宽，通常是结果向量每次都被重新分配或者没有复用。

5. 避坑清单：五次让我拿到假数据的现场

这一节按现象、原因、解决三段式记，都是我在不同机器上真实撞过的。基准程序跑出漂亮数字的时候，先别高兴，把这五条挨个排除一遍。

5.1 第一次调用慢十倍，后面又正常

现象：脚本跑的第一组数据耗时是后续几组的 5～10 倍，单独重跑那一组又恢复正常。

原因：cuBLAS 的 handle 创建、内核选择、以及部分版本的 JIT 编译都发生在首次调用，这段开销被算进了计时。另外首次分配显存也可能触发页表建立，属于一次性成本。

解决：在开始计时前固定跑 5 次同参数的调用并cudaDeviceSynchronize()。如果参数组合很多，可以在每组开始前都做一次预热，代价是总时长增加约 10%，换来的是数据干净。这个习惯一旦养成，后面所有基准都不会再被首调用污染。

5.2 数值对得上，吞吐只有别人的一半

现象：小尺寸校验通过，大尺寸结果也对，但 TFLOPS 明显低于同型号卡的公开数据。

原因：多半是行主序转列主序时用了CUBLAS_OP_T而不是调整参数顺序。加了转置标志之后，cuBLAS 要么走一个额外的转置内核，要么用非合并的访问模式读矩阵，吞吐直接腰斩。

解决：回到 2.2 节的写法，不改数据只改参数，两个操作标志都传CUBLAS_OP_N，靠交换 A/B 的位置和调整m/n/k的顺序来表达行主序语义。改完之后如果吞吐仍然上不去，再检查lda是否等于矩阵行数。

5.3 FP32 跑出了超过理论峰值的数字

现象：一张卡的 FP32 非张量峰值大概几十 TFLOPS，基准程序却报出一百多。

原因：数学模式没显式设置，驱动或库版本默认启用了 TF32 张量核心路径。输入还是 FP32，内部降精度到 10 位尾数计算，吞吐翻几倍，精度悄悄掉了。

解决：基准程序启动时立刻调用cublasSetMathMode(handle, CUBLAS_DEFAULT_MATH)，把这个默认行为钉死。需要测 TF32 就跑单独一组，结果文件名带上tf32后缀。任何时候都不要让同一个 CSV 里混着两种精度的数。

5.4 GEMV 曲线出现断崖，被当成算力问题

现象：K 增大到某个值之后 GB/s 直接掉一半，有人据此判断"这块卡不适合大 K 的 GEMV"。

原因：矩阵超出了 L2 缓存，从缓存命中变成纯显存访问。这不是卡的缺陷，是任何卡都会有的拐点，只是位置不同。

解决：把输出指标从 TFLOPS 换成 GB/s，同时打印算术强度。强度低于 10 的结果一律归类为带宽受限，和算力数据分开放。拐点位置记下来，作为这块卡的缓存特征，后续做 shape 选型时能直接参考。

5.5 同一张卡两次运行差 20%

现象：完全相同的命令，隔十分钟跑一次，结果差 20% 以上。

原因：两个可能。一是温度或功耗墙，长时间跑大尺寸 GEMM 会让核心降频；二是机器上还有别的进程占着同一张卡，哪怕只占了 10% 的 SM 也会明显影响大尺寸内核。

解决：跑之前用nvidia-smi确认没有其他进程，把--iters控制在总时长 30 秒以内避免过热。有条件的话锁频再跑，锁频前后各测一次，两次结果一致才说明测量本身是稳的。这条是最容易忽略的，因为大部分人只跑一次就下结论。

6. 进阶验证：把基准数字落到 roofline 和选型决策上

有了干净的 GEMM 和 GEMV 数据，下一步不是画柱状图，而是把它们放进同一张 roofline 里看位置。做法很简单：横轴算术强度，纵轴实测吞吐，再画两条斜线——一条是带宽上限（GB/s × 强度），一条是算力上限（峰值 TFLOPS）。

以常见的 4096³ GEMM 为例，搬运量约为3 × 4096² × 4字节，约 201MB；运算量2 × 4096³，约 137 GFLOP，算术强度在 680 FLOP/Byte 量级，远在算力上限那条水平线上，属于深度计算受限。同样 4096×4096 的 GEMV，强度只有 0.5 FLOP/Byte，稳稳落在带宽斜线的低端。两条曲线之间的交点，就是这块卡从带宽受限切换到算力受限的分界，数值上大致等于峰值FLOPS / 峰值带宽，常见数据中心卡在 20～40 之间。

真正有用的判断发生在中间地带。比如M=4096, N=4096, K=64这种长条形 GEMM，强度掉到 30 左右，正好卡在交点附近，这时候光看 TFLOPS 会误判，得同时报 GB/s 和 TFLOPS，看哪个先撞墙。我一般会在扫描脚本里直接算好强度列，跑完先按强度排序，把 10 以下和 100 以上的分开，中间的单独拎出来人工看。

数值验证这一步也别省。基准程序跑得快不代表算得对，尤其是开了 TF32 或者 FP16 之后：

import numpy as np # benchmark 用 --dump 把 C 以行主序 FP32 写到 c.bin ref = a.astype(np.float32) @ b.astype(np.float32) # CPU 参考 c = np.fromfile("c.bin", dtype=np.float32).reshape(512, 512) rel = np.abs(c - ref) / (np.abs(ref) + 1e-6) print("max_rel=%.3e mean_rel=%.3e" % (rel.max(), rel.mean()))

阈值按精度模式卡：纯 FP32 的max_rel应该小于 1e-5，TF32 小于 1e-2，FP16 输入加 FP32 累加小于 1e-3。超过阈值不是"误差大一点无所谓"，而是说明布局或者数学模式搞错了，这时候的性能数字一律作废。

最后一个技巧：把 cuBLASLt 的 heuristic 查询也接进基准。cublasLtMatmul允许传一个 workspace 上限和若干算法偏好，让库自己挑内核，同时支持 epilogue 融合（比如 bias、GELU）。同一组 shape 下，cuBLASLt 选出的算法有时候比默认路径快 10%～20%，差异在大 batch 小矩阵上尤其明显。基准程序里加一个--use-lt开关，把两条路径的数字并排存，选型时就有据可依，而不是拍脑袋。

从那以后我每次跑基准，都强制走一遍"预热 → 数值校验 → 带宽交叉验证 → 锁频复测"这四步，少一步都觉得手上的数字不踏实。希望帮到你。

本文还有配套的精品资源，点击获取