cuBLAS GEMM/GEMV 基准测试：计时、缓存与带宽解析

简介：cublas_benchmarks 是一份面向 CUDA/C++ 开发者与高性能计算学习者的 CuBLAS 性能测试工程，围绕 GEMM、GEMV 两类基础线性代数算子在 GPU 上的执行表现展开，可用于科学计算、机器学习与深度学习场景下的性能评估和调优参考。压缩包共 10 个文件，约 13KB，以 4 个 .cu CUDA 源文件、4 个 .h 头文件为主，并配有 Makefile 与 README，分别承载测试逻辑实现、数据结构与函数声明、编译构建和使用说明，目录组织简洁，便于按功能模块阅读与二次扩展。已有 1095 人学习下载。代码通过设置不同矩阵与向量规模调用 CuBLAS 函数，配合 GPU 同步与耗时记录，帮助读者观察不同参数、CUDA 版本和硬件配置下的性能差异，理解基准测试的完整流程，并为进一步优化计算策略、选择合适算子实现提供可复用参考，对需要复现实验、对比实现或排查性能瓶颈的读者尤具价值。

1. 同一张卡跑出三倍差距：cuBLAS 的 GEMM 与 GEMV 基准到底在测什么

同一张卡、同一组方阵，两份 cublas_benchmarks 跑出来的 GEMM 有效算力能差出两三倍，差别往往不在硬件，而在计时方式、L2 缓存命中情况和参数传法。这类基准做的事很具体：把 cuBLAS 里的 GEMM（矩阵乘矩阵）和 GEMV（矩阵乘向量）按不同维度、不同精度、不同转置组合各跑一遍，记录耗时并换算成算力（TFLOPS）和带宽（GB/s），最终得到一条“这张卡在当前形状下能跑多快”的曲线。它对应两类高频需求：自研 kernel 到底比 cuBLAS 快还是慢，以及硬件选型时瓶颈在算力侧还是显存带宽侧。写自定义算子要基线、调推理吞吐要找瓶颈、做验收要一个可复现数字的人，都适合按这套流程跑一遍，但必须先把测量规则定对。

2. 搭起基准骨架：依赖、目录与一个不误导人的计时函数

基准代码本身不难，难的是让它报出来的数字可被信任。我见过太多版本是用std::chrono包住cublasSgemm调用，然后再cudaDeviceSynchronize()，看起来合理，实际把主机端提交耗时也算进去了；矩阵小的时候，误差能到百分之几十。所以整个骨架要从计时函数开始搭，而不是从 kernel 开始搭。

2.1 依赖清单与目录结构

需要的东西不多，但每一样都要能自证存在，避免写到一半才发现链接不上。

目录按“源码 / 脚本 / 结果”三分，结果目录一律带时间戳，避免两次实验互相覆盖：

cublas_benchmarks/ ├── src/ │ ├── timer.cuh # 计时与错误检查宏 │ ├── bench_gemm.cu # GEMM 扫描 │ └── bench_gemv.cu # GEMV 扫描 ├── scripts/ │ ├── run_sweep.sh # 批量跑不同形状 │ └── plot.py # 读 csv 出曲线 └── results/ └── 20240101_1200/ ├── gemm.csv └── gemv.csv

目录本身不影响性能，但它决定了你能不能在一个月后复现当时那条曲线。血泪经验是：形状、精度、迭代次数、flush 开关这些参数必须一起落进 csv，只存耗时的话，回头根本解释不了那个尖峰是怎么来的。

2.2 用 cudaEvent 写一个能复用的计时函数

核心思路是把“测什么”和“怎么测”分开：被测体是一个可调用对象，计时器只负责 warmup、迭代、事件打点，以及把每次迭代前的缓存清理动作单独计价再扣掉。

// src/timer.cuh #pragma once #include <cuda_runtime.h> #include <cublas_v2.h> #include <cstdio> #include <cstdlib> #include <functional> #define CK(x) do { cudaError_t _e = (x); \ if (_e != cudaSuccess) { std::fprintf(stderr, "CUDA %s @%s:%d\n", \ cudaGetErrorString(_e), __FILE__, __LINE__); std::exit(1); } } while (0) #define CB(x) do { cublasStatus_t _s = (x); \ if (_s != CUBLAS_STATUS_SUCCESS) { std::fprintf(stderr, "cuBLAS %d @%s:%d\n", \ (int)_s, __FILE__, __LINE__); std::exit(1); } } while (0) // body: 被测调用；pre: 每次迭代前的 L2 清理，不做就传空 lambda inline float timeit(const std::function<void()>& body, const std::function<void()>& pre, int warmup, int iters) { for (int i = 0; i < warmup; ++i) { pre(); body(); } CK(cudaDeviceSynchronize()); // warmup 结果全部丢弃 cudaEvent_t e0, e1; CK(cudaEventCreate(&e0)); CK(cudaEventCreate(&e1)); // 先单独量一遍 pre 的成本，稍后从总时间里扣掉 CK(cudaEventRecord(e0)); for (int i = 0; i < iters; ++i) pre(); CK(cudaEventRecord(e1)); CK(cudaEventSynchronize(e1)); float pre_ms = 0.f; CK(cudaEventElapsedTime(&pre_ms, e0, e1)); CK(cudaEventRecord(e0)); for (int i = 0; i < iters; ++i) { pre(); body(); } CK(cudaEventRecord(e1)); CK(cudaEventSynchronize(e1)); float all_ms = 0.f; CK(cudaEventElapsedTime(&all_ms, e0, e1)); cudaEventDestroy(e0); cudaEventDestroy(e1); return (all_ms - pre_ms) / iters; // 单次迭代毫秒数 }

逻辑说明：cudaEvent打在流上，记录的是设备侧时间点，天然排除了主机提交开销，这是它比std::chrono靠谱的根本原因。用std::function包被测体，是为了让同一个计时器同时服务 GEMM、GEMV、batched 三种调用，不必写三份。

参数说明：warmup建议不低于 5，cuBLAS 句柄在首次调用时要做内部初始化，同时 GPU 时钟也需要一段时间的负载才能稳定频率，这两件事都会让第一次迭代明显偏慢；iters不能照抄，应该按“单次耗时 × 迭代数 ≥ 200 ms”反推，否则事件分辨率和时钟抖动会主导结果；pre的扣减是必须的，因为清理缓存的 kernel 与 GEMM 在同一个流上串行执行，会直接计入总量。

清理缓存的实现不必复杂，写一遍就行，关键是容量要对齐设备属性：

// 按设备的 L2 实际容量分配，不要写死一个数字 cudaDeviceProp prop; CK(cudaGetDeviceProperties(&prop, 0)); size_t flush_bytes = (size_t)prop.l2CacheSize * 2; float* d_flush; CK(cudaMalloc(&d_flush, flush_bytes)); auto flush = [&]{ CK(cudaMemsetAsync(d_flush, 0, flush_bytes)); };

提示：l2CacheSize在不同架构上差别很大，写死 40 MB 或 60 MB 都会在某类卡上失效。翻车方式很隐蔽——缓存没挤干净，测出来的算力反而“超过峰值”，还以为是代码写得好。

2.3 编译参数与第一次 baseline 跑通

编译只需要-lcublas，但几个开关值得固定下来：

nvcc -O3 -std=c++14 -arch=native -lineinfo \ src/bench_gemm.cu -lcublas -o build/bench_gemm ./build/bench_gemm --m 4096 --n 4096 --k 4096 --dtype fp32 --iters 50

-O3影响的是宿主侧循环，对 device 代码影响有限但没坏处；-arch=native让编译器按当前设备生成代码，省得手填架构号填错导致跑在兼容模式上；-lineinfo保留行号信息，配合 profiler 能直接定位到具体调用点。

命令行参数我一般这么定义，原因是它决定了后面 csv 的列：

第一次跑通的目标不是拿数字，而是确认三件事：程序能正常退出、CB()宏没有触发、耗时量级和手算的 FLOPs 对得上。如果 4096 的三次方跑出来只有几 GFLOPS，先别怀疑硬件，去看是不是把beta传成了未初始化的值，或者lda传成了k而不是m。

3. GEMM 基准的账要对得上：column-major、leading dimension 与 TFLOPS 规则

GEMM 是这套基准里最容易“数字好看但没意义”的部分。原因有三个：cuBLAS 只认列主序，而绝大多数人手里的矩阵是行主序；leading dimension 传小了会越界读，程序照样跑得飞快，因为它在读别人的内存；算力的分子分母如果口径不统一，和别人的数据根本没法比。

3.1 行主序的矩阵怎么交给列主序的 cuBLAS

cuBLAS 的约定是列主序，内存里连续的那一维是行下标。行主序的A(m×k)在列主序视角下等价于A_cm(k×m) = A_rm^T。想算C_rm = A_rm * B_rm，把它整体转置一次：C_rm^T = B_rm^T * A_rm^T，左右两边换位。落到调用上，就是把 M 和 N 互换、A 和 B 的位置互换：

// 内存按行主序存放：A 是 m×k，B 是 k×n，C 是 m×n // 目标：C_rm = alpha * A_rm * B_rm + beta * C_rm // 列主序视角下 C_cm = B_cm * A_cm，所以 m/n 互换、A/B 换位 const float alpha = 1.0f, beta = 0.0f; CB(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, // 注意顺序是 N, M, K &alpha, dB, N, // 第一个操作数现在是 B，lda = n dA, K, // 第二个操作数是 A，ldb = k &beta, dC, N)); // ldc = n（行主序下 C 的行宽）

逻辑说明：这个变换不改变任何数值结果，只是让“内存里的行主序布局”和“cuBLAS 眼里的列主序布局”对齐，省掉了显式转置的拷贝开销。对 4096 规模的方阵，一次显式转置的拷贝代价足以让测出来的性能掉三成以上，属于典型的自己给自己制造瓶颈。

参数说明：lda是第一个操作数在主维度上的行数。当transa = CUBLAS_OP_N时，A按列主序是M × K，因此lda >= max(1, M)；当transa = CUBLAS_OP_T时lda >= max(1, K)。上面的调用里 A、B 换位后，各自的ld就是K和N。这三个ld只要有一个传小，cuBLAS 不会报错，它会老老实实按你给的步长去读，读到相邻矩阵的数据，结果错得毫无征兆，性能却正常。

3.2 TFLOPS 的分子分母：2MNK、缓存清理与迭代次数

一次M×N×K的 GEMM 要做M*N次点积，每次点积K次乘加，每次乘加算两个浮点操作：

double flops = 2.0 * (double)M * N * K; double tflops = flops / (ms * 1e-3) / 1e12; // ms 是 timeit 返回的均值 // 对照用的等效带宽（口径写清楚，别和 GEMV 的混用） double bytes = 4.0 * ((double)M * K + (double)K * N + (double)M * N); double gbps = bytes / (ms * 1e-3) / 1e9;

分子上的2是固定的，M*N*K也是固定的，真正会因为选择不同而产生分歧的是分母：算进去的是不是只有 kernel 时间、包不包含alpha/beta的额外乘法、有没有把C的读写算进去。我一般的做法是明确声明“只统计单次 cuBLAS 调用在设备侧的耗时，迭代取均值，剔除缓存清理成本”，并在 csv 里用一列记下flush=on/off。

L2 清理到底影响多大，取决于矩阵工作集和 L2 容量的比值。A、B、C三块加起来小于 L2 时，第二轮迭代开始数据就在片上，测出来的“算力”会显著高于真实场景。这类数字拿去汇报会被质疑，而且质疑得对。判据很简单：把flush打开再跑一遍，两次差得离谱，说明之前的数字就是缓存红利。

注意：打开flush后小矩阵的数字会掉得很明显，这不是退步，而是回到了“数据必须从显存重新搬一次”的真实情况。

3.3 三种精度入口：SGEMM、HGEMM 与 TF32 快路

同一套扫描脚本要能跑三种精度，接口层必须分开写，因为alpha/beta的类型跟着计算精度走。

// 1) FP32 基线 const float a_f = 1.0f, b_f = 0.0f; CB(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &a_f, dB, N, dA, K, &b_f, dC, N)); // 2) FP32 输入走 TF32 快速路径：先切换 math mode，仍然调 sgemm CB(cublasSetMathMode(handle, CUBLAS_TF32_TENSOR_OP_MATH)); CB(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &a_f, dB, N, dA, K, &b_f, dC, N)); // 3) FP16：alpha/beta 必须是 half 类型 __half a_h = __float2half(1.0f), b_h = __float2half(0.0f); CB(cublasHgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &a_h, dB16, N, dA16, K, &b_h, dC16, N)); // 4) 显式指定计算精度：cublasGemmEx，便于做 A/B 对照 CB(cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &a_f, dB, CUDA_R_32F, N, dA, CUDA_R_32F, K, &b_f, dC, CUDA_R_32F, N, CUBLAS_COMPUTE_32F_FAST_TF32, CUBLAS_GEMM_DEFAULT));

逻辑说明：第 2 种和第 4 种在数学上都是“FP32 输入、低精度累加或低精度乘法”，区别在于第 2 种是会话级开关，一次设置对后续所有sgemm生效，容易污染其他实验；第 4 种是调用级，作用域干净，做对比实验时我更倾向用它。第 3 种的坑在alpha/beta的类型，传成float指针会读出错误的位模式。

参数说明：CUBLAS_COMPUTE_32F是纯 FP32 累加，CUBLAS_COMPUTE_32F_FAST_TF32才允许走低精度乘法单元；CUBLAS_GEMM_DEFAULT是让库自己选内核，不要在这里手填算法号，除非你在做算法级对比。

三种精度的对照关系建议固定成一张表，长期放在结果目录里：

同一组 4096 方阵，FP32 基线和 TF32 路径的差距通常在数倍量级。如果跑出来只差几个百分点，先别高兴，去查 math mode 有没有真的生效——这是第 5 章第 3 条要展开的事。

4. GEMV 是带宽题不是算力题：算术强度、batched 调用与验收线

GEMM 的优化目标是逼近算力峰值，GEMV 的目标完全不同。把 GEMV 的 TFLOPS 拿去和 GEMM 比，结论一定是“慢得离谱”，但这个结论没有信息量。GEMV 的算术强度低到几乎只受显存带宽约束，衡量它的指标应该是 GB/s，而不是 TFLOPS。

4.1 为什么 GEMV 永远无法达到峰值算力

算术强度定义为浮点操作数除以搬运字节数。对y = A * x，A是m×n，一次计算要做2mn次浮点操作，至少要搬4mn字节的矩阵（FP32），于是：

算术强度 ≈ 2mn / (4 * (mn + m + n)) ≈ 0.5 FLOP/Byte

这个值几乎不随矩阵规模变化，因为m、n变大时分子分母同比放大。而硬件的“拐点”等于峰值算力除以峰值带宽，现代 GPU 上这个比值通常在几十 FLOP/Byte 的量级。0.5 和几十之间差了两个数量级，意味着 GEMV 无论怎么调，都不可能落在算力受限区。

所以 GEMV 的基准只有一件事值得做：把实测带宽和理论带宽的比值压到尽可能高。理论带宽可以用 deviceQuery 给出的显存位宽乘以等效频率再除以 8 估算，不需要额外工具。

4.2 cublasSgemv 的调用与 batched 版本写法

单条 GEMV 的接口很直接，坑集中在lda和incx/incy上：

// y = alpha * A * x + beta * y // A 为 m×n 列主序；OP_N 时 lda >= max(1, m) const float alpha = 1.0f, beta = 0.0f; CB(cublasSgemv(handle, CUBLAS_OP_N, M, N, &alpha, dA, M, // lda = m dx, 1, // incx &beta, dy, 1)); // incy

逻辑说明：cublasSgemv在内部会按m的大小走不同的 kernel，m较小、n较大时倾向于按行切分，反之按列切分。这个选择你控制不了，也不需要控制，但要知道它存在——同一段代码在m=64, n=65536和m=65536, n=64两种形状下性能曲线完全不同，扫描时两个方向都要跑。

参数说明：incx、incy是向量元素的步长，单位是元素个数，不是字节数；传 0 会被拒绝，传大于 1 的值会走非合并访问路径，性能会明显下降，除非你确实在做跨步取数。

当一个模型里同时有几百条 GEMV 时，逐条调用会被启动开销主导。两种替代写法，实测差异不小：

// 写法一：把 n=1 的 GEMM 用 strided batched 一次发出去 CB(cublasSgemmStridedBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, M, 1, N, // n 固定为 1 &alpha, dA, M, (long long)M * N, // 每个 batch 的 A 起点间隔 dB, N, (long long)N, // 每个 batch 的 x，长度 n &beta, dC, M, (long long)M, // 每个 batch 的 y，长度 m batch_count)); // 写法二：指针数组版，适合每个 batch 的数据不在连续显存里 CB(cublasSgemvBatched(handle, CUBLAS_OP_N, M, N, &alpha, dA_array, M, dx_array, 1, &beta, dy_array, 1, batch_count));

逻辑说明：SgemmStridedBatched要求每个 batch 的数据在显存里等距排列，用一个stride就能定位，这种布局在推理框架里最常见；SgemvBatched接受设备端的指针数组，灵活但多一次指针解引用。前者对连续布局更快，后者对离散布局是唯一选择。

参数说明：strideA传M * N是按“每个 batch 一个完整矩阵”算的，如果各 batch 共享同一个矩阵（多向量乘同一权重），strideA要传 0；这一点很容易漏，漏了就会读到越界地址。

4.3 用实测带宽占理论带宽的比例当验收线

GEMV 的结果表里，我固定只留三列：形状、耗时、实测带宽，外加一列“占理论带宽比例”。比例低于 60% 就说明还有优化空间，通常是向量化没打开或者 batch 太小；比例超过 100% 则说明数据被缓存兜住了，这个数字不能对外用。

批量扫描用 shell 就够，不必上框架：

for n in 1024 2048 4096 8192 16384; do ./build/bench_gemv --m 4096 --n $n --iters 200 --csv results/gemv.csv done

这里--iters给到 200 是因为 GEMV 单次耗时短，迭代少的话事件分辨率不够。扫描完成后用一段脚本把 csv 画成“带宽 vs 规模”的曲线，形状通常是一条先升后平的线：上升段是启动开销被摊薄，平台段就是真实带宽上限。

5. 排查清单：这类基准最容易出问题的 5 个地方

跑了几十轮之后，我把遇到的问题归成两类：一类是计时的问题，数字不可信但代码不报错；一类是参数的问题，数字好看但结果不对。两类都不会给你任何提示，只能靠固定流程去撞。

5.1 计时与同步相关的三个问题

现象一：小矩阵跑出来的算力超过硬件理论峰值。原因几乎总是缓存。工作集小于 L2 时，第二轮迭代开始A、B、C全在片上，搬运成本被抹掉，测出来的就不再是显存路径。解决方式是按l2CacheSize的两倍分配 flush 缓冲，每次迭代前写一遍，并把 flush 成本单独测出来扣掉；如果懒得做，至少要在结果里标注矩阵规模和“flush=off”，让别人知道这个数字的边界。

现象二：第一次迭代特别慢，后面突然就稳了。两个原因叠加：cuBLAS 句柄在首次调用时需要做内部初始化，GPU 时钟也需要持续的负载才能进入高频率状态。解决方式是把 warmup 次数显式提上来（我一般给 10 次），warmup 结果全部丢弃，并且在正式计时前插一次cudaDeviceSynchronize()，确保没有残留任务干扰起点。

现象三：同一份二进制连跑两次，结果差 15% 以上。常见原因是设备被其他进程共享、机箱散热导致降频、或者迭代总时长太短被抖动主导。解决方式是把单次迭代的总时长拉到 200 ms 以上，正式结果取三轮的中位数而不是单轮均值，并记录运行时的时钟频率；如果条件允许，独占设备跑，别一边编译一边测。

5.2 参数与规则相关的两个问题

现象四：FP32 和 TF32 只差几个百分点，明显不合理。原因是低精度路径没真的生效。用cublasSgemm时忘记cublasSetMathMode(handle, CUBLAS_TF32_TENSOR_OP_MATH)，或者用了cublasGemmEx但 compute type 填了CUBLAS_COMPUTE_32F，都会静默回落到纯 FP32。解决方式是把 compute type 打印到日志里，并固定一个对照实验：同一个形状分别跑COMPUTE_32F和COMPUTE_32F_FAST_TF32，两个数字应该明显分层，分不开就是没生效。

现象五：性能正常但数值不对，CPU 参考实现一比就错。最常见的是lda/ldb/ldc传小，或者行主序与列主序搞混。越界读不会触发任何错误，程序反而因为读到了更“近”的内存而跑得更快。解决方式是在扫描之前先做小尺寸校验：取m=7, n=5, k=3这种不规整的形状，用三重循环的 CPU 版本算出参考值，逐元素比对相对误差，误差在 1e-5 以内才允许进入大规模扫描。这个校验代码只要二十行，但能拦住后面几天的无效数据。

注意：这五个问题里，前三个影响数字可信度，后两个影响结论正确性。前三个靠流程解决，后两个只能靠交叉校验，没有捷径。

6. 让数字站得住：交叉验证、cublasLt 试探与结果表的写法

一套基准跑完，最尴尬的局面是有人问“你这个数字准吗”，而你只能回答“我觉得准”。要避免这个局面，我会在出报告前做两件事。

第一件是内部交叉验证。同一个矩阵乘，用cublasSgemm（纯 FP32 路径）和cublasGemmEx且 compute type 设为CUBLAS_COMPUTE_32F各跑一遍，两者在数学上等价，耗时差异应该落在 5% 以内；如果差出一截，说明其中一条路径选了不同的内核，这个差异本身值得记录。数值上再用 CPU 参考实现对小尺寸逐元素比对，避免“跑得快但算错了”这种最难查的情况。

第二件是 adjudicate 一下自己写 kernel 的必要性。cuBLAS 默认给的是通用启发式选择，cuBLAS 之上还有一层cublasLt接口，可以拿到候选算法列表和它们各自需要的 workspace 大小，逐个试一遍再取最优：

cublasLtMatmulHeuristicResult_t heur; int returned = 0; CB(cublasLtMatmulAlgoGetHeuristic(lt_handle, lt_op, &layout_a, &layout_b, &layout_c, &layout_d, &pref, &heur, 1, &returned)); // returned > 0 时 heur.algo 才是有效候选，别跳过这个判断

这段代码的价值不在能提速多少，而在于它能告诉你“通用启发式离最优解还有多远”。如果差距只有 3%，自研 kernel 的投入基本不值得；如果差 30% 以上，说明这个形状确实没被库覆盖好，那才是值得动手的信号。翻车点很明确：returned为 0 时heur是未初始化的，直接拿去调cublasLtMatmul会返回参数错误，一定要先判。

最后是结果表的列设计，我固定成下面这样，多一列都嫌累赘：

“占估算峰值”这一列是给自己看的，峰值用 SM 数乘以每 SM 每时钟的 FMA 数再乘频率估算就行，不必追求精确——它的作用是让你一眼看出离理论上限还差多远，而不是拿去当宣传数字。

我自己踩得最深的一次，是拿一份没开 flush 的小矩阵数据去论证某条推理路径“已经接近硬件极限”，后来被人用同样代码加上缓存清理重跑，数字掉了一半多。从那以后我的习惯是：任何小于 L2 容量的形状，结果表里必须带flush=on，否则数字一律不写进结论。基准代码谁都能写，难的是让报出来的每个数都能被另一个人用另一台机器复现，希望帮到你。

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