1.1 main 控制流程#

1
try { return run(parse_options(argc, argv)); }

先看 main()，这里先执行 parse_options，成功返回配置后才调用 run。

大致的主控制流程如下：

1
命令行文本
2
    ↓ parse_options：转换类型，检查参数约束
3
Options 配置
4
    ↓ run：选择路径并组织实验
5
    ├─ help / list / plan：输出信息后返回
6
    └─ GPU 实验：准备数据 → 启动核函数 → 校验 → 测量或采集

1.2 匿名命名空间#

作用：把实现细节留在本翻译单元

1
namespace
2
{
3
constexpr int kInvalidArguments = 2;
4
constexpr int kNotImplemented = 3;
5
constexpr int kValidationFailed = 4;
6
// Options、解析函数、validate、run 等也在其中。
7
}

没有名字的 namespace 即匿名命名空间。这里用它收纳本翻译单元（一个源文件经过预处理、展开头文件后交给编译器的内容）的内部实现，避免把通用名字（如 run、Options）作为跨文件接口暴露出去。

constexpr int 在这里声明编译期常量；用 kValidationFailed 表达返回原因，这样可以避免到处散落 return 4 ，增加可读性。

1.3 Options#

作用：把一次实验的配置组织成一个值，类似于配置模式。

1
struct Options
2
{
3
    int mode = 0;
4
    size_t m = 512, n = 512, k = 512;
5
    int warmup = 10;
6
    int iterations = 100;
7
    double atol = 1e-4, rtol = 1e-4;
8
    bool help = false, list = false, plan_only = false;
9
    bool check_only = false, profile = false;
10
};

与 class 类似，执行 Options options; 时，这里的各成员会使用相应默认值，不需要额外写一个只负责赋值的构造函数。

实验配置组合成 Options 后，解析器负责构造配置，执行函数通过 run(const Options& options) 读取配置。

注意：Options 本身没有强制所有约束。例如仍可手动写入 options.iterations = 0。当前程序依靠 parse_options 检查其负责的输入约束，run 再检查设备相关限制。这里其实体现了本次设计的核心：解析器负责构造配置，但是不保证所有的配置都是合法的，具体是否合法需要按照各自的业务在后续判断。

2.1 检测命令行字符串合规性#

1
size_t parse_unsigned(const std::string& text, const char* name)
2
{
3
    if (text.empty() || text.find_first_not_of("0123456789") != std::string::npos)
4
        throw std::invalid_argument(std::string(name) + " must be an unsigned integer");
5

6
    const auto value = std::stoull(text);
7
    if (value > std::numeric_limits<size_t>::max())
8
        throw std::invalid_argument(std::string(name) + " is too large");
9

10
    return static_cast<size_t>(value);
11
}

1
parse_unsigned("512", "M");
2
parse_unsigned("abc", "warmup"); // 报错会带上 warmup

1
const auto value = std::stoull(text);
2
// 此处相当于：
3
// const unsigned long long value = std::stoull(text);

1
std::numeric_limits<size_t>::max()

1
throw std::invalid_argument(std::string(name) + " is too large");

1
// main()...
2
catch (const std::exception& error)
3
{
4
    std::fprintf(stderr, "Error: %s\nUse --help for usage.\n", error.what());
5
    return kInvalidArguments;
6
}

1
int parse_int(const std::string& text, const char* name)
2
{
3
    const size_t value = parse_unsigned(text, name);
4
    if (value > static_cast<size_t>(std::numeric_limits<int>::max()))
5
        throw std::invalid_argument(std::string(name) + " is too large");
6
    return static_cast<int>(value);
7
}

1
double parse_tolerance(const std::string& text)
2
{
3
    size_t used = 0;
4
    const double value = std::stod(text, &used);
5
    if (used != text.size() || !std::isfinite(value) || value < 0)
6
        throw std::invalid_argument("tolerances must be finite and nonnegative");
7
    return value;
8
}

1
size_t used = 0;
2
double value = std::stod("0.01abc", &used);
3
// value 为 0.01，used 为 4。

2.2 矩阵大小限制#

1
size_t matrix_elements(size_t rows, size_t columns, size_t element_bytes)
2
{
3
    const size_t limit = std::numeric_limits<size_t>::max();
4
    if (rows == 0 || columns == 0 || rows > limit / columns)
5
        throw std::invalid_argument("matrix dimensions must be positive and fit size_t");
6
    const size_t count = rows * columns;
7
    if (count > limit / element_bytes)
8
        throw std::invalid_argument("matrix byte size does not fit size_t");
9
    return count;
10
}

WARNING

不能先乘再检查!

rows、columns 都是 size_t，乘积也按该无符号类型计算。超出范围时结果按模回绕，比较看到的已经不是数学上的真实乘积。

1
// 错误思路：
2
if (rows * columns > limit) { /* 报错 */ }

因此，无符号乘积不会通过 > size_t最大值 来暴露问题 标准草案：无符号整数运算

• 把乘法约束改写成除法约束

对正整数 columns，希望满足：rows × columns ≤ limit，该表达式等价于 rows ≤ floor(limit / columns)

因此，代码先检查 rows > limit / columns，通过后才计算乘积。整数除法自然向下取整，正好给出允许的最大行数。

• 短路求值保护除法

1
rows == 0 || columns == 0 || rows > limit / columns

这里条件的排列具有安全意义：当 columns == 0 时，内建 || 不再计算后面的除法，从而避免除零。检查前置条件必须发生在使用这个条件的运算之前。

1
matrix_elements(options.m, options.n, sizeof(double));

解析阶段对 C 的形状使用 sizeof(double)，对应 CPU double 参考结果的大小；GPU 的 float 输出需要的字节数更少。sizeof(T) 返回类型 T 占用的字节数，类型形式不需要创建对象。

NOTE

当前调用处都传入 sizeof(float) 或 sizeof(double)，保证 element_bytes > 0。函数自身没有检查第三个参数是否为零，因此不能假设它能安全接受任意 size_t 参数组合。

检查通过仅表示大小的算术表达没有超出 size_t 范围，并不保证主机内存、显存或容器容量足够，实际分配仍可能失败。

2.3 parse_options：从零散文本组装实验配置#

前面的解析函数检查一个数值，这个函数则负责识别参数角色，并检查一次实验的整体配置。

1
Options 默认值 + argv 文本序列
2
    ↓ 第一阶段：扫描参数
3
    ├─ 无值开关：设置布尔成员
4
    ├─ 带值选项：消费下一个参数，转换后设置成员
5
    └─ 位置参数：按原顺序暂存到 positional
6
    ↓ 第二阶段：解释位置参数并检查整体约束
7
    mode / M / N / K → 检查模式、矩阵大小、迭代次数、选项互斥
8
    ↓
9
按值返回 Options

1
Options options;
2
std::vector<std::string> positional;
3
for (int i = 1; i < argc; ++i)
4
{
5
    const std::string arg = argv[i];
6
    // 按参数角色处理 arg。
7
}

1
else if (arg == "--warmup" || arg == "--iterations" ||
2
         arg == "--atol" || arg == "--rtol")
3
{
4
    if (++i >= argc) throw std::invalid_argument(arg + " needs a value");
5
    if (arg == "--warmup") options.warmup = parse_int(argv[i], "warmup");
6
    else if (arg == "--iterations") options.iterations = parse_int(argv[i], "iterations");
7
    else if (arg == "--atol") options.atol = parse_tolerance(argv[i]);
8
    else options.rtol = parse_tolerance(argv[i]);
9
}

1
下标：  0       1    2       3        4    5    6         7
2
参数：sgemm     1   128   --warmup    2   64   32   --check-only

1
else if (arg.rfind("--", 0) == 0)
2
    throw std::invalid_argument("unknown option: " + arg);
3
else positional.push_back(arg);

1
if (positional.size() > 4)
2
    throw std::invalid_argument("expected: mode M N K");
3
if (!positional.empty()) options.mode = parse_int(positional[0], "mode");
4
gemm_variant(options.mode);
5
if (positional.size() > 1) options.m = parse_unsigned(positional[1], "M");
6
if (positional.size() > 2) options.n = parse_unsigned(positional[2], "N");
7
if (positional.size() > 3) options.k = parse_unsigned(positional[3], "K");

1
return options;

2.4 print_help：相邻字符串字面量的拼接#

1
std::puts(
2
    "Usage: sgemm ...\n"
3
    "Options:\n"
4
    "  --help  Show help\n");

这里相邻的普通字符串字面量会在翻译阶段合并成一个字符串；没有逗号，因此 puts 仍只接收一个参数。这种写法方便把长帮助文本按行排版，不需要运行时通过 std::string 的 + 拼接。标准草案：字符串字面量

源码换行本身不会成为输出换行，文本里的 \n 才负责换行；puts 还会在输出末尾额外写入一个换行。字面量相邻拼接只适用于这种语法，两个字符串变量不能直接并排书写来拼接。

本节待自测：若收到 sgemm 2 --warmup 0 64 --iterations 5 32 16，positional 最终有哪些字符串？哪些 Options 成员被覆盖？暂不要求运行 GPU，先按读取游标推导。

3. validate：数值判定与错误诊断

作用：比较已经位于主机内存的 GPU 输出与 CPU double 参考结果，统计误差并返回是否通过。数据拷贝与后续流程由调用者负责。

3.1 接口：借用输入，返回判定#

1
bool validate(const std::vector<float>& actual, const std::vector<double>& reference,
2
              const Options& options, const char* phase)

actual、reference 和 options 都通过 const 引用读取，沿用前文的借用方式。phase 是本次校验的阶段标签，例如 "before"、"after"，用于让日志说明结果来自哪次检查。

这里有三类输出：bool 返回是否通过；汇总日志报告错误规模和最严重程度；首个错误日志提供定位线索。数值不匹配时返回 false，由 run 决定是否停止实验并返回校验失败退出码。

这种分工让同一个校验函数可以在实验的不同阶段复用。

3.2 绝对误差与混合容差#

1
const double error = std::isfinite(actual[i])
2
    ? std::abs(static_cast<double>(actual[i]) - reference[i])
3
    : std::numeric_limits<double>::infinity();
4
const double allowed = options.atol + options.rtol * std::abs(reference[i]);

对有限的实际结果，单元素通过条件为：

1
error = |actual - reference|
2
allowed = atol + rtol × |reference|
3
error ≤ allowed

atol 是绝对容差，提供接近零时仍可使用的误差额度；rtol 是相对容差，使允许误差随参考值的量级增大。当前公式将两者相加。

以默认 atol=rtol=0.0001 为例，下面按公式计算：

相同的绝对误差，在不同量级上可能得到不同判定。参考值为零时，相对容差项也为零，绝对容差继续提供误差额度。

浮点累加的顺序和精度会影响结果，CPU 与 GPU 也可能采用不同的运算组织方式，因此需要根据算法需求选择合适的容差。NVIDIA：参考结果比较、浮点数精度

把实际值提升为 double，是在 double 精度下计算对照误差；它无法恢复 GPU float 计算已经舍入掉的信息。CPU double 结果提高了对照计算的精度，也仍属于有限精度计算。

3.3 条件运算符：构造诊断值#

1
条件 ? 条件为真时的表达式 : 条件为假时的表达式

?: 是条件运算符，会产生一个值，只有选中的分支会被求值。因此，它适合在声明 const 变量时，根据条件选择初始化值。标准草案：条件运算符

在 error 的初始化中，有限 actual 使用误差计算；非有限 actual 统一映射成正无穷，作为失败诊断值。std::numeric_limits<double>::infinity() 取得浮点类型的正无穷值。

为什么要显式处理 NaN？NaN 参与这里的 > 等有序比较时结果为 false，若只靠 error > allowed，可能漏掉非法结果。因此，后面的失败条件也显式包含 !std::isfinite(actual[i])。

3.4 scaled：误差占允许额度的多少倍#

1
const double scaled = !std::isfinite(error) ? std::numeric_limits<double>::infinity() :
2
    (allowed > 0 ? error / allowed :
3
     (error == 0 ? 0 : std::numeric_limits<double>::infinity()));

通常情况下，scaled = error / allowed。它是一个以容差为基准的无量纲比值：0.5 表示用掉一半允许额度，1 表示恰好到达边界，2 表示误差达到允许值的两倍。

这个指标用到了 atol 和 rtol，不能把它简单记为“相对误差”。

嵌套表达式可以按下面的分支理解。以下为讲解用等价展开，未修改源码；inf 代表正无穷：

1
const double inf = std::numeric_limits<double>::infinity();
2
double scaled;
3
if (!std::isfinite(error))
4
    scaled = inf;
5
else if (allowed > 0)
6
    scaled = error / allowed;
7
else if (error == 0)
8
    scaled = 0;
9
else
10
    scaled = inf;

允许 atol 和 rtol 同时为零，就需要明确规定零容差下的诊断行为：误差为零时记为 0；误差大于零时记为无穷。条件分支避免直接计算 0 / 0。

设计思维：除了定义通过条件，还应为异常值和边界输入定义可解释的诊断结果。

3.5 一次遍历，收集不同层次的信息#

1
max_absolute = std::max(max_absolute, error);
2
max_scaled = std::max(max_scaled, scaled);
3
if (!std::isfinite(actual[i]) || error > allowed)
4
{
5
    if (failed == 0) first_bad = i;
6
    ++failed;
7
}

两个最大值可能来自不同元素；first_bad 也不一定对应最大误差。它们分别服务于误差概况和快速定位。

if (failed == 0) 只在首次失败时记录下标，随后递增失败计数。循环继续扫描所有元素，才能获得完整数量和最大误差，而不是在首次失败后立即结束。

first_bad 初始为 0 并不表示“没有错误”；是否存在错误由 failed 决定。否则，下标 0 真正出错时就容易与“没找到”混淆。

这里的比较扫描时间复杂度为 O(M×N)，额外统计空间为 O(1)；构造 CPU GEMM 参考结果的 O(M×N×K) 开销属于 run 中的另一个步骤。

3.6 从线性下标定位到矩阵坐标#

1
first_bad / options.n  // row
2
first_bad % options.n  // col

对于行主序 C[M,N]：i = row * N + col，因此除以行宽 N 得到行号，取余得到列号。

例如 N=4、first_bad=6，对应零起始坐标 row=1、col=2。这让诊断信息可以直接关联到核函数的输出索引。

3.7 日志格式也属于工程接口#

1
Check | phase=... | status=... | checked=full | elements=... | failed=...
2
Mismatch | row=... | col=... | expected=... | actual=...

稳定的记录类型和 key=value 字段，同时方便人阅读和脚本提取。成功时输出汇总，失败时补充首个错误详情，避免一次打印整个错误矩阵。

源码中的格式说明符可以这样读：

g 根据数值选择普通或科学计数形式，精度表示有效数字位数，并非固定小数位数；默认还会去掉不必要的尾随零。printf 的格式说明与实参类型需要匹配，不能因为某些类型在当前平台同宽就随意混用。Microsoft：printf 格式说明

3.8 当前函数依赖的输入契约#

调用者需要保证两份数组形状匹配、长度为 M×N、N>0，并提供有限的参考结果及合理的有限非负容差。当前 run 按相同 count_c 分配 actual 和 reference，再产生参考数据。

validate 本身没有重新检查数组长度或 reference 的有限性。若将它提取到其他工程中复用，需要重新决定这些前提由调用者保证，还是由函数自行检查。

PASS 表示本次输入在指定容差下通过了完整输出比较，不等于证明所有输入都正确，也不替代越界、竞争和同步依赖方面的检查。

本节待自测：reference=0、atol=0、rtol=0 时，actual=0 与 actual=0.125 各自得到怎样的 error、allowed、scaled 和判定？先按分支推导，不要求运行 GPU。

4.1 用接口区分描述、计划与执行#

1
const auto& variant = gemm_variant(options.mode);
2
const auto launch = make_gemm_launch(options.mode, options.m, options.n);

variant 借用已有的 GemmVariant 描述，launch 则持有按当前形状计算得到的 GemmLaunch 值。这里的引用或值选择与被调用函数的返回接口相对应，不必把所有 auto 都写成引用。

将元数据、启动配置与真正执行拆开后，--plan 可以只计算并展示 grid、block 和 shared 大小；实际执行前再检查设备的线程数、网格大小和 shared 容量限制。昂贵的参考计算和资源分配安排在这些检查之后。

其中 cudaDeviceProp properties{}; 用空花括号把这个结构体初始化到零等初始状态，随后由 cudaGetDeviceProperties 填充。

4.2 准备可复现的输入与独立参考#

输入使用固定公式生成，A、B 的公式不同，数值包含正负号，并乘以 0.03125f，即可以用二进制精确表示的 1/32。这样可以复现实验输入，也避免统一常量输入掩盖部分索引错误。

CPU 参考计算采用 row → inner → col 的顺序：最内层 col 连续访问 B 的一行与 reference 的一行，提高局部性；double 用于参考累加。这一过程位于事件计时区间外。

run 通过 CudaBuffer<float> 持有显存，通过 .get() 借用设备指针交给 CUDA API。主机 vector 的 .data() 同样提供其连续存储的指针。这些调用提供底层接口所需的地址，资源所有权仍属于对应容器或封装对象，具体 RAII 实现留到头文件部分。

4.3 lambda：把已经绑定上下文的操作交给计时器#

1
auto gemm = [&]() {
2
    launch_gemm(options.mode, device_a.get(), device_b.get(), device_c.get(),
3
                options.m, options.n, options.k, launch);
4
};
5
auto check_result = [&](const char* phase) {
6
    CUDA_CHECK(cudaMemcpy(actual.data(), device_c.get(), count_c * sizeof(float), cudaMemcpyDeviceToHost));
7
    return validate(actual, reference, options, phase);
8
};

lambda 表达式产生一个具有编译器生成类型的闭包对象，auto 接住这个对象；gemm() 则调用它。标准草案：闭包类型

[&] 使这里需要捕获的外部局部变量按引用捕获；lambda 可以使用现有配置和缓冲区，无需把它们重新列为每次调用的参数。check_result 的 phase 是每次调用时传入的参数，它的 bool 返回类型由 return 表达式推导。标准草案：lambda 捕获

这里选择引用捕获也与所有权设计一致：CudaBuffer 禁止复制，直接将 gemm 改成 [=] 会尝试复制捕获的缓冲区对象，无法编译。引用捕获借用资源，不转移所有权，也不延长被捕获对象的生命周期；当前两个 lambda 都只在 run 内使用。

1
const float ms = timer.benchmark(gemm, options.warmup, options.iterations);

这个接口把职责拆开：run 定义“一次工作是什么”，计时器负责“预热几次、执行几次、怎样测量”。捕获型 lambda 携带上下文，计时器通过模板接收可调用对象，具体模板和引用语义在 cuda_timer.cuh 一节继续展开。

当前 GEMM 接口要求每次覆盖 C，结果只依赖 A、B。重复调用因此使用相同工作负载，不会因为持续累加旧 C 而改变数值含义。

4.4 输出哨兵：使漏写能够显现#

1
CUDA_CHECK(cudaMemset(device_c.get(), 0xff, count_c * sizeof(float)));

cudaMemset 按字节填充；每个 float 的四个字节都写成 0xff，形成位模式 0xffffffff，在 CUDA 使用的 IEEE 754 binary32 表示中是 NaN。它不是把每个 float 赋成数值 255。NVIDIA：cudaMemset

如果核函数遗漏某个输出，保留下来的 NaN 会被 validate 拒绝。若预先清零，遗漏位置恰好应为零时可能误通过。这体现了主动设计测试初始状态的思路。

当前代码只在首次执行前填充一次哨兵；后续 after 校验检查最终结果，不单独证明每一轮都完整写过输出。

4.5 将正确性检查放在测量边界两侧#

首次执行后先同步并检查 before，避免对当前输入下结果不合格的实现继续报告性能；普通计时和 profile 路径结束后再检查 after。

gemm 调用完成只表示主机端的启动操作返回，GPU 核函数可能尚未完成。当前代码在读取首次结果和关闭采集范围前使用 cudaDeviceSynchronize，在事件计时路径中由计时器等待 stop 事件完成。这些同步让资源使用、读回结果与测量边界对齐。

4.6 吞吐率：明确计算量与时间单位#

代码按 GEMM 的常用计数约定使用约 2×M×N×K FLOP：一次乘加计两次浮点运算。ms 是平均毫秒数，因此：

1
GFLOP/s = FLOP / (ms × 10^-3) / 10^9
2
        = FLOP / (ms × 10^6)

这解释了 flops / (static_cast<double>(ms) * 1e6) 的单位换算。结果在 after 校验通过、时间有限且为正后才输出。

本节待思考：若把 cudaMemcpy 与 validate 一起移进传给 benchmark 的 lambda，得到的时间还对应现在的工作负载吗？重点说明测量范围，而不是背语法。

1. 先读懂一条实验命令#

PowerShell 是解释并执行命令的 shell，也是一门脚本语言。终端是输入命令、显示输出的界面；.ps1 文件则保存 PowerShell 代码。

本节先掌握命令调用规则，脚本内部的管道、重定向和参数展开后续再讲。

1
PowerShell
2
    ├─ 直接启动 sgemm.exe → 程序内 parse_options 解释参数
3
    └─ 运行 run_benchmark.ps1 → PowerShell 绑定脚本参数
4
                                  ↓
5
                             脚本逐次启动 sgemm.exe

1
.\build\labs\sgemm\Release\sgemm.exe 0 127 193 65 --plan

1
Write-Output 2+3     # 输出文本 2+3
2
Write-Output (2+3)   # 输出数值 5

1
Get-Location
2
Set-Location 'D:\CG\Kernel_Optimization\cuda-labs'

1
.\labs\sgemm\scripts\run_benchmark.ps1 -Modes 0,2 -Shapes '3x5x7','17x19x23' -Repeat 1 -CheckOnly

1
$taskSize = 128
2
'M=$taskSize'     # 输出 M=$taskSize
3
"M=$taskSize"     # 输出 M=128

1
.\labs\sgemm\scripts\run_benchmark.ps1 -BuildDir 'build debug' -Modes 0 -CheckOnly

1
$taskExe = '.\build\labs\sgemm\Release\sgemm.exe'
2
& $taskExe 0 127 193 65 --plan

1
& '.\build\labs\sgemm\Release\sgemm.exe' 0 127 193 65 --plan

1
.\labs\sgemm\scripts\run_benchmark.ps1 `
2
    -Modes 0,2 `
3
    -Shapes '3x5x7','17x19x23' `
4
    -Repeat 1 `
5
    -CheckOnly

1
Get-Help .\labs\sgemm\scripts\run_benchmark.ps1 -Parameter *
2
.\build\labs\sgemm\Release\sgemm.exe --help