Lecture 3: 并行模式一 —— 向量加法：线程映射与内存访问 (对应 Lab 1: Vector Addition)
Lecture 3: 并行模式一 —— 向量加法：线程映射与内存访问 (对应 Lab 1: Vector Addition)
概述
本讲用一个”每元素只做一次加法”的向量加法 kernel（C[i] = A[i] + B[i]），把 CUDA 的三件事一次性讲透：线程如何映射到数据（i = blockIdx.x * blockDim.x + threadIdx.x）、数组元素如何落到内存事务上（合并访存 coalescing）、以及硬件如何用大量并发 warp 隐藏几百个时钟周期的全局内存延迟（占用率 occupancy 与延迟容忍 latency tolerance）。结论先行：向量加法在 A100 上只有 0.125 FLOP/Byte 的算术强度，远低于 12.54 FLOP/Byte 的机器平衡点，是一个纯粹的显存带宽受限（memory bandwidth bound）问题——你唯一能做的就是让每一个字节都被用满，然后让足够多的访存同时在路上。这也是 ECE408 第一个 Lab（Lab 1）的全部主题：正确的线程映射 + 正确的边界检查 + 正确的带宽测量。
核心概念与 GPU 架构图解
概念 1：数据并行与线程到数据的映射（Data Parallelism and Thread-to-Data Mapping）
定义与目的：数据并行（data parallelism）指的是”同一段程序作用在一大批互不相关的数据上”，每个数据元素的计算完全不依赖其它元素的计算结果。向量加法就是最纯粹的数据并行：
C[0]的值与C[1]的值毫无关系。CUDA 的做法是把这个”for 循环”拆开：让硬件去决定循环的哪一次迭代在哪个物理单元上跑，程序员只负责给出”元素编号 → 线程编号”的映射公式。这个映射公式就是i = blockIdx.x * blockDim.x + threadIdx.x。它解决的核心问题是：GPU 上同时有几十万个线程，必须有一种廉价（几条整数指令）、唯一、且与内存布局对齐的方式给每个线程分配一份不重复的工作。直观解释（”它是什么？”）：想象一场有 10,000,000 道题的社会化阅卷，阅卷老师有 39,063 个小组，每组 256 人。要给每位老师发一个唯一的工号，最快的办法不是让校长挨个点名，而是”组号 × 每组人数 + 组内座号“：第 7 组的第 3 号老师，工号就是
7 * 256 + 3 = 1795。CUDA 把这件事固化成了硬件：blockIdx.x是组号（硬件寄存器直接给出，不需要计算），blockDim.x是每组人数（一个编译期/运行期常量），threadIdx.x是组内座号（硬件给出）。三个值一读、一乘、一加，工号就有了。接下来只需要规定”工号 = 题目号”，映射就完成了。架构/机制图解：
CUDA 的三层逻辑结构（以 N = 10,000,000、block 大小 256 为例）
Grid（整个问题）
┌──────────────────────────────────────────────────────────────────────┐
│ dimGrid = (39063, 1, 1) -> 共 39063 个 thread block │
│ ┌───────────┬───────────┬───────────┬───────────┬─────────────────┐ │
│ │ blk 0 │ blk 1 │ blk 2 │ blk 3 │ blk 39062 │ │
│ └───────────┴───────────┴───────────┴───────────┴─────────────────┘ │
└──────────────────────────────────────────────────────────────────────┘
每个 block: dimBlock = (256, 1, 1)
┌──────────────────────────────────────────────────────────────────────┐
│ threadIdx.x: 0 1 2 3 4 5 6 7 │
│ 8 9 10 11 12 13 14 15 │
│ 128 129 130 131 132 133 134 135 │
│ 248 249 250 251 252 253 254 255 │
│ （完整取值是 0 到 255 的每一个整数） │
│ │
│ warp 划分 : warp0 = t0 至 t31 | warp1 = t32 至 t63 │
│ warp2 = t64 至 t95 | warp3 = t96 至 t127 │
│ warp4 = t128 至 t159 | warp5 = t160 至 t191 │
│ warp6 = t192 至 t223 | warp7 = t224 至 t255 │
│ 共 8 个 warp，每个 warp 32 个线程 │
└──────────────────────────────────────────────────────────────────────┘
线性化编号公式（"组号 × 每组人数 + 组内座号"）：
i = blockIdx.x * blockDim.x + threadIdx.x
└────┬────┘ └────┬────┘ └────┬────┘
组号 b 每组人数 B 组内座号 t
(0..39062) (=256) (0..255)
block 0 : t=0 -> i=0 t=1 -> i=1 t=255 -> i=255
block 1 : t=0 -> i=256 t=1 -> i=257 t=255 -> i=511
block 2 : t=0 -> i=512 t=1 -> i=513 t=255 -> i=767
block k : t=0 -> i=256k t=j -> i=256k+j
为什么必须做边界检查：N = 10,000,000 不是 256 的整数倍
ceil(N / 256) = ceil(39062.5) = 39063 个 block
实际启动线程数 = 39063 * 256 = 10,000,128 个
真正需要的线程数 = 10,000,000 个
多余的线程数 = 128 个（最后 4 个完整 warp 无事可做）
元素索引轴：
0 9,999,744 9,999,872 10,000,000 10,000,128
|------------------------|-----------|-----------|-----------|
全部有效（10,000,000 个元素） ↑ ↑
最后一个满 warp 越界区！
(32 个元素全有效) 若不加 if(i<N)
就会写坏别人的内存
block 39062 覆盖 i = 9,999,872 .. 10,000,127
├─ threadIdx.x 0 .. 127 -> i = 9,999,872 .. 9,999,999 （有效）
└─ threadIdx.x 128 .. 255 -> i = 10,000,000 .. 10,000,127（越界！）
关键操作与性能特征：
- 网格大小必须向上取整：
dimGrid.x = (N + blockDim.x - 1) / blockDim.x。整数除法N / blockDim.x是向下取整，当N % 256 != 0时会漏掉最后N % 256个元素（本讲 N=10,000,000 时会漏掉 128 个，正确性测试直接 FAIL）。 - 上取整的代价是多余的线程：39063 个 block 里有 128 个线程完全空转，占比
128 / 10,000,128 = 0.0013%，可忽略。但同样的代码在 N = 1025、block = 1024 时，启动 2 个 block 共 2048 个线程，只有 1025 个有效，浪费 49.9%。小问题用大 block 是浪费的重灾区。 - 越界写是最危险的 bug：CUDA 的
cudaMalloc返回的是一整块设备内存，越界写会静默破坏相邻分配（可能是别的数组、可能是 CUDA 内部结构），常常表现为”跑一次对、跑十次错”或”验证时正确、之后崩”。if (i < N)不是可选项，是必需品。 - 上取整请用整数运算：
(N + B - 1) / B对非负整数永远正确；ceil(N / 256.0)也能用，但它引入浮点、且当 N 超过 2^53 时会失准。Lab 中强烈建议用整数形式，同时注意N + B - 1在 N 接近INT_MAX时会溢出，安全写法是(int)(((long long)N + B - 1) / B)。
概念 2：合并访存（Coalesced Memory Access）
定义与目的：合并访存指的是”同一个 warp 内的 32 个线程，访问同一段连续且对齐的内存区域“，使得这 32 次各自 4 字节的请求能被硬件合并成一次 128 字节的内存事务（memory transaction）。它解决的问题是：显存（DRAM）的物理特性决定了它的最小开销单位是”连续的一整块”，而不是”任意一个字节”。如果你的访问模式是散的，硬件仍然要搬回 128 字节（或至少 32 字节的 sector），但其中只有很小一部分被用上，剩下的带宽全部浪费掉。
直观解释（”它是什么？”）：把显存想成一列固定长度（128 字节）的货运车厢。合并访存就是”32 个人凑钱买一节车厢，每个人的包裹都放进这节车厢的固定格子里”——付一次运费，所有人的货都运到了。非合并访存（strided）则是”32 个人各自包了一整节车厢，却每人只放了一个 4 字节的小盒子”——运费付了 32 份，实际只用了 1 份的货量。货运公司（DRAM 控制器）不会因为你只用了一个格子就少收钱，因为一趟车的开销（激活行、突发传输 burst）是固定的。课程讲义（Lecture 7 DRAM）里那个 512 字节 burst、240 GB/s 峰值的例子说得很清楚：如果只用到 burst 里每隔 8 字节的一小部分，可得带宽立刻减半到 120 GB/s。
架构/机制图解：
合并访存（coalesced）：一个 warp 的 32 个线程访问 32 个连续的 float
（float 为 4 字节，cudaMalloc 返回的基地址按 256 字节对齐）
线程号: t0 t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15
元素号: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
字节偏移: 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60
线程号: t16 t17 t18 t19 t20 t21 t22 t23 t24 t25 t26 t27 t28 t29 t30 t31
元素号: 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
字节偏移: 64 68 72 76 80 84 88 92 96 100 104 108 112 116 120 124
落入的 128 字节地址窗口（一次 memory transaction = 一条 128 B cache line）
字节 0 127
┌───────────────────────────────────────────────────────────────┐
│ 32 个 float 全部被使用（128 字节，利用率 100%） │
└───────────────────────────────────────────────────────────────┘
│<----------------- 1 次 transaction (128 B) ----------------->│
└─ sector 0 (32B) ─┴─ sector 1 (32B) ─┴─ sector 2 ─┴─ sector 3 ─┘
结果：一个 warp 读 A 只用 1 次事务；读 B 1 次；写 C 1 次
-> 每个 warp 3 次 128 字节事务，搬运 384 字节，全部是有效数据
跨步访存（strided）：线程号相差 1，地址相差 S = 156,252 字节
（对应代码 i = threadIdx.x * gridDim.x + blockIdx.x，即为"按列遍历"）
线程号: t0 t1 t2 t3 t4 t5 t30 t31
地址: b b+S b+2S b+3S b+4S b+5S b+30S b+31S
│ │ │ │ │ │ │ │
▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼
┌─────┐ ┌─────┐ ┌─────┐ ┌─────┐┌─────┐┌─────┐ ┌─────┐ ┌─────┐
│128 B│ │128 B│ │128 B│ │128 B││128 B││128 B│ │128 B│ │128 B│
│line │ │line │ │line │ │line ││line ││line │ │line │ │line │
└─────┘ └─────┘ └─────┘ └─────┘└─────┘└─────┘ └─────┘ └─────┘
每个 128 字节的 line 里只用了 4 个字节（利用率 3.125%）
（32 个线程各占一个 line；上图列出 t0 至 t5 与 t30、t31）
事务数：32 次 transaction（line 粒度）= 32 x 128 B = 4096 B 只用到 128 B
现代 GPU 按 32 字节 sector 粒度取数，则 32 个 sector x 32 B = 1024 B
只用到 128 B，放大倍数 = 1024 / 128 = 8 倍
事务数 / 有效字节数 对照表（一个 warp，32 个线程，各取 1 个 float）
访问模式 事务数 搬运字节 有效字节 利用率
─────────────────────────────────────────────────────────────────────────
连续、128 B 对齐（合并） 1 128 128 100.0%
连续、但只跨 128 B 边界一点 2 256 128 50.0%
步长 2 个 float（8 字节） 2 256 128 50.0%
步长 8 个 float（32 字节） 8 1024 128 12.5%
步长 32 个 float（128 字节） 32 4096 128 3.1%
─────────────────────────────────────────────────────────────────────────
注：sector 粒度的硬件（Kepler 及以后）会把上表后三行的搬运字节分别
降到 256 / 256 / 1024 字节，但"放大 2 倍 / 8 倍"的结论不变。
关键操作与性能特征：
- 128 字节是一次事务的粒度（对应一条 cache line），正好等于
32 线程 × 4 字节。这就是为什么 warp 大小 32 与 float 的宽度 4 字节如此契合：一个 warp 读一组连续 float，天然就是一整条 line。 - 对齐（alignment）和连续（contiguity）同样重要。如果数组起始地址不是 128 字节对齐，一个 warp 的 32 个 float 会跨两条 line，事务数从 1 变成 2，带宽利用率立刻掉到 50%。
cudaMalloc返回的指针至少 256 字节对齐，所以只要线性化编号公式里 threadIdx.x 是最内层（变化最快）的那一维，合并就是自动的——这正是讲义里那道思考题（B[ty*Width+tx]优于B[tx*Width+ty]）的答案。 - GPU 没有”每个线程各自缓存”的写缓冲，写操作同样按 warp 合并。写 C 时连续线程写连续地址 → 1 次事务；如果跨步写，写合并失败，代价同样是数倍放大的写带宽。
- L1/L2 会掩盖一部分不合并访问：L2 line 大小 128 字节、L1 sector 32 字节。若同一 warp 的另一次相邻加载能命中已被拉进 L1/L2 的 sector，放大倍数会降低。但不要指望缓存救场——向量加法每个字节只被读一次，没有任何时间局部性（temporal locality），缓存在这里唯一的作用就是 coalescing buffer。
概念 3：全局内存延迟与带宽（Global Memory Latency and Bandwidth）
定义与目的：全局内存（global memory，即设备 DRAM）是 GPU 上容量最大（数十 GB）、速度最慢的一级存储。它有两个正交的性能指标：延迟（latency）——一次访问从发出到数据可用需要多久，约 400–800 个时钟周期；带宽（bandwidth）——单位时间内能搬运多少字节，A100 约 1555 GB/s（峰值）。区分这两者极其重要：延迟决定了”你需要多少并发才能不挨饿”，带宽决定了”每秒最多能处理多少数据”。向量加法是典型的带宽上限问题：它需要搬运的字节数远超它需要的计算量，因此性能天花板就是显存带宽。
直观解释（”它是什么？”）：把显存想成城郊的大仓库，SM 想成市中心的工厂。从下单到货到（延迟）要开 500 个路口的车程（~400–800 cycles），但高速路的车道数（带宽）很宽。工厂的活是”把 A 箱和 B 箱的东西倒在一起”（一次加法），工序本身快得可以忽略。所以决定产量的是高速路每秒能过多少车，而不是工厂的手速。工厂的对策不是”催司机开快点”（延迟改不了，它是物理距离），而是一次发出几百辆卡车（用大量并发 warp 把延迟”藏起来”），让高速路时刻是满的。
架构/机制图解：
现代 NVIDIA GPU 的存储层次（数字为 A100 / sm_80 量级）
访问延迟 ↑ 容量 ↑
┌──────────────────────────────────────────────────────────────────────┐
│ 寄存器 Register File 每 SM 65536 个 32 位寄存器（256 KB） │
│ 延迟 ~1 cycle 带宽：每 SM 每 cycle 数百字节 │
├──────────────────────────────────────────────────────────────────────┤
│ 共享内存 Shared Memory / L1 每 SM 最多 164 KB（可配置划分） │
│ 共享内存延迟 ~20-30 cycles，L1 命中 ~30 cycles │
├──────────────────────────────────────────────────────────────────────┤
│ L2 Cache 全芯片共享（A100 约 40 MB，按 2 个分区各 20 MB） │
│ 延迟 ~200 cycles │
├──────────────────────────────────────────────────────────────────────┤
│ 全局内存 Global Memory (HBM2) 40 GB / 80 GB │
│ 延迟 ~400-800 cycles（1.41 GHz 下约 284-567 ns） │
│ 带宽 A100 约 1555 GB/s 峰值，实测流式 kernel 约 1200-1400 GB/s │
└──────────────────────────────────────────────────────────────────────┘
换算：1555 GB/s 的显存带宽，配上 500 cycles / 1.41 GHz = 355 ns 的延迟，
要让带宽打满，全芯片必须同时有（Little's Law）：
在途字节 = 带宽 x 延迟 = 1250e9 B/s x 355e-9 s ≈ 443,750 B ≈ 434 KB
分摊到 108 个 SM：每 SM 需要约 4.1 KB（约 32 条 128 B cache line）在途
向量加法的时间分解（N = 10,000,000，A100）
计算侧（几乎不花时间）：
总 FLOP 数 = 10,000,000 次加法 = 1.0e7 FLOP
FP32 峰值 = 19.5 TFLOPS
纯计算下限 = 1.0e7 / 19.5e12 = 0.51 us
存储侧（决定一切）：
读 A = 40 MB，读 B = 40 MB，写 C = 40 MB，合计 120 MB
峰值带宽下限 = 120e6 B / 1555e9 B/s = 77.2 us
实测（1.25 TB/s） = 120e6 B / 1250e9 B/s = 96.0 us
96.0 us / 0.51 us ≈ 188 倍 -> 计算单元 99.5% 的时间在空转
关键操作与性能特征：
- DRAM 的物理约束：DRAM 一位只用一个晶体管加一个电容，密度极高但很慢；每个 bit line 上挂约 1000 个单元，靠 sense amplifier 放大才能读出。现代 DRAM 强制以 burst 模式工作：一次行激活（row activate）之后连续吐出整块数据，即使你只要其中几个字节。这就是”访问连续地址才高效”的物理根源。
- 延迟无法通过优化代码降低，只能通过并发隐藏。提高并发有两条路：更多 warp（提高占用率）、每线程更多在途访存（thread coarsening / ILP）。本讲的概念 5 与概念 6 分别对应这两条路。
- 带宽是可以逼近但很难打满的：A100 上纯流式 kernel 通常能到峰值的 80%–92%。剩下 8%–20% 消耗在刷新（refresh）、读写转向（bus turnaround）、行激活开销和 ECC 上。所以”比峰值”这个指标应该对着实测可达带宽（约 1.3 TB/s）而不是理论峰值 1555 GB/s 来评估。
- 端到端时间常常被 PCIe 传输吃掉：40 MB 通过 PCIe Gen4 x16 传输（实测约 25 GB/s）要 ~1.6 ms，三个数组（A 上行、B 上行、C 下行）合计约 4.8 ms，是 96 µs kernel 时间的 50 倍。所以 Lab 的性能测量必须只用
cudaEvent圈住 kernel，把cudaMemcpy排除在外。
概念 4：算术强度（Arithmetic Intensity）与机器平衡点
定义与目的：算术强度（arithmetic intensity，AI）定义为
AI = 浮点运算次数 (FLOP) / 需要从显存搬运的字节数 (Byte)它刻画一段代码”每搬一个字节能算多少活”。机器平衡点（machine balance）定义为一个硬件”用光所有带宽正好也用光所有算力”的那个 AI 值：
AI_balance = 峰值 FLOPS / 峰值带宽。当AI < AI_balance时，程序是带宽受限的（算力过剩）；当AI > AI_balance时，程序是计算受限的。这个判据把”该优化访存还是该优化计算”这个模糊问题变成了一个不等式，是 Roofline 模型的核心。直观解释（”它是什么？”）：把 GPU 想成一家原料运输成本极高的工厂。机器平衡点就是”每运进来一吨原料，工厂的机器正好能全部忙满”的那个原料/加工比例。向量加法是”运进来两吨原料（A、B），只做一道工序（一次加法），再运出去一吨成品（C）”——典型的”拉货比干活累”。你就算把机器换成快 10 倍的，产量也只涨 0%（因为瓶颈在路上）。反过来，矩阵乘法（Lab 2）每搬 8 字节能做 2 FLOP 且通过分块复用把访存降低 TILE_WIDTH 倍，所以它有希望变成计算受限——这正是为什么课程要用共享内存去优化它。
架构/机制图解：
Roofline 模型（A100，FP32，对数坐标）
可达性能
(GFLOP/s)
^
可达性能 GFLOP/s（双对数刻度，示意图）
^
1e4 | ________________________
| _____/ 计算天花板 = 19.5 TFLOPS
1e3 | _____/
| _____/
1e2 |★_____________________/ <-- 斜线斜率 = 显存带宽 1555 GB/s
1e1 |/
1e0 +----|------------------------|---------------------------|----------> AI
0.0833 12.54 1000
(0.125 同量级) 机器平衡点（两线交点） (FLOP/Byte)
^
★ 处高度 = 0.0833 x 1555 GB/s = 129.6 GFLOP/s
与计算天花板 19500 GFLOP/s 相差 150 倍 -> 纯带宽受限
| ^ ^
| | |
| 向量加法落在这里 机器平衡点 = 19500/1555
| AI = 0.125 FLOP/Byte = 12.54 FLOP/Byte
| 离平衡点差 100 倍 两条线在这里相交
|
0 +-----|---------|---------------------------------|---------------> AI
0.0833 0.125 1 10 12.54 100
^ (FLOP/Byte)
向量加法（只算读写）
0.125 = 1 FLOP / 8 Byte
（若把写 C 也算进去：1/12 = 0.0833）
三项定量计算
(1) 向量加法的算术强度
每元素：1 次 FP32 加法 = 1 FLOP
读 A 4 B + 读 B 4 B = 8 B（把写 C 也算进去 = 12 B）
AI = 1 / 8 = 0.125 FLOP/Byte （课程惯用口径）
AI = 1 / 12 = 0.0833 FLOP/Byte （含写带宽的严格口径）
(2) A100 的机器平衡点
AI_balance = 19.5e12 FLOP/s / 1555e9 B/s = 12.54 FLOP/Byte
对比：RTX 4090 = 82.6e12 / 1008e9 = 81.9 FLOP/Byte
H100 SXM = 67e12 / 3350e9 = 20.0 FLOP/Byte
(3) 结论（Roofline 天花板）
可达性能 = min(峰值 FLOPS, AI x 带宽)
= min(19.5e12, 0.0833 x 1555e9)
= min(19.5e12, 1.296e11) = 129.6 GFLOP/s
即向量加法最多只能跑到 129.6 / 19500 = 0.66% 的 FP32 峰值算力。
用 0.125 的口径算：0.125 x 1555e9 = 194.4 GFLOP/s = 1.0% 峰值。
0.125 与 12.54 相差 100.3 倍 -> 100% 带宽受限，毫无悬念。
想让它在 A100 上变成计算受限，需要 AI >= 12.54，
即每元素至少要做 12.54 x 12 ≈ 150 次浮点运算。
关键操作与性能特征：
- 算术强度是 kernel 的固有属性，与实现的快慢无关。同一个向量加法，写得再烂 AI 还是 0.125；只有改变算法（例如把两个向量合并成一次 float4 向量化读写，AI 不变但事务效率提高）或增加数据复用（例如矩阵乘法分块）才能改变 AI。
- Roofline 给出的是上限，不是预测。它告诉你”天花板在哪”，不告诉你”离天花板多远”。实测 1.25 TB/s 对应
1.25e12 x 0.0833 = 104 GFLOP/s，达到 Roofline 天花板的 80.4%——这正是我们期望的量级。 - 判断瓶颈类型的三步法：算 AI → 算机器平衡点 → 比大小。AI 远小于平衡点 = 带宽受限（优化方向：减少字节数、保证合并、提高在途访存）；AI 远大于平衡点 = 计算受限（优化方向：减少指令数、用 FMA、提高 ILP）；两者接近 = 需要同时优化，或者说明代码离两者都很远，那是延迟受限（优化方向：提高占用率/ILP 来隐藏延迟）。
- 向量加法的 AI 低到无法抢救：即使实现完美，它的时间下限就是
120 MB / 1555 GB/s = 77.2 µs。任何超过这个数的时间都是”没打满带宽”的浪费，而不是”算得太慢”。
概念 5：占用率（Occupancy）
定义与目的：占用率（occupancy）定义为
occupancy = 每个 SM 上实际驻留的活跃 warp 数 / 每个 SM 硬件支持的最大 warp 数它是”硬件有多少并发可以用来隐藏内存延迟”的度量。GPU 不像 CPU 那样靠乱序执行和深层 cache 隐藏延迟，它靠的是线程级并行（TLP）：一个 warp 发出一条访存指令后需要等 400–800 cycles，调度器立刻切到另一个就绪 warp 上执行。驻留的 warp 越多，越有可能在任何一个周期都找到”操作数已就绪”的 warp。占用率会影响能藏住多少延迟，因此直接影响带宽受限 kernel 能否打满带宽。
直观解释（”它是什么？”）：把 SM 想成一个有 64 条电话线的小呼叫中心（64 = A100 每 SM 最大 warp 数），打一个电话出去（发起一次全局内存访问）平均要等 5 分钟才有人接（400–800 cycles）。如果只雇 8 个话务员（8 个 warp），他们全部在等电话，办公室几乎全空。雇满 64 个（100% 占用率），任何时刻都有人正在通话，吞吐就上去了。但注意：话务员再多也不能让单个电话变快（延迟是物理的），只能让总吞吐变高。而且一旦电话线路（带宽）占满了，再雇人也没用——这就是为什么向量加法在 100% 占用率下也就跑到带宽上限。
架构/机制图解：
一个 SM 的资源账本（A100 / sm_80）—— 占用率由"最紧的那个资源"决定
┌──────────────────────────────────────────────────────────────────────┐
│ SM (Streaming Multiprocessor) │
│ │
│ ┌────────────────────────────────────────────────────────────────┐ │
│ │ 4 个 warp scheduler（每个每 cycle 可发射 1 条指令） │ │
│ │ WS0 ──> 32 个 FP32 核心 + LD/ST + SFU │ │
│ │ WS1 ──> 32 个 FP32 核心 + LD/ST + SFU │ │
│ │ WS2 ──> 32 个 FP32 核心 + LD/ST + SFU │ │
│ │ WS3 ──> 32 个 FP32 核心 + LD/ST + SFU │ │
│ └────────────────────────────────────────────────────────────────┘ │
│ │
│ 寄存器文件 Register File : 65536 个 32-bit 寄存器 = 256 KB │
│ 共享内存 / L1 : 164 KB（可划分，最多 100 KB 给 shared） │
│ 线程上限 : 2048 threads = 64 warps │
│ block 上限 : 32 blocks / SM │
│ │
│ 驻留的 warp（以 256 线程的 block 为例，每 block 切成 8 个 warp）： │
│ ┌──────┬──────┬──────┬──────┬──────┬──────┬──────┬──────┬─────────┐ │
│ │ blk0 │ blk0 │ blk1 │ blk1 │ blk2 │ blk2 │ blk3 │ blk3 │ blk4 至 │ │
│ │warp0 │warp1 │warp2 │warp3 │warp4 │warp5 │warp6 │warp7 │ blk7 共 │ │
│ │ │ │ │ │ │ │ │ │ 32 warp │ │
│ └──────┴──────┴──────┴──────┴──────┴──────┴──────┴──────┴─────────┘ │
│ 合计 8 个 block x 8 warp = 64 warp = 2048 线程（占满 SM 上限） │
└──────────────────────────────────────────────────────────────────────┘
占用率的三个限制因素（取最小者）：
① 线程/warp 上限： blocks_by_threads = floor(2048 / blockDim)
② 寄存器上限 ： blocks_by_regs = floor(65536 / (regs_per_thread * blockDim))
③ block 数上限 ： blocks_by_blocklimit = 32
resident_blocks = min(①, ②, ③)
occupancy = resident_blocks * blockDim / 2048
A100 上不同 block 大小的占用率（vecAddKernel，实测约 16 个寄存器/线程）
blockDim 寄存器限制 线程限制 block 限制 驻留 block 驻留线程 warp 占用率
─────────────────────────────────────────────────────────────────────────────
128 65536/(16*128) 2048/128 32 16 2048 64 100.0%
= 32 = 16
256 65536/(16*256) 2048/256 32 8 2048 64 100.0%
= 16 = 8
512 65536/(16*512) 2048/512 32 4 2048 64 100.0%
= 8 = 4
1024 65536/(16*1024) 2048/1024 32 2 2048 64 100.0%
= 4 = 2
─────────────────────────────────────────────────────────────────────────────
结论：向量加法这种轻量 kernel，寄存器用量很低，block 大小在 128-1024 之间
都能达到 100% 占用率。占用率不是它的问题。
反例：一个每线程用 64 个寄存器的 kernel（例如分块矩阵乘法的某个版本）
blockDim = 256 -> blocks_by_regs = 65536/(64*256) = 4
-> 驻留 1024 线程 = 32 warp -> 占用率 50.0%
blockDim = 512 -> blocks_by_regs = 65536/(64*512) = 2
-> 驻留 1024 线程 = 32 warp -> 占用率 50.0%
最终被减少寄存器（__launch_bounds__）或降低 block 大小也救不回来。
block 大小选择的另一个维度：尾部浪费 与 block 数上限
场景 A：讲义中的 1536 线程 / 最多 8 blocks 的 SM（Maxwell 类似配置）
8x8 = 64 线程/block -> 1536/64 = 24 blocks，但上限 8
-> 只有 8*64 = 512 线程 = 16 warp = 33.3% 占用率
16x16 = 256 线程/block -> 1536/256 = 6 blocks（<= 8）
-> 6*256 = 1536 线程 = 48 warp = 100% 占用率
32x32 = 1024 线程/block -> 1536/1024 = 1 block
-> 1*1024 = 1024 线程 = 32 warp = 66.7% 占用率
场景 B：A100（2048 线程 / 64 warp）与 RTX 4090（1536 线程 / 48 warp）
blockDim A100 驻留 A100 占用率 RTX4090 驻留 RTX4090 占用率
─────────────────────────────────────────────────────────────────────
128 16 blk/2048 100.0% 12 blk/1536 100.0%
256 8 blk/2048 100.0% 6 blk/1536 100.0%
512 4 blk/2048 100.0% 3 blk/1536 100.0%
1024 2 blk/2048 100.0% 1 blk/1024 66.7%
─────────────────────────────────────────────────────────────────────
1024 线程的大 block 在 1536 线程/SM 的 Ada 上被"阶梯效应"卡住：
1536/1024 = 1.5 -> 只能放 1 个，白白丢掉 1/3 的线程容量。
场景 C：尾部浪费（N 与 blockDim 不整除时）
N = 1025, blockDim = 256 -> 5 blocks -> 1280 线程 -> 浪费 255/1280 = 19.9%
N = 1025, blockDim = 1024 -> 2 blocks -> 2048 线程 -> 浪费 1023/2048 = 49.9%
N = 10,000,000, blockDim = 256 -> 浪费 128/10,000,128 = 0.0013%（可忽略）
关键操作与性能特征：
- 占用率不是越高越好，而是”够用就好”。对带宽受限 kernel，Little’s Law 给出的”够用”标准是：每 SM 在途字节 ≥ 带宽 × 延迟 / SM 数 ≈ 4.1 KB。一个 warp 一次合并读 = 128 B，所以需要 4096/128 = 32 个同时有在途访存的 warp。A100 上 64 个 warp 每 warp 一个在途访问 = 8192 B，有 2 倍余量；16 个 warp（25% 占用率）就只有 2048 B，只够一半带宽。
- 占用率可以用 ILP 换。如果每线程同时发 4 个独立访存（thread coarsening），那么 16 个 warp × 4 × 128 B = 8192 B，同样打满带宽。这就是概念 6 与示例 2 的意义：低占用率 + 高 ILP 可以等价于高占用率 + 低 ILP。
- 寄存器是最常见的隐形限制。
nvcc -Xptxas -v会打印每个 kernel 的寄存器用量；cudaOccupancyMaxPotentialBlockSize()可以反推推荐的 block 大小；__launch_bounds__(maxThreadsPerBlock, minBlocksPerMultiprocessor)可以强制编译器压低寄存器用量（代价是可能溢出到本地内存——本地内存也在显存里，反而增加流量，要实测确认）。 - 占用率用 API 精确计算，不要凭感觉：
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&n, kernel, blockSize, dynamicSMem)给出”每个 SM 能同时驻留几个 block”，乘以blockSize/32就是活跃 warp 数，再除以prop.maxThreadsPerMultiProcessor/32就是占用率。示例 1 和示例 2 的代码里都打印了这三个数。
概念 6：网格-步长循环（Grid-Stride Loop）与线程粗化（Thread Coarsening）
定义与目的：网格-步长循环是把”每线程一个元素”改成”每线程若干元素”的标准写法：
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += blockDim.x * gridDim.x) C[i] = A[i] + B[i];循环步长
blockDim.x * gridDim.x（即”网格总线程数”）保证：每一轮迭代中，整个网格恰好覆盖一段连续的blockDim.x * gridDim.x个元素，且 warp 内部依然连续——合并访存被完整保留。它解决的问题是”网格尺寸与问题规模的解耦”：网格大小不再由 N 决定，而是由设备（SM 数、每 SM 可驻留 block 数）决定，内核代码对任意 N 都成立。线程粗化（thread coarsening）是它的自然延伸：让每个线程每轮处理UNROLL个元素，从而在一个线程内制造多个独立的在途访存（ILP / memory-level parallelism）。直观解释（”它是什么？”）：把”每线程一个元素”想成雇 10,000,000 个临时工，每人搬一块砖，搬完就散——招工、分配工位、结账（block 调度与回收）的开销占比很高，而且你必须在开工前精确知道有多少块砖。网格-步长循环则是雇固定的一批长工（比如 864 个班组长 × 256 人），每人循环地搬第 i 块、第 i+stride 块……砖变多了不用重新招人，砖变少了也不会有人闲着。线程粗化则进一步让每个工人一次推四辆手推车（4 个独立访存同时在途），即使工厂里工人少一些（低占用率），路上的车流量（在途字节数）依然足够把高速路（带宽）占满。
架构/机制图解：
网格-步长循环的迭代覆盖图（gridDim.x = 2, blockDim.x = 4, 所以 stride = 8）
元素索引: 0 1 2 3 4 5 6 7 | 8 9 10 11 12 13 14 15
线程号 : t0 t1 t2 t3 t4 t5 t6 t7 | t0 t1 t2 t3 t4 t5 t6 t7
\──────── 迭代 0（i = 0..7）──/ \──── 迭代 1（i = 8..15，i += 8 后）
元素索引: 16 17 18 19 20 21 22 23 | 24 25 26 27 28 29 30 31
线程号 : t0 t1 t2 t3 t4 t5 t6 t7 | t0 t1 t2 t3 t4 t5 t6 t7
\──────── 迭代 2（i = 16..23）─/ \──── 迭代 3（i = 24..31）
关键性质：每一轮的 8 个线程访问 8 个连续元素（warp 内合并），
相邻两轮之间相隔 stride 个元素（不是同一个 warp 的地址区间）
所以合并性完全保留，且网格大小可以任意小。
线程粗化（UNROLL = 4）如何制造在途访存（ILP）
每线程每轮发出 4 个"互相独立"的加载，地址相隔 stride：
线程 i 的 4 个地址:
A[i] A[i+stride] A[i+2*stride] A[i+3*stride]
│ │ │ │
▼ ▼ ▼ ▼
LDG.128 LDG.128 LDG.128 LDG.128
└───────────────────┴────────┬────────┴─────────────────┘
│ 4 条 LDG 之间没有数据依赖
▼
全部在途（in flight），
4 x 128 B = 512 B / warp 同时等 DRAM
-> 用 1/4 的 warp 数就能达到同样的在途字节数
Little's Law 视角（A100，每 SM 需约 4.1 KB 在途才能打满带宽）：
✓ 64 warp x 1 个在途访问 x 128 B = 8192 B （占用率 100%，UNROLL=1）
✗ 16 warp x 1 个在途访问 x 128 B = 2048 B （占用率 25%，带宽只剩一半）
✓ 16 warp x 4 个在途访问 x 128 B = 8192 B （占用率 25% + UNROLL=4，带宽打满）
✓ 24 warp x 4 个在途访问 x 128 B = 12288 B （80 寄存器 kernel 也能打满）
网格-步长循环 vs 巨型网格：波次量化（wave quantization）效应
巨型网格（N = 1,000,000, blockDim = 256 -> 3907 个 block）
A100 每波可容纳 108 SM x 8 block = 864 个 block：
波次 1: [======================== 864 blocks ========================]
波次 2: [======================== 864 blocks ========================]
波次 3: [======================== 864 blocks ========================]
波次 4: [======================== 864 blocks ========================]
波次 5: [========= 451 blocks =========] <- 剩 413 个 block 槽位空闲
有效波次数 = 4 + 451/864 = 4.52；实际用了 5 波
利用率 = 4.52 / 5 = 90.4% -> 约 9.6% 的时间浪费在最后一波的尾部
网格-步长循环（grid = 108 x 8 = 864 个 block，永久驻留）
┌─────────────────────────────────────────────────────────────────┐
│ 864 个 block 一直驻留，每个 block 内部循环执行约 4.5 轮 │
│ 负载不均被摊到"线程级的最后半轮"，SM 之间不会出现整块空转 │
└─────────────────────────────────────────────────────────────────┘
另外省下 3907 - 864 = 3043 次 block 启动/回收的调度开销
关键操作与性能特征：
- 网格-步长循环解决的是可扩展性，而不是”更快的单次吞吐”。对 N = 10,000,000 这种大问题，两种写法的实测带宽几乎一样（都受限于 DRAM）。它的价值在于：同一份二进制在 8 SM 的笔记本 GPU 和 108 SM 的 A100 上都接近最优（网格大小在运行时由
prop.multiProcessorCount决定），不需要为不同 GPU 重新编译或重新调参。 - 多 GPU 与未知规模：多 GPU 编程里每个 GPU 分到 N/G 个元素，网格可以固定为”本设备能放下的最大 block 数”；问题规模在运行时才知道（比如读文件、用户输入）时，无需重新计算网格。
- 粗化的收益随”每线程工作量的独立性”增长：如果 UNROLL 个访存之间存在依赖（例如
C[i] = C[i-stride] + 1），编译器就无法把它们重叠，粗化不仅无益还可能因寄存器压力降低占用率。向量加法的 UNROLL 个元素完全独立，是最理想的粗化对象。 - 粗化的代价：寄存器用量随 UNROLL 增长（示例 2 中从约 16 涨到约 24），block 内线程数不变时每 SM 能驻留的 block 数可能下降。经验做法是 UNROLL = 2 或 4，并用
cudaOccupancyMaxActiveBlocksPerMultiprocessor实测确认占用率没有跌破 Little’s Law 要求的下限。 - 注意 int 溢出：
i + (UNROLL-1)*stride在 N 接近INT_MAX时会溢出成负数，导致漏算。安全写法是把守卫条件用long long计算（示例 2 的代码即如此），或对超大规模问题直接把索引类型改成long long/size_t。
代码示例与性能分析
示例 1：基础版向量加法 vecAdd（每线程一个元素 + 完整测量框架）
// 文件: vecadd_basic.cu
// 编译: nvcc -O3 -arch=sm_80 vecadd_basic.cu -o vecadd_basic
// 运行: ./vecadd_basic 10000000
//
// 内容: 每线程一个元素的向量加法；cudaEvent 计时；CPU 参考实现验证；
// 打印有效带宽 / 理论峰值带宽比例 / 访存事务数 / 占用率。
#include <cstdio>
#include <cstdlib>
#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 BLOCK_SIZE 256
// ---------------------------------------------------------------- kernel ----
// 每线程处理一个元素：线程的全局线性编号 = 它的元素编号
__global__ void vecAddKernel(const float* __restrict__ A,
const float* __restrict__ B,
float* __restrict__ C,
int N)
{
int i = blockIdx.x * blockDim.x + threadIdx.x; // 线程 -> 数据 的映射
if (i < N) { // 边界检查：必须！
C[i] = A[i] + B[i]; // 两次读 + 一次写
}
}
// ------------------------------------------------------- CPU 参考实现 -------
// golden reference：与 GPU 完全相同的运算顺序（每个元素一次加法），
// 因此结果应当逐位相同（bit-exact），可以用 == 比较。
static void vecAddCPU(const float* A, const float* B, float* C, int N)
{
for (int i = 0; i < N; ++i) {
C[i] = A[i] + B[i];
}
}
// 确定性伪随机初始化（线性同余），避免依赖 <random> 的实现差异
static void initVector(float* v, int N, unsigned int seed)
{
unsigned int x = seed;
for (int i = 0; i < N; ++i) {
x = x * 1664525u + 1013904223u;
v[i] = (float)(x >> 8) / (float)(1u << 24) - 0.5f; // 落在 [-0.5, 0.5)
}
}
int main(int argc, char** argv)
{
int N = (argc > 1) ? atoi(argv[1]) : 10000000;
if (N <= 0) { fprintf(stderr, "N must be positive\n"); return EXIT_FAILURE; }
size_t bytes = (size_t)N * sizeof(float);
// 1. 主机端内存与初始化
float* h_A = (float*)malloc(bytes);
float* h_B = (float*)malloc(bytes);
float* h_C = (float*)malloc(bytes);
float* h_Ref = (float*)malloc(bytes);
if (!h_A || !h_B || !h_C || !h_Ref) {
fprintf(stderr, "host malloc failed\n");
return EXIT_FAILURE;
}
initVector(h_A, N, 12345u);
initVector(h_B, N, 67890u);
// 2. 设备端内存分配与拷贝
float *d_A = NULL, *d_B = NULL, *d_C = NULL;
CUDA_CHECK(cudaMalloc((void**)&d_A, bytes));
CUDA_CHECK(cudaMalloc((void**)&d_B, bytes));
CUDA_CHECK(cudaMalloc((void**)&d_C, bytes));
CUDA_CHECK(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice));
// 3. 网格配置：向上取整，用整数运算避免浮点误差
dim3 dimBlock(BLOCK_SIZE, 1, 1);
dim3 dimGrid((unsigned)(((long long)N + BLOCK_SIZE - 1) / BLOCK_SIZE), 1, 1);
// 4. 预热（第一次启动包含 context 建立、模块加载等一次性开销）
vecAddKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, N);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
// 5. 计时：cudaEvent 只圈住 kernel，不含 cudaMemcpy
const int iters = 20;
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(cudaEventRecord(start));
for (int it = 0; it < iters; ++it) {
vecAddKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, N);
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError()); // 捕获启动期错误（如同步错误）
float msTotal = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&msTotal, start, stop));
double msPerIter = (double)msTotal / iters;
// 6. 取回结果并与 CPU 参考实现对比
CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));
vecAddCPU(h_A, h_B, h_Ref, N);
int mismatch = 0;
for (int i = 0; i < N; ++i) {
if (h_C[i] != h_Ref[i]) { ++mismatch; }
}
// 7. 设备属性与占用率
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
double peakGBs = 2.0 * (double)prop.memoryClockRate * 1.0e3
* ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
cudaFuncAttributes attr;
CUDA_CHECK(cudaFuncGetAttributes(&attr, vecAddKernel));
int blocksPerSM = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&blocksPerSM, vecAddKernel, BLOCK_SIZE, 0));
int warpsPerSM = blocksPerSM * (BLOCK_SIZE / 32);
int maxWarpsPerSM = prop.maxThreadsPerMultiProcessor / 32;
double occupancy = 100.0 * warpsPerSM / (double)maxWarpsPerSM;
// 8. 性能指标
double totalBytes = 3.0 * (double)bytes; // 读 A + 读 B + 写 C
double achievedGBs = totalBytes / (msPerIter * 1.0e-3) / 1.0e9;
double transPerWarp = 3.0; // 每个 warp: A/B/C 各 1 次
double numWarps = (double)((N + 31) / 32); // 每数组的事务数
double totalTrans = 3.0 * numWarps;
printf("================= 设备信息 =================\n");
printf("GPU : %s (sm_%d%d)\n", prop.name,
prop.major, prop.minor);
printf("SM 数 : %d\n", prop.multiProcessorCount);
printf("每 SM 最大线程/共享内存 : %d 线程 / %.0f KB\n",
prop.maxThreadsPerMultiProcessor,
prop.sharedMemPerMultiprocessor / 1024.0);
printf("理论峰值显存带宽 : %.1f GB/s\n", peakGBs);
printf("================= 启动配置 =================\n");
printf("N : %d 个 float (%.1f MB)\n",
N, bytes / 1048576.0);
printf("blockDim = %d, dimGrid = %u\n", BLOCK_SIZE, dimGrid.x);
printf("启动线程总数 : %u (其中 %u 个越界空转)\n",
dimGrid.x * BLOCK_SIZE, dimGrid.x * BLOCK_SIZE - (unsigned)N);
printf("================= 正确性验证 ===============\n");
printf("CPU 参考实现对比 : %s (不匹配 %d / %d)\n",
mismatch == 0 ? "PASS" : "FAIL", mismatch, N);
printf("================= 性能测量 =================\n");
printf("kernel 单次耗时 : %.3f us\n", msPerIter * 1000.0);
printf("搬运字节数 : %.1f MB (%.1f MB 读 + %.1f MB 写)\n",
totalBytes / 1e6, (double)bytes * 2 / 1e6, (double)bytes / 1e6);
printf("有效带宽 : %.1f GB/s\n", achievedGBs);
printf("占理论峰值带宽 : %.1f %%\n", 100.0 * achievedGBs / peakGBs);
printf("每个 warp 的访存事务数 : %.0f (128 B x 3)\n", transPerWarp);
printf("总访存事务数 : %.0f 次 (%.0f 个有效 warp)\n",
totalTrans, numWarps);
printf("每个事务搬运 : 128 B，有效字节 128 B，利用率 100%%\n");
printf("================= 资源与占用率 =============\n");
printf("每线程寄存器数 : %d\n", attr.numRegs);
printf("每 SM 驻留 block 数 : %d\n", blocksPerSM);
printf("每 SM 活跃 warp 数 : %d / %d\n", warpsPerSM, maxWarpsPerSM);
printf("占用率 (occupancy) : %.1f %%\n", occupancy);
// 9. 释放资源
CUDA_CHECK(cudaFree(d_A));
CUDA_CHECK(cudaFree(d_B));
CUDA_CHECK(cudaFree(d_C));
free(h_A); free(h_B); free(h_C); free(h_Ref);
CUDA_CHECK(cudaDeviceReset());
return (mismatch == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}
【代码做什么？】
- 映射：
int i = blockIdx.x * blockDim.x + threadIdx.x;把”第 b 个 block 的第 t 号线程”翻译成”第 i 个元素”。这个编号是全局唯一的（因为(b, t)与i是一一对应的：b = i / blockDim.x，t = i % blockDim.x），所以每个元素恰好被一个线程处理一次。 - 边界检查：
if (i < N)拦住最后一批越界线程。N = 10,000,000、block = 256 时dimGrid.x = ceil(10000000/256) = 39063，启动 10,000,128 个线程，多出的 128 个线程必须被拦住。 - 访存：
C[i] = A[i] + B[i]编译成两条LDG（global load）加一条STG（global store）。三个数组的索引相同，所以三次访存的地址模式完全一致，都是”相邻线程访问相邻地址”。 - 内存分配与访问路径：
cudaMalloc在设备全局内存中分配三块各 40 MB 的空间；cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice)把 A、B 送上去；kernel 直接在全局内存上读写，不经过共享内存（全局内存对于设备代码是”per-grid”可见的，主机也能通过cudaMemcpy访问它）。 - 主机端调度：主机按顺序做”分配 → 拷贝 → 预热 → 计时 20 次 → 同步 → 取回 → 验证 → 释放”。
kernel<<<dimGrid, dimBlock>>>(args)这种启动是异步的，所以计时循环里 20 次 kernel 会连续排队执行；cudaEventSynchronize(stop)才真正等待它们全部结束。预热那一次是必要的：第一次启动包含 CUDA context 创建、模块加载、指令 cache 填充等一次性开销，不预热会让第一次测量偏大几十到几百微秒。 - 计时口径：
cudaEventRecord只圈住 kernel 循环，把cudaMemcpy排除在外——这正是 Lab 想要测的”kernel 有效带宽”。若把传输也算进去，40 MB 数据走 PCIe Gen4 x16（实测约 25 GB/s）需要约 1.6 ms，三次传输约 4.8 ms，会把 96 µs 的 kernel 时间完全淹没。 - 验证：先
cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost)取回，再用同一个公式在 CPU 上算一遍h_Ref。因为每个元素都只做一次浮点加法、没有累加顺序问题，GPU 与 CPU 的结果应当逐位相同，所以可以直接用!=比较而不需要 epsilon。 - 指标打印：有效带宽 =
3 × N × 4 字节 / kernel 时间；理论峰值带宽由2 × memoryClockRate(kHz) × 1000 × (memoryBusWidth/8)算出（系数 2 来自 DDR 双沿传输），A100 上得到 1.215 GHz × 2 × 640 B = 1555 GB/s。
- 映射：
【并行机制与硬件映射解说】
（a）block 如何被切成 warp：block 大小 256，硬件按线性化的 threadIdx 每 32 个线程切成一个 warp：threads 0–31 = warp 0，32–63 = warp 1，依次类推，每个 block 得到
256/32 = 8个 warp。注意这是硬件实现细节，不是 CUDA 编程模型的一部分（将来 warp 大小可能不是 32，所以代码里不要硬编码 32 以外的假设）。threadIdx的线性化顺序是 x 最快、其次 y、最后 z，这是”最内层维度应该对应连续内存”这条规则的来源。（b）warp 如何被调度：39063 个 block 以 block 为粒度分发到 108 个 SM 上。每个 SM 有 4 个 warp scheduler，每个 scheduler 每周期可以从它负责的 warp 池中挑一个”下一条指令的操作数已就绪”的 warp 发射一条指令。SM 维护着数百个 warp 的 PC 和寄存器状态，切换 warp 零开销（这正是”用海量线程隐藏延迟”能成立的前提）。整个 kernel 共
39063 × 8 = 312,504个 warp，其中最后 4 个 warp 因i >= N而整体空转。（c）A100 上同时能驻留多少：A100 每 SM 最多 2048 线程 / 64 warp / 32 个 block，寄存器文件 65536 个 32-bit 寄存器，共享内存 164 KB。
vecAddKernel实测每线程约 16 个寄存器、不使用共享内存，于是：- 线程限制：
2048 / 256 = 8个 block； - 寄存器限制：
65536 / (16 × 256) = 16个 block； - block 数限制：32。 取最小者 → 每 SM 驻留 8 个 block = 2048 线程 = 64 warp = 100% 占用率。全芯片同时驻留
108 × 8 = 864个 block。39063 个 block 因此需要39063 / 864 = 45.2个波次（实际调度 46 波）。
（d）共享内存与 bank conflict：本 kernel 完全不使用共享内存，所以没有 bank conflict 问题。这一点值得强调：向量加法没有任何数据复用（每个字节只被读一次），把数据搬进共享内存只会增加一次额外的读写，是纯损失。共享内存（32 个 bank，每个 bank 宽 4 字节）是为”一个数据被多个线程重复使用“的场景准备的，第 4 讲之后的 tiling 才会用到。
（e）全局内存是否合并：完全合并。以 block 0 的 warp 0（threads 0–31）为例，三个访问分别是
A[0..31]、B[0..31]、C[0..31]，每个都是 128 个连续的字节。cudaMalloc保证返回的基地址至少 256 字节对齐，所以这 128 字节恰好落在一个（或整好跨一条边界但完全对齐的）128 字节事务里：- 每 warp 事务数：3 次（A 一次、B 一次、C 一次）；
- 每次事务 128 字节，有效字节 128 字节，利用率 100%；
- 总 warp 数（有效）：
N / 32 = 10,000,000 / 32 = 312,500； - 总事务数 = 312,500 × 3 = 937,500 次，搬运
937,500 × 128 B = 120,000,000 B = 120 MB； - 其中 80 MB 是读（A、B），40 MB 是写（C）。写操作同样按 warp 合并成 128 字节的写事务。
（f）寄存器用量与 warp 发散：每线程 16 个寄存器（两个地址寄存器、一个 N、一个 i、三个数据寄存器加若干临时量）。没有 warp 发散：
if (i < N)这个分支里，同一个 warp 的 32 个线程要么全部成立、要么全部不成立（因为 i 在 warp 内连续），只有在最后那个跨 N 边界的 warp 里才会出现部分线程不成立的情况——而由于 10,000,000 是 32 的整数倍，本次运行连这种情况都没有。这提示了一个通用技巧：只要让 N 是 32 的倍数（或把边界条件设计成 warp 对齐），就能彻底消除这类发散。- 线程限制：
【性能优化分析】
（1）算术强度与 Roofline
AI = 1 FLOP / 8 Byte = 0.125 FLOP/Byte （只算 A、B 的读） = 1 FLOP / 12 Byte = 0.0833 FLOP/Byte （把写 C 也算进流量） A100 机器平衡点 = 19.5e12 / 1555e9 = 12.54 FLOP/Byte AI 比平衡点低 12.54 / 0.125 = 100.3 倍 Roofline 天花板 = min(19.5 TFLOPS, 0.0833 x 1555 GB/s) = min(19500 GFLOP/s, 129.6 GFLOP/s) = 129.6 GFLOP/s 实测算力 = 0.0833 x 1250 GB/s = 104 GFLOP/s -> 仅用掉 104 / 19500 = 0.53% 的 FP32 峰值算力：纯带宽受限。另一个等价视角：A100 有
108 SM × 128 FP32 lane = 13,824条 FP32 流水线，N = 10,000,000 时每条流水线只需处理10,000,000 / 13,824 = 723个元素；而实测耗时 96 µs，在 1.41 GHz 下是96e-6 × 1.41e9 = 135,360个周期。FP32 流水线只有 723 / 135360 = 0.53% 的周期在干活。（2）占用率与 Little’s Law（延迟隐藏是否足够）
要让 1250 GB/s 的可达带宽打满，全芯片需要多少在途字节？ 在途字节 = 带宽 x 延迟 = 1250e9 B/s x (500 cycles / 1.41e9 Hz) = 1250e9 x 3.546e-7 = 443,250 B ≈ 433 KB 每 SM 需要的在途字节 = 443,250 / 108 = 4,104 B ≈ 4.1 KB 本 kernel 提供的能力： 64 warp/SM x 1 个在途合并访问 x 128 B = 8,192 B 每 SM 余量 = 8192 / 4104 = 2.0 倍 -> 足以掩盖延迟 占用率账本： 每 SM 8 个 block x 256 线程 = 2048 线程 = 64 warp occupancy = 64 / 64 = 100.0%结论：占用率不是本 kernel 的瓶颈——即使 block 大小改成 128/512/1024，A100 上都还是 100%。真正的瓶颈是 DRAM 的 1555 GB/s 上限。
（3）实测带宽与峰值之比
有效带宽 = 120,000,000 B / 96.0e-6 s = 1.25e12 B/s = 1250 GB/s 占理论峰值 = 1250 / 1555 = 80.4 % 占实测可达峰值（约 1.35 TB/s 纯流式）= 1250 / 1350 = 92.6 %80% 是 A100 上这个 kernel 的正常水平。剩下的 20% 来自：DRAM 刷新、读写转向（读 80 MB 写 40 MB 会有一部分总线转向开销）、行激活、ECC 以及 kernel 的头尾（launch 与 tail）。
（4）瓶颈判定与优化方向
瓶颈类型 = 内存带宽受限（AI 0.125 ≪ 12.54），且延迟隐藏充分（2 倍余量）。这意味着：
- 提高占用率 / 增加 ILP：无效，带宽已经接近打满。
- 减少指令数（比如用
float4向量化）：略微有效——它减少的是指令发射开销，不是字节数；在本 kernel 里能带来 1%–3% 的收益，主要价值是减少索引计算和提升每线程的在途字节数。 - 减少字节数：唯一真正有效的手段。例如
C = A + B之后再读 C 做别的运算是完全可以避免的；或者如果问题只需要读不需要写，流量可以从 12 B/元素降到 8 B/元素，理论上限时间从 77.2 µs 降到 51.5 µs。 - 用
float4（128 位）访存：这是本讲最值得记住的优化。让每个线程处理 4 个连续 float（128 位对齐），则一个 warp 依然是一次 128 字节事务，但每线程的指令数降到 1/4，在途字节数变成 4 倍（每 warp 每个LDG.128就是 512 字节），可以显著降低对占用率的要求。这就是示例 2 的思路。
示例 2：网格-步长循环 + 线程粗化（vecAddCoarsened）
// 文件: vecadd_gridstride.cu
// 编译: nvcc -O3 -arch=sm_80 vecadd_gridstride.cu -o vecadd_gridstride
// 运行: ./vecadd_gridstride 10000000
//
// 内容: 三个版本横向对比
// (1) vecAddBasicKernel 每线程 1 个元素，网格 = ceil(N/blockDim)
// (2) vecAddGridStrideKernel 网格-步长循环，网格由设备规模决定
// (3) vecAddCoarsenedKernel 网格-步长 + 每线程 UNROLL 个元素（粗化）
// 并打印每个版本的耗时、带宽、寄存器数、驻留 block 数与占用率。
#include <cstdio>
#include <cstdlib>
#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 UNROLL 4
// ---------------------------------------------- (1) 基线：每线程一个元素 ----
__global__ void vecAddBasicKernel(const float* __restrict__ A,
const float* __restrict__ B,
float* __restrict__ C,
int N)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < N) {
C[i] = A[i] + B[i];
}
}
// ------------------------------------- (2) 网格-步长循环（每轮一个元素） ----
// 步长 = 网格总线程数；每轮整个网格覆盖一段连续内存，合并性完好。
__global__ void vecAddGridStrideKernel(const float* __restrict__ A,
const float* __restrict__ B,
float* __restrict__ C,
int N)
{
int stride = blockDim.x * gridDim.x;
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < N; i += stride) {
C[i] = A[i] + B[i];
}
}
// ------------------------- (3) 网格-步长 + 线程粗化（每线程 UNROLL 个） ----
// 同一轮里的 UNROLL 个访存互相独立（地址相隔 stride），
// 4 条 LDG 同时在途 -> 每 warp 512 B / 每数组，制造 memory-level parallelism。
// 注意：间隔取 stride 而不是 blockDim.x，才能保持每个访问各自合并。
__global__ void vecAddCoarsenedKernel(const float* __restrict__ A,
const float* __restrict__ B,
float* __restrict__ C,
int N)
{
const int stride = blockDim.x * gridDim.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;
// 主循环：保证 UNROLL 个访问全部在界内（用 long long 避免 int 溢出）
for (; (long long)i + (long long)(UNROLL - 1) * stride < (long long)N;
i += UNROLL * stride) {
float a0 = A[i];
float a1 = A[i + (long long)1 * stride];
float a2 = A[i + (long long)2 * stride];
float a3 = A[i + (long long)3 * stride];
float b0 = B[i];
float b1 = B[i + (long long)1 * stride];
float b2 = B[i + (long long)2 * stride];
float b3 = B[i + (long long)3 * stride];
C[i] = a0 + b0;
C[i + (long long)1 * stride] = a1 + b1;
C[i + (long long)2 * stride] = a2 + b2;
C[i + (long long)3 * stride] = a3 + b3;
}
// 尾部：剩余元素按 stride 逐个处理
for (; i < N; i += stride) {
C[i] = A[i] + B[i];
}
}
// --------------------------------------------------------------- 工具函数 ---
#define TIME_LAUNCH(kernelName, grid, block, iters, outMs) \
do { \
kernelName<<<(grid), (block)>>>(d_A, d_B, d_C, N); \
CUDA_CHECK(cudaGetLastError()); \
CUDA_CHECK(cudaDeviceSynchronize()); \
cudaEvent_t s__, e__; \
CUDA_CHECK(cudaEventCreate(&s__)); \
CUDA_CHECK(cudaEventCreate(&e__)); \
CUDA_CHECK(cudaEventRecord(s__)); \
for (int it__ = 0; it__ < (iters); ++it__) { \
kernelName<<<(grid), (block)>>>(d_A, d_B, d_C, N); \
} \
CUDA_CHECK(cudaEventRecord(e__)); \
CUDA_CHECK(cudaEventSynchronize(e__)); \
CUDA_CHECK(cudaGetLastError()); \
float ms__ = 0.0f; \
CUDA_CHECK(cudaEventElapsedTime(&ms__, s__, e__)); \
(outMs) = (double)ms__ / (iters); \
CUDA_CHECK(cudaEventDestroy(s__)); \
CUDA_CHECK(cudaEventDestroy(e__)); \
} while (0)
static void initVector(float* v, int N, unsigned int seed)
{
unsigned int x = seed;
for (int i = 0; i < N; ++i) {
x = x * 1664525u + 1013904223u;
v[i] = (float)(x >> 8) / (float)(1u << 24) - 0.5f;
}
}
static void vecAddCPU(const float* A, const float* B, float* C, int N)
{
for (int i = 0; i < N; ++i) { C[i] = A[i] + B[i]; }
}
static int verify(const float* h_got, const float* h_ref, int N, const char* name)
{
int bad = 0;
for (int i = 0; i < N; ++i) {
if (h_got[i] != h_ref[i]) {
if (bad < 3) {
printf(" [%s] 首个不匹配 i=%d: got %.9g, want %.9g\n",
name, i, h_got[i], h_ref[i]);
}
++bad;
}
}
printf(" %-26s %s (不匹配 %d / %d)\n", name,
bad == 0 ? "PASS" : "FAIL", bad, N);
return bad;
}
int main(int argc, char** argv)
{
int N = (argc > 1) ? atoi(argv[1]) : 10000000;
if (N <= 0) { fprintf(stderr, "N must be positive\n"); return EXIT_FAILURE; }
size_t bytes = (size_t)N * sizeof(float);
float* h_A = (float*)malloc(bytes);
float* h_B = (float*)malloc(bytes);
float* h_C = (float*)malloc(bytes);
float* h_Ref = (float*)malloc(bytes);
if (!h_A || !h_B || !h_C || !h_Ref) {
fprintf(stderr, "host malloc failed\n");
return EXIT_FAILURE;
}
initVector(h_A, N, 12345u);
initVector(h_B, N, 67890u);
float *d_A = NULL, *d_B = NULL, *d_C = NULL;
CUDA_CHECK(cudaMalloc((void**)&d_A, bytes));
CUDA_CHECK(cudaMalloc((void**)&d_B, bytes));
CUDA_CHECK(cudaMalloc((void**)&d_C, bytes));
CUDA_CHECK(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice));
// ---- 让 CUDA 用占用率模型推荐一个 block 大小 -------------------------
int suggestedBlock = 0, minGridSize = 0;
CUDA_CHECK(cudaOccupancyMaxPotentialBlockSize(
&minGridSize, &suggestedBlock, vecAddCoarsenedKernel, 0, 0));
int blockSize = (suggestedBlock >= 128 && suggestedBlock <= 512)
? suggestedBlock : 256;
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
int numSMs = prop.multiProcessorCount;
double peakGBs = 2.0 * (double)prop.memoryClockRate * 1.0e3
* ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
// ---- 三个版本的网格配置 ---------------------------------------------
int gridBasic = (int)(((long long)N + blockSize - 1) / blockSize);
int blocksPerSM = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&blocksPerSM, vecAddGridStrideKernel, blockSize, 0));
int gridStride = numSMs * blocksPerSM;
int needed = (int)(((long long)N + blockSize - 1) / blockSize);
if (gridStride > needed) { gridStride = needed; }
int blocksPerSMC = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&blocksPerSMC, vecAddCoarsenedKernel, blockSize, 0));
int gridCoarse = numSMs * blocksPerSMC;
if (gridCoarse > needed) { gridCoarse = needed; }
printf("================= 设备与配置 =================\n");
printf("GPU : %s (sm_%d%d), %d 个 SM, 峰值带宽 %.1f GB/s\n",
prop.name, prop.major, prop.minor, numSMs, peakGBs);
printf("推荐 block 大小 : %d (cudaOccupancyMaxPotentialBlockSize)\n",
suggestedBlock);
printf("实际使用 block : %d\n", blockSize);
printf("N : %d (%.1f MB/数组)\n", N, bytes / 1048576.0);
printf("版本 1 网格 : %d 个 block (每线程 1 个元素)\n", gridBasic);
printf("版本 2 网格 : %d 个 block (网格-步长, 每线程多轮)\n", gridStride);
printf("版本 3 网格 : %d 个 block (网格-步长 + UNROLL=%d)\n",
gridCoarse, UNROLL);
const int iters = 20;
double msBasic = 0.0, msStride = 0.0, msCoarse = 0.0;
double totalBytes = 3.0 * (double)bytes;
printf("================= 正确性验证 =================\n");
vecAddCPU(h_A, h_B, h_Ref, N);
vecAddBasicKernel<<<gridBasic, blockSize>>>(d_A, d_B, d_C, N);
CUDA_CHECK(cudaGetLastError()); CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));
int bad1 = verify(h_C, h_Ref, N, "vecAddBasicKernel");
vecAddGridStrideKernel<<<gridStride, blockSize>>>(d_A, d_B, d_C, N);
CUDA_CHECK(cudaGetLastError()); CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));
int bad2 = verify(h_C, h_Ref, N, "vecAddGridStrideKernel");
vecAddCoarsenedKernel<<<gridCoarse, blockSize>>>(d_A, d_B, d_C, N);
CUDA_CHECK(cudaGetLastError()); CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));
int bad3 = verify(h_C, h_Ref, N, "vecAddCoarsenedKernel");
printf("================= 性能测量 =================\n");
TIME_LAUNCH(vecAddBasicKernel, gridBasic, blockSize, iters, msBasic);
TIME_LAUNCH(vecAddGridStrideKernel, gridStride, blockSize, iters, msStride);
TIME_LAUNCH(vecAddCoarsenedKernel, gridCoarse, blockSize, iters, msCoarse);
cudaFuncAttributes a1, a2, a3;
CUDA_CHECK(cudaFuncGetAttributes(&a1, vecAddBasicKernel));
CUDA_CHECK(cudaFuncGetAttributes(&a2, vecAddGridStrideKernel));
CUDA_CHECK(cudaFuncGetAttributes(&a3, vecAddCoarsenedKernel));
int maxWarpsPerSM = prop.maxThreadsPerMultiProcessor / 32;
int w1 = blocksPerSM * (blockSize / 32);
int w3 = blocksPerSMC * (blockSize / 32);
printf(" %-26s %10s %10s %9s %9s %9s\n",
"版本", "耗时(us)", "带宽(GB/s)", "峰值占比", "寄存器", "占用率");
printf(" %-26s %10.2f %10.1f %8.1f%% %9d %8.1f%%\n",
"1 基础(1 元素/线程)", msBasic * 1000.0,
totalBytes / (msBasic * 1.0e-3) / 1e9,
100.0 * (totalBytes / (msBasic * 1.0e-3) / 1e9) / peakGBs,
a1.numRegs, 100.0 * w1 / maxWarpsPerSM);
printf(" %-26s %10.2f %10.1f %8.1f%% %9d %8.1f%%\n",
"2 网格-步长(1 元素)", msStride * 1000.0,
totalBytes / (msStride * 1.0e-3) / 1e9,
100.0 * (totalBytes / (msStride * 1.0e-3) / 1e9) / peakGBs,
a2.numRegs, 100.0 * w1 / maxWarpsPerSM);
printf(" %-26s %10.2f %10.1f %8.1f%% %9d %8.1f%%\n",
"3 网格-步长+粗化(4 元素)", msCoarse * 1000.0,
totalBytes / (msCoarse * 1.0e-3) / 1e9,
100.0 * (totalBytes / (msCoarse * 1.0e-3) / 1e9) / peakGBs,
a3.numRegs, 100.0 * w3 / maxWarpsPerSM);
printf(" (总流量 = %.1f MB = 读 A + 读 B + 写 C)\n", totalBytes / 1e6);
CUDA_CHECK(cudaFree(d_A));
CUDA_CHECK(cudaFree(d_B));
CUDA_CHECK(cudaFree(d_C));
free(h_A); free(h_B); free(h_C); free(h_Ref);
CUDA_CHECK(cudaDeviceReset());
int bad = bad1 + bad2 + bad3;
printf("总体验证结果: %s\n", bad == 0 ? "ALL PASS" : "FAIL");
return (bad == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}
【代码做什么？】
- 三个 kernel 并存，共用同一套主机端框架，因此它们的耗时、带宽、寄存器、占用率可以在同一台机器上直接横向比较——这是做性能实验的正确方式：控制变量（相同 N、相同 block 大小、相同数据、相同迭代次数、相同计时方式）。
- 版本 1 与示例 1 相同：
dimGrid = ceil(N/256) = 39063，每线程一个元素，多出的 128 个线程被if (i < N)拦住。 - 版本 2 的循环：
for (int i = tid; i < N; i += stride)，其中stride = blockDim.x * gridDim.x。网格大小改为numSMs × blocksPerSM（A100 上是108 × 8 = 864），每个线程平均执行10,000,000 / (864 × 256) = 45.2轮。循环条件i < N同时承担了边界检查的职责——网格-步长循环天然是边界安全的，因为它不会启动越界线程，只在处理越界索引时退出。 - 版本 3 在主循环里一次处理
UNROLL = 4个元素，地址分别相隔1×stride、2×stride、3×stride。关键细节：间隔取stride而不是blockDim.x。若用blockDim.x作为间隔，那么同一线程的 4 个元素在内存里相邻，但相邻线程的访问会互相错开，覆盖率会出现漏洞（例如 grid=2、block=2、UNROLL=2 时会漏掉索引 6 和 7），同时”一个 block 的 4 个块”也会与相邻 block 的区间重叠。用stride作间隔则保证：每一轮里每个访问各自都是 warp 级连续的，而且全体线程的访问集合恰好覆盖[0, N)一次。 - 尾部处理：主循环的守卫条件
i + (UNROLL-1)*stride < N保证 4 个地址全部在界内；一旦不成立就退出主循环，由后面的for (; i < N; i += stride)用单元素方式收尾。这两段循环的覆盖集合是不重叠的：主循环覆盖的是”整组的 UNROLL 个 stride 间隔”，尾部循环从同一个 i 出发继续以 stride 步进，正好接上。 - 网格大小由设备决定：
cudaOccupancyMaxActiveBlocksPerMultiprocessor返回”每个 SM 能同时驻留几个这个 kernel 的 block”，乘以prop.multiProcessorCount就是”一次性能装满整个 GPU”的网格大小。cudaOccupancyMaxPotentialBlockSize则反向推荐一个能最大化占用率的 block 大小（本例中若推荐值不在 128–512 的合理区间就回退到 256，因为 1024 这类大 block 会带来尾部浪费和阶梯效应）。 __restrict__的作用：三个指针都标注了__restrict__，向编译器承诺”这三块内存不会互相重叠”。这让编译器可以放心地把A[i]、B[i]的加载提前并重排（不需要担心C[i]的写入会改变A[i+1]的值），是让 UNROLL 的四条LDG真正同时在途的必要条件。没有__restrict__时，编译器必须保守地假设写 C 可能影响后面的读 A，从而把访存串行化。verify()用逐位比较：每个输出元素都是一次独立的 FP32 加法，没有累加顺序问题，所以三个版本的结果都必须与 CPU==相等。任何不匹配都说明映射或边界处理有 bug。
【并行机制与硬件映射解说】
（a）三个版本的 warp 组织完全相同：block 大小都是 256 → 每 block 8 个 warp；warp 内线程仍然按
threadIdx.x连续划分，所以每个单独的访存指令依然是 32 个连续地址 = 1 次 128 字节事务。粗化并没有破坏合并性——它只是让每个 warp 每轮发出 4 条这样的LDG，而不是 1 条。（b）粗化如何提高访存并行度（MLP）：版本 3 的 4 条
LDG之间没有任何数据依赖（4 个地址由i和常量偏移算出，互不相关），所以硬件可以在等待第一条数据返回之前就把 4 条都发出去。每条LDG携带 128 字节（一个 warp），于是每个 warp 每轮有 512 字节在途。整个 warp 能同时容纳的LDG数量受”每 SM 的 LSU（load/store unit）队列深度”限制（Ampere 上每 SM 在途访存请求上限是数百个），4 条远未触顶。（c）占用率与资源：A100 每 SM 2048 线程 / 64 warp / 65536 寄存器 / 32 个 block 上限。
- 版本 1、2：约 16 寄存器/线程 →
65536/(16×256) = 16个 block（寄存器不限制），线程限制2048/256 = 8个 block → 8 block = 64 warp = 100%。 - 版本 3：8 个活跃 float（
a0..a3、b0..b3）加索引，实测约 24 寄存器/线程 →65536/(24×256) = 10个 block（寄存器不限制），线程限制 8 → 8 block = 64 warp = 100%。 三个版本的占用率都是 100%，所以本次对比不是在做”占用率 vs 带宽”的实验，而是在验证”在占用率相同的前提下，粗化是否还能带来收益”。答案是：在这个 N 和这台 GPU 上，三者带宽基本相同（都接近 DRAM 上限），说明内存系统已经完全饱和。
（d）共享内存 / bank conflict：三个 kernel 都不使用共享内存，无 bank conflict。
（e）寄存器与 warp 发散：版本 3 的尾部循环
for (; i < N; i += stride)在最后一个 warp 里可能出现部分线程退出、部分线程继续的情况（因为 N 不一定被 stride 整除），这会造成轻微的 warp 发散。但由于它只在最后一轮发生一次，代价可以忽略。真正需要注意的是：粗化把”边界处理”从”最多一个 warp 发散”扩大到了”最多 UNROLL 个 warp 发散”，仍然是一次性开销。- 版本 1、2：约 16 寄存器/线程 →
【性能优化分析】
（1）实测对比（A100，N = 10,000,000，block = 256，20 次迭代取平均）
┌──────────────────────────┬──────────┬────────────┬─────────┬────────┬─────────┐ │ 版本 │ 耗时(us) │ 带宽(GB/s) │ 峰值占比│ 寄存器 │ 占用率 │ ├──────────────────────────┼──────────┼────────────┼─────────┼────────┼─────────┤ │ 1 基础(1 元素/线程) │ 96.5 │ 1244 │ 80.0% │ 16 │ 100.0% │ │ 2 网格-步长(1 元素) │ 96.0 │ 1250 │ 80.4% │ 16 │ 100.0% │ │ 3 网格-步长+粗化(4 元素) │ 95.2 │ 1261 │ 81.1% │ 24 │ 100.0% │ └──────────────────────────┴──────────┴────────────┴─────────┴────────┴─────────┘ 差距在 1.4% 以内 —— 三者都已打满显存带宽，剩下的 19% 是 DRAM 物理开销。这个结果本身就是本讲最重要的实验结论之一：当一个 kernel 已经是纯带宽受限时，任何”提高并行度”的优化都不会有可观测的收益。你能做的只有减少字节数。
（2）粗化的真正价值：用 ILP 换占用率
粗化的收益在寄存器压力大、占用率被迫降低的 kernel 上才会显现。构造一个具体场景：假设 kernel 每线程用 80 个寄存器（例如每线程要维护更多中间量）。
版本 2（UNROLL = 1）在 80 寄存器、block = 256 时： 每 SM 可驻留 block 数 = floor(65536 / (80 x 256)) = floor(3.2) = 3 驻留线程 = 3 x 256 = 768 = 24 warp 占用率 = 24 / 64 = 37.5% 每 SM 在途字节 = 24 warp x 1 个访问 x 128 B = 3,072 B 打满带宽需要每 SM 约 4,104 B（Little's Law）： 可得带宽 ≈ 1250 x (3072 / 4104) = 936 GB/s -> 只有可达带宽的 75% 版本 3（UNROLL = 4）在同样的 80+8 寄存器假设下： 每 SM 在途字节 = 24 warp x 4 个访问 x 128 B = 12,288 B（需要的 3 倍） -> 带宽重新回到 1250 GB/s 附近这就是”低占用率 + 高 ILP ≡ 高占用率 + 低 ILP“的定量版本，也是现代 GPU kernel（FlashAttention 一类）普遍采用”少量线程、大量寄存器、强 ILP”风格的原因。
（3）算术强度与 Roofline（三版本相同）
AI = 1 FLOP / 8 B = 0.125，或含写为 1/12 = 0.0833 FLOP/Byte A100 平衡点 = 12.54 FLOP/Byte -> 带宽受限 Roofline 天花板 = 0.0833 x 1555 = 129.6 GFLOP/s 版本 3 实测算力 = 0.0833 x 1261 = 105.0 GFLOP/s = 天花板的 81.0% 向量加法的时间下限（不可突破）：120 MB / 1555 GB/s = 77.2 us（4）瓶颈判定与优化方向
三个版本都是带宽受限 + 延迟隐藏充分。因此优化方向只剩：
- 减少字节数（唯一有大收益的方向）：把
float换成更窄的类型（如果精度允许）、避免不必要的中间数组、融合多个 pass。 - 使用向量化访存
float4：每线程处理 4 个相邻float（而不是相隔 stride 的 4 个），一条LDG.128就搬运 16 字节/线程 = 512 字节/warp（4 次事务）。它把索引计算指令数降到 1/4，在算术强度极低的 kernel 上通常能带来 1%–3% 的提升（本 kernel 收益很小，因为指令发射不是瓶颈）。 - 不做：提高占用率（已 100%）、增加 unroll（已饱和）、使用共享内存（无复用，纯增开销）。
- 减少字节数（唯一有大收益的方向）：把
示例 3：二维网格矩阵加法 matrixAdd（二维线程映射）
// 文件: matrixadd_2d.cu
// 编译: nvcc -O3 -arch=sm_80 matrixadd_2d.cu -o matrixadd_2d
// 运行: ./matrixadd_2d 4096 4096
//
// 内容: 用二维 grid / 二维 block 处理二维数据（C = A + B，矩阵加法）。
// 展示 row-major 线性化、二维边界检查、以及"哪个维度该当 threadIdx.x"。
#include <cstdio>
#include <cstdlib>
#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 TILE 32 // block 形状 32 x 32 = 1024 线程
// 二维线程映射：threadIdx.x 对应"列"（row-major 中最快变化的维度），
// 这样同一个 warp 的 32 个线程访问同一行里连续的 32 个 float = 128 字节。
__global__ void matrixAddKernel(const float* __restrict__ A,
const float* __restrict__ B,
float* __restrict__ C,
int width, int height)
{
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (col < width && row < height) { // 两个维度都要检查
int idx = row * width + col; // row-major 线性化
C[idx] = A[idx] + B[idx];
}
}
static void initMatrix(float* m, long n, unsigned int seed)
{
unsigned int x = seed;
for (long i = 0; i < n; ++i) {
x = x * 1664525u + 1013904223u;
m[i] = (float)(x >> 8) / (float)(1u << 24) - 0.5f;
}
}
int main(int argc, char** argv)
{
int width = (argc > 1) ? atoi(argv[1]) : 4096;
int height = (argc > 2) ? atoi(argv[2]) : 4096;
if (width <= 0 || height <= 0) {
fprintf(stderr, "width/height must be positive\n");
return EXIT_FAILURE;
}
long N = (long)width * (long)height;
size_t bytes = (size_t)N * sizeof(float);
float* h_A = (float*)malloc(bytes);
float* h_B = (float*)malloc(bytes);
float* h_C = (float*)malloc(bytes);
float* h_Ref = (float*)malloc(bytes);
if (!h_A || !h_B || !h_C || !h_Ref) {
fprintf(stderr, "host malloc failed\n");
return EXIT_FAILURE;
}
initMatrix(h_A, N, 12345u);
initMatrix(h_B, N, 67890u);
float *d_A = NULL, *d_B = NULL, *d_C = NULL;
CUDA_CHECK(cudaMalloc((void**)&d_A, bytes));
CUDA_CHECK(cudaMalloc((void**)&d_B, bytes));
CUDA_CHECK(cudaMalloc((void**)&d_C, bytes));
CUDA_CHECK(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice));
// 二维网格：每个维度都要向上取整
dim3 dimBlock(TILE, TILE, 1);
dim3 dimGrid((unsigned)((width + TILE - 1) / TILE),
(unsigned)((height + TILE - 1) / TILE), 1);
matrixAddKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, width, height);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
const int iters = 20;
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(cudaEventRecord(start));
for (int it = 0; it < iters; ++it) {
matrixAddKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, width, height);
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float msTotal = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&msTotal, start, stop));
double msPerIter = (double)msTotal / iters;
CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));
for (long i = 0; i < N; ++i) { h_Ref[i] = h_A[i] + h_B[i]; }
long mismatch = 0;
for (long i = 0; i < N; ++i) { if (h_C[i] != h_Ref[i]) { ++mismatch; } }
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
double peakGBs = 2.0 * (double)prop.memoryClockRate * 1.0e3
* ((double)prop.memoryBusWidth / 8.0) / 1.0e9;
int blocksPerSM = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&blocksPerSM, matrixAddKernel, TILE * TILE, 0));
double totalBytes = 3.0 * (double)bytes;
double achievedGBs = totalBytes / (msPerIter * 1.0e-3) / 1.0e9;
printf("================= 矩阵加法 %d x %d =================\n", width, height);
printf("GPU : %s, 峰值带宽 %.1f GB/s\n", prop.name, peakGBs);
printf("blockDim : (%d, %d), dimGrid = (%u, %u)\n",
TILE, TILE, dimGrid.x, dimGrid.y);
printf("总线程数 : %u (每个线程一个元素)\n",
dimGrid.x * dimGrid.y * TILE * TILE);
printf("正确性验证 : %s (不匹配 %ld / %ld)\n",
mismatch == 0 ? "PASS" : "FAIL", mismatch, N);
printf("kernel 单次耗时 : %.3f us\n", msPerIter * 1000.0);
printf("有效带宽 : %.1f GB/s (%.1f%% of %.1f GB/s 峰值)\n",
achievedGBs, 100.0 * achievedGBs / peakGBs, peakGBs);
printf("每 SM 驻留 block 数 : %d (%d 线程/SM)\n",
blocksPerSM, blocksPerSM * TILE * TILE);
printf("占用率 : %.1f %%\n",
100.0 * blocksPerSM * TILE * TILE / prop.maxThreadsPerMultiProcessor);
printf("每 warp 事务数 : 3 (A/B/C 各 1 次 128 B 事务)\n");
CUDA_CHECK(cudaFree(d_A));
CUDA_CHECK(cudaFree(d_B));
CUDA_CHECK(cudaFree(d_C));
free(h_A); free(h_B); free(h_C); free(h_Ref);
CUDA_CHECK(cudaDeviceReset());
return (mismatch == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}
【代码做什么？】
- 二维映射：
col = blockIdx.x * blockDim.x + threadIdx.x、row = blockIdx.y * blockDim.y + threadIdx.y。这是把一维公式分别应用到两个维度上，两个维度互相独立。 - row-major 线性化：C/C++ 的二维数组在内存里按行连续存放，元素
(row, col)的线性地址是row * width + col（讲义中的例子：M[2][1]在 4 列矩阵里编号为2*4+1 = 9）。因此”行”是慢变维度、”列”是快变维度。 - 哪个维度当
threadIdx.x：必须让快变（列）维度对应threadIdx.x，因为硬件按threadIdx.x最快线性化来划分 warp。这样同一个 warp 的 32 个线程访问的是同一行里连续的 32 个 float = 128 连续字节 = 1 次事务。如果反过来（col对应threadIdx.y），warp 内 32 个线程会跨越 32 行，每行取 4 字节，步长width*4 = 16384字节，事务数从 1 变成 32，带宽利用率掉到 3.1%。 - 二维边界检查：
if (col < width && row < height)。维度必须分别检查——只有一维需要检查时另一维可以省，但两个维度都写了更安全。讲义中的例子（76×62 图像用 16×16 的 block）需要ceil(76/16)=5列 ×ceil(62/16)=4行 = 20 个 block，覆盖 80×64，多出的区域全靠这两个条件拦住。 - 网格配置：
dimGrid = (ceil(width/TILE), ceil(height/TILE), 1)，两个维度都向上取整；dimBlock = (TILE, TILE, 1)，本程序 TILE = 32，即每个 block 1024 线程。 - 主机端流程与示例 1 相同（分配 → 拷贝 → 预热 → 计时 → 取回 → 用 CPU 参考实现逐位验证 → 打印带宽与占用率 → 释放）。
- 二维映射：
【并行机制与硬件映射解说】
（a）warp 划分：block 形状 32×32，
threadIdx的线性化顺序是 x 最快、y 次之、z 最后。所以 warp 0 =threadIdx.y = 0, threadIdx.x = 0..31（第 0 行），warp 1 =threadIdx.y = 1, threadIdx.x = 0..31（第 1 行），依此类推，一个 block 恰好 32 个 warp，每个 warp 对应一整行。（b）A100 上的驻留量：每 SM 2048 线程 / 64 warp。block 大小 1024 →
2048/1024 = 2个 block（64 warp）→ 100% 占用率；寄存器约 16 个/线程 →65536/(16×1024) = 4个 block，不构成限制。但在 RTX 4090（1536 线程/SM）上，1536/1024 = 1个 block → 1024 线程 = 32 warp → 只有 66.7% 占用率。这是”大 block 的阶梯效应”的典型例子：同样是 1024 线程的 block，在不同 GPU 上占用率相差 33%。稳健的选择是 16×16 = 256 线程的 block，它在 1536、2048 两种 SM 配置下都能整除并达到 100%。（c）全局内存访问是否合并：以矩阵
4096×4096为例，width = 4096。block (0,0) 的 warp 0 处理row = 0, col = 0..31，线性索引idx = 0*4096 + 0..31 = 0..31，字节地址 0..127。基地址 256 字节对齐 → 1 次 128 字节事务，利用率 100%。warp 1 处理row = 1, col = 0..31，idx = 4096..4127，字节地址 16384..16511，同样是 1 次对齐的事务。每个 warp 3 次事务（A、B、C 各一次），与一维向量加法完全一致——二维映射本身不影响合并性，影响合并性的是”哪个维度最快变化”。（d）共享内存 / bank conflict：不使用共享内存，无 bank conflict。矩阵加法同样没有数据复用（每个元素只被读一次），共享内存是纯开销。
（e）寄存器与发散：约 16 个寄存器/线程（两个坐标、一个 width、一个 idx、三个数据）。发散很轻微：只有落在
col >= width或row >= height区域的 warp 会发散。4096 是 32 的整数倍，所以本次运行完全没有发散；若换成4096×4094，则每行的最后一个 warp 会有 30 个线程空转（因为 y 方向有 2 行是越界的，整个 warp 都不活跃）。用 16×16 的 block 时，越界区域是”半个 warp 的粒度”（因为一个 warp 只覆盖一行 16 个线程的两行），发散开销比 32×32 更细。【性能优化分析】
（1）算术强度与 Roofline
每个元素：1 FLOP（一次加法），访存 读 4 B + 读 4 B + 写 4 B = 12 B AI = 1 / 12 = 0.0833 FLOP/Byte（含写）；只算读时 = 1/8 = 0.125 FLOP/Byte A100 平衡点 = 19.5e12 / 1555e9 = 12.54 FLOP/Byte AI 比平衡点低 100 倍以上 -> 纯带宽受限，与向量加法结论完全一致 数据规模：4096 x 4096 x 4 B = 64 MiB/矩阵，总流量 192 MiB = 201.3 MB 时间下限 = 201.3e6 / 1555e9 = 129.5 us 实测（按 1250 GB/s）= 201.3e6 / 1250e9 = 161.0 us 对应算力 = 0.0833 x 1250e9 = 104 GFLOP/s = FP32 峰值的 0.53%注意：把矩阵加法当成一个”更大的向量加法”是完全正确的——
row * width + col只是一个仿射变换，访问模式依然是纯顺序的。矩阵加法比向量加法多出来的唯一东西是”索引计算的指令数”和”两个维度各自的边界检查”，这些都属于指令开销，不是带宽开销。（2）占用率
block = 32 x 32 = 1024 线程 A100 : floor(2048/1024) = 2 block/SM -> 2048 线程 = 64 warp -> 100.0% RTX4090: floor(1536/1024) = 1 block/SM -> 1024 线程 = 32 warp -> 66.7% 改用 16 x 16 = 256 线程： A100 : floor(2048/256) = 8 block/SM -> 2048 线程 -> 100.0% RTX4090: floor(1536/256) = 6 block/SM -> 1536 线程 -> 100.0%这就是讲义中那个经典分析（”1536 线程 / 最多 8 block 的 SM 上，8×8、16×16、32×32 哪个最好”）的推广：16×16 是跨硬件最稳健的二维 block 形状，8×8 会被”每 SM block 数上限”卡住（例如上限 8 时只能放 8×64 = 512 线程 = 25% 的可用线程容量），32×32 会被”每 SM 线程数上限”卡住。
（3）瓶颈判定与优化方向 瓶颈 = 带宽受限（AI 0.0833 ≪ 12.54），延迟隐藏充分（A100 上 100% 占用率，或在 4090 上用 16×16 时 100%）。
每 SM 在途字节（A100, 32x32 block, 2 block/SM = 64 warp）： 64 warp x 1 个在途访问 x 128 B = 8,192 B >= 需要的 4,104 B（2 倍余量） 这是否足够？ 是。所以矩阵加法在 A100 上同样能达到约 80% 的峰值带宽。优化方向：
- 二维 kernel 里把
width、height传成int是够的（4096 不会溢出），但如果规模可能超过 2^31 个元素，row * width + col必须用long long。这是二维 kernel 最常见的溢出点。 - 如果不需要二维语义（比如整个矩阵统一加一个常数），直接用一维映射：更少的索引指令、更简单的边界检查，实测通常快 1%–2%。
- 只有在需要邻域访问（模糊、卷积、stencil）时，二维 block 的形状才真正重要——那时需要共享内存来做 halo 交换，block 形状会影响 halo 占比和 bank conflict。
- 二维 kernel 里把
Lab 1 的工程习惯与评分口径
（一）三条必须养成的工程习惯
① CPU 端 golden reference —— "没有参考实现的性能数字毫无意义"
- 在 CPU 上用最朴素的双重循环/单重循环算一遍，作为标准答案。
- 对欧拉式 kernel（每元素独立、无累加），结果应当逐位相同（bit-exact），
可以直接用 == 比较；对累加/reduction 类 kernel 才需要 epsilon。
- 在 GPU 上跑出"更快"的结果但没验证 = 得到一个错误的数字，Lab 直接 0 分。
- 建议把验证做成函数 verify()，在任何性能测量之前先跑一遍。
② 边界检查 —— "最后一个 block 的越界线程一定会捅娄子"
- 网格向上取整：dimGrid = (N + blockDim - 1) / blockDim（整数运算）。
- kernel 里：if (i < N) { C[i] = A[i] + B[i]; }，二维时两个维度分别检查。
- 测试用例要专门包含 N 不是 blockDim 整数倍的情况
（例如 N = 1000、N = 1025、N = 1），这些正是能暴露 bug 的输入。
③ CUDA 错误检查 —— "静默失败是 CUDA 的默认行为"
- 每一个 CUDA API 调用都用 CUDA_CHECK 宏包起来。
- kernel 启动后用 cudaGetLastError() 检查启动期错误（非法配置、越界等）。
- 计时循环之后必须 cudaDeviceSynchronize()（或 cudaEventSynchronize），
否则测到的是"启动 kernel 的时间"而不是"kernel 执行的时间"。
- 强烈建议在 debug 构建里跑一次 cuda-memcheck（或 compute-sanitizer），
它能直接指出越界访问的地址和线程号。
开发循环（推荐的调试顺序）
写 kernel
-> 用 N = 1、N = 31、N = 33、N = 1024、N = 1025 五个小规模跑通（验证边界）
-> 与 CPU 参考实现逐位对比，全部 PASS
-> 再放大到 N = 10^6 / 10^7 测性能
-> 用 cudaEvent 只圈住 kernel，迭代 >= 20 次取平均（消除抖动）
-> 先预热 1 次，丢弃第一次的测量
-> 报告"达到峰值带宽的百分比"，而不是裸的毫秒数
切忌：先调性能，后补正确性。这样做的学生最后往往两者都没有。
（二）评分口径（课程 Lab 的一般形式）
┌─────────────────────────────────────────────────────────────────────┐
│ 第一部分：正确性测试（通常占大头，不正确则性能分无效） │
│ - 用课程提供的 harness 把学生写的 .cu 与参考答案对比 │
│ - 多个规模：小规模（暴露边界 bug）、中等规模、大规模（暴露溢出） │
│ - 多个数据分布：随机值、全零、极端值（验证数值稳定性） │
│ - 容差：整数/单次加法类要求精确；归约类允许相对误差 1e-5 量级 │
│ - 常见的"部分分"设计：只处理 N 为 blockDim 整数倍的情况得部分分 │
├─────────────────────────────────────────────────────────────────────┤
│ 第二部分：性能测试（在指定 GPU 上测，通常是 NCSA Delta 的 A100） │
│ - 指标一般是有效带宽（GB/s）或相对基线 kernel 的加速比 │
│ - 常见门槛：达到理论峰值带宽的 70%-85%（向量加法类） │
│ - 计时只圈 kernel，不含 cudaMemcpy；多次迭代取平均 │
│ - 有的版本会给出不同 N 的曲线，考察"随规模变化的趋势"的解释 │
├─────────────────────────────────────────────────────────────────────┤
│ 第三部分：实验报告（不少版本要求提交） │
│ - 说明线程映射公式、网格大小计算、边界处理的正确性论证 │
│ - 给出带宽测量方法（如何排除传输、如何预热、迭代次数） │
│ - 用算术强度 / Roofline / 占用率解释"为什么这个 kernel 是带宽受限" │
│ - 讨论：为什么继续优化（更多 unroll、更大 block）收益很小 │
└─────────────────────────────────────────────────────────────────────┘
注：以上是课程 Lab 的一贯口径的概括，具体门槛以当学期 handout 为准。
性能优化技巧总结
- 用
i = blockIdx.x * blockDim.x + threadIdx.x做最简映射：它只需一条 IMAD 指令，且天然使相邻线程访问相邻地址——映射公式本身就是合并访存的前提。 - 网格大小一定向上取整：
(N + blockDim.x - 1) / blockDim.x。向下取整会漏掉尾部元素，是 Lab 1 最常见的丢分点。 - kernel 内必须做边界检查：越界写会静默破坏其它内存；二维问题两个维度都要检查。
- 让
threadIdx.x对应内存中变化最快的维度：warp 按threadIdx.x最快线性化，这样 32 个线程的访问才是 128 连续字节 = 1 次事务。 - 用
__restrict__标注指针：承诺”这些指针不重叠”，编译器才敢把多条访存提前并重排，粗化/UNROLL 才能真正形成在途并行。 - block 大小选 128/256/512，别默认用 1024：1024 会在 1536 线程/SM 的 GPU 上因阶梯效应只能放 1 个 block（占用率掉到 66.7%），且尾部浪费比例更大。
- 用
cudaOccupancyMaxActiveBlocksPerMultiprocessor实测占用率，不要凭寄存器数猜——寄存器、共享内存、block 数上限三个约束取最小值才是真相。 - 用”每 SM 在途字节 = 带宽 × 延迟 / SM 数”检查延迟隐藏是否足够（A100 约 4.1 KB/SM）：不够就提高占用率，或每线程多发几个独立访存（thread coarsening）。
- 粗化时让同线程的多个访问地址相隔
stride（网格总线程数），而不是相隔blockDim.x：这样每一个访问各自仍然合并，且覆盖率不会出现漏洞。 - 只圈 kernel 计时：
cudaEvent包住 kernel 循环，排除cudaMemcpy；预热一次；迭代 20 次以上取平均。 - 报告”峰值带宽百分比”而不是绝对毫秒：绝对时间随 GPU 型号变化，百分比才说明你的 kernel 有没有打满带宽。
- 减少字节数永远比提高并行度有效：当 kernel 已带宽受限时，任何占用率/ILP 优化都不会有可观测收益，唯一的出路是少搬数据（或提升数据复用——那是共享内存的领域）。
关键要点
- 线程映射公式
i = blockIdx.x * blockDim.x + threadIdx.x是数据并行的全部骨架：它把”第 b 个 block 的第 t 号线程”唯一地映射到一个数据元素，使硬件可以完全自主地调度 block（block 之间的执行顺序是任意的、不确定的，所以绝不能让 block 之间互相依赖）。 - 一个 warp 的 32 次 4 字节访问，如果地址连续且 128 字节对齐，会被合并成 1 次 128 字节事务（利用率 100%）；如果跨步到 32 个不同的 cache line，就是 32 次事务（利用率 3.1%，sector 粒度下也是 8 倍放大）。合并与否直接决定带宽的有效利用率，是最重要的性能分水岭。
- 向量加法的算术强度只有 0.125 FLOP/Byte（含写为 0.0833），而 A100 的机器平衡点是 12.54 FLOP/Byte，低了 100 倍：它必然、也只能是纯带宽受限。它的时间下限是
3 × N × 4 字节 / 峰值带宽，N = 10⁷ 时是 77.2 µs；实测 96 µs（1.25 TB/s，峰值的 80.4%）已经接近极限，FP32 流水线只有约 0.5% 的周期在工作。 - 占用率用”活跃 warp 数 / SM 最大 warp 数”计算，由线程数上限、寄存器上限、block 数上限三者取最小值决定：A100 上 2048 线程 / 64 warp / 65536 寄存器 / 32 block。轻量 kernel（16 寄存器）在 128–1024 的任何 block 大小下都能到 100%，但寄存器用量大的 kernel 会被寄存器卡住（80 寄存器 + 256 线程 → 只有 3 个 block → 37.5%）。
- 占用率的意义是”藏住 400–800 周期的全局内存延迟”：Little’s Law 给出量化标准——每 SM 需要约 4.1 KB 在途字节（A100）。64 warp × 1 个在途合并访问 = 8 KB（2 倍余量，够用）；16 warp × 1 = 2 KB（只够一半带宽）；而 16 warp × 4 个独立访存 = 8 KB（又够了）——这就是线程粗化用 ILP 换占用率的定量依据。
- 网格-步长循环把”网格大小”与”问题规模”解耦：网格由设备（SM 数 × 每 SM 可驻留 block 数）决定，每线程循环处理多个元素，因而同一份代码能自适应任意 GPU、任意 N、多 GPU 划分和运行时才知道的规模，同时消除波次量化（wave quantization）的尾部浪费和大量 block 的调度开销。
常见陷阱与注意事项
- 整数除法取整导致网格覆盖不足：
dimGrid = N / blockDim.x在 N 不是 blockDim.x 整数倍时漏算最后的元素（N = 10,000,000 时漏 128 个，正确性测试直接失败）→ 改用(N + blockDim.x - 1) / blockDim.x，并用long long计算分子避免N + B - 1溢出。 - 忘记边界检查（或只检查一个维度）导致越界访问：
if (i < N)漏写时，最后一个 block 的多余线程会读写自己数组之外的内存（可能写坏相邻分配甚至 CUDA 内部结构），表现为”时对时错”或延迟崩溃，用compute-sanitizer才能定位；二维 kernel 里只写if (col < width)而遗漏row < height，在高度不是 blockDim.y 整数倍时同样越界 → 网格向上取整与边界检查是一对，不能只写一半；二维问题的越界区域在行、列两个方向都存在，必须分别检查。附带一点：if (i < N)本身会在跨界的那一个 warp 里造成 warp 发散（branch divergence）——warp 内一部分线程走 then 分支、一部分走 else 分支，硬件必须串行执行两条路径——但由于它只在每个 kernel 的最后一轮发生一次，代价可以忽略；真正要避免的是”用threadIdx做条件去划分大量工作”（例如if (threadIdx.x > 2) { 大量代码 } else { 大量代码 }），正确做法是让分支粒度对齐 warp 大小，写成if (threadIdx.x / 32 > 2)，这样任何 warp 都只走一条路径。 - 忘记
cudaDeviceSynchronize()/ 未同步就测量：kernel 启动是异步的，在 kernel 后立刻读主机内存或用 CPU 时钟计时，测到的是启动开销而不是执行时间，还会读到未完成的结果 → 用cudaEvent计时并在cudaEventSynchronize之后读cudaEventElapsedTime；验证结果前也必须同步。 - 忘记检查 CUDA API 返回值与 kernel 启动错误：
cudaMalloc失败（显存不足）、<<<>>>配置非法（block 超过 1024、共享内存超出）都会返回错误码但不中断程序，后续代码在一片未初始化的内存上运行 → 所有 API 用CUDA_CHECK包裹；kernel 启动后立刻cudaGetLastError()。 - host/device 指针混用与
cudaMemcpy方向写反：把主机malloc得到的指针直接传给 kernel 会触发非法内存访问并报 “an illegal memory access was encountered”（一旦出现这个错误，整个 CUDA context 都会被污染，后续所有调用都失败）；把cudaMemcpyHostToDevice与cudaMemcpyDeviceToHost弄混则会把设备内存的垃圾拷到主机 → 设备端内存必须用cudaMalloc（或cudaMallocManaged），主机端不能解引用cudaMalloc返回的指针；记住cudaMemcpy的参数顺序是(dst, src, count, kind)，kind描述的是”数据从哪里来”。 - 跨步（strided）访问破坏合并性：为了”看起来更整齐”而写
i = threadIdx.x * gridDim.x + blockIdx.x（按列遍历），会让一个 warp 的 32 个线程访问相隔 156 KB 的地址，事务数从 1 涨到 32（sector 粒度下流量放大 8 倍），实测带宽掉到 1/8 → 始终让threadIdx.x（线性化最快的维度）对应连续内存；__restrict__与向量化都不能弥补这个问题。 - 在已饱和的带宽受限 kernel 上继续优化并行度：把 unroll 从 4 加到 8、把 block 从 256 改到 1024、加
__launch_bounds__提高占用率，在向量加法上都不会有可观测收益（三版本实测差距在 1.4% 以内）→ 先用算术强度和 Roofline 判断瓶颈类型，带宽受限时只做”减少字节数”和”保证合并”这两件事。 - 粗化时用
blockDim.x而不是stride作为线程内多元素的间隔：会让不同线程的覆盖区间互相重叠或出现空洞（测试 N 不是 stride 整数倍时结果错误且难以察觉）→ 同一个线程的多个元素地址必须相隔blockDim.x * gridDim.x，这样每个访问各自合并、整体恰好无重复无遗漏地覆盖[0, N)。
思考题（带答案）
Q1. 一个 1D 数组有 N = 1,000,000 个 float，用 blockDim = 256、每线程一个元素的向量加法处理。请给出网格大小、启动的线程总数、多出的空转线程数、一个 warp 的访存事务数、总事务数，以及在 A100 上打满峰值带宽所需的最短时间。
网格大小 = ceil(1,000,000 / 256) = 3907 个 block；启动线程总数 = 3907 × 256 = 1,000,192，多出 192 个空转线程（占 0.019%）。每个 warp（32 个连续 float）对 A、B、C 各有 1 次 128 字节事务，即 3 次/warp；有效 warp 数 = 1,000,000 / 32 = 31,250，总事务数 = 31,250 × 3 = 93,750 次，共搬运 93,750 × 128 B = 12,000,000 B = 12 MB。A100 峰值带宽 1555 GB/s，最短时间 = 12e6 / 1555e9 = 7.72 µs。注意这个规模太小，kernel 的 launch 开销（约 3–5 µs）已经与之同量级，这也是为什么性能实验要用 N = 10⁷ 以上。
Q2. 某 kernel 中线程 t 访问 A[(t % 32) * 4096 + (t / 32)]（float 数组，t 为 warp 内线程号 0..31）。请分析其访存效率；如果这是 Lab 1 的向量加法 kernel，实际带宽会变成多少（设 A100 可达带宽 1250 GB/s）？
warp 内相邻线程（t 和 t+1）的地址相差 4096 个 float = 16,384 字节。32 个线程访问的地址就是 4096 × t（t 取 0 到 31），即元素下标 0、4096、8192、12288 一直到 126976，也就是 32 个彼此相隔 16 KB 的地址。每个地址落在不同的 128 字节 cache line（也落在不同的 32 字节 sector），所以事务数 = 32 次/warp，用 sector 粒度算是 32 × 32 B = 1024 B 搬回 128 B 有用数据，流量放大 8 倍。实际带宽 ≈ 1250 / 8 = 156 GB/s，只有原来的 12.5%，if (i < N) 也拦不住这种错误（它是地址映射错误，不是越界）。
Q3. A100 的 FP32 峰值是 19.5 TFLOPS、显存带宽 1555 GB/s。有人提出”把向量加法改写成一个每元素做 16 次浮点运算的 kernel（例如同时做一次多项式求值）”，声称这样能”更好地利用 GPU 的算力”。请判断这个说法，并算出这个 kernel 在新算术强度下是带宽受限还是计算受限。
原来的 AI = 1 FLOP / 8 B = 0.125 FLOP/Byte。新 kernel 每元素 16 FLOP，流量仍是 8 B（读 A、B）或 12 B（含写 C）：AI = 16 / 12 = 1.333 FLOP/Byte。A100 的机器平衡点 = 19.5e12 / 1555e9 = 12.54 FLOP/Byte。1.333 < 12.54，所以它仍然是带宽受限的，Roofline 天花板 = 1.333 × 1555e9 = 2.07 TFLOP/s，仅占 FP32 峰值的 10.6%。要想让它在 A100 上变成计算受限，需要 AI ≥ 12.54，即每元素至少 12.54 × 12 ≈ 150 次浮点运算。这个例子说明：“多算一点”几乎不可能让一个访存密集的 kernel 变成计算密集的——要跨越 100 倍的算术强度鸿沟，必须改变数据复用方式（分块、共享内存），而不是增加每元素的运算次数。