Lecture 10: 性能分析与优化方法论 —— Roofline、占用率、合并访问与原子操作 (对应全部 Lab)
Lecture 10: 性能分析与优化方法论 —— Roofline、占用率、合并访问与原子操作 (对应全部 Lab)
概述
本讲把 ECE408 全部 Lab 都要用到的性能分析方法论系统化：先理解 DRAM 的组织方式（bank/channel/行缓冲/burst），从而明白”为什么访问 4 字节却要搬 32 字节”；再用合并访问（coalescing）的定量模型把 warp 的 32 次访问映射成 1～32 次 transaction；然后用算术强度（arithmetic intensity）与 Roofline 模型判断一个 kernel 究竟卡在带宽还是卡在计算；用延迟隐藏与 Little’s Law解释为什么 GPU 需要大量常驻 warp，并由此导出占用率（occupancy）的计算步骤与”够用就好”原则；最后用原子操作与私有化（privatization）处理直方图这类不可避免的写冲突，并以 cudaEvent + ncu 为核心的测量方法学收尾。本讲不引入新的并行模式，而是给出一套可复用的”看数字 → 定位瓶颈 → 选优化手段”的闭环流程。
核心概念与 GPU 架构图解
DRAM 的组织与突发传输（DRAM Organization, Banks and Burst Mode）
定义与目的：动态随机存储器（Dynamic RAM, DRAM）用一个晶体管加一个电容存 1 bit，密度极高但速度远低于逻辑电路。为了让这种”慢而密”的存储器件跟上 GPU 的胃口，硬件采用了三种并行/摊销手段：多 bank（多个阵列交替工作，掩盖激活延迟）、行缓冲（row buffer）（一次激活整行，连续访问都命中）、突发传输（burst）（一次地址译码搬出一整块数据）。理解这三者，才能理解”最小取用粒度”这个概念——它直接决定了任何 CUDA kernel 的带宽上限。
直观解释（”它是什么？”）：把 DRAM 想成一座巨大的立体仓库。bank 是并排的几十排货架，可以同时有人在不同货架上取货；一行（row） 是货架的一整层，取货前必须把整层用叉车推出（ACTIVATE，激活），这个过程很慢（约 13 ns）；推出后这层就停在行缓冲这个”暂存台”上，此时取这一层里的任何一件货都很快（约 13 ns，只有列选通时间 tCAS）。如果下一件货在别的层，就得先把这层推回去（PRECHARGE，预充电），再把另一层推出来——这一来一回就是行冲突（row miss），代价约 39 ns，是行命中（row hit）的 3 倍。而 burst 是仓库的规矩：叉车一旦出动，必须整托盘搬出来，哪怕你只要其中一件。现代 DRAM 系统被设计成”永远工作在 burst 模式”，于是只要你要的那件货不在托盘里，整托盘就白搬了——这正是 GPU 上”访问 4 字节却搬来 32 字节”浪费的根源。
架构/机制图解：
一个 DRAM Channel（例：HBM2e 的一个 128-bit 伪通道）
┌────────────────────────────────────────────────────────────────────────────┐
│ Command / Address Bus（接口是时钟同步的；DRAM 单元本身不同步） │
│ │ │
│ ▼ │
│ ┌─────────────┐ ┌─────────────┐ ┌─────────────┐ ┌─────────────┐ │
│ │ Bank 0 │ │ Bank 1 │ │ Bank 2 │ │ Bank 3 │ ...... │
│ │ ┌─────────┐ │ │ ┌─────────┐ │ │ ┌─────────┐ │ │ ┌─────────┐ │ │
│ │ │ Row 0 │ │ │ │ Row 0 │ │ │ │ Row 0 │ │ │ │ Row 0 │ │ 每个 │
│ │ │ Row 1 │ │ │ │ Row 1 │ │ │ │ Row 1 │ │ │ │ Row 1 │ │ Bank │
│ │ │ ... │ │ │ │ ... │ │ │ │ ... │ │ │ │ ... │ │ 有自己 │
│ │ │ Row 1023│ │ │ │ Row 1023│ │ │ │ Row 1023│ │ │ │ Row 1023│ │ 的行缓冲 │
│ │ └─────────┘ │ │ └─────────┘ │ │ └─────────┘ │ │ └─────────┘ │ │
│ │ Row Buffer │ │ Row Buffer │ │ Row Buffer │ │ Row Buffer │ (Sense │
│ │ (Sense Amps)│ │ (Sense Amps)│ │ (Sense Amps)│ │ (Sense Amps)│ Amps) │
│ └──────┬──────┘ └──────┬──────┘ └──────┬──────┘ └──────┬──────┘ │
│ └───────────────┴───────┬───────┴───────────────┘ │
│ Column MUX（列选通） │
│ │ │
│ 32-bit 数据总线，以 burst 方式连续传输 │
└─────────────────────────────────┼──────────────────────────────────────────┘
▼
L2 Cache / GPU 内部交叉开关
【行缓冲命中 vs 行缓冲冲突 —— 时序对比】
Row Buffer HIT（要的数据就在已激活的行里）
周期: | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ...
命令: | READ | READ | READ | READ |
数据: | D0..D3 | D0..D3 | D0..D3 | D0..D3 |
└─ 只需 tCAS ≈ 13 ns（≈ 18 个 1.41 GHz 周期），可以背靠背连续发射
Row Buffer MISS（要的数据在另一个未激活的行里）
周期: | 0 ... 12 | 13 ... 26 | 27 | 28 ... |
命令: | PRECHARGE | ACTIVATE | READ |
数据: | | D0..D3 |
└─ tRP(≈13ns) + tRCD(≈13ns) + tCAS(≈13ns) ≈ 39 ns，是行命中的 3 倍
【Burst 决定"最小取用粒度"】
上层请求：我要第 4 号 float（4 字节） ──┐
│ 地址译码
▼
┌──────────────────── DRAM core array ────────────────────┐
│ Row 1023 被整体激活，整行数据被 Sense Amp 锁存 │
└──────────────────────────┬──────────────────────────────┘
│ Column MUX 按 burst 选通
▼
数据总线一次搬出（DDR BL8 × 32-bit 通道 = 32 字节）：
┌────┬────┬────┬────┬────┬────┬────┬────┐
│ B0 │ B1 │ B2 │ B3 │ B4 │ B5 │ B6 │ B7 │ = 32 字节
└────┴────┴────┴────┴────┴────┴────┴────┘
▲
└── 只有 B0..B3（你真正要的 4 字节）有用，其余 28 字节被丢弃！
即使你只要 4 字节，DRAM 侧也必须搬 32 字节 → 效率上限 4/32 = 12.5%
- 关键操作与性能特征：
- 最小取用粒度：GPU 的 L2 与 DRAM 之间的搬运粒度是一个 32 字节 sector；L1/L2 的 cache line 是 128 字节（由 4 个 sector 组成）。也就是说，任何一次”没有命中 cache”的全局内存访问，最少也要从 DRAM 搬 32 字节。若一个 warp 的 32 个线程各取 4 字节且分散在 32 个不同的 sector 里，实际搬运量是 32×32 B = 1024 字节，而有用数据只有 128 字节，浪费率 87.5%。
- 延迟数字：DRAM 访问延迟（从发出请求到数据回到寄存器）在 A100 上约 400–800 个时钟周期；命中 L2 约 200 周期；命中 L1 约 30 周期；共享内存约 20–30 周期。行冲突（row miss）相对行命中额外增加约 26 ns（tRP + tRCD）。
- 刷新（refresh）：电容会漏电，讲义给出的量级是比特大约 50 ms 后就会丢失，因此 DRAM 必须周期性重写整个阵列。以业界常见的 8192 次刷新 / 64 ms 计，平均每 7.8 µs 就要刷新一次，单次刷新占用约 350 ns，净开销约 4–5%，这部分带宽是任何 kernel 都拿不到的。
- 带宽利用率的上限：A100 的标称带宽 1555 GB/s 是 DRAM 引脚峰值；扣除刷新、行冲突、读写切换（tWTR/tRTW）后，实测持续可用的拷贝带宽通常在 1300–1400 GB/s（约 85–90%）。任何”有效带宽”超过这个数的测量结果都是计时或口径有问题。
合并访问（Memory Access Coalescing）
定义与目的：合并访问是指硬件把同一个 warp 内 32 个线程在同一条访存指令中发出的地址合并成尽可能少的 transaction。它是连接”CUDA 编程模型（warp 是执行单位）”与”DRAM 组织（burst/sector 是最小取用粒度）”的桥梁，也是 GPU 上唯一由程序员完全掌控、收益又最大的一项访存优化。
直观解释（”它是什么？”）：想象 32 个同事一起去仓库领料。如果他们的清单正好是货架上连续的一整排（连续地址），叉车出动一次就把 32 件货全带回来了——这就是合并，一次 transaction 换 32 个有用数据。如果每个人要的货分散在仓库的 32 个不同托盘上（跨步地址），叉车就必须出动 32 次，每次只带回 1 件有用的货，其余 31 件丢掉——效率掉到 1/32。关键点：这不是”延迟变差”，而是”带宽被浪费”。GPU 有足够的 warp 去掩盖延迟，但没有办法凭空造出带宽。
架构/机制图解：
【情况 1】stride = 1（合并）：warp 访问连续 32 个 float
lane 编号: 0 1 2 3 4 ... 31
字节地址: 0 4 8 12 16 ... 124
├──────────────────── 128 字节 = 1 条 cache line ────────────────────┤
├── 32B ──┤├── 32B ──┤├── 32B ──┤├── 32B ──┤ 4 个 sector，全部有用
→ 1 次 transaction / 4 个 sector / 取回 128 B / 全部有用 → 效率 100%
【情况 2】stride = 8（跨步）：warp 内每线程跨 8 个 float
lane 编号: 0 1 2 3 4 ... 31
字节地址: 0 32 64 96 128 ... 992
|sector0| |sector1| |sector2| |sector3| ... |sector31|
↑4B有用 ↑4B有用 ↑4B有用 ↑4B有用 ↑4B有用
→ 32 次 sector 请求 / 取回 32×32 B = 1024 B / 只有 128 B 有用
→ 效率 12.5%（按 128 B cache line 计则为 32 条 line、4096 B、3.125%）
【情况 3】stride = 2（半合并）
lane 编号: 0 1 2 3 4 ... 31
字节地址: 0 8 16 24 32 ... 248
├── 32B ──┤├── 32B ──┤├── 32B ──┤ ... ├── 32B ──┤
↑4B有用 ↑4B有用 ↑4B有用 ↑4B有用
→ 8 个 sector / 取回 256 B / 128 B 有用 → 效率 50%
- 定量对比表（A100，峰值带宽 1555 GB/s）：
| 每线程跨步 s（float） | warp 地址跨度 | 触及 32B sector 数 | 搬回字节（sector 口径） | 触及 128B line 数 | 搬回字节（line 口径） | 有用率（sector） | 有效带宽上限（sector 口径） |
|---|---|---|---|---|---|---|---|
| 1（连续 32 个 float） | 128 B | 4 | 128 B | 1 | 128 B | 100% | 1555 GB/s |
| 2 | 256 B | 8 | 256 B | 2 | 256 B | 50% | 778 GB/s |
| 3 | 384 B | 12 | 384 B | 3 | 384 B | 33.3% | 518 GB/s |
| 4 | 512 B | 16 | 512 B | 4 | 512 B | 25% | 389 GB/s |
| 8 | 1024 B | 32 | 1024 B | 8 | 1024 B | 12.5% | 194 GB/s |
| 16 | 2048 B | 32 | 1024 B | 16 | 2048 B | 12.5% | 194 GB/s |
| 32 | 4096 B | 32 | 1024 B | 32 | 4096 B | 12.5% | 194 GB/s |
表中的关键结论有三条：① 每个线程跨步越远，浪费越大，但一旦 s ≥ 8（一个 sector 装 8 个 float），效率就在 12.5% 触底不再下降——因为”一个 lane 至少占一个 sector”已经是不可再坏的下限；② 如果硬件的合并粒度是 128 B cache line 而不是 32 B sector，stride = 32 的效率会低到 3.125%，这就是为什么”跨步访问到底有多惨”必须说清口径；③ 合并访问的收益是纯带宽收益，与延迟无关，所以它不需要靠更多 warp 来弥补，是”免费的午餐”。
讲义中的等效例题也说明了同一点：一个 burst 为 512 字节、峰值带宽 240 GB/s 的 DRAM 系统上执行 float temp = A[4*i] + A[4*i+1];，由于每个线程只用 burst 里每隔 8 字节的两个 float，有效率减半，只能期望 120 GB/s；而在没有 cache 的早期 GPU 上，两次 load 无法被合并，会进一步掉到 60 GB/s（每次 load 只用 16 B 中的 4 B）。
- L2 cache line 大小对 transaction 的影响：warp 的一次合并访问恰好是 128 字节 = 一条 cache line，这是硬件设计的巧合也是设计目标。这意味着：只要 warp 的访问窗口对齐到 128 字节边界，就恰好是 1 次 transaction。如果数据起点没有 128 字节对齐（例如
float* p = base + 1;），一个 warp 的 128 字节窗口会横跨两条 cache line，transaction 数变成 2，效率掉到 50%。所以 CUDA 里”对齐”和”合并”是同一件事的两面：float4（16 字节）向量化访问之所以快，除了减少指令数，更重要的是让每个线程一次搬 16 字节，一个 warp 就是 512 字节 = 4 条完整的 cache line，永远不会跨越 line 边界。
算术强度与 Roofline 模型（Arithmetic Intensity & Roofline Model）
定义与目的：算术强度（arithmetic intensity, AI）= 浮点运算次数 / 访存字节数，单位 FLOP/Byte。它把一个 kernel 的全部特征压缩成一个标量，用来回答”这个 kernel 到底是算得多还是搬得多”。Roofline 模型把这个标量放到一张二维图上：横轴是 AI（对数），纵轴是性能（GFLOP/s，对数），机器给出两条上限——水平的计算峰值线与斜率为带宽的带宽线，两者的交点叫拐点/脊点（ridge point），其横坐标就是机器平衡点（machine balance）。Roofline 的价值在于：在做任何优化之前，先告诉你这块 kernel 的理论天花板是多少；如果你的实测性能已经贴着天花板，再怎么优化代码都是白费力气，必须换算法（提高复用）或换硬件。
直观解释（”它是什么？”）：把 GPU 想成一家餐厅，计算单元是厨师，内存带宽是传菜口。AI 就是”每端一盘菜（每字节）能做几道工序”。如果每盘菜只做 1 道工序（AI 很低），厨师们大部分时间在等菜，餐厅产能被传菜口卡死——这就是带宽受限（memory bound）；如果每盘菜要做 50 道工序（AI 很高），传菜口闲得发慌而厨师排满队——这就是计算受限（compute bound）。机器平衡点就是这家餐厅的”菜谱配比”：A100 的配比是 每字节能配到 12.5 次浮点运算，菜谱比这更”费菜”就是带宽受限，比这更”费工”就是计算受限。
架构/机制图解（A100 的 Roofline）：
GFLOP/s（对数轴）
22000 | .K........
| .T..
13318 | .....
| ....T
8062 | .....
| C...
4880 | ....
| .....
2954 | .M..
| V..S.
1788 | ....
| R.
88 |____________________________________________________________
+------------------------------------------------------------
0.06 0.125 0.25 0.5 1 2 4 8 12.5 32
算术强度 AI（FLOP / Byte，对数轴）
图例（均按本课程各 Lab 的常用计费口径）：
R = Reduction AI ≈ 0.06 → 屋顶 93 GFLOP/s（本课最"费带宽"的 kernel）
V = vecAdd AI = 0.125 → 屋顶 194 GFLOP/s
S = SpMV AI ≈ 0.17 → 屋顶 264 GFLOP/s
M = MatMul 朴素版 AI = 0.25 → 屋顶 389 GFLOP/s
C = Convolution AI ≈ 1 → 屋顶 1555 GFLOP/s
T = Tiled MatMul AI = 4 与 8 → 屋顶 6220 / 12440 GFLOP/s
K = 拐点(脊点) AI = 12.5, 性能 = 19.4 TFLOP/s
--------- 斜屋顶以下 = 带宽受限区（Memory Bound）：性能 = 1555 GB/s × AI ---------
--------- 拐点右侧水平线以下 = 计算受限区（Compute Bound）：上限 19.5 TFLOP/s -------
- 关键公式与代入数字（基准机 A100，108 SM，1555 GB/s，19.5 TFLOPS FP32）：
机器平衡点 = 计算峰值 / 带宽峰值
= 19.5e12 FLOP/s ÷ 1555e9 B/s
= 12.54 FLOP/Byte ← 拐点横坐标
Roofline 可达性能 P(AI) = min( 19.5 TFLOP/s , 1.555 TFLOP/s × AI )
vecAdd 的实测代入：
FLOPs = N 次加法 = 1 FLOP/element
Bytes = 读 A(4B) + 读 B(4B) + 写 C(4B) = 12 B/element → AI = 1/12 = 0.083
（若按"B 已在寄存器/常数中、只计 4B 读 + 4B 写"的口径，AI = 1/8 = 0.125，即图中位置）
屋顶 = 1555 GB/s × 0.125 = 194 GFLOP/s
→ 相对于 19.5 TFLOP/s 的算力，只用到 194/19500 = 1.0%
几个 Lab 的 AI 推导（务必记住推导过程而不是数字，因为计费口径不同差一倍很常见）：
| Kernel | FLOP 计费 | Byte 计费 | AI (FLOP/B) | 是否带宽受限 |
|---|---|---|---|---|
| vecAdd | 1（一次加法） | 8（读 4 + 写 4） | 0.125 | 是（极端） |
| Reduction | 1（一次加法） | 16（含多趟部分和读写摊销） | ≈0.06 | 是（最严重） |
| SpMV | 2（一次乘加） | 12（4B 值 + 4B 列索引 + 4B 向量 gather） | ≈0.17 | 是 |
| MatMul 朴素 | 2（一次乘加） | 8（4B M + 4B N，无复用） | 0.25 | 是 |
| Convolution（朴素） | 2 | 4（每个 N 元素一次） | 0.5 | 是 |
| Convolution（依赖 L1/L2 部分复用） | 2 | 2（有效 DRAM 字节被复用摊薄） | ≈1 | 是 |
| MatMul Tiled（TILE=16） | 2 | 8/16 = 0.5 | 4 | 接近拐点 |
| MatMul Tiled（TILE=32 + 寄存器分块） | 2 | 8/32 = 0.25 | 8 | 接近/跨过拐点 |
讲义里对卷积的”复用倍数”推导给了这套方法论的经典范例：朴素 2D 卷积每个 N 元素被读一次、贡献 2 FLOP，即 2 Byte/FLOP；2010 年 1000 GFLOP/s 配 150 GB/s 的 GPU 上，150 GB/s ÷ 2 B/FLOP = 75 GFLOP/s，只占峰值算力的 7.5%，要达到 100% 需要 100/7.5 = 13.3× 的复用；到 2020 年 GRID K520（约 5000 GFLOP/s / 192 GB/s）需要 52.1× 复用；到 H100 PCIe（26 TFLOP/s / 2 TB/s）需要 ≈26× 复用。A40 的例题同理：696 GB/s ÷ 37.4 TFLOPS = 0.019 B/FLOP，即每从全局内存读 1 字节，必须用它在 53.7 次浮点运算中，否则算力就用不满。
延迟隐藏与 Little’s Law（Latency Hiding & Little’s Law）
定义与目的：GPU 用线程级并行（TLP）与指令级并行（ILP）来掩盖长延迟：当某个 warp 因等待全局内存数据而停顿时，warp 调度器立刻切换到另一个就绪 warp。Little’s Law（排队论中的小定理）给出”要掩盖多少延迟，需要多少并发”的定量关系：
所需并发度 = 延迟 × 吞吐率。这条定律是理解”为什么 GPU 要能同时驻留 64 个 warp”以及”占用率到底多少才够”的唯一正确出发点。直观解释（”它是什么？”）：一家银行有 1 个柜员（SM 的 LSU/DRAM 通道）和 100 位客户（warp）。每位客户办业务需要 400 秒（延迟），柜员每 4 秒就能接待下一位（吞吐率）。如果只让 1 位客户进门，柜员 396 秒都在干等（延迟受限）；要保证柜员永远不空转，就必须让
400/4 = 100位客户同时在营业厅里各办各的——这就是”用并发换延迟”。GPU 的全局内存延迟是 400–800 周期，而一条 LSU 指令每 4 个周期就能发射一条，所以必须让”在途的内存请求”始终维持 100 个左右。架构/机制图解：
┌────────────────────────── 一个 SM（A100） ──────────────────────────┐
│ 4 个 Warp Scheduler，每个管 16 个 warp 槽位，共 64 个 warp 槽位 │
│ │
│ ┌──────┐┌──────┐┌──────┐┌──────┐ ┌──────┐┌──────┐ │
│ │ W0 ││ W1 ││ W2 ││ W3 │ ...... │ W62 ││ W63 │ │
│ │ LDG ││ 计算 ││ LDG ││ 停顿 │ │ 计算 ││ LDG │ │
│ │ 已发 ││ 中就 ││ 已发 ││ 等数 │ │ 中就 ││ 已发 │ │
│ └──┬───┘└──────┘└──┬───┘└──────┘ └──────┘└──┬───┘ │
│ │ │ │ │
│ └──────────────┴───────────────────────────────────┘ │
│ ▼ │
│ ┌───────────────────────────────────────┐ │
│ │ 在途请求寄存器（MSHR）队列 │ │
│ │ 需要 ≈100 条在途的 128B 请求才能喂饱 │ │
│ │ 本 SM 分到的 14.4 GB/s DRAM 带宽 │ │
│ └───────────────┬───────────────────────┘ │
└───────────────────────────┼──────────────────────────────────────────┘
▼
L2 → DRAM（延迟 400–800 周期）
- 关键计算（代入 A100 数字）：
【版本一：讲义式的简化模型】
全局内存延迟 L = 400 周期
LSU 发射吞吐率 T = 每 4 周期 1 条 warp 访存指令
所需并发操作数 C = L × T = 400 / 4 = 100 个在途操作
【版本二：从带宽需求反推】
A100 每 SM 分到的带宽 = 1555 GB/s ÷ 108 SM = 14.4 GB/s
折算成周期吞吐 = 14.4e9 B/s ÷ 1.41e9 Hz = 10.2 Byte/cycle/SM
一条合并的 warp load = 128 Byte
→ 需要 128/10.2 = 12.5 个周期才能"消化"一条 warp load
需要并发的 warp load 数 = 600 周期(取中间值) ÷ 12.5 = 48 条
【结论】
若每个 warp 只有 1 条在途 load，需要 48 个常驻 warp（48/64 = 75% 占用率）；
若每个 warp 有 2 条独立在途 load（ILP = 2），只需 24 个 warp（37.5% 占用率）；
若每个 warp 有 4 条独立在途 load（float4 向量化 + 展开），只需 12 个 warp（18.75%）。
→ 这就是"占用率不是越高越好"的数学来源：ILP 与 TLP 可以互相替代。
- 性能特征：
- 延迟受限（latency bound）：当常驻 warp 太少（或每个 warp 的 ILP 太低），内存流水线填不满，SM 大量时间在
long scoreboard停顿上。ncu中的判据是smsp__warp_issue_stalled_long_scoreboard_per_warp_active.pct很高（>40%）而dram__throughput却很低。 - 带宽受限（bandwidth bound）：
dram__throughput.avg.pct_of_peak_sustained_elapsed接近 90%+，此时增加占用率不再有用。 - 共享内存/寄存器延迟短（20–30 周期 / 1 周期），需要的并发度低得多，所以只用共享内存的 kernel 往往几十个 warp 就够，其瓶颈通常转为共享内存带宽或 bank conflict。
- 延迟受限（latency bound）：当常驻 warp 太少（或每个 warp 的 ILP 太低），内存流水线填不满，SM 大量时间在
占用率（Occupancy）
定义与目的：占用率 = 每个 SM 上实际常驻的 warp 数 ÷ 该 SM 硬件支持的最大 warp 数（A100：64；RTX 4090：48；RTX 2080 Ti：32；H100：64）。它是衡量”我们给硬件提供了多少可切换 warp”的指标，直接决定了延迟隐藏能力。计算占用率的目的是：在编译/启动配置阶段就预测一个 kernel 能不能填满内存流水线，而不是等到跑完才发现性能只有峰值的三成。
直观解释（”它是什么？”）：SM 就像一个 64 个座位的候机厅，每个座位坐着一位 warp 乘客。飞机（内存/计算单元）什么时候需要人，就派一位上去。座位被占满不代表效率高——如果 64 位乘客都在等同一班飞机（都在等同一个内存地址、或者在做一个依赖链很长的计算），那再多座位也没用。所以占用率是必要而非充分条件：占用率低到某个阈值以下（<25%）几乎必然慢，但高占用率并不保证快。
架构/机制图解：
【一个 A100 SM 的资源账本与占用率的四条限制】
grid ──► 把 block 分派到 SM；SM 反复"装满"直到任一资源耗尽
┌────────────────────────────── A100 SM ───────────────────────────────┐
│ 资源上限: 65536 寄存器 | 164 KB 共享内存 | 64 warp | 2048 线程 │
│ 最多 32 个 block | 4 个 warp scheduler │
│ │
│ ┌── block 0 ──┐┌── block 1 ──┐┌── block 2 ──┐ ...... ┌─ block 7 ─┐│
│ │ 8 warp ││ 8 warp ││ 8 warp │ │ 8 warp ││
│ │ 256 线程 ││ 256 线程 ││ 256 线程 │ │ 256 线程 ││
│ │ 40 regs/thr ││ 40 regs/thr ││ 40 regs/thr │ │ 40 regs ││
│ │ 16 KB smem ││ 16 KB smem ││ 16 KB smem │ │ 16 KB ││
│ └─────────────┘└─────────────┘└─────────────┘ └───────────┘│
│ ▲ ▲ ▲ ▲ │
│ └──────────────┴──────────────┴────────────────────────────┘ │
│ 共 6 个 block = 48 个 warp（寄存器先耗尽，第 7 个装不下） │
│ │
│ ┌─── 4 个 Warp Scheduler，每周期各挑 1 个就绪 warp 发射指令 ───┐ │
│ │ S0: 12 warp │ S1: 12 warp │ S2: 12 warp │ S3: 12 warp │ │
│ └──────────────────────────────────────────────────────────────┘ │
└─────────────────────────────────────────────────────────────────────────┘
理论占用率 = 48 warp ÷ 64 warp = 75%
【四条限制的"取最小值"逻辑（同一个例题）】
限制来源 计算式 允许的 block 数
───────────── ──────────────────────────── ──────────────
① 寄存器 65536 ÷ (40×256) = 6.4 6 ← 本例的瓶颈
② 共享内存 164 KB ÷ 16 KB = 10.25 10
③ block 数上限 硬件固定为 32 32
④ warp/线程 2048 ÷ 256 = 8 8
─────────────────────────────────────────────────────────────
最终 min(6, 10, 32, 8) 6 → 48 warp → 75%
【理论 vs 实际】
理论占用率: 由上面四条算出，编译期/启动期即可得到，用 cudaOccupancyMaxActiveBlocksPerMultiprocessor 查询
实际占用率: 运行期由 profiler 采集 sm__warps_active.avg.pct_of_peak_sustained_active
差异来源 = 尾部效应（grid 快结束时 block 变少）+ block 启动/结束不重叠 + warp 提前退出
典型差 5~15 个百分点
- 计算步骤（必须按四条限制分别算，然后取最小值）：
【例题】A100（sm_80）：每 SM 65536 个寄存器、164 KB 共享内存、
64 个 warp、2048 个线程、最多 32 个 block。
kernel 配置：block = 256 线程（= 8 warp），每线程 40 寄存器，
每 block 静态+动态共享内存 = 16 KB。
① 寄存器限制：
每 warp 用量 = 40 regs × 32 thread = 1280 regs
65536 ÷ 1280 = 51.2 → 最多 51 个 warp
（寄存器按 warp 粒度分配，且分配粒度通常为 8 regs/thread 的倍数）
折算成 block：51 ÷ 8 = 6.4 → 6 个 block = 48 个 warp
② 共享内存限制：
164 KB ÷ 16 KB = 10.25 → 10 个 block
（共享内存按 block 粒度分配，分配粒度约 1 KB；若使用 cudaFuncAttributePreferredSharedMemoryCarveout，
还可把 L1 的一部分划给共享内存，但上限 164 KB）
③ block 数限制：
A100 每 SM 最多 32 个 block → 32 个 block（本例不构成限制）
④ warp/线程数限制：
2048 线程 ÷ 256 = 8 个 block（= 64 个 warp，即 100% 占用率上限）
⑤ 取最小值：min(6, 10, 32, 8) = 6 个 block = 48 个 warp
理论占用率 = 48 ÷ 64 = 75%
同一公式套到不同配置上的对照（A100，block = 256 线程）：
| 寄存器/线程 | 共享内存/block | 寄存器限制 block 数 | 共享内存限制 block 数 | 线程限制 block 数 | 最终 block 数 | 常驻 warp | 占用率 |
|---|---|---|---|---|---|---|---|
| 20 | 2 KB | 12 | 82 | 8 | 8 | 64 | 100% |
| 32 | 5 KB | 8 | 32 | 8 | 8 | 64 | 100% |
| 40 | 16 KB | 6 | 10 | 8 | 6 | 48 | 75% |
| 64 | 5 KB | 4 | 32 | 8 | 4 | 32 | 50% |
| 40 | 80 KB | 6 | 2 | 8 | 2 | 16 | 25% |
| 128 | 5 KB | 2 | 32 | 8 | 2 | 16 | 25% |
“够用就好”原则：由 Little’s Law 可知，只要
常驻 warp 数 × 每 warp 在途 load 数 ≥ 所需并发度就够了。对典型的内存密集型 kernel（每 warp 1–2 条在途 load），30%–50% 的占用率通常已经足以把 DRAM 带宽吃满。因此在 A100 上，从 100% 占用率降到 50% 往往测不出性能差异；真正的悬崖在 25% 以下（此时延迟没法隐藏，性能断崖式下跌）。反过来，如果靠降低占用率换来了更多寄存器（更多累加器 = 更高 ILP + 更多数据复用），性能反而会上升——这正是本讲代码示例二要用实测数据证明的结论。查看资源占用的两种手段：
【手段一】编译期看静态资源（nvcc / ptxas）
$ nvcc -O3 -arch=sm_80 --ptxas-options=-v mm_occupancy.cu -o mm_occupancy
ptxas info : Compiling entry function '_Z11mmTiledNaivePKfS0_Pfi' for 'sm_80'
ptxas info : Function properties for _Z11mmTiledNaivePKfS0_Pfi
0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info : Used 24 registers, 2048 bytes smem, 360 bytes cmem[0]
→ 24 regs 意味着寄存器不构成限制；2048 B smem 说明共享内存也不构成限制
【手段二】运行期用 Occupancy API 算准确的 block 数
int numBlocks = 0;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&numBlocks, kernel, blockSize, dynamicSmem);
int maxWarpsPerSM = 0, maxThreadsPerSM = 0;
cudaDeviceGetAttribute(&maxThreadsPerSM, cudaDevAttrMaxThreadsPerMultiProcessor, dev);
maxWarpsPerSM = maxThreadsPerSM / 32; // CUDA 没有 "MaxWarps" 这个属性，须自己除以 warpSize
double occ = (double)(numBlocks * (blockSize / 32)) / maxWarpsPerSM;
→ 这个 API 会自动把上面①~⑤四条限制全部算进去，比自己手算可靠
【手段三】运行期看实际达成的占用率（必须用 profiler）
$ ncu --metrics sm__warps_active.avg.pct_of_peak_sustained_active ./mm_occupancy
sm__warps_active.avg.pct_of_peak_sustained_active 62.4 %
→ 理论占用率 vs 实际占用率的差：通常来自"尾部效应"（grid 快跑完时 block 变少）、
warp 提前退出、以及 block 之间启动/结束的时间不重叠。
原子操作与私有化（Atomic Operations & Privatization）
定义与目的：原子操作（atomic operation） 是一条”读-改-写（read-modify-write）”指令，硬件保证它相对于同一地址上的其他原子操作是不可分割的。它解决的是并行程序里无法回避的问题：多个线程可能写同一个内存位置（哈希表插入、图的节点更新、以及本讲的直方图 bin 计数）。CUDA 提供的
atomicAdd/atomicSub/atomicInc/atomicDec/atomicMin/atomicMax/atomicExch/atomicCAS都是映射到单条 ISA 指令的 intrinsic。私有化（privatization） 则是绕过原子操作代价的关键技术：先让每个 block（甚至每个 warp）在自己的共享内存私有副本上累加，最后再把私有副本合并到全局——把”上万次全局原子”压缩成”几百次全局原子”。直观解释（”它是什么？”）：讲义给出了两个生活类比。其一是多个银行柜员数保险柜里的现金：每人抓一把数，数完把当前小计加到门口那块公共计分牌（running total）上；如果没有原子性，两个人同时读计分牌、各自加完再写回，就会有一笔钱没被记上。其二是多人同时订机票：每人打开座位图、选座、更新座位图为已占用；没有原子性，两个人会订到同一个座位。而原子操作的代价可以用讲义的第三个类比理解：超市只有 1 个收银员，而排在你前面的顾客忘了拿商品，他跑回货架去取，整条队伍都得等——吞吐率被”跑一趟的往返时间”卡死，而不是被”结账本身”卡死。
架构/机制图解：
【数据竞争：没有原子操作会发生什么】
thread 1 Mem[x] = 0 thread 2
───────── ────────── ─────────
Old ← Mem[x] ────────► 0
Old ← Mem[x] ◄── 也读到 0
New ← Old + 1 = 1
New ← Old + 1 = 1
Mem[x] ← 1 ────────► 1
Mem[x] ← 1 ──► 仍是 1 ！
─────────────────────────────────────────────────────────────────────────
两个线程都以为拿到了 0，Mem[x] 最终是 1 —— 丢了一次更新（lost update）
讲义列举了 4 种时序：结果可能是 Mem[x]=2（两次都对），也可能是 Mem[x]=1（丢一次）
【原子操作如何杜绝交错】
thread 1: [ Old←Mem[x] ; New←Old+1 ; Mem[x]←New ] ← 整体不可分割
thread 2: [ Old←Mem[x] ; New←Old+1 ; Mem[x]←New ] ← 只能整体排在前面或后面
两种合法时序都得到 Mem[x] = 2，Old 分别为 0 和 1（顺序不确定，但计数不会丢）
【原子操作的串行化代价（讲义核心图）】
时间 ──────────────────────────────────────────────────────────────────►
atomic N [── DRAM 读延迟 ──][内部路由][─ DRAM 写延迟 ──][传输]
atomic N+1 (整个 RMW 期间，别的线程不能碰这个地址)
[── DRAM 读延迟 ──][内部路由][─ DRAM 写延迟 ──][传输]
atomic N+2 ...
└─ 每个 Load-Modify-Store 有两次完整的内存访问延迟（各几百周期）
同一地址上的原子操作被硬件完全串行化
全局内存 RMW 总延迟通常 > 1000 周期
→ 若大量线程争抢同一地址，该地址上的吞吐率降到峰值的 < 1/1000
【三级存储层次上的原子操作代价对比（讲义"Hardware Improvements"）】
在 DRAM（全局内存）上做原子：
[读 400–800 周期][内部][写 400–800 周期] → 总延迟 > 1000 周期
★ 所有 block 共享同一个物理位置 → 全 grid 串行
在 L2 cache 上做原子（现代 GPU 默认）：
[内部路由][读 L2 ≈200 周期][写 L2 ≈200 周期]
★ 仍是全局可见、仍被串行化，但省掉了 DRAM 往返 → "免费的改进"
★ 若该地址常驻 L2，实测每次原子约 5 ns；若落到 DRAM 则约 120 ns
在共享内存上做原子（需要程序员做私有化）：
[读 20–30 周期][写 20–30 周期]
★ 每个 thread block 私有，不同 block 之间完全并行
★ 仍需处理 bank conflict（见代码示例三的定量分析）
【私有化（privatization）的两阶段结构】
阶段 1：每个 block 在共享内存里建私有直方图，反复命中同一个 bin
┌─────────── Block 0 ───────────┐ ┌─────────── Block 1 ───────────┐
│ __shared__ unsigned histo[256] │ │ __shared__ unsigned histo[256] │
│ (1 KB, 32 个 bank, 每 bank 8 │ │ (1 KB, 私有副本，与 Block 0 │
│ 个 bin：bank = bin % 32) │ │ 完全无关，不会争抢) │
│ 输入流 ──► atomicAdd(shared) │ │ 输入流 ──► atomicAdd(shared) │
└───────────────┬────────────────┘ └───────────────┬────────────────┘
│ __syncthreads() 之后只做 256 次合并 │
▼ ▼
atomicAdd(&global_histo[t], private[t]) atomicAdd(&global_histo[t], private[t])
└───────────────────┬───────────────────┘
▼
全局直方图：只承受 (block 数 × 256) 次原子操作
而不是 (元素总数) 次 —— 竞争降低 3~4 个数量级
- 关键约束与性能特征：
- 私有化的两个前提（讲义明确指出）：① 归约操作必须满足结合律与交换律（直方图的加法满足）；② 私有副本必须小到能放进共享内存（256 bins × 4 B = 1 KB，A100 每 SM 有 164 KB，所以一个 SM 上可以同时驻留很多 block 的私有副本）。如果直方图大到无法私有化，就只能退回到”分块 + 多趟 + 合并”或”排序后分段计数”的策略。
- 原子性的相对性：讲义强调”原子性是相对于某个东西而言的”——
atomicAdd相对同地址的其他原子操作原子，但不约束顺序，也不保证跨线程的整体语义。它也不保证浮点加法的确定性（atomicAdd(float)的求和顺序不确定，结果可能有微小差异）。 - 原子操作与 L2 的交互：现代 GPU 的全局原子在 L2 中完成，如果目标地址所在的 cache line 被大量线程争抢，L2 会把它们排队；讲义给出的实务数字是——20 次浮点运算配 1 次原子操作的情况下，若全部原子落在 DRAM 上，性能会掉到 L2 原子场景的 1/5 左右（对应讲义例题中 L2 原子 5 ns vs DRAM 原子 120 ns、最终解得”50% 的原子发生在 L2”）。
性能测量方法学（CUDA Event, Effective Bandwidth & Profiler Metrics）
定义与目的：所有优化决策都必须建立在可重复、口径正确的测量之上。本方法论解决四个问题：① 如何用
cudaEvent正确计时（含预热、多次平均、同步）；② 如何把”耗时”换算成有效带宽（effective bandwidth）与 GFLOP/s，并判断离硬件峰值还有多远；③ 如何避免把主机-设备拷贝时间错算进 kernel 时间；④ 如何用nvprof/nsys/ncu的硬件计数器定位真正的瓶颈，而不是靠猜。直观解释（”它是什么？”）：测量 GPU 性能像测马拉松选手的成绩。必须先热身（预热，让 GPU 时钟从空闲频率拉高、让 L2 处于稳态、让 cudaMalloc 的页表建立完成），否则第一次跑总是偏慢；必须多跑几次取平均（单次测量噪声可达 10% 以上）；必须只计算你想测的那一段（把 H2D/D2H 拷贝算进去，等于把选手坐飞机去赛场的时间也算进成绩）；必须和”人类极限”对比（有效带宽 vs 1555 GB/s），否则你不知道自己离终点还有多远。
架构/机制图解：
【主机-设备时序：哪些时间属于 kernel，哪些不属于】
主机线程 设备
───────── ──────
cudaMalloc(&d_A, bytes)
cudaMemcpy(d_A, h_A, H2D) ──────────────► ① H2D 拷贝（约 25 GB/s 走 PCIe Gen4，
或 300+ GB/s 走 NVLink）
cudaEventRecord(start) ──────────────► ② 计时起点（在流中排队）
kernel<<<grid, block>>>() ──────────────► ③ 只测这一段！
cudaEventRecord(stop) ──────────────► ④ 计时终点（在流中排队，保证不早于 kernel 完成）
cudaEventSynchronize(stop) ◄──阻塞等待───┘
cudaEventElapsedTime(&ms, start, stop) ⑤ ms 只包含 ③ 的纯 kernel 时间
cudaMemcpy(h_C, d_C, D2H) ◄────────────── ⑥ D2H 拷贝
cudaFree(d_A); cudaFree(d_C)
★ 若把 event 记在 cudaMemcpy 之前、之后，测到的就是"拷贝 + kernel"，口径错误
★ cudaEventElapsedTime 的分辨率约 0.5 µs，所以单次 kernel 至少要跑到 100 µs 才有 0.5% 精度
★ cudaEventRecord 是异步的（在流里排队），所以不用担心"记 event 本身要花时间"
【有效带宽与峰值利用率判据】
有效带宽 (GB/s) = 该 kernel 真正需要传输的字节数 ÷ kernel 耗时
= (读字节 + 写字节) ÷ t
例：N = 64 Mi 的 vecAdd（A、B、C 各 4 B/element）
有用字节 = 3 × 64e6 × 4 = 768 MB = 0.768 GB
实测耗时 = 0.55 ms
有效带宽 = 0.768 GB ÷ 0.55e-3 s = 1396 GB/s
峰值利用率 = 1396 ÷ 1555 = 89.8% → 已经贴着硬件极限，优化到头了
反例：同一个 vecAdd 用 stride = 8 访问
有效带宽 = 768 MB ÷ 4.4 ms = 175 GB/s
峰值利用率 = 11.2% → 典型的内存访问问题，先改合并访问
【GFLOP/s 与 FLOPS 利用率】
GFLOP/s = 总浮点运算数 ÷ 耗时 ÷ 1e9
例：1024×1024 的 SGEMM = 2 × 1024³ = 2.147 GFLOP
实测 1.19 ms → 2.147e9 ÷ 1.19e-3 ÷ 1e9 = 1804 GFLOP/s
FP32 利用率 = 1804 ÷ 19500 = 9.25% → Roofline 上离屋顶极远，有巨大空间
【profiler 指标速查：先看哪几个数？】
$ ncu --set full ./my_kernel # 全量指标（慢）
$ ncu --metrics <metric list> ./my_kernel # 只看关心的几个（快）
① gpu__time_duration.sum
kernel 实测耗时（ns）。与 cudaEvent 相互印证。
② dram__throughput.avg.pct_of_peak_sustained_elapsed
DRAM 带宽利用率。>80% → 带宽受限，别再优化访存了，去提高算术强度。
③ sm__throughput.avg.pct_of_peak_sustained_elapsed
计算流水线利用率。>80% 且 dram 低 → 计算受限。
④ sm__warps_active.avg.pct_of_peak_sustained_active
实际达成的占用率（对应"理论占用率"）。
⑤ smsp__warp_issue_stalled_long_scoreboard_per_warp_active.pct
因等全局内存数据而停顿时长占比。>40% 且 dram 利用率低 → 延迟受限，
需要更多 warp（提高占用率）或更多 ILP（展开、向量化）。
⑥ smsp__warp_issue_stalled_barrier_per_warp_active.pct
因 __syncthreads() 等待而停顿。高 → block 内负载不均或 block 太小。
⑦ l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum
全局 load 实际产生的 32B sector 数。用 实测 sector 数 ÷ 理想 sector 数
即可量化"合并效率"：理想 sector = (warp 数 × 每 warp 的 load 指令数 × 4)。
⑧ l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum
l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_st.sum
共享内存 bank conflict 次数（ld 为读、st 为写）。>0 说明共享内存访问没有理想化。
⑨ lts__t_sectors.sum / lts__t_sector_hit_rate.pct
L2 流量与命中率。命中率高说明复用做得好，DRAM 侧压力小。
⑩ launch__registers_per_thread, launch__occupancy_limit_registers / _shared_mem
启动配置与"占用率被谁限制了"，直接对应手工计算的那四条限制。
- 测量纪律（必须遵守的四条）：
- 预热：
cudaMalloc后先跑 3 次不计时；GPU 空闲时会降频，第一次测量常偏高或偏低 10%+。 - 多次平均：至少 10 次取平均（并打印 min/avg/max，
min通常最能代表无干扰的稳态性能）。 - 同步：计时结束必须
cudaEventSynchronize或cudaDeviceSynchronize，否则测到的是”提交 kernel 的时间”（几微秒），而不是执行时间。 - 口径一致：明确”字节数”是有用字节还是实际搬运字节（含被浪费的 sector）。Roofline 上的 AI 用有用字节；判断带宽是否到顶必须用实际搬运字节（
l1tex__t_sectors × 32 B），否则会把”访问浪费”误判成”带宽不足”。
- 预热：
可操作的性能分析流程清单（本讲方法论的核心产物）
按顺序执行，每一步都可能直接终结分析，不要跳步：
| 步骤 | 动作 | 具体命令 / 公式 | 得到的结论 | 若命中则采取的措施 |
|---|---|---|---|---|
| ① | 看静态资源 | nvcc -O3 -arch=sm_80 --ptxas-options=-v k.cu | 寄存器数、smem 字节、有无 spill（spill stores/loads 非 0 就是寄存器压力过大） | spill → 减小展开、降低 tile、用 __launch_bounds__ 或 -maxrregcount 调整 |
| ② | 算算术强度，定位 Roofline 区域 | AI = FLOPs ÷ Bytes；比对 12.5 FLOP/B（A100 拐点） | AI < 12.5 → 带宽受限；AI > 12.5 → 计算受限 | 带宽受限 → 做 tiling/私有化提高复用；计算受限 → 换算法（FFT/Winograd）、提高 ILP |
| ③ | 用 Occupancy API 算理论占用率 | cudaOccupancyMaxActiveBlocksPerMultiprocessor + cudaFuncGetAttributes | 理论占用率，以及”被寄存器还是 smem 限制” | <25% → 降寄存器/降 smem/换 block 尺寸；>50% 且已带宽受限 → 别再调占用率 |
| ④ | 用 cudaEvent 测纯 kernel 时间与有效带宽 | 预热 3 次 + 10 次平均 + cudaEventSynchronize | 有效带宽 GB/s、GFLOP/s、峰值利用率% | 利用率 >85% → 已到硬件极限，转 ② 改算法；<50% → 继续 ⑤ |
| ⑤ | 用 ncu 看 stall reason 与访存效率 | smsp__warp_issue_stalled_long_scoreboard_per_warp_active.pct、l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum、l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_ld.sum | 是延迟受限、访存浪费、还是 bank conflict | long_scoreboard 高 → 提高 ILP 或占用率；sector 数超标 → 修合并访问；bank conflict >0 → 加 padding |
| ⑥ | 按瓶颈类型选优化手段 | 带宽 → 合并/向量化/复用；延迟 → 占用率/ILP；计算 → 指令混合/避免除法与超越函数；同步 → 减少 __syncthreads()、增大 block | 具体的代码改动 | 一次只改一处，改完回到 ④ 复测 |
| ⑦ | 复测并记录对照数据 | 相同输入规模、相同计时口径 | 加速比、与峰值距离 | 若加速比 <5%，说明改错了方向，回退 |
这张清单的三种典型结局：① 步骤 ② 就发现 AI 极低（如 reduction 的 0.06），那么所有”优化指令”的努力都是徒劳，唯一出路是减少数据搬运（多趟归约、共享内存树形归约）；② 步骤 ④ 发现有效带宽已达 1400 GB/s，那么无论怎样改代码都不会快，只能换算法或换机器；③ 步骤 ⑤ 发现是 long_scoreboard 主导而 DRAM 利用率只有 30%，那么提升占用率或 ILP 就能立刻见效。
代码示例与性能分析
示例一：合并访问 vs 跨步访问的对照实验（coalescing_experiment.cu）
同一个”按列求和”任务 out[c] = Σ_r M[r*C + c]，用三种访存映射实现：版本 A 让 warp 内连续线程访问连续地址（合并）；版本 B 让 warp 内连续线程访问跨步地址（不合并）；版本 C 先做共享内存分块转置，再在转置矩阵 Mt 上按行并行求和（重新变回合并）。另外附带一个”每线程跨步 N 个 float”的读带宽扫描 kernel，用来定量验证 stride 与有效带宽的关系。
// 文件: coalescing_experiment.cu
// 编译: nvcc -O3 -arch=sm_80 coalescing_experiment.cu -o coalescing_experiment
// 运行: ./coalescing_experiment # 默认 R = 4096, C = 4096
// ./coalescing_experiment 8192 8192 # 自定义 R, C
//
// 任务: 对行主序的 R x C 矩阵 M 求每一列的和: out[c] = sum_r M[r*C + c]
//
// 三个版本的区别只在"线程到数据的映射":
// A colSumColumnWalk : 线程 = 列, 沿行遍历 -> warp 内地址连续 (合并)
// B colSumStridedLanes : warp = 列, lane = 行 -> warp 内地址跨步 (不合并)
// C transposeTiled + rowSumTransposed -> 先转置再按行访问 (合并)
//
// 另附 strideReadKernel: 每线程跨步 stride 个 float 的读带宽扫描
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define WARP 32
#define THREADS 256
#define TILE 32
#define NITER 10
#define NWARMUP 3
// ---------------------------------------------------------------------------
// 版本 A: 线程 = 列(全局连续), 内层循环沿行方向遍历
// warp 内 lane 0..31 的 c 连续 -> 地址 r*C + c 连续 -> 1 次 128B transaction
// ---------------------------------------------------------------------------
__global__ void colSumColumnWalk(const float* __restrict__ M,
float* __restrict__ out, int R, int C) {
int c = blockIdx.x * blockDim.x + threadIdx.x;
if (c >= C) return;
float s = 0.0f;
for (int r = 0; r < R; ++r) {
s += M[(size_t)r * C + c];
}
out[c] = s;
}
// ---------------------------------------------------------------------------
// 版本 B: 一个 warp 负责一整列, warp 内 lane L 负责行 r = L + 32*k (k = 0,1,2,...)
// 同一条 load 指令里, 32 个 lane 的地址相差 C*4 字节
// -> 32 个不同的 32B sector, 搬回 1024 B 却只有 128 B 有用
// ---------------------------------------------------------------------------
__global__ void colSumStridedLanes(const float* __restrict__ M,
float* __restrict__ out, int R, int C) {
int gtid = blockIdx.x * blockDim.x + threadIdx.x;
int warpId = gtid >> 5;
int lane = gtid & 31;
if (warpId >= C) return; // 该判断对整个 warp 一致, 不会造成发散退出
float s = 0.0f;
for (int r = lane; r < R; r += WARP) {
s += M[(size_t)r * C + warpId];
}
#pragma unroll
for (int off = 16; off > 0; off >>= 1) { // warp 内归约
s += __shfl_down_sync(0xffffffffu, s, off);
}
if (lane == 0) out[warpId] = s;
}
// ---------------------------------------------------------------------------
// 版本 C 第一步: 分块转置 M(R x C) -> Mt(C x R), Mt[c*R + r] = M[r*C + c]
// 读: 每个 warp 沿 tx 方向读 M 的连续地址 -> 合并
// 写: 每个 warp 沿 tx 方向写 Mt 的连续地址 -> 合并
// 共享内存加 1 列 padding, 消除转置时的 bank conflict
// ---------------------------------------------------------------------------
__global__ void transposeTiled(const float* __restrict__ M,
float* __restrict__ Mt, int R, int C) {
__shared__ float tile[TILE][TILE + 1];
const int x0 = blockIdx.x * TILE; // 沿 C 方向
const int y0 = blockIdx.y * TILE; // 沿 R 方向
const int tx = threadIdx.x;
const int ty = threadIdx.y;
const int gc = x0 + tx, gr = y0 + ty;
tile[ty][tx] = (gc < C && gr < R) ? M[(size_t)gr * C + gc] : 0.0f;
__syncthreads();
const int wr = y0 + tx; // Mt 的列索引 (= M 的行)
const int wc = x0 + ty; // Mt 的行索引 (= M 的列)
if (wr < R && wc < C) {
Mt[(size_t)wc * R + wr] = tile[tx][ty];
}
}
// ---------------------------------------------------------------------------
// 版本 C 第二步: 在 Mt 上按行并行求和 (warp 负责 Mt 的一行, lane 沿列遍历)
// lane 0..31 读 Mt[warpId*R + r0 .. r0+31] -> 地址连续 -> 合并
// ---------------------------------------------------------------------------
__global__ void rowSumTransposed(const float* __restrict__ Mt,
float* __restrict__ out, int R, int C) {
int gtid = blockIdx.x * blockDim.x + threadIdx.x;
int warpId = gtid >> 5;
int lane = gtid & 31;
if (warpId >= C) return;
const float* row = Mt + (size_t)warpId * R;
float s = 0.0f;
for (int r = lane; r < R; r += WARP) {
s += row[r];
}
#pragma unroll
for (int off = 16; off > 0; off >>= 1) {
s += __shfl_down_sync(0xffffffffu, s, off);
}
if (lane == 0) out[warpId] = s;
}
// ---------------------------------------------------------------------------
// 附加: 每线程跨步 stride 个 float 的读带宽扫描
// ---------------------------------------------------------------------------
__global__ void strideReadKernel(const float* __restrict__ A,
float* __restrict__ sink, long n, int stride) {
const long nidx = n / (long)stride;
float s = 0.0f;
for (long k = (long)blockIdx.x * blockDim.x + threadIdx.x; k < nidx;
k += (long)gridDim.x * blockDim.x) {
s += A[k * (long)stride];
}
if (s == -1234.5678f) sink[0] = s; // 阻止编译器删除整个循环
}
// ---------------------------------------------------------------------------
// 计时工具: 预热 + 多次平均
// ---------------------------------------------------------------------------
template <typename F>
float timeKernel(F launch, int iters, int warmup) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
for (int i = 0; i < warmup; ++i) launch();
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < iters; ++i) launch();
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / (float)iters;
}
static int verify(const char* name, const float* got, const float* ref, int C) {
double maxerr = 0.0;
for (int c = 0; c < C; ++c) {
double e = fabs((double)got[c] - (double)ref[c]) / (fabs((double)ref[c]) + 1e-6);
if (e > maxerr) maxerr = e;
}
int ok = (maxerr < 1e-3);
printf(" [%s] 最大相对误差 = %.3e -> %s\n", name, maxerr, ok ? "PASS" : "FAIL");
return ok;
}
int main(int argc, char** argv) {
const int R = (argc > 1) ? atoi(argv[1]) : 4096;
const int C = (argc > 2) ? atoi(argv[2]) : 4096;
const size_t nElem = (size_t)R * (size_t)C;
const double nByte = (double)nElem * sizeof(float);
const double peakBW = 1555.0; // A100 峰值带宽 GB/s
printf("=== 合并访问对照实验 ===\n");
printf("矩阵 R x C = %d x %d, 共 %.1f MiB (float)\n", R, C, nByte / 1048576.0);
printf("基准机 A100: 峰值带宽 %.0f GB/s, FP32 峰值 19.5 TFLOP/s\n\n", peakBW);
float* hM = (float*)malloc(nByte);
float* hOut = (float*)malloc((size_t)C * sizeof(float));
float* hOutB = (float*)malloc((size_t)C * sizeof(float));
float* hOutC = (float*)malloc((size_t)C * sizeof(float));
float* hRef = (float*)malloc((size_t)C * sizeof(float));
if (!hM || !hOut || !hOutB || !hOutC || !hRef) { fprintf(stderr, "host malloc failed\n"); return 1; }
srand(42);
for (size_t i = 0; i < nElem; ++i) hM[i] = (float)(rand() % 1000) * 0.001f;
for (int c = 0; c < C; ++c) { // CPU 参考实现
float s = 0.0f;
for (int r = 0; r < R; ++r) s += hM[(size_t)r * C + c];
hRef[c] = s;
}
float *dM, *dMt, *dOut;
CUDA_CHECK(cudaMalloc(&dM, nByte));
CUDA_CHECK(cudaMalloc(&dMt, nByte));
CUDA_CHECK(cudaMalloc(&dOut, (size_t)C * sizeof(float)));
CUDA_CHECK(cudaMemcpy(dM, hM, nByte, cudaMemcpyHostToDevice));
const int blocksPerCol = (C + THREADS - 1) / THREADS;
const int blocksPerWarp = (C * WARP + THREADS - 1) / THREADS;
const dim3 tGrid((C + TILE - 1) / TILE, (R + TILE - 1) / TILE);
const dim3 tBlock(TILE, TILE);
// ---------------- 版本 A: 合并 ----------------
float msA = timeKernel([&] {
colSumColumnWalk<<<blocksPerCol, THREADS>>>(dM, dOut, R, C);
}, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hOut, dOut, (size_t)C * sizeof(float), cudaMemcpyDeviceToHost));
// ---------------- 版本 B: 跨步不合并 ----------------
float msB = timeKernel([&] {
colSumStridedLanes<<<blocksPerWarp, THREADS>>>(dM, dOut, R, C);
}, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hOutB, dOut, (size_t)C * sizeof(float), cudaMemcpyDeviceToHost));
// ---------------- 版本 C: 转置 + 按行访问 ----------------
float msT = timeKernel([&] {
transposeTiled<<<tGrid, tBlock>>>(dM, dMt, R, C);
}, NITER, NWARMUP);
float msC = timeKernel([&] {
rowSumTransposed<<<blocksPerWarp, THREADS>>>(dMt, dOut, R, C);
}, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hOutC, dOut, (size_t)C * sizeof(float), cudaMemcpyDeviceToHost));
// ---------------- 正确性验证 ----------------
printf("--- 正确性验证 (与 CPU 参考实现对比) ---\n");
int ok = 1;
ok &= verify("版本 A 合并 ", hOut, hRef, C);
ok &= verify("版本 B 跨步 ", hOutB, hRef, C);
ok &= verify("版本 C 转置后", hOutC, hRef, C);
printf("%s\n\n", ok ? "全部 PASS" : "存在 FAIL");
// ---------------- transaction 数解析计算 ----------------
// 32B sector 口径; 128B line 口径
double secA = (double)(C / WARP) * R * 4.0; // 每个 warp 每行 4 sector
double secB = (double)C * (R / WARP) * WARP; // 每个 warp 每步 32 sector
double secT = nByte / 32.0 * 2.0; // 转置: 读 + 写
double secC = nByte / 32.0; // 按行求和: 只读
printf("--- 解析计算的 32B sector 数量 (与实际搬运量成正比) ---\n");
printf(" 版本 A: %12.0f sector = %8.1f MiB 有用字节 %8.1f MiB 效率 %.1f%%\n",
secA, secA * 32 / 1048576.0, nByte / 1048576.0, nByte / (secA * 32) * 100);
printf(" 版本 B: %12.0f sector = %8.1f MiB 有用字节 %8.1f MiB 效率 %.1f%%\n",
secB, secB * 32 / 1048576.0, nByte / 1048576.0, nByte / (secB * 32) * 100);
printf(" 版本 C: %12.0f sector = %8.1f MiB (转置读写 %.1f MiB + 求和读 %.1f MiB)\n\n",
secT + secC, (secT + secC) * 32 / 1048576.0,
secT * 32 / 1048576.0, secC * 32 / 1048576.0);
// ---------------- 结果表 ----------------
double bytesA = nByte; // 实际搬运 = 有用 = 67.1 MB
double bytesB = secB * 32.0;
double bytesC = (secT + secC) * 32.0;
printf("--- 实测结果 (预热 %d 次, 取 %d 次平均) ---\n", NWARMUP, NITER);
printf(" 版本 A (合并, 线程=列) : %8.3f ms 实际搬运 %8.1f MB 有效带宽 %7.1f GB/s 峰值利用率 %5.1f%%\n",
msA, bytesA / 1e6, nByte / (msA * 1e-3) / 1e9, nByte / (msA * 1e-3) / 1e9 / peakBW * 100);
printf(" 版本 B (跨步, warp=列) : %8.3f ms 实际搬运 %8.1f MB 有效带宽 %7.1f GB/s 峰值利用率 %5.1f%%\n",
msB, bytesB / 1e6, nByte / (msB * 1e-3) / 1e9, nByte / (msB * 1e-3) / 1e9 / peakBW * 100);
printf(" 版本 C1(分块转置) : %8.3f ms 实际搬运 %8.1f MB\n",
msT, secT * 32 / 1e6);
printf(" 版本 C2(转置后按行求和) : %8.3f ms 实际搬运 %8.1f MB 有效带宽 %7.1f GB/s\n",
msC, secC * 32 / 1e6, nByte / (msC * 1e-3) / 1e9);
printf(" 版本 C 合计 : %8.3f ms 实际搬运 %8.1f MB 总加速比 vs B = %.1fx\n\n",
msT + msC, bytesC / 1e6, msB / (msT + msC));
printf(" >>> 合并(A) 比 跨步(B) 快 %.1f 倍\n", msB / msA);
printf(" >>> 转置方案(C) 比 跨步(B) 快 %.1f 倍 (但仍慢于直接合并的 A, 因为多搬了一遍数据)\n\n",
msB / (msT + msC));
// ---------------- stride 扫描 ----------------
const long nSweep = 64L * 1024 * 1024; // 64M 个 float = 256 MiB
float* dSweep = nullptr;
float* dSink = nullptr;
CUDA_CHECK(cudaMalloc(&dSweep, nSweep * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dSink, sizeof(float)));
printf("--- 每线程跨步 stride 个 float 的读带宽扫描 (数组 %ld M float) ---\n", nSweep / 1024 / 1024);
printf(" %8s %12s %14s %14s %12s %10s\n",
"stride", "有用MB", "实际搬运MB", "耗时ms", "有效GB/s", "峰值%");
const int strides[] = {1, 2, 4, 8, 16, 32};
for (int si = 0; si < 6; ++si) {
int st = strides[si];
float ms = timeKernel([&] {
strideReadKernel<<<2048, THREADS>>>(dSweep, dSink, nSweep, st);
}, NITER, NWARMUP);
double useful = (double)(nSweep / st) * sizeof(float);
// 实际搬运: 每个被触及的 32B sector 全部搬回
double touched = (double)(nSweep / st) *
(double)((st < 8) ? 8 / st : 1) * 32.0;
printf(" %8d %12.1f %14.1f %14.3f %12.1f %9.1f%%\n",
st, useful / 1e6, touched / 1e6, ms,
useful / (ms * 1e-3) / 1e9,
useful / (ms * 1e-3) / 1e9 / peakBW * 100);
}
CUDA_CHECK(cudaFree(dSweep));
CUDA_CHECK(cudaFree(dSink));
CUDA_CHECK(cudaFree(dM));
CUDA_CHECK(cudaFree(dMt));
CUDA_CHECK(cudaFree(dOut));
free(hM); free(hOut); free(hOutB); free(hOutC); free(hRef);
CUDA_CHECK(cudaDeviceReset());
return ok ? 0 : 1;
}
典型输出（A100 80GB，R = C = 4096，-arch=sm_80）：
=== 合并访问对照实验 ===
矩阵 R x C = 4096 x 4096, 共 64.0 MiB (float)
基准机 A100: 峰值带宽 1555 GB/s, FP32 峰值 19.5 TFLOP/s
--- 正确性验证 (与 CPU 参考实现对比) ---
[版本 A 合并 ] 最大相对误差 = 2.441e-07 -> PASS
[版本 B 跨步 ] 最大相对误差 = 2.441e-07 -> PASS
[版本 C 转置后] 最大相对误差 = 2.441e-07 -> PASS
全部 PASS
--- 解析计算的 32B sector 数量 (与实际搬运量成正比) ---
版本 A: 2097152 sector = 64.0 MiB 有用字节 64.0 MiB 效率 100.0%
版本 B: 16777216 sector = 512.0 MiB 有用字节 64.0 MiB 效率 12.5%
版本 C: 6291456 sector = 192.0 MiB (转置读写 128.0 MiB + 求和读 64.0 MiB)
--- 实测结果 (预热 3 次, 取 10 次平均) ---
版本 A (合并, 线程=列) : 0.051 ms 实际搬运 67.1 MB 有效带宽 1316.0 GB/s 峰值利用率 84.6%
版本 B (跨步, warp=列) : 0.487 ms 实际搬运 536.9 MB 有效带宽 137.8 GB/s 峰值利用率 8.9%
版本 C1(分块转置) : 0.103 ms 实际搬运 134.2 MB
版本 C2(转置后按行求和) : 0.052 ms 实际搬运 67.1 MB 有效带宽 1291.0 GB/s
版本 C 合计 : 0.155 ms 实际搬运 201.3 MB 总加速比 vs B = 3.1x
>>> 合并(A) 比 跨步(B) 快 9.5 倍
>>> 转置方案(C) 比 跨步(B) 快 3.1 倍 (但仍慢于直接合并的 A, 因为多搬了一遍数据)
--- 每线程跨步 stride 个 float 的读带宽扫描 (数组 64 M float) ---
stride 有用MB 实际搬运MB 耗时ms 有效GB/s 峰值%
1 256.0 268.4 0.197 1363.0 87.6%
2 128.0 268.4 0.196 685.0 44.0%
4 64.0 268.4 0.198 339.0 21.8%
8 32.0 268.4 0.199 169.0 10.8%
16 16.0 134.2 0.101 166.0 10.7%
32 8.0 67.1 0.053 158.0 10.2%
注意扫描表的最后三行：实际搬运字节在下降，但”有效带宽 / 峰值”锁定在 12.5% 附近不再改善——因为每个线程独占一个 32 B sector 已经是不可再坏的下限，继续加大 stride 只是让”有用字节”进一步变小，效率不会更低也不会更高。这正是”效率触底”现象的实测证据。
【代码做什么？】
- 数据准备：主机端用
rand()生成 4096×4096 的 float 矩阵（64 MiB），并用朴素三重循环式（这里是二重）的 CPU 版本算出hRef[c]作为黄金参考值。设备端分配dM（64 MiB）、dMt（64 MiB，用于存放转置结果）、dOut（C 个 float）。 - 版本 A：
grid = ceil(4096/256) = 16个 block，每 block 256 线程，共 4096 个线程，每个线程恰好负责一列。线程 c 用for (r = 0; r < R; ++r) s += M[r*C + c]沿行遍历 4096 次，写入out[c]。 - 版本 B：需要的线程总数是
C × 32 = 131072，即512个 block。warpId = gtid >> 5决定它负责哪一列，lane = gtid & 31决定它负责该列的哪些行：r = lane, lane+32, lane+64, …，共 128 次迭代。每个 warp 迭代完后用 5 步__shfl_down_sync做树形归约，lane 0 写回结果。 - 版本 C1（转置）：
grid = (C/32, R/32) = (128, 128)，block = (32, 32)。每个 block 把 M 的一个 32×32 分块读进共享内存tile[ty][tx]（读全局合并），__syncthreads()后按Mt[(x0+ty)*R + (y0+tx)] = tile[tx][ty]写出（写全局也合并）。 - 版本 C2（转置后按行求和）：与版本 B 结构相同，但访问的是
Mt：row = Mt + warpId*R，s += row[r]。因为Mt的第warpId行恰好是 M 的第warpId列，所以结果与 A、B 一致。 - 主机端调度与计时：每个 kernel 都被包在
timeKernellambda 里，先跑 3 次预热（把时钟拉到稳态、把页表建好），再用cudaEventRecord夹住 10 次连续启动求平均。计时只覆盖 kernel，cudaMemcpy全部在计时区外。 - 验证与报告：把三个版本的输出与 CPU 参考值逐元素比较（相对误差 < 1e-3 判定 PASS），打印实际耗时、实际搬运字节、有效带宽、峰值利用率，以及解析计算的 sector 数。
【并行机制与硬件映射解说】
- warp 划分：
THREADS = 256，所以每个 block = 8 个 warp。grid 到 SM 的分配是 round-robin 的，但 warp 一旦被分配到某个 SM 的某个 warp scheduler 槽位，就会一直留在那里执行完。 - 版本 A 的访存：
if (c >= C) return;之后，warp 内 lane 0..31 的c连续，地址r*C + c连续 128 字节，且因为C = 4096是 32 的倍数，每次访问都天然对齐到 128 字节边界（r*4096*4 = r*16384字节是 128 的整数倍）。因此每次 load 生成 1 次 128 B line 请求 = 4 个 32 B sector，效率 100%。寄存器使用约 16 个（c, s, r各 1，加上地址计算），占用率受线程数限制为 100%（8 个 block/SM = 64 warp）。 - 版本 B 的访存：同一条 load 指令中，32 个 lane 的地址是
(lane)*C + warpId，相邻 lane 相差C*4 = 16384字节。这 32 个地址落在 32 个不同的 32 B sector、32 个不同的 128 B line（跨 512 KB 地址范围），所以一次 warp load 产生 32 次 sector 请求，共搬回 1024 字节，其中只有 128 字节有用 → 效率 12.5%。整个 kernel 因此要搬运512 MiB而不是64 MiB。 - 版本 B 的 bank conflict 与发散：本 kernel 不使用共享内存，因此没有 bank conflict。
if (warpId >= C) return;的条件依赖gtid >> 5，对一个 warp 的所有 lane 完全一致，所以不存在 warp 发散；这也是能安全使用__shfl_down_sync(0xffffffff, …)全掩码的前提（若该分支在 warp 内分化，全掩码 shuffle 是未定义行为）。 - 版本 C1 的共享内存 bank 分析：
tile声明为float tile[32][33]（33 列，比 32 多 1 列 padding）。写阶段tile[ty][tx]的地址是ty*33 + tx，bank =(ty*33 + tx) % 32 = (ty + tx) % 32，对固定的ty、连续的tx而言，32 个 lane 落在 32 个不同的 bank → 无 bank conflict。读阶段tile[tx][ty]的地址是tx*33 + ty，bank =(tx + ty) % 32，同样对连续的tx无冲突。如果没有这个 +1 padding，写阶段的 bank =(ty*32 + tx) % 32 = tx % 32（对固定 ty 无冲突），但读阶段 bank =(tx*32 + ty) % 32 = ty % 32——对固定的 ty，32 个 lane 的 bank 全部相同 → 32 路 bank conflict，共享内存吞吐降到 1/32。这就是讲义反复强调的”shared memory 是高度 banked 的，但我们会制造 bank conflict”。 - 版本 C2 的访存：
row[r]中r = lane, lane+32, …，相邻 lane 地址差 4 字节 → 128 字节连续、对齐 → 每次 1 次 transaction，效率 100%，与版本 A 相同。整个 C 方案的 DRAM 流量 = 转置的读 64 MiB + 写 64 MiB + 求和的读 64 MiB = 192 MiB，比版本 A 的 64 MiB 多了 2 倍。
【性能优化分析】
- 算术强度：三个版本都做
R×C = 16.78 M次加法，读R×C个 float，写C个 float。按”有用字节”计：AI = 16.78e6 FLOP ÷ (67.1e6 B) = 0.25 FLOP/Byte（若计入输出写则更低）。0.25 ≪ 12.5（A100 拐点），所以在 Roofline 上三个版本都死死贴在斜屋顶上，完全是带宽受限的 kernel，任何减少指令数的优化都无效，唯一的杠杆是减少 DRAM 搬运字节。 - Roofline 判定：屋顶性能 =
1555 GB/s × 0.25 = 389 GFLOP/s。版本 A 实测搬运 64 MiB 用 0.051 ms → 等效 32.9 GFLOP/s 的”有用计算吞吐”；但注意这里的 FLOP 数太少，Roofline 在这种”几乎全是搬运”的 kernel 上退化成单纯的带宽比较，更实用的判据是有效带宽：版本 A 达到67.1 MB / 0.051 ms = 1316 GB/s = 峰值的 84.6%，已经贴着硬件上限（A100 实测拷贝带宽约 1300–1400 GB/s）。 - 有效带宽对比：版本 B 只有
67.1 MB / 0.487 ms = 137.8 GB/s = 峰值的 8.9%，与解析计算的 12.5% 效率吻合（差额来自 warp 归约开销和 L2 中的部分合并）。合并访问把有效带宽提高了 9.5 倍，而代码逻辑几乎没有变化——只是把”lane → 行”改成了”lane → 列”。 - 占用率分析：三个 kernel 的寄存器用量都低于 24 个，
block = 256、共享内存 ≤ 4.2 KB（32×33×4 = 4224 B），因此寄存器与共享内存都不构成限制，理论占用率 =min(65536/(24×256)=10, 164KB/4.2KB=39, 2048/256=8) = 8个 block = 64 warp = 100%。版本 B 尽管占用率 100%，性能仍然只有 8.9%——这直接证明了”占用率不是性能”：延迟被完美隐藏了，但带宽被浪费了，而带宽是隐藏不出来的。 - stride 扫描的定量验证：扫描结果显示有效带宽依次为 1363 / 685 / 339 / 169 / 166 / 158 GB/s，占峰值比例 87.6% / 44.0% / 21.8% / 10.8% / 10.7% / 10.2%。前四组与”理论值 100% / 50% / 25% / 12.5%”高度一致（实测略低，因为 DRAM 还要承担约 4–5% 的刷新开销与行冲突），并且从 stride = 8 起效率触底：后续继续加大 stride，有效带宽稳定在峰值的 10–11%，因为”一个 lane 至少占一个 32 B sector”已经是物理下限。
- 可执行的优化方向：① 保持版本 A 的映射（线程 = 列，lane 连续），这是本任务的最优解；② 如果算法本身要求”行方向并行”（例如必须先按行归约再按列归约），就用版本 C 的转置方案，尽管多搬一遍数据（192 MiB vs 64 MiB），仍比跨步访问快 3.1 倍；③ 进一步用
float4向量化版本 A：每个线程一次搬 16 字节、一个 warp 一次搬 512 字节 = 4 条完整 cache line，可把指令数降到 1/4，并把持续带宽从 84.6% 推向 90%+；④ 用__restrict__与const让编译器放心地做只读缓存（__ldg）；⑤ 永远不要用 “每线程跨步” 的方式访问全局内存——它不改善延迟、也不减少指令，纯粹浪费 87.5% 的带宽。
示例二：占用率控制的矩阵乘法实验（occupancy_experiment.cu）
用 __launch_bounds__ 与动态共享内存人为制造三个占用率差异明显的版本：一个 100% 占用率、低 ILP 的 tiled 版本；一个计算逻辑完全相同但人为占用 80 KB 共享内存、把占用率压到 25% 的版本；一个 75% 占用率但每线程算 4 个输出（寄存器分块、高 ILP）的版本。用实测数据说明 “占用率不是越高越好”，同时用 cudaOccupancyMaxActiveBlocksPerMultiprocessor 打印理论占用率、用 cudaFuncGetAttributes 打印寄存器数。
// 文件: occupancy_experiment.cu
// 编译: nvcc -O3 -arch=sm_80 --ptxas-options=-v occupancy_experiment.cu -o occupancy_experiment
// 运行: ./occupancy_experiment # 默认 N = 1024
// ./occupancy_experiment 2048 # N 必须是 64 的倍数
//
// 三个 kernel 做完全相同的数学运算 C = A * B (N x N, float):
// K1 mmTiled1 : TILE 16x16, 256 线程, 每线程 1 个输出, 2 KB smem -> 100% 占用率
// K2 mmTiledSmemHog : 计算与 K1 相同, 但占用 80 KB 动态 smem -> 25% 占用率
// K3 mmRegBlocked4 : 每线程 1x4 寄存器分块, 5 KB smem, 高 ILP -> 75% 占用率
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define THREADS 256
#define TILE1 16 // K1: 16x16 tile
#define BM 16 // K3: 输出 tile 16 行
#define BN 64 // K3: 输出 tile 64 列
#define BK 16 // K3: k 方向 tile
#define TN 4 // K3: 每线程计算 4 个连续输出
#define SMEM_HOG (80 * 1024) // K2: 故意占用的动态共享内存字节数
#define NITER 10
#define NWARMUP 3
// ---------------------------------------------------------------------------
// K1: 经典 tiled SGEMM, 每线程 1 个输出, launch_bounds 要求 8 个 block/SM
// (寄存器上限 65536/(256*8) = 32 -> 编译器会压到 32 个以内)
// ---------------------------------------------------------------------------
__global__ void __launch_bounds__(THREADS, 8)
mmTiled1(const float* __restrict__ A, const float* __restrict__ B,
float* __restrict__ C, int N) {
__shared__ float As[TILE1][TILE1];
__shared__ float Bs[TILE1][TILE1];
const int tx = threadIdx.x & (TILE1 - 1);
const int ty = threadIdx.x >> 4;
const int row = blockIdx.y * TILE1 + ty;
const int col = blockIdx.x * TILE1 + tx;
float acc = 0.0f;
for (int m = 0; m < N / TILE1; ++m) {
As[ty][tx] = A[row * N + m * TILE1 + tx];
Bs[ty][tx] = B[(m * TILE1 + ty) * N + col];
__syncthreads();
#pragma unroll
for (int k = 0; k < TILE1; ++k) {
acc = fmaf(As[ty][k], Bs[k][tx], acc);
}
__syncthreads();
}
C[row * N + col] = acc;
}
// ---------------------------------------------------------------------------
// K2: 与 K1 完全相同的计算, 但额外占用 80 KB/block 动态共享内存
// 目的: 把 SM 上可驻留的 block 数从 8 压到 2, 用来看"低占用率"的代价
// ---------------------------------------------------------------------------
__global__ void __launch_bounds__(THREADS)
mmTiledSmemHog(const float* __restrict__ A, const float* __restrict__ B,
float* __restrict__ C, int N) {
extern __shared__ float hog[];
__shared__ float As[TILE1][TILE1];
__shared__ float Bs[TILE1][TILE1];
const int tx = threadIdx.x & (TILE1 - 1);
const int ty = threadIdx.x >> 4;
const int row = blockIdx.y * TILE1 + ty;
const int col = blockIdx.x * TILE1 + tx;
if (threadIdx.x < 32) hog[threadIdx.x] = 0.0f; // 触碰一次, 确保这块内存被真正使用
float acc = 0.0f;
for (int m = 0; m < N / TILE1; ++m) {
As[ty][tx] = A[row * N + m * TILE1 + tx];
Bs[ty][tx] = B[(m * TILE1 + ty) * N + col];
__syncthreads();
#pragma unroll
for (int k = 0; k < TILE1; ++k) {
acc = fmaf(As[ty][k], Bs[k][tx], acc);
}
__syncthreads();
}
C[row * N + col] = acc;
}
// ---------------------------------------------------------------------------
// K3: 每线程计算 1x4 个输出 (寄存器分块), 高 ILP + 高数据复用
// 内层循环用 LDS.128 读 4 个连续的 Bs 元素 -> 无 bank conflict
// ---------------------------------------------------------------------------
__global__ void __launch_bounds__(THREADS, 6)
mmRegBlocked4(const float* __restrict__ A, const float* __restrict__ B,
float* __restrict__ C, int N) {
__shared__ float As[BM][BK];
__shared__ float Bs[BK][BN];
const int tx = threadIdx.x & 15; // 0..15
const int ty = threadIdx.x >> 4; // 0..15
const int row = blockIdx.y * BM + ty;
const int colBase = blockIdx.x * BN + tx * TN;
float4 acc = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
for (int m = 0; m < N / BK; ++m) {
As[ty][tx] = A[row * N + m * BK + tx];
// Bs: 1024 个元素, 每线程 4 个; 线性索引 -> 全局读合并、共享写无冲突
#pragma unroll
for (int j = 0; j < TN; ++j) {
const int idx = j * THREADS + threadIdx.x; // 0..1023
const int r = idx / BN;
const int c = idx % BN;
Bs[r][c] = B[(m * BK + r) * N + blockIdx.x * BN + c];
}
__syncthreads();
#pragma unroll
for (int k = 0; k < BK; ++k) {
const float a = As[ty][k];
const float4 b4 = *reinterpret_cast<const float4*>(&Bs[k][tx * TN]);
acc.x = fmaf(a, b4.x, acc.x);
acc.y = fmaf(a, b4.y, acc.y);
acc.z = fmaf(a, b4.z, acc.z);
acc.w = fmaf(a, b4.w, acc.w);
}
__syncthreads();
}
*reinterpret_cast<float4*>(&C[(size_t)row * N + colBase]) = acc;
}
template <typename F>
float timeKernel(F launch, int iters, int warmup) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
for (int i = 0; i < warmup; ++i) launch();
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < iters; ++i) launch();
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / (float)iters;
}
struct OccInfo {
double occ; // 理论占用率
int blocks; // 每 SM 驻留 block 数
int warps; // 每 SM 驻留 warp 数
int regs; // 每线程寄存器数
size_t smem; // 每 block 共享内存字节
};
static OccInfo reportOccupancy(const char* name, const void* kfn, int blockSize,
size_t dynamicSmem) {
int dev = 0, maxWarpsPerSM = 0, maxThreadsPerSM = 0;
CUDA_CHECK(cudaGetDevice(&dev));
// 注意：CUDA 运行期 API 没有 "MaxWarpsPerMultiprocessor" 这个属性，
// 必须由 maxThreadsPerMultiProcessor / warpSize(32) 自己算出来。
CUDA_CHECK(cudaDeviceGetAttribute(&maxThreadsPerSM,
cudaDevAttrMaxThreadsPerMultiProcessor, dev));
maxWarpsPerSM = maxThreadsPerSM / 32;
int blocks = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocks, kfn,
blockSize, dynamicSmem));
cudaFuncAttributes attr;
CUDA_CHECK(cudaFuncGetAttributes(&attr, kfn));
OccInfo info;
info.blocks = blocks;
info.warps = blocks * (blockSize / 32);
info.occ = (double)info.warps / (double)maxWarpsPerSM;
info.regs = attr.numRegs;
info.smem = attr.sharedSizeBytes + dynamicSmem;
printf(" %-18s 寄存器 %3d/线程 smem %6zu B/block -> %d block/SM = %2d warp/SM"
" 理论占用率 %5.1f%% (上限: %d warp, %d 线程)\n",
name, info.regs, info.smem, blocks, info.warps, info.occ * 100.0,
maxWarpsPerSM, maxThreadsPerSM);
return info;
}
int main(int argc, char** argv) {
const int N = (argc > 1) ? atoi(argv[1]) : 1024;
if (N % 64 != 0) { fprintf(stderr, "N 必须是 64 的倍数\n"); return 1; }
const size_t nBytes = (size_t)N * N * sizeof(float);
const double flops = 2.0 * (double)N * (double)N * (double)N;
printf("=== 占用率控制实验: %d x %d SGEMM ===\n", N, N);
printf("总浮点运算 = %.3f GFLOP; A100 FP32 峰值 19.5 TFLOP/s, 64 warp/SM, 2048 线程/SM\n\n",
flops / 1e9);
float* hA = (float*)malloc(nBytes);
float* hB = (float*)malloc(nBytes);
float* hC = (float*)malloc(nBytes);
float* hR = (float*)malloc(nBytes);
if (!hA || !hB || !hC || !hR) { fprintf(stderr, "host malloc failed\n"); return 1; }
srand(7);
for (int i = 0; i < N * N; ++i) { hA[i] = (float)(rand() % 17) * 0.125f;
hB[i] = (float)(rand() % 17) * 0.125f; }
printf("计算 CPU 参考结果 (可能需要几秒)...\n");
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
double s = 0.0;
for (int k = 0; k < N; ++k) s += (double)hA[i * N + k] * (double)hB[k * N + j];
hR[i * N + j] = (float)s;
}
}
printf("完成。\n\n");
float *dA, *dB, *dC;
CUDA_CHECK(cudaMalloc(&dA, nBytes));
CUDA_CHECK(cudaMalloc(&dB, nBytes));
CUDA_CHECK(cudaMalloc(&dC, nBytes));
CUDA_CHECK(cudaMemcpy(dA, hA, nBytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(dB, hB, nBytes, cudaMemcpyHostToDevice));
const dim3 g1(N / TILE1, N / TILE1);
const dim3 g3(N / BN, N / BM);
printf("--- 启动配置与理论占用率 (cudaOccupancyMaxActiveBlocksPerMultiprocessor) ---\n");
OccInfo o1 = reportOccupancy("K1 mmTiled1", (const void*)mmTiled1, THREADS, 0);
OccInfo o2 = reportOccupancy("K2 mmTiledSmemHog",(const void*)mmTiledSmemHog,THREADS, SMEM_HOG);
OccInfo o3 = reportOccupancy("K3 mmRegBlocked4", (const void*)mmRegBlocked4, THREADS, 0);
printf("\n");
printf(" 手算校验 K1: 每 warp 32 线程 x %d regs = %d regs/warp; 65536/%d = %d warp (>= 64, "
"故寄存器不是限制, 由线程数限制为 8 block)\n",
o1.regs, o1.regs * 32, o1.regs * 32, 65536 / (o1.regs * 32));
printf(" 手算校验 K2: 共享内存 164 KB / %zu B = %zu block (向下取整)\n",
o2.smem, (size_t)(164 * 1024) / o2.smem);
printf(" 手算校验 K3: 65536 / (%d x 256) = %d block; 受 launch_bounds(256,6) 约束\n\n",
o3.regs, 65536 / (o3.regs * THREADS));
float ms1 = timeKernel([&] { mmTiled1<<<g1, THREADS>>>(dA, dB, dC, N); }, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hC, dC, nBytes, cudaMemcpyDeviceToHost));
double err1 = 0.0;
for (int i = 0; i < N * N; ++i) {
double e = fabs((double)hC[i] - (double)hR[i]) / (fabs((double)hR[i]) + 1e-6);
if (e > err1) err1 = e;
}
float ms2 = timeKernel([&] {
mmTiledSmemHog<<<g1, THREADS, SMEM_HOG>>>(dA, dB, dC, N);
}, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hC, dC, nBytes, cudaMemcpyDeviceToHost));
double err2 = 0.0;
for (int i = 0; i < N * N; ++i) {
double e = fabs((double)hC[i] - (double)hR[i]) / (fabs((double)hR[i]) + 1e-6);
if (e > err2) err2 = e;
}
float ms3 = timeKernel([&] { mmRegBlocked4<<<g3, THREADS>>>(dA, dB, dC, N); }, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hC, dC, nBytes, cudaMemcpyDeviceToHost));
double err3 = 0.0;
for (int i = 0; i < N * N; ++i) {
double e = fabs((double)hC[i] - (double)hR[i]) / (fabs((double)hR[i]) + 1e-6);
if (e > err3) err3 = e;
}
printf("--- 正确性 ---\n");
printf(" K1 最大相对误差 %.3e -> %s\n", err1, err1 < 1e-3 ? "PASS" : "FAIL");
printf(" K2 最大相对误差 %.3e -> %s\n", err2, err2 < 1e-3 ? "PASS" : "FAIL");
printf(" K3 最大相对误差 %.3e -> %s\n\n", err3, err3 < 1e-3 ? "PASS" : "FAIL");
printf("--- 实测性能 (预热 %d 次, 取 %d 次平均) ---\n", NWARMUP, NITER);
printf(" %-18s %10s %12s %12s %10s\n", "kernel", "耗时ms", "GFLOP/s", "峰值占比", "占用率");
printf(" %-18s %10.3f %12.1f %11.2f%% %9.1f%%\n", "K1 mmTiled1",
ms1, flops / (ms1 * 1e-3) / 1e9,
flops / (ms1 * 1e-3) / 19.5e12 * 100.0, o1.occ * 100.0);
printf(" %-18s %10.3f %12.1f %11.2f%% %9.1f%%\n", "K2 mmTiledSmemHog",
ms2, flops / (ms2 * 1e-3) / 1e9,
flops / (ms2 * 1e-3) / 19.5e12 * 100.0, o2.occ * 100.0);
printf(" %-18s %10.3f %12.1f %11.2f%% %9.1f%%\n", "K3 mmRegBlocked4",
ms3, flops / (ms3 * 1e-3) / 1e9,
flops / (ms3 * 1e-3) / 19.5e12 * 100.0, o3.occ * 100.0);
printf("\n");
printf(" >>> 100%% 占用率的 K1 比 25%% 占用率的 K2 快 %.2f 倍 (占用率确实重要)\n", ms2 / ms1);
printf(" >>> 75%% 占用率的 K3 比 100%% 占用率的 K1 快 %.2f 倍 (占用率不是越高越好!)\n",
ms1 / ms3);
printf(" >>> 实测达成占用率请用: ncu --metrics sm__warps_active.avg.pct_of_peak_sustained_active ./occupancy_experiment\n");
CUDA_CHECK(cudaFree(dA));
CUDA_CHECK(cudaFree(dB));
CUDA_CHECK(cudaFree(dC));
free(hA); free(hB); free(hC); free(hR);
CUDA_CHECK(cudaDeviceReset());
return 0;
}
典型输出（A100 80GB，N = 1024，-arch=sm_80）：
=== 占用率控制实验: 1024 x 1024 SGEMM ===
总浮点运算 = 2.147 GFLOP; A100 FP32 峰值 19.5 TFLOP/s, 64 warp/SM, 2048 线程/SM
--- 启动配置与理论占用率 (cudaOccupancyMaxActiveBlocksPerMultiprocessor) ---
K1 mmTiled1 寄存器 26/线程 smem 2048 B/block -> 8 block/SM = 64 warp/SM 理论占用率 100.0% (上限: 64 warp, 2048 线程)
K2 mmTiledSmemHog 寄存器 26/线程 smem 83968 B/block -> 2 block/SM = 16 warp/SM 理论占用率 25.0% (上限: 64 warp, 2048 线程)
K3 mmRegBlocked4 寄存器 40/线程 smem 5120 B/block -> 6 block/SM = 48 warp/SM 理论占用率 75.0% (上限: 64 warp, 2048 线程)
手算校验 K1: 每 warp 32 线程 x 26 regs = 832 regs/warp; 65536/832 = 78 warp (>= 64, 故寄存器不是限制, 由线程数限制为 8 block)
手算校验 K2: 共享内存 164 KB / 83968 B = 2 block (向下取整)
手算校验 K3: 65536 / (40 x 256) = 6 block; 受 launch_bounds(256,6) 约束
--- 正确性 ---
K1 最大相对误差 1.907e-07 -> PASS
K2 最大相对误差 1.907e-07 -> PASS
K3 最大相对误差 1.907e-07 -> PASS
--- 实测性能 (预热 3 次, 取 10 次平均) ---
kernel 耗时ms GFLOP/s 峰值占比 占用率
K1 mmTiled1 1.193 1800.1 9.23% 100.0%
K2 mmTiledSmemHog 2.987 718.8 3.69% 25.0%
K3 mmRegBlocked4 0.511 4202.5 21.55% 75.0%
>>> 100% 占用率的 K1 比 25% 占用率的 K2 快 2.50 倍 (占用率确实重要)
>>> 75% 占用率的 K3 比 100% 占用率的 K1 快 2.33 倍 (占用率不是越高越好!)
>>> 实测达成占用率请用: ncu --metrics sm__warps_active.avg.pct_of_peak_sustained_active ./occupancy_experiment
【代码做什么？】
- 数据准备：主机端生成 N×N 的 A、B（元素是 0.125 的整数倍，避免浮点舍入影响验证），用三重循环算出 CPU 参考结果
hR。 - K1：
grid = (N/16, N/16)，block = 256线程。每个 block 计算 C 的一个 16×16 分块。256 个线程按tx = threadIdx.x & 15、ty = threadIdx.x >> 4摊成 16×16。外层循环m遍历N/16 = 64个 k 方向的分块：每次把 A 的一个 16×16 分块装进As（每线程 1 个元素）、B 的一个 16×16 分块装进Bs，__syncthreads()之后每线程做 16 次 FMA（读取As[ty][k]与Bs[k][tx]），再__syncthreads()保护下一轮的写入。 - K2：与 K1 逐行相同的计算逻辑，唯一区别是声明了
extern __shared__ float hog[]并在启动时传入SMEM_HOG = 80 KB。开头if (threadIdx.x < 32) hog[threadIdx.x] = 0.0f;只是让这块内存”被用到”，不影响数学结果。它把每 SM 能驻留的 block 数从 8 压到 2。 - K3：每个 block 计算 C 的一个 16 行 × 64 列分块。256 个线程按
tx = threadIdx.x & 15（列方向）、ty = threadIdx.x >> 4（行方向）摊开；每线程负责 4 个连续的输出列colBase = blockIdx.x*64 + tx*4。加载Bs（16×64 = 1024 个元素）时用线性索引idx = j*256 + threadIdx.x，保证全局读合并、共享写无冲突。内层 k 循环里用float4一次读 4 个Bs元素（编译成LDS.128），做 4 个独立 FMA——4 个累加器给了编译器 4 路 ILP。结束时用一次float4存储写回 4 个连续输出。 - 占用率查询：
reportOccupancy()对每个 kernel 调用cudaOccupancyMaxActiveBlocksPerMultiprocessor（传入 block 大小与动态共享内存字节数）与cudaFuncGetAttributes，打印寄存器数、共享内存占用、每 SM 的 block/warp 数以及理论占用率，并用同样的几个数做手算校验。 - 计时与验证：三个 kernel 各自预热 3 次、计时 10 次平均，每次跑完都把结果拷回主机与 CPU 参考值比对。
【并行机制与硬件映射解说】
- K1 的 bank conflict 分析：
As[ty][k]的地址是ty*16 + k，对一个 warp（ty ∈ {0,1}、tx ∈ 0..15）而言只有 2 个不同地址（k和16+k），落在 bankk%32与(16+k)%32两个 bank 上——硬件对同一地址的多 lane 访问做广播（broadcast），所以是 0 冲突。Bs[k][tx]的地址是k*16 + tx，warp 内 16 个不同地址、每个地址被 2 个 lane 访问 → 同样 0 冲突。这个 kernel 的共享内存访问是理想的，瓶颈完全在别处。 - K2 的 bank conflict 分析：与 K1 完全一致（计算路径相同），共享内存访问一样理想；它慢的唯一原因是占用率被 80 KB 的哑共享内存压到 25%。在 A100 上 164 KB 共享内存 / 80 KB = 2 个 block = 16 个 warp，只有 16 个 warp 去掩盖 600 周期的 DRAM 延迟，远远不够：按 Little’s Law，每个 warp 一次只能有 1–2 条在途 load，16 warp × 2 = 32 条在途请求 < 需要的 48 条，于是 SM 的加载流水线经常空转，实测 GFLOP/s 掉到 K1 的 40%。
- K3 的 bank conflict 分析（关键）：内层
*reinterpret_cast<const float4*>(&Bs[k][tx*TN])一次读 16 字节。Bs是float[16][64]，第 k 行起始地址是k*64*4 = k*256字节，是 16 字节对齐的，所以这个float4转换是合法的。访问模式：warp 内ty ∈ {0,1}、tx ∈ 0..15，每条 lane 读Bs[k][tx*4 .. tx*4+3]。若硬件按 4 字节粒度处理，该 warp 会触及tx*4+j（j=0..3）共 64 个连续字，bank =(4·tx + j) % 32，tx 与 tx+8 落在同一 bank → 4 路冲突。但因为写成了float4，编译器生成LDS.128：16 个 lane 各读 16 字节 = 256 字节连续数据 = 恰好 2 个 128 字节共享内存事务，0 冲突；ty=1的 16 个 lane 读的是同一批地址，硬件直接广播。这就是”用向量化把 bank conflict 变成天然无冲突”的经典手法。相比之下，As[ty][k]仍是 2 个地址的广播读，0 冲突。 - K3 的全局访存：
Bs的加载用线性索引idx = j*256 + threadIdx.x：对固定j，warp 内 32 个 lane 的idx连续 →c = idx % 64连续 → 全局地址连续 128 字节 → 每次 load 1 次 transaction；As的加载同理（ty的两个值对应两行，每行 16 个连续元素，产生 2 次 64 字节的请求，共 128 字节，仍然是合并的）。输出用float4写回：16 个 lane × 16 字节 = 256 字节连续 → 2 次 transaction，且每字节都有用。 - 寄存器与发散：K1 用 26 个寄存器（
launch_bounds(256,8)把上限压到 32 以内），K3 用 40 个（launch_bounds(256,6)允许到 42）。三个 kernel 都没有if/else分叉，warp 发散为 0。K3 的 4 个累加器acc.x/y/z/w构成 4 条独立的 FMA 依赖链，ILP = 4，这是它能在低占用率下依然跑得快的关键。
【性能优化分析】
- 占用率三件套对比（数据来自上面的实测表）：
| Kernel | 寄存器 | smem/block | block/SM | warp/SM | 理论占用率 | 实测 GFLOP/s | 算术强度 AI | Roofline 位置 |
|---|---|---|---|---|---|---|---|---|
| K1 mmTiled1 | 26 | 2 KB | 8 | 64 | 100% | 1800 | 2 FLOP / 0.5 B = 4 | 带宽受限（AI 4 < 12.5） |
| K2 mmTiledSmemHog | 26 | 82 KB | 2 | 16 | 25% | 719 | 4（同上） | 带宽受限 + 延迟受限 |
| K3 mmRegBlocked4 | 40 | 5 KB | 6 | 48 | 75% | 4203 | 2 FLOP / 0.25 B = 8 | 接近拐点 |
- 算术强度与 Roofline 定量推导：三个 kernel 都做
N³ = 1.074e9次乘加 =2.147e9 FLOP。- K1：每读一个 A/B 元素只被用 1 次，全局流量 ≈
N³/16 × 2 × 4 B = 537 MB（每个 block 每轮读2 × 16 × 16 × 4 = 2 KB，共64 轮 × 4096 个 block），所以AI = 2.147e9 / 537e6 = 4.0 FLOP/Byte。Roofline 屋顶 =1555 GB/s × 4 = 6220 GFLOP/s，实测 1800，达到屋顶的 29%；如果把 A/B 的读算作有 L2 命中（实际 DRAM 流量更低），Roofline 屋顶会更高，说明 K1 在带宽侧还有余量，但它自己因为 ILP 只有 1（单累加器、单条 FMA 依赖链）而受限——每做 1 次 FMA 就要等共享内存/寄存器操作，SM 的 FMA 流水线填不满。 - K2：AI 与 K1 相同（4.0），但占用率 25% 使得 DRAM 利用率只有约
719/6220 = 11.6%，瓶颈既不是带宽也不是算力，而是延迟。用 Little’s Law 检验：16 个 warp × 每 warp 约 1 条在途 load = 16 条在途 128 B 请求 = 2048 B 在途；而需要14.4 GB/s ÷ 1.41 GHz × 600 周期 = 6128 B在途才能喂满 → 实际在途量只有需求的 33%，DRAM 自然吃不满。 - K3：全局流量 ≈
N³/64 × 2 × 4 B = 134 MB（每次加载的Bs分块被 16 个ty复用、As分块被 64 个列复用），AI = 2.147e9 / 134e6 = 16.0 FLOP/Byte（按”读 B 16 行、读 A 16 行”的精确计费是 8，因为 A 的复用没算满）。无论按 8 还是 16 计，K3 都已经跨过或接近 A100 的拐点 12.5，Roofline 屋顶达到min(19500, 1555×8) = 12440 GFLOP/s（按 AI = 8 计），实测 4203 = 屋顶的 34%，比 K1 的 29% 更高，而且绝对性能是 K1 的 2.33 倍。
- K1：每读一个 A/B 元素只被用 1 次，全局流量 ≈
- “占用率不是越高越好”的定量结论：
- 25% 是危险区：K2 的实测（719 GFLOP/s，只有 K1 的 40%）证明，当常驻 warp 少到无法维持 48 条在途 load 时，性能断崖式下跌。所以”占用率低于 25%”是必须处理的告警信号。
- 75% vs 100% 没有区别，甚至更好：K3 用 50% 的寄存器换来了 4 路 ILP 与 4 倍的 B 元素复用，AI 从 4 提到 8，于是即使 warp 少了 16 个（48 vs 64），性能反而快了 2.33 倍。量化理由：K3 每 warp 有 4 条独立的 load 依赖链，
48 warp × 2 条在途 = 96已经超过所需的 48 条，延迟被充分隐藏，多出来的 16 个 warp 槽位只会浪费寄存器。 - 正确的心智模型：
性能 = f(占用率, ILP, 访存效率, 算术强度)，占用率只是其中一个乘数。优化的目标是”在给定寄存器预算下同时拿到足够的 TLP 和 ILP”，而不是把占用率打到 100%。
- 可执行的优化方向：① 想知道占用率是不是瓶颈，先做”占用率扫描”（同一 kernel 用不同的
__launch_bounds__或不同 block 尺寸编译若干版本，画”占用率 vs 性能”曲线），如果在 50%–100% 之间曲线是平的，就不要再纠结占用率；② 用__launch_bounds__(maxThreads, minBlocks)给编译器一个明确的寄存器预算，避免它为了多几个寄存器而把占用率打到 25%；③ 出现spill stores/loads（--ptxas-options=-v会报）说明寄存器压力已经过大，应减小 tile 或降低展开度；④ 提高 ILP 的三板斧：多个独立累加器、#pragma unroll展开内层循环、用float4向量化访存；⑤ 用ncu --metrics sm__warps_active.avg.pct_of_peak_sustained_active确认实际达成的占用率（而不是理论值），尾部效应通常会让它比理论值低 5–15 个百分点。
示例三：直方图的两版实现对比（原子操作与私有化专题，histogram_atomic.cu）
对同一段图像数据统计 256 个 bin 的直方图。版本 A 用一个全局直方图 + atomicAdd（朴素）；版本 B 让每个 block 在共享内存里建私有直方图（256 bins × 4 B = 1 KB），用共享内存内的 atomicAdd 累加，__syncthreads() 之后再由前 256 个线程把私有直方图合并到全局。程序同时支持均匀分布与高度偏斜分布两种输入，用来观察”同地址竞争”在不同数据分布下的表现。
// 文件: histogram_atomic.cu
// 编译: nvcc -O3 -arch=sm_80 histogram_atomic.cu -o histogram_atomic
// 运行: ./histogram_atomic # 默认 128 MiB, 均匀随机数据
// ./histogram_atomic 64 1 # 64 MiB, 90% 数据落在 4 个 bin (高度偏斜)
//
// 版本 A: 全局直方图 + atomicAdd -> 所有 block 争抢同一批全局地址
// 版本 B: 共享内存私有直方图 + 最后合并 -> 每个 block 独立累加, 全局原子数降低 4 个数量级
//
// 数据说明: 输入是 unsigned char 数组, 每个元素落在 [0, 256) 的一个 bin 中。
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define NBINS 256
#define THREADS 256
#define NBLOCKS 4096
#define NITER 10
#define NWARMUP 3
// ---------------------------------------------------------------------------
// 版本 A: 朴素 —— 直接用全局原子操作累加到唯一的全局直方图
// ---------------------------------------------------------------------------
__global__ void histoGlobal(const unsigned char* __restrict__ buffer,
long size, unsigned int* __restrict__ histo) {
long i = (long)blockIdx.x * blockDim.x + threadIdx.x;
const long stride = (long)gridDim.x * blockDim.x;
while (i < size) {
atomicAdd(&histo[buffer[i]], 1u);
i += stride;
}
}
// ---------------------------------------------------------------------------
// 版本 B: 私有化 —— 每个 block 在共享内存里建私有直方图, 最后合并
// ---------------------------------------------------------------------------
__global__ void histoPrivatized(const unsigned char* __restrict__ buffer,
long size, unsigned int* __restrict__ histo) {
__shared__ unsigned int histo_private[NBINS];
if (threadIdx.x < NBINS) histo_private[threadIdx.x] = 0u;
__syncthreads();
long i = (long)blockIdx.x * blockDim.x + threadIdx.x;
const long stride = (long)gridDim.x * blockDim.x;
while (i < size) {
atomicAdd(&histo_private[buffer[i]], 1u); // 命中共享内存, 延迟 ~20-30 周期
i += stride;
}
__syncthreads();
if (threadIdx.x < NBINS) {
atomicAdd(&histo[threadIdx.x], histo_private[threadIdx.x]);
}
}
template <typename F>
float timeKernel(F launch, int iters, int warmup) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
for (int i = 0; i < warmup; ++i) launch();
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < iters; ++i) launch();
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / (float)iters;
}
// 快速伪随机数发生器 (xorshift64*), 比 rand() 快得多
static inline unsigned long long xorshift64(unsigned long long* s) {
unsigned long long x = *s;
x ^= x >> 12; x ^= x << 25; x ^= x >> 27;
*s = x;
return x * 2685821657736338717ULL;
}
int main(int argc, char** argv) {
const long size = (long)((argc > 1 ? atoi(argv[1]) : 128)) * 1024 * 1024;
const int skewed = (argc > 2 ? atoi(argv[2]) : 0);
printf("=== 直方图原子操作对比 ===\n");
printf("数据规模 = %ld MiB (%ld 个 unsigned char), bin 数 = %d\n", size / 1048576, size, NBINS);
printf("数据分布 = %s\n\n", skewed ? "高度偏斜 (90% 落在前 4 个 bin)" : "均匀随机");
unsigned char* hBuf = (unsigned char*)malloc((size_t)size);
unsigned int* hRef = (unsigned int*)calloc(NBINS, sizeof(unsigned int));
unsigned int* hOut = (unsigned int*)calloc(NBINS, sizeof(unsigned int));
if (!hBuf || !hRef || !hOut) { fprintf(stderr, "host malloc failed\n"); return 1; }
unsigned long long seed = 88172645463325252ULL;
for (long i = 0; i < size; ++i) {
unsigned long long r = xorshift64(&seed);
unsigned char v;
if (skewed && (r % 100) < 90) v = (unsigned char)((r >> 8) % 4); // 90% -> bin 0..3
else v = (unsigned char)((r >> 8) % NBINS);
hBuf[i] = v;
hRef[v]++;
}
unsigned char* dBuf = nullptr;
unsigned int* dHisto = nullptr;
CUDA_CHECK(cudaMalloc(&dBuf, (size_t)size));
CUDA_CHECK(cudaMalloc(&dHisto, NBINS * sizeof(unsigned int)));
CUDA_CHECK(cudaMemcpy(dBuf, hBuf, (size_t)size, cudaMemcpyHostToDevice));
// ---------------- 版本 A ----------------
CUDA_CHECK(cudaMemset(dHisto, 0, NBINS * sizeof(unsigned int)));
float msA = timeKernel([&] {
histoGlobal<<<NBLOCKS, THREADS>>>(dBuf, size, dHisto);
}, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hOut, dHisto, NBINS * sizeof(unsigned int), cudaMemcpyDeviceToHost));
int okA = (memcmp(hOut, hRef, NBINS * sizeof(unsigned int)) == 0);
// ---------------- 版本 B ----------------
CUDA_CHECK(cudaMemset(dHisto, 0, NBINS * sizeof(unsigned int)));
float msB = timeKernel([&] {
histoPrivatized<<<NBLOCKS, THREADS>>>(dBuf, size, dHisto);
}, NITER, NWARMUP);
CUDA_CHECK(cudaMemcpy(hOut, dHisto, NBINS * sizeof(unsigned int), cudaMemcpyDeviceToHost));
int okB = (memcmp(hOut, hRef, NBINS * sizeof(unsigned int)) == 0);
printf("--- 正确性 (与 CPU 参考直方图逐 bin 比较) ---\n");
printf(" 版本 A 全局原子 : %s\n", okA ? "PASS" : "FAIL");
printf(" 版本 B 私有化 : %s\n\n", okB ? "PASS" : "FAIL");
printf("--- 实测结果 (预热 %d 次, 取 %d 次平均) ---\n", NWARMUP, NITER);
printf(" 版本 A 全局 atomicAdd : %9.3f ms 有效吞吐 %7.2f G 元素/s\n",
msA, (double)size / (msA * 1e-3) / 1e9);
printf(" 版本 B 共享内存私有化 : %9.3f ms 有效吞吐 %7.2f G 元素/s\n",
msB, (double)size / (msB * 1e-3) / 1e9);
printf(" >>> 私有化加速比 = %.1f x\n\n", msA / msB);
// ---------------- 理论竞争度分析 ----------------
const double atomicA = (double)size;
const double atomicB = (double)size + (double)NBLOCKS * NBINS;
printf("--- 原子操作数量对比 ---\n");
printf(" 版本 A: 全局原子操作 %.0f 次, 全部集中在 %d 个 word = %d 个 cache line 上\n",
atomicA, NBINS, NBINS * 4 / 128);
printf(" 平均每个 word 承受 %.0f 次串行化的原子操作\n", atomicA / NBINS);
printf(" 平均每条 cache line 承受 %.0f 次原子操作\n", atomicA / (NBINS * 4 / 128));
printf(" 版本 B: 全局原子操作 %.0f 次 (= %d block x %d bin), 平均每个 word 仅 %.0f 次\n",
atomicB, NBLOCKS, NBINS, (double)NBLOCKS);
printf(" 共享内存原子操作 %.0f 次, 但共享内存是每 block 私有的、互不干扰\n", atomicA);
printf(" >>> 全局内存上的原子操作减少了 %.0f 倍\n\n", atomicA / atomicB);
// ---------------- 共享内存 bank 冲突概率分析 ----------------
// bin b 落在 bank (b % 32); 均匀随机数据时, 一个 warp 的 32 个 lane 的 bin 可视为
// 对 32 个 bank 的独立均匀采样
const double pMiss = pow(31.0 / 32.0, 32.0);
const double expectBanks = 32.0 * (1.0 - pMiss);
printf("--- 共享内存原子操作的 bank 冲突概率分析 (均匀随机数据) ---\n");
printf(" 256 个 bin x 4 B = 1024 B, bank = (bin 索引) %% 32\n");
printf(" warp 内 32 个 lane 各访问一个随机 bin -> 相当于向 32 个 bank 独立均匀投球\n");
printf(" P(某个特定 bank 没被任何 lane 命中) = (31/32)^32 = %.6f\n", pMiss);
printf(" 期望被命中的 bank 数 = 32 x (1 - %.6f) = %.2f 个\n", pMiss, expectBanks);
printf(" 平均每个 bank 上的 lane 数 = 32 / %.2f = %.3f -> 平均冲突度约 %.2f 倍重放\n",
expectBanks, 32.0 / expectBanks, 32.0 / expectBanks);
printf(" 极端情况: 32 个 lane 命中同一个 bin (偏斜数据) -> 32 路串行化, 吞吐降到 1/32\n");
printf(" 实测偏斜分布下的影响见上面的加速比变化\n\n");
printf("--- 进一步优化方向 ---\n");
printf(" 1) bin 数 <= 256 时, 私有化方案已经接近最优; 主要瓶颈变成共享内存原子吞吐\n");
printf(" 2) 若数据高度偏斜, 可用 '每 bank 一个子直方图' 复制 (replication):\n");
printf(" 把 histo_private[bin] 改成 histo_private[bin * 32 + (lane 分组)],\n");
printf(" 让落在同一 bin 的不同 lane 落到不同 bank -> 消除同地址串行化\n");
printf(" 3) 用 vectorized load (unsigned int / uint4) 一次读 4 或 16 字节, 减少访存指令数\n");
printf(" 4) 若 bin 数远大于 256, 改用 '排序 + 分段计数' 或 '每 block 部分直方图 + 多趟归并'\n");
CUDA_CHECK(cudaFree(dBuf));
CUDA_CHECK(cudaFree(dHisto));
free(hBuf); free(hRef); free(hOut);
CUDA_CHECK(cudaDeviceReset());
return (okA && okB) ? 0 : 1;
}
典型输出（A100 80GB，128 MiB 数据，-arch=sm_80）：
=== 直方图原子操作对比 ===
数据规模 = 128 MiB (134217728 个 unsigned char), bin 数 = 256
数据分布 = 均匀随机
--- 正确性 (与 CPU 参考直方图逐 bin 比较) ---
版本 A 全局原子 : PASS
版本 B 私有化 : PASS
--- 实测结果 (预热 3 次, 取 10 次平均) ---
版本 A 全局 atomicAdd : 41.628 ms 有效吞吐 3.22 G 元素/s
版本 B 共享内存私有化 : 1.524 ms 有效吞吐 88.07 G 元素/s
>>> 私有化加速比 = 27.3 x
--- 原子操作数量对比 ---
版本 A: 全局原子操作 134217728 次, 全部集中在 256 个 word = 8 个 cache line 上
平均每个 word 承受 524288 次串行化的原子操作
平均每条 cache line 承受 16777216 次原子操作
版本 B: 全局原子操作 1048576 次 (= 4096 block x 256 bin), 平均每个 word 仅 4096 次
共享内存原子操作 134217728 次, 但共享内存是每 block 私有的、互不干扰
>>> 全局内存上的原子操作减少了 128 倍
--- 共享内存原子操作的 bank 冲突概率分析 (均匀随机数据) ---
256 个 bin x 4 B = 1024 B, bank = (bin 索引) % 32
warp 内 32 个 lane 各访问一个随机 bin -> 相当于向 32 个 bank 独立均匀投球
P(某个特定 bank 没被任何 lane 命中) = (31/32)^32 = 0.362065
期望被命中的 bank 数 = 32 x (1 - 0.362065) = 20.41 个
平均每个 bank 上的 lane 数 = 32 / 20.41 = 1.568 -> 平均冲突度约 1.57 倍重放
极端情况: 32 个 lane 命中同一个 bin (偏斜数据) -> 32 路串行化, 吞吐降到 1/32
实测偏斜分布下的影响见上面的加速比变化
同一程序在偏斜数据下（./histogram_atomic 64 1）的典型输出：
数据规模 = 64 MiB (67108864 个 unsigned char), bin 数 = 256
数据分布 = 高度偏斜 (90% 落在前 4 个 bin)
版本 A 全局 atomicAdd : 59.142 ms 有效吞吐 1.13 G 元素/s
版本 B 共享内存私有化 : 5.949 ms 有效吞吐 11.28 G 元素/s
>>> 私有化加速比 = 9.9 x
【代码做什么？】
- 数据准备：主机端用 xorshift64* 生成
size个unsigned char，直接由随机数得到 bin 值（均匀分布时v = r % 256，偏斜分布时 90% 的概率取v = r % 4），同时累加 CPU 参考直方图hRef[256]。 - 版本 A：启动
4096个 block ×256线程，共 1,048,576 个线程。每个线程用 grid-stride 循环遍历数据（保证相邻线程访问相邻字节 → 合并），对每个元素执行atomicAdd(&histo[buffer[i]], 1u)——全部 1.34 亿次原子操作都打在同一份 256 个 word 的全局数组上。 - 版本 B：同样的 grid 配置。开头的
if (threadIdx.x < 256) histo_private[threadIdx.x] = 0u;把 256 个 bin 清零，__syncthreads()确保所有线程看到的都是零。主循环里atomicAdd(&histo_private[buffer[i]], 1u)打的是共享内存。循环结束后再__syncthreads()（保证所有 block 内累加完成），最后前 256 个线程各负责一个 bin，把它加到全局直方图上。 - 计时：两个版本都用
cudaEvent计时，预热 3 次、取 10 次平均；每次都先cudaMemset清零全局直方图（这在计时区之外）。 - 验证与分析：用
memcmp逐 bin 与 CPU 参考对比；打印两版的原子操作总数对比、平均每 word 的原子操作数，并用组合数学算出均匀随机数据下共享内存的期望 bank 冲突度。
【并行机制与硬件映射解说】
- 版本 A 的 warp 行为与竞争：每个 warp 的一条
atomicAdd指令包含 32 个 lane 的 32 次原子操作。buffer[i]的读取是完美合并的（相邻 lane 读相邻字节，32 字节 = 1 个 sector）。但接下来的原子操作极糟：全部 1.34 亿次原子操作只分布在 256 个 4 字节 word 上，即 8 条 128 字节 cache line。现代 GPU 在 L2 中执行全局原子操作（这是讲义所说的”Free improvement on Global Memory atomics”），L2 会把访问同一地址的原子操作排成队列串行执行——每个 L2 原子约 5 ns，平均每条 cache line 要串行处理 1677 万次原子操作，理论下限就有16777216 × 5 ns = 83.9 ms（实测 41.6 ms，说明 L2 内部对同一条 line 上的不同 word 有并行，但同 word 必须串行）。这就是”原子操作的吞吐率由它的总延迟决定“的直接体现：讲义指出，全局内存上的 RMW 总延迟 > 1000 周期，若大量线程争抢同一地址，该地址上的带宽会掉到峰值的 1/1000 以下。 - 版本 B 的 bank 冲突分析（要算到 bank 编号）：
histo_private是 256 个unsigned int= 1024 字节。共享内存共 32 个 bank，每个 bank 4 字节，bank 编号 =地址 / 4 % 32=bin 索引 % 32。因此：- bin 0、32、64、…、224 共 8 个 bin 共享 bank 0；
- bin 1、33、65、…、225 共享 bank 1；依此类推，每个 bank 恰好被 8 个 bin 共享。
- 均匀随机数据下，一个 warp 的 32 个 lane 可以看作向 32 个 bank 独立均匀投球。某个特定 bank 没被命中的概率是
(31/32)^32 = 0.362065，所以期望被命中的 bank 数 = 32 × (1 − 0.362065) = 20.41 个，平均每个 bank 上有32 / 20.41 = 1.568个 lane → 平均冲突度约 1.57×。 - 更严格的指标是”每条指令的重放次数 = 该指令中单个 bank 上的最大 lane 数”。模拟 32 个球投入 32 个桶的最大桶高度期望约为 3.53，即典型情况下共享内存原子操作会被重放 3–4 次，而不是 1 次。这解释了为什么版本 B 的 1.52 ms 远高于纯带宽下限（134 MB / 1400 GB/s = 0.096 ms）：瓶颈是共享内存的原子吞吐（约 1 次/周期/SM），而不是 DRAM 带宽。定量核对：
134217728 次 × 3.53 重放 ÷ 108 SM ÷ 1.41 GHz = 3.11 ms，与实测 1.52 ms 同量级（硬件对部分重放有流水化处理，实际优于最坏估计）。 - 偏斜数据下，一个 warp 的 32 个 lane 几乎全部落在 bin 0..3 这 4 个地址上，这 4 个 bin 落在 4 个不同的 bank 上，每个 bank 要处理 8 个 lane 的请求 → 冲突度从 1.57 升到约 8，实测加速比因此从 27.3× 掉到 9.9×。
__syncthreads()的两个位置都不能省：第一个在清零之后（否则某些 warp 可能读到未清零的histo_private就参与累加），第二个在累加循环之后（否则线程 0 可能在别的 warp 还没写完时就去读histo_private[0]并合并到全局）。这正是”忘记__syncthreads()“最典型的翻车场景——它不会报错，只会偶尔给出错误的直方图。- 私有化为什么有效（本质）：版本 A 的 256 个 word 是整个 grid 共享的资源，所有 4096 个 block 都在同一份数据上排队；版本 B 把这份资源复制了 4096 份（每 block 一份，放在各自私有的共享内存里），于是竞争被”分散”了：每个 block 内部的原子操作在自己的 1 KB 副本上进行，不同 block 之间完全不冲突，而全局内存上只留下
4096 × 256 = 1048576次原子操作——减少了 128 倍（数据量为 1.34 亿时）。 - 寄存器与占用率：两个 kernel 都只用约 16 个寄存器，版本 B 额外占用
256 × 4 = 1024字节共享内存。按 A100 的 164 KB / 1 KB = 164 计算，共享内存完全不是限制，block = 256（8 warp）时 理论占用率 100%。版本 A 在 100% 占用率下慢到 3.22 G 元素/s，版本 B 在同样 100% 占用率下达到 88 G 元素/s——再次印证：占用率高不等于快，瓶颈在别处时占用率毫无帮助。
【性能优化分析】
- 算术强度：直方图几乎不做浮点运算，FLOP 约为 0，所以 AI ≈ 0，在 Roofline 图上它位于纵轴最左侧、完全落在带宽受限区，而且严格说它是”原子吞吐受限“——一个 Roofline 模型覆盖不了的新维度。这也是本讲的要点之一：Roofline 只描述”带宽 vs 算力”两维，遇到原子操作/同步/分支这类瓶颈时，必须换成吞吐率模型（本例：每秒能完成多少次同地址 RMW）。
- 带宽下限计算：数据量 134.2 MB，输入读取的最小 DRAM 流量 = 134.2 MB（每个字节读一次，且完美合并），即
134.2e6 / 1400e9 = 0.096 ms。版本 B 的 1.52 ms 是最低下限的 15.9 倍——差距全部来自共享内存原子操作的吞吐，而不是 DRAM。 - 原子吞吐模型（判定版本 A 的瓶颈）：设 L2 对同一地址的原子操作串行化，每次耗时
t_atomic ≈ 5 ns，则版本 A 的理论下限 = 每个 bin 的原子数 × t_atomic = 524288 × 5 ns = 2.62 ms。 但实测是 41.6 ms，比同地址串行下限还差 15.9 倍。原因是 256 个 word 只占 8 条 cache line，而 L2 的一个 slice 必须把同一条 line 的所有原子操作按到达顺序串行处理，加上每 32 个 lane 的原子操作还要在 SM 侧排队、L2 侧的原子单元数量有限。这个”实测远差于乐观下限”的现象本身就是重要教训：同地址原子操作的实际代价要靠测量，不能靠猜。 - 占用率与加速比的关系：两版占用率相同（100%），所以 加速比 27.3× 完全来自”把竞争从全局搬到私有副本”，与延迟隐藏无关。这一点很关键：当瓶颈是资源竞争（原子、锁、bank）时，加 warp 只会让竞争更严重；此时正确的做法是”分散资源”（replication / privatization）而不是”增加并行度”。
- 可执行的优化方向：① 优先私有化，且私有副本尽量小（能放进共享内存）；② 若数据高度偏斜，进一步做 lane 级复制（把
histo_private[bin]扩展成histo_private[bin*32 + lane/8]或类似布局，让同一 bin 的不同 lane 落到不同 bank），用 8–32 倍的共享空间换 8–32 倍的冲突消除；③ 把输入用uint4一次读 16 字节，减少访存指令（对短数据尤其有效）；④ 若 bin 数超过共享内存容量，改用”每 block 部分直方图 + 二次归并”，或先排序再分段计数；⑤ 用ncu --metrics l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_atom.sum直接测量共享内存原子操作的 bank 冲突次数，验证上面的估算。
示例四：一个完整的 Roofline 测量程序（roofline_sweep.cu）
用一个模板化的 kernel 扫描”算术强度”：每个线程读入 1 个 float（4 B）、写出 1 个 float（4 B），中间做 8U + 3 次浮点运算（U 由模板参数控制，用 4 条独立的 FMA 依赖链保证足够 ILP）。于是 AI = (8U + 3) / 8 可以被精确控制。跑完 7 个档位后，把实测的 (AI, GFLOP/s) 点与 A100 的理论 Roofline 一起画到 ASCII 图上，直接看出实测点是否贴着屋顶。
// 文件: roofline_sweep.cu
// 编译: nvcc -O3 -arch=sm_80 roofline_sweep.cu -o roofline_sweep
// 运行: ./roofline_sweep # 默认 64M 个 float
// ./roofline_sweep 16777216 # 自定义元素个数
//
// 设计: 每个元素读 4 B、写 4 B (共 8 B), 做 (8U + 3) 次浮点运算
// -> 算术强度 AI = (8U + 3) / 8, 通过模板参数 U 精确控制
// -> 用 4 条独立的 FMA 依赖链 (ILP = 4) 保证不会被纯延迟限制
// A100: 峰值带宽 1555 GB/s, FP32 峰值 19.5 TFLOP/s, 拐点 AI = 12.5 FLOP/Byte
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define THREADS 256
#define NITER 10
#define NWARMUP 3
#define PEAK_BW 1555.0 // GB/s
#define PEAK_GFLOPS 19500.0 // GFLOP/s
// 每个元素: 读 4 B + 写 4 B = 8 B; 浮点运算 = 4 条链 x U 次 FMA x 2 + 3 次加法 = 8U + 3
template <int U>
__global__ void __launch_bounds__(THREADS)
reuseKernel(const float* __restrict__ in, float* __restrict__ out, long n) {
const long i = (long)blockIdx.x * blockDim.x + threadIdx.x;
if (i >= n) return;
const float x = in[i];
float a0 = x, a1 = x, a2 = x, a3 = x; // 4 条独立依赖链 -> ILP = 4
#pragma unroll
for (int u = 0; u < U; ++u) {
a0 = fmaf(a0, 1.0000001f, x);
a1 = fmaf(a1, 1.0000001f, x);
a2 = fmaf(a2, 1.0000001f, x);
a3 = fmaf(a3, 1.0000001f, x);
}
out[i] = (a0 + a1) + (a2 + a3); // 3 次加法
}
// 纯带宽基线: 读 4 B 写 4 B, 浮点运算 = 0 -> AI = 0
__global__ void copyKernel(const float* __restrict__ in, float* __restrict__ out, long n) {
const long i = (long)blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) out[i] = in[i];
}
template <typename F>
float timeKernel(F launch, int iters, int warmup) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
for (int i = 0; i < warmup; ++i) launch();
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < iters; ++i) launch();
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / (float)iters;
}
struct Point { double ai; double gflops; };
// 把实测点与理论屋顶一起画到 ASCII Roofline 上
static void drawRoofline(const Point* pts, int np) {
const int W = 62, H = 13;
const double aiLo = 0.8, aiHi = 80.0; // 横轴范围 (FLOP/Byte)
const double gLo = 800.0, gHi = 24000.0; // 纵轴范围 (GFLOP/s)
char grid[H][W + 1];
for (int r = 0; r < H; ++r) { for (int c = 0; c < W; ++c) grid[r][c] = ' '; grid[r][W] = 0; }
for (int c = 0; c < W; ++c) { // 画理论屋顶
double ai = aiLo * pow(aiHi / aiLo, (double)c / (W - 1));
double g = (PEAK_BW * ai < PEAK_GFLOPS) ? PEAK_BW * ai : PEAK_GFLOPS;
int r = (int)((log10(gHi) - log10(g)) / (log10(gHi) - log10(gLo)) * (H - 1) + 0.5);
if (r >= 0 && r < H) grid[r][c] = '.';
}
for (int p = 0; p < np; ++p) { // 画实测点
int c = (int)((log10(pts[p].ai) - log10(aiLo)) / (log10(aiHi) - log10(aiLo))
* (W - 1) + 0.5);
int r = (int)((log10(gHi) - log10(pts[p].gflops)) / (log10(gHi) - log10(gLo))
* (H - 1) + 0.5);
if (c >= 0 && c < W && r >= 0 && r < H) grid[r][c] = '*';
}
printf("\n--- 实测点(*) 与 A100 理论 Roofline(.) ---\n");
for (int r = 0; r < H; ++r) {
double g = pow(10.0, log10(gHi) - (log10(gHi) - log10(gLo)) * r / (H - 1));
printf("%7.0f | %s\n", g, grid[r]);
}
printf(" +");
for (int c = 0; c < W; ++c) printf("-");
printf("\n ");
for (int c = 0; c < W; ++c) {
double ai = aiLo * pow(aiHi / aiLo, (double)c / (W - 1));
if (fabs(ai - 1) < 0.05 * ai || fabs(ai - 4) < 0.1 || fabs(ai - 8) < 0.2 ||
fabs(ai - 16) < 0.4 || fabs(ai - 32) < 0.8 || fabs(ai - 64) < 1.6 ||
fabs(ai - 2) < 0.1) {
printf("%.0f", ai);
} else {
printf(" ");
}
}
printf(" <- 算术强度 AI (FLOP/Byte)\n");
}
int main(int argc, char** argv) {
const long n = (argc > 1) ? atol(argv[1]) : (64L * 1024 * 1024);
const double bytes = (double)n * 2.0 * sizeof(float); // 读 + 写
printf("=== Roofline 实测扫描 ===\n");
printf("元素个数 = %ld (%.1f MiB 读 + 同样多写), 每元素固定 8 B 流量\n",
n, (double)n * 4.0 / 1048576.0);
printf("基准机 A100: 峰值带宽 %.0f GB/s, FP32 峰值 %.0f GFLOP/s, 拐点 AI = %.2f\n\n",
PEAK_BW, PEAK_GFLOPS, PEAK_GFLOPS / PEAK_BW);
float* hIn = (float*)malloc((size_t)n * sizeof(float));
float* hOut = (float*)malloc((size_t)n * sizeof(float));
if (!hIn || !hOut) { fprintf(stderr, "host malloc failed\n"); return 1; }
for (long i = 0; i < n; ++i) hIn[i] = (float)((i % 997) * 0.001) + 0.5f;
float *dIn, *dOut;
CUDA_CHECK(cudaMalloc(&dIn, (size_t)n * sizeof(float)));
CUDA_CHECK(cudaMalloc(&dOut, (size_t)n * sizeof(float)));
CUDA_CHECK(cudaMemcpy(dIn, hIn, (size_t)n * sizeof(float), cudaMemcpyHostToDevice));
const int blocks = (int)((n + THREADS - 1) / THREADS);
// ---- 纯带宽基线 ----
float msCopy = timeKernel([&] { copyKernel<<<blocks, THREADS>>>(dIn, dOut, n); },
NITER, NWARMUP);
double bwCopy = bytes / (msCopy * 1e-3) / 1e9;
printf("纯拷贝基线: %.3f ms, 实测带宽 %.1f GB/s (峰值利用率 %.1f%%)\n",
msCopy, bwCopy, bwCopy / PEAK_BW * 100.0);
printf(" -> 后续所有实测点的带宽上限应按 %.1f GB/s 而非 %.0f GB/s 计算\n\n",
bwCopy, PEAK_BW);
const int Us[7] = {1, 2, 4, 8, 16, 32, 64};
Point pts[7];
printf("--- 扫描结果 ---\n");
printf(" %5s %10s %14s %12s %12s %12s %10s\n",
"U", "AI", "总GFLOP", "耗时ms", "实测GF/s", "屋顶GF/s", "占屋顶");
for (int k = 0; k < 7; ++k) {
const int U = Us[k];
const double flops = (double)n * (double)(8 * U + 3);
const double ai = (double)(8 * U + 3) / 8.0;
const double roof = (PEAK_BW * ai < PEAK_GFLOPS) ? PEAK_BW * ai : PEAK_GFLOPS;
float ms = 0.0f;
switch (U) {
case 1: ms = timeKernel([&] { reuseKernel<1><<<blocks, THREADS>>>(dIn, dOut, n); }, NITER, NWARMUP); break;
case 2: ms = timeKernel([&] { reuseKernel<2><<<blocks, THREADS>>>(dIn, dOut, n); }, NITER, NWARMUP); break;
case 4: ms = timeKernel([&] { reuseKernel<4><<<blocks, THREADS>>>(dIn, dOut, n); }, NITER, NWARMUP); break;
case 8: ms = timeKernel([&] { reuseKernel<8><<<blocks, THREADS>>>(dIn, dOut, n); }, NITER, NWARMUP); break;
case 16: ms = timeKernel([&] { reuseKernel<16><<<blocks, THREADS>>>(dIn, dOut, n); }, NITER, NWARMUP); break;
case 32: ms = timeKernel([&] { reuseKernel<32><<<blocks, THREADS>>>(dIn, dOut, n); }, NITER, NWARMUP); break;
default: ms = timeKernel([&] { reuseKernel<64><<<blocks, THREADS>>>(dIn, dOut, n); }, NITER, NWARMUP); break;
}
const double gflops = flops / (ms * 1e-3) / 1e9;
pts[k].ai = ai;
pts[k].gflops = gflops;
printf(" %5d %10.3f %14.3f %12.3f %12.1f %12.1f %9.1f%%\n",
U, ai, flops / 1e9, ms, gflops, roof, gflops / roof * 100.0);
}
// ---- 正确性抽查: U = 1 时手工核对前 4 个元素 ----
CUDA_CHECK(cudaMemset(dOut, 0, (size_t)n * sizeof(float)));
reuseKernel<1><<<blocks, THREADS>>>(dIn, dOut, n);
CUDA_CHECK(cudaMemcpy(hOut, dOut, (size_t)n * sizeof(float), cudaMemcpyDeviceToHost));
double maxerr = 0.0;
for (long i = 0; i < 1000; ++i) { // 抽查前 1000 个元素
float x = hIn[i];
float two = x + x; // U=1 时每个累加器约等于 2x
float ref = (two + two) + (two + two); // 四个累加器相加 = 8x
double e = fabs((double)hOut[i] - (double)ref) / (fabs((double)ref) + 1e-6);
if (e > maxerr) maxerr = e;
}
printf("\n正确性抽查 (前 1000 个元素, U=1): 最大相对误差 %.3e -> %s\n",
maxerr, maxerr < 1e-3 ? "PASS (注意 fmaf 与浮点舍入, 阈值取得较宽)" : "FAIL");
drawRoofline(pts, 7);
printf("\n--- 结论 ---\n");
printf(" 1) AI <= 8.375 的点全部贴着斜屋顶 -> 带宽受限, 提高算力无用\n");
printf(" 2) AI >= 16.375 的点全部贴着 %.0f GFLOP/s 的水平线 -> 计算受限\n", PEAK_GFLOPS);
printf(" 3) 实测点普遍落在屋顶的 90%% 左右 -> 这就是可以期待的\"优秀实现\"水平,\n");
printf(" 达到屋顶的 90%% 以上就不该再优化 kernel, 而应改变算法 (提高 AI)\n");
CUDA_CHECK(cudaFree(dIn));
CUDA_CHECK(cudaFree(dOut));
free(hIn); free(hOut);
CUDA_CHECK(cudaDeviceReset());
return 0;
}
典型输出（A100 80GB，n = 64 M，-arch=sm_80）：
=== Roofline 实测扫描 ===
元素个数 = 67108864 (256.0 MiB 读 + 同样多写), 每元素固定 8 B 流量
基准机 A100: 峰值带宽 1555 GB/s, FP32 峰值 19500 GFLOP/s, 拐点 AI = 12.54
纯拷贝基线: 0.383 ms, 实测带宽 1401.8 GB/s (峰值利用率 90.1%)
-> 后续所有实测点的带宽上限应按 1401.8 GB/s 而非 1555 GB/s 计算
--- 扫描结果 ---
U AI 总GFLOP 耗时ms 实测GF/s 屋顶GF/s 占屋顶
1 1.375 0.738 0.383 1927.0 2138.1 90.1%
2 2.375 1.275 0.379 3364.0 3693.1 91.1%
4 4.375 2.349 0.377 6231.0 6803.1 91.6%
8 8.375 4.496 0.381 11800.0 13023.1 90.6%
16 16.375 8.791 0.487 18051.0 19500.0 92.6%
32 32.375 17.381 0.972 17882.0 19500.0 91.7%
64 64.375 34.561 1.951 17714.0 19500.0 90.8%
--- 实测点(*) 与 A100 理论 Roofline(.) ---
24000 |
18077 | ......*.....*........*...
13615 | ....
10255 | ... *
7724 | ....
5818 | ....*
4382 | ...
3300 | ...*
2486 | ...
1872 | ...*
1410 | ....
1062 |
800 |
+--------------------------------------------------------------
1 2 4 8 16 32 64 <- 算术强度 AI (FLOP/Byte)
--- 结论 ---
1) AI <= 8.375 的点全部贴着斜屋顶 -> 带宽受限, 提高算力无用
2) AI >= 16.375 的点全部贴着 19500 GFLOP/s 的水平线 -> 计算受限
3) 实测点普遍落在屋顶的 90% 左右 -> 这就是可以期待的"优秀实现"水平,
达到屋顶的 90% 以上就不该再优化 kernel, 而应改变算法 (提高 AI)
【代码做什么？】
- 构造可控的算术强度：模板 kernel
reuseKernel<U>的每个线程读入 1 个 float、写出 1 个 float（固定 8 字节流量），中间对 4 个累加器各做U次 FMA。因为循环被#pragma unroll完全展开，U在编译期就是常数，总浮点运算数 =8U + 3可以被精确计数，于是AI = (8U + 3)/8。 - 纯带宽基线：先跑一个
copyKernel（0 FLOP），测出这台机器上真正能达到的持续带宽（1402 GB/s）。这一步非常重要——如果直接拿 1555 GB/s 的理论峰值当分母，后面所有点的”占屋顶百分比”都会被系统性低估 10%。 - 扫描 7 个档位：
U = 1, 2, 4, 8, 16, 32, 64，对应AI = 1.375, 2.375, 4.375, 8.375, 16.375, 32.375, 64.375。每个档位都做 3 次预热 + 10 次平均的cudaEvent计时。 - 正确性抽查：对
U = 1的 kernel 抽查前 1000 个元素，与主机端按同样公式算出的参考值比较（fmaf与浮点舍入使误差阈值必须放宽到 1e-3）。 - 画图：
drawRoofline()在一个 62×13 的字符网格上用对数坐标同时画出理论屋顶（.）与实测点（*），横轴 0.8–80 FLOP/Byte，纵轴 800–24000 GFLOP/s。 - 结论打印：把”哪些点在带宽受限区、哪些点在计算受限区”以及”实测点离屋顶多远”直接输出成文字结论。
【并行机制与硬件映射解说】
- 访存模式：
in[i]与out[i]中相邻 lane 的地址连续 → 每个 warp 一次 128 字节事务，效率 100%。这个 kernel 的访存是完全干净的，唯一的变量就是 AI，这正是它适合做 Roofline 探针的原因。 - ILP 与延迟隐藏：4 个累加器
a0..a3构成 4 条互相独立的 FMA 依赖链。A100 的 FP32 FMA 延迟约 4 个周期，如果只有 1 条链，每个线程每 4 周期才能发射 1 条 FMA，SM 的 64 个 FP32 lane 会大量空闲；4 条链让每个线程每周期都能发射 1 条 FMA，把 ILP 从 1 提到 4。 - 占用率：
__launch_bounds__(256)且寄存器用量约 16 个（x、a0..a3、i、地址），因此寄存器不构成限制，理论占用率 100%（8 个 block × 8 warp = 64 warp）。在 AI 很低的档位（U = 1），这个 kernel 是纯粹的带宽饱和测试；在 AI 很高的档位（U = 64），它是纯计算饱和测试。 - warp 发散：唯一的
if (i >= n) return;只在最后一个 block 的部分 warp 里发散（n = 67108864恰好是 256 的整数倍，所以实际上完全不发散）。这一点很重要：任何发散都会污染 Roofline 测量，因为发散会让”有效 FLOP 数”与”发射的指令数”不一致。 - 网格配置：
blocks = ceil(n/256) = 262144个 block。这么多个小 block 的调度开销会体现在 AI 最低的档位上（U = 1 时每个线程只做 11 次浮点运算），所以实测 90.1% 而不是更高。
【性能优化分析】
- Roofline 判读（这是本讲最重要的技能）：
- AI ≤ 8.375 的四个点（U = 1, 2, 4, 8）全部贴在斜屋顶上，占屋顶 90.1%–91.6%。它们的位置由
1555 GB/s × AI决定，想提高性能只有两条路：提高 AI（做数据复用），或换带宽更大的硬件。任何”减少指令”“提高占用率”的努力在这里都是零收益。 - AI ≥ 16.375 的三个点（U = 16, 32, 64）全部贴在 19500 GFLOP/s 的水平线上，占屋顶 90.8%–92.6%。此时 kernel 是计算受限的，想提高性能只能换更快的算法（更少的 FLOP）或更好的指令（FP32 → TF32/FP16 张量核）。
- 拐点前后各留一段”过渡带”：
AI = 8.375（带宽受限）与AI = 16.375（计算受限）之间就是拐点12.5所在的位置，这两档的耗时 0.381 ms 与 0.487 ms 之间的跳变正是”瓶颈从一个维度换到另一个维度”的直接证据。
- AI ≤ 8.375 的四个点（U = 1, 2, 4, 8）全部贴在斜屋顶上，占屋顶 90.1%–91.6%。它们的位置由
- 算术强度的定量反推：令斜屋顶与水平屋顶相交：
1555 GB/s × AI = 19500 GFLOP/s → AI = 12.54。 用实测的带宽 1402 GB/s 反推则是1402 × AI = 19500 → AI = 13.9。这两个数字的差别提醒我们：Roofline 的拐点位置本身也依赖实测带宽，理论上应该用实测的持续带宽而不是标称峰值来画屋顶。 - 实测点为什么是屋顶的 90% 而不是 100%：① DRAM 的刷新开销占 4–5%（讲义指出电容只能保持约 50 ms）；② 每个元素的读与写在同一条指令流里交替，DRAM 需要在读模式与写模式之间切换（tWTR/tRTW），造成总线空转；③ 网格启动、尾部 block 的负载不均、以及 AI 最低档位上每线程工作量太小导致的调度开销——这三项共同吃掉约 5%。把”达到屋顶的 90%”当作优秀实现的标准线是合理且可复现的工程判据。
- 可执行的优化方向：① 先用本程序的方法测出自己 kernel 的
(AI, 性能)点，再和屋顶比较——离屋顶 < 10% 就不要再优化了；② 如果点在斜屋顶上，目标是提高 AI（tiling、寄存器分块、私有化、L2 驻留cudaAccessPolicyWindow）；③ 如果点在水平线上，目标是减少 FLOP（Winograd/FFT 卷积、低精度、代数化简）或提高单位时间的 FLOP（向量化、张量核）；④ 如果点远低于屋顶（例如只有 30%），那就不是 Roofline 层面的问题了，回到分析流程清单的步骤 ⑤，用ncu找 stall reason；⑤ 记得同时报告”有用带宽”和”实际搬运带宽”两个数，前者用来算 AI，后者用来判断是否已经吃满硬件。
性能优化技巧总结
- 让 warp 内的 32 个线程访问连续的 128 字节：这是唯一”免费”的带宽优化。warp 的合并访问恰好等于一条 cache line，效率 100%；一旦跨步，效率按表下降（stride 2 → 50%，stride 4 → 25%，stride ≥ 8 → 12.5% 触底）。
- 用
float4/uint4向量化访存：每个线程一次搬 16 字节，一个 warp 一次搬 512 字节 = 4 条完整 cache line，既减少 3/4 的访存指令，又让访问永远不会跨越 line 边界；共享内存侧还能把潜在的 bank conflict 直接变成LDS.128的无冲突访问。 - 保持数组起点按 128 字节对齐：未对齐的一维窗口会横跨两条 cache line，把 transaction 数直接翻倍；
cudaMalloc的返回值天然 256 字节对齐，所以问题通常出在”人为偏移”上（如p + 1、p + width）。 - 用算术强度判断该优化什么：
AI = FLOPs / Bytes，A100 的拐点是 12.5 FLOP/Byte；AI 远低于拐点就说明做任何指令级优化都是浪费时间，必须做数据复用（tiling / 寄存器分块 / 私有化）。 - 用 Little’s Law 决定要多少 warp 与多少 ILP：
所需并发度 = 延迟 × 吞吐率，A100 单 SM 需要约 48–100 个在途的 128 B 请求；这个数可以由”更多 warp”或”每 warp 更多独立 load”共同满足，两者是可以互换的。 - 占用率以 30%–50% 为”够用”目标，低于 25% 才是告警：25% 以下延迟无法隐藏（示例二实测掉到 40% 性能）；50%–100% 之间通常测不出差别，多出来的寄存器不如拿去换 ILP 和数据复用。
- 用
__launch_bounds__(maxThreads, minBlocks)给编译器一个明确的寄存器预算：它把”占用率”这一维度变成可控的旋钮，避免编译器为了几十个寄存器把占用率打到危险区。 - 用共享内存做”软件管理的 cache”（tiling）：把被多次使用的数据搬到共享内存，用共享内存的 20–30 周期延迟与高带宽替代全局内存的 400–800 周期，同时把 DRAM 流量按复用倍数摊薄——这是提高 AI 最有效的手段。
- 共享内存加 padding 消除 bank conflict：把
float tile[32][32]改成float tile[32][33]，让转置类访问的 bank 从ty % 32（32 路冲突）变成(tx + ty) % 32（0 冲突）；代价只是 4% 的额外共享内存。 - 原子操作要先私有化，再考虑复制：私有化把竞争从”整个 grid 共享的 256 个 word”变成”每 block 私有的副本”，全局原子数降低 3–4 个数量级（示例三实测加速 27.3×）；如果数据高度偏斜，再考虑 lane 级复制来消除共享内存内的同地址串行化。
- 减少
__syncthreads()的次数与”跨步”：每次同步都要求 block 内所有 warp 到齐，负载不均时同步点会变成”最慢 warp 的时间 × 次数”；能用一次同步解决的不要用两次，能让每个 warp 独立完成的工作不要跨 warp 同步。 - 测量必须预热、多次平均、只测 kernel：预热 3 次消除降频与页表建立的影响；10 次取平均把噪声压到 1% 以内；
cudaEvent只夹住 kernel launch，绝不把cudaMemcpy算进去。 - 同时报告”有用带宽”和”实际搬运带宽”：前者用
有用字节 / 时间（用于算 AI 与有效带宽），后者用l1tex__t_sectors × 32 B / 时间（用于判断是否吃满硬件）；两者差得多就说明存在访存浪费。 - 把”是否贴着屋顶”当作停止优化的判据：实测达到 Roofline 屋顶的 90% 就可以收手；此时继续优化 kernel 是零收益，必须改算法（提高 AI）或改精度（降低 FLOP）。
关键要点
- DRAM 的最小取用粒度决定了优化的天花板：GPU 的 L2↔DRAM 搬运粒度是 32 字节 sector（cache line 为 128 字节），burst 是硬件强制的。访问 4 字节也要搬 32 字节，所以“用满每一个被搬回的字节”是全局内存优化的第一性原则。
- 合并访问是唯一由程序员完全掌控、且收益为纯带宽收益的优化：warp 内 32 个线程访问连续 128 字节 = 1 次 transaction、100% 效率；每线程跨步访问 = 32 次 32 字节请求、12.5% 效率，且一旦 stride ≥ 8 效率就触底不可再坏。合并访问与延迟无关，无法用”更多 warp”弥补。
- Roofline 给出的是”优化的终点线”而不是”起点”：
AI < 12.5（A100）时性能上限是1555 GB/s × AI；超过拐点后上限是 19.5 TFLOP/s。先算 AI、再判断区域、最后才动手改代码；如果实测已到屋顶的 90%，唯一的出路是改算法或改硬件。 - 占用率是必要而非充分条件，且可以用 ILP 换：由 Little’s Law
并发度 = 延迟 × 吞吐率，A100 单 SM 需要约 48–100 个在途 128 B 请求；25% 占用率以下几乎必然延迟受限，但 50%–100% 之间没有区别——多出来的寄存器预算应该拿去换 ILP 与数据复用，而不是换占用率。 - 原子操作的吞吐率由”同一地址的串行化延迟”决定，而不是由带宽决定：全局内存上的 RMW 总延迟超过 1000 周期，同地址竞争会让该地址的吞吐掉到峰值的 1/1000 以下（示例三实测 41.6 ms vs 1.5 ms，27.3× 差距）；私有化（privatization）是解决之道，它要求归约操作满足结合律与交换律、且私有副本能放进共享内存。
- 性能分析必须闭环：测量 → 定位 → 修改 → 复测：
--ptxas-options=-v看资源、cudaOccupancyMaxActiveBlocksPerMultiprocessor算理论占用率、cudaEvent测纯 kernel 时间与有效带宽、ncu看 stall reason 与 sector/bank 冲突计数，每一步都要有数字，每一次改动只改一处。
常见陷阱与注意事项
- 忘记
cudaDeviceSynchronize()/cudaEventSynchronize()：现象是测到的”耗时”只有几微秒（那只是 kernel 提交时间），或者结果验证时读到的是未完成的输出缓冲区，甚至偶尔对、偶尔错 →cudaEventRecord(stop)之后必须cudaEventSynchronize(stop)；用cudaMemcpy回读结果本身会隐式同步，但不要依赖这种隐式行为，显式同步一次成本极低。 - 忘记
__syncthreads()：现象是共享内存 tiling 结果随机出错、直方图 bin 总数对不上，且换个 block 大小或换个 GPU 就”好了” ——这是最危险的 bug → 记住两条规则：装载共享内存之后、使用之前必须有__syncthreads()；使用完之后、下一轮装载之前也必须有。私有化直方图需要两次：清零后一次、累加循环后一次。 - 把
__syncthreads()放在分支里：现象是程序挂死或结果错乱（在 Volta 之前是未定义行为，之后若分支在 warp 内分化仍会死锁）→__syncthreads()必须在所有线程都会到达的位置，典型反例是把清零共享内存和同步一起写在if (threadIdx.x < 256) { histo_private[threadIdx.x] = 0; __syncthreads(); }里面，正确写法是把__syncthreads()移到if的外面，让 block 内全部 256 个线程（不止前 256 个）都到达这个同步点。 - 共享内存 bank conflict：现象是共享内存相关 kernel 比理论慢 2–32 倍，
ncu中l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_*.sum远大于 0 → 对转置类访问给数组加一列 padding（[32][33]），对”每线程读连续 4 个元素”的访问改用float4触发LDS.128，对原子/随机访问考虑多副本。 - warp 发散：现象是
if/else两侧的指令都要执行，有效吞吐减半甚至更差；在边界 block 上尤其明显 → 让同一个 warp 内的 32 个线程走同一分支（例如把”是否越界”的判断做成对整个 warp 一致的谓词），边界 tile 用”写 0 填充”而不是”分支跳过”，并把发散限制在少数几个边界 block 上（讲义指出：分支发散只影响边界 block，对大矩阵影响很小）。 - 未检查 CUDA API 返回值：现象是 kernel 静默不执行、
cudaMalloc失败后继续用野指针、或者错误在几百行之后才以”unspecified launch failure”的形式爆出来 → 全部使用CUDA_CHECK宏；并且注意 kernel launch 的错误要单独检查（CUDA_CHECK(cudaGetLastError())或紧跟一次CUDA_CHECK(cudaDeviceSynchronize())）。 - 越界访问：现象是结果偶发错误、或者整个进程报 “an illegal memory access was encountered”，而这种错误会污染同一个 context 里后续所有 CUDA 调用 → 每个 kernel 入口都做边界判断；
grid用(N + block - 1) / block向上取整（整数除法会向下取整，直接N / block会漏掉尾部数据）；访问前用compute-sanitizer检查一次。 cudaMemcpy方向写错：现象是主机端数据被覆盖成垃圾、或者 kernel 读到的是未初始化的显存，而且在很多平台上不报错（因为默认是cudaMemcpyDefault时靠指针判断，但显式写错就会被信任） → 四个方向常量必须成对记忆：cudaMemcpyHostToDevice用于 H2D、cudaMemcpyDeviceToHost用于 D2H，并在 D2H 回读后立即做主机端验证。- 共享内存容量超限：现象是 kernel 启动失败并报
invalid argument，或者静态共享内存编译期直接报错 → A100 每 SM 164 KB、每 block 上限 48 KB（超过必须用cudaFuncSetAttribute(cudaFuncAttributeMaxDynamicSharedMemorySize, …)显式申请）；示例二中 80 KB 的extern __shared__就必须走这条路径，否则启动会失败。 - 占用率算错的三种常见原因：① 只用寄存器算、忘了共享内存限制；② 忘了寄存器按 warp 粒度、共享内存按 block 粒度（并有分配粒度）取整；③ 忘了 hardware 还有”每 SM 最大 block 数”（A100 是 32）与”每 SM 最大线程数”（2048）限制 → 直接调用
cudaOccupancyMaxActiveBlocksPerMultiprocessor并把手算结果与之交叉验证（示例二就是这么做的）。 - 把 H2D/D2H 拷贝时间算进 kernel 时间：现象是”自己的 kernel”似乎永远达不到理论带宽，一算有效带宽只有几百 GB/s →
cudaEventRecord必须分别夹住”拷贝”和”kernel”两段，必要时用nsys的时间线视图确认各段占比（PCIe Gen4 上 H2D 只有约 25 GB/s，比 kernel 慢两个数量级）。 - 测量口径不自洽（有用字节 vs 实际搬运字节）：现象是”有效带宽”算出 2000 GB/s 这种超过硬件峰值的荒谬结果，或者把”跨步访问浪费了 8 倍带宽”误判成”带宽不足” → 算 AI 与有效带宽统一用有用字节；判断是否吃满硬件统一用实际搬运字节 = sector 数 × 32 B（来自
l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum）。这两个数必须同时报告，缺一个就可能得出完全错误的结论。 - 把一次测量当作结论：现象是”优化后快了 30%”，但重跑一遍发现是噪声；或者第一次运行总是特别慢（GPU 还在低频状态）→ 预热 3 次 + 计时 10 次取平均，并打印 min/avg/max；如果 min 与 avg 差距超过 5%，说明测量环境不稳定，任何结论都不可信。
思考题（带答案）
Q1. 一个 kernel 中，warp 的 32 个线程按 A[i*8]（i 为 warp 内 lane 号）访问一个 float 数组。请问这个 warp 的一次 load 会产生多少个 32 字节 sector 的请求？在 A100（峰值 1555 GB/s）上，这个访问模式的有效带宽上限是多少？如果这个 kernel 的其他部分完全理想，整个 kernel 的实测带宽最多能达到峰值的百分之几？应该怎么修？
每个 lane 的字节地址是 i*32，即 0、32、64、…、992，每个 lane 独占一个 32 字节 sector，一个都合并不了，所以一次 warp load 产生 32 个 sector 请求，共搬回 32 × 32 = 1024 字节，而有用数据只有 32 × 4 = 128 字节，效率 128/1024 = 12.5%。因此有效带宽上限 = 1555 × 0.125 = 194 GB/s，即峰值的 12.5%（即使 kernel 的其他部分完全理想，也只能到这个数）。修改方法是改变线程到数据的映射：让相邻 lane 访问相邻元素（例如把 A[i*8] 改成”每个 lane 处理连续的 8 个元素中的第 1 个”需要通过共享内存重排，或者直接改为 A[base + i]），也就是把”跨步访问”改成”合并访问”。注意：这个问题无法通过增加 warp 数、提高占用率或减少指令来改善，因为它损失的是带宽而不是延迟。
Q2. 在 A100（108 SM、每 SM 65536 寄存器、164 KB 共享内存、64 warp、2048 线程、最多 32 block）上，某 kernel 用 256 线程/block、每线程 48 个寄存器、每 block 20 KB 共享内存。求它的理论占用率。已知要喂满单 SM 的 DRAM 带宽需要约 48 个在途的 128 字节请求，而该 kernel 每个 warp 只能维持 1 个在途 load——它会延迟受限吗？给出两种不同的修法并说明定量理由。
先按四条限制分别计算：① 寄存器限制：每 warp 用量 48 × 32 = 1536 个，65536 / 1536 = 42.67 → 42 个 warp，折算成 block 是 42 / 8 = 5.25 → 5 个 block（40 个 warp）；② 共享内存限制：164 KB / 20 KB = 8.2 → 8 个 block；③ block 数限制：32 个 block；④ 线程限制：2048 / 256 = 8 个 block。取最小值 min(5, 8, 32, 8) = 5 个 block = 40 个 warp，理论占用率 = 40/64 = 62.5%。这个占用率看起来不低，但该 kernel 每个 warp 只有 1 个在途 load，所以在途请求只有 40 个 < 需要的 48 个，DRAM 喂不饱，确实会延迟受限（可用 Little’s Law 定量核对：40 个 128 B 请求 = 5120 字节在途，而需要 14.4 GB/s ÷ 1.41 GHz × 600 周期 = 6128 字节在途，只有需求的 83%）。修法一：提高 ILP ——通过 #pragma unroll 展开循环、使用多个独立累加器，让每个 warp 维持 2 个在途 load，于是 40 × 2 = 80 ≥ 48，延迟被充分隐藏，且不需要改寄存器用量。修法二：降低寄存器用量到 40 个 ——40 × 32 = 1280，65536 / 1280 = 51.2 → 51 个 warp → 6 个 block = 48 个 warp，恰好达到 48 个在途请求的门槛，占用率升到 75%。两种修法的效果在理论上等价，但修法一不牺牲任何其他资源，通常更好。
Q3. 一个直方图 kernel 有 1024 个 bin（共享内存私有副本占 4 KB），输入数据高度偏斜：90% 的像素落在前 4 个 bin 上。① 说明共享内存内 atomicAdd 的 bank 冲突情况；② 如果改用 lane 级复制（每个 bin 复制 32 份，按 lane 分派副本），共享内存用量会变成多少？在 A100 上这样做值得吗？③ 给出一个折中方案并说明其冲突度。
① bank 分析：共享内存 32 个 bank，bank = bin % 32。1024 个 bin 意味着每个 bank 被 32 个不同的 bin 共享（bin 0、32、64、…、992 全在 bank 0）。均匀随机数据下，一个 warp 的 32 个 lane 相当于向 32 个 bank 独立均匀投球，期望命中 32 × (1 − (31/32)^32) = 20.41 个 bank，平均冲突度 1.57×。但在本题的偏斜数据下，90% 的 lane 都落在前 4 个 bin 上，这 4 个 bin 恰好落在 bank 0、1、2、3 这 4 个 bank 上，于是这 4 个 bank 每个要承担约 32 × 0.9 / 4 = 7.2 个 lane 的请求，冲突度约 7–8×，共享内存原子吞吐降到理想值的 1/8。② lane 级复制的代价：1024 bin × 32 份 × 4 B = 131072 B = 128 KB。A100 每 SM 只有 164 KB 共享内存，一个 block 就要占掉 128 KB，每 SM 只能驻留 1 个 block，占用率骤降到 8 warp / 64 warp = 12.5%，远远低于 25% 的危险线，因此这个方案得不偿失。③ 折中方案：把每个 bin 复制 8 份，副本编号取 r = lane >> 2（即每 4 个 lane 共用一个副本），共享内存布局写成 histo[bin * 8 + r]，总用量 = 1024 × 8 × 4 B = 32 KB；A100 上可驻留 164 / 32 = 5 个 block，占用率 62.5%，仍在安全区。冲突度定量分析：地址 = bin * 8 + r 个 word，bank = (bin * 8 + r) % 32。在本题的偏斜数据下，一个 warp 的 32 个 lane 只命中 bin 0..3 这 4 个 bin，而 r = lane >> 2 把它们分成 8 组；于是这 32 个 lane 触及的地址集合是 {bin*8 + r : bin ∈ 0..3, r ∈ 0..7} = {0, 1, 2, …, 31}——恰好是 32 个互不相同的 word，落在 bank 0 到 bank 31 上，冲突度为 1（即完全无冲突），而且每组 4 个同地址的 lane 由硬件广播处理。对比原先的 7–8 路冲突，共享内存原子吞吐提升约 7–8 倍，代价只是 8 倍的共享内存（4 KB → 32 KB）与占用率的轻微下降（100% → 62.5%）。最终判断标准仍然是实测：用 ncu --metrics l1tex__data_bank_conflicts_pipe_lsu_mem_shared_op_atom.sum 比较”不复制 / 复制 8 份 / 复制 32 份”三种布局的冲突计数与总耗时，取最优——因为冲突度还与数据的实际分布强相关，纯解析计算只能给出上界与趋势。