Lecture 6: GPU Architecture and CUDA Programming (continued)

目录 · ← l5 · l7 →

Lecture 6: GPU Architecture and CUDA Programming (continued)

1. 章节标题与概述

Lecture 6: GPU Architecture and CUDA Programming（continued）（GPU 体系结构与 CUDA 编程（续））

• 本讲核心问题：上一讲回答了”CUDA 的抽象是什么、GPU 硬件大致长什么样”；本讲要回答的是“为什么同一个 CUDA 程序能差十倍，以及怎样把它写快”。讲义把答案压缩成三条可以直接动手的优化判据：coherent warps（warp 内不发散）、coalesced memory access（合并访存）、高数据复用（shared memory tiling + register tiling），并用两个完整 case study——矩阵乘（Case Study 1）与并行归约（Case Study 2）——把这三条判据从”口号”推到”每一行代码、每一个数字”。

• 涉及的主要硬件/软件机制：
  • 硬件侧：SM（Streaming Multiprocessor）内部的 sub-core / warp selector / fetch-decode / SIMD 功能单元（fp32、int、fp64、LSU、Tensor Core）与寄存器堆；片上 shared memory + L1；L2 与 HBM（High Bandwidth Memory，高带宽显存，讲义标注 ~1 TB/s 量级）；GPU work scheduler（线程块调度器，按资源需求动态把 block 映射到 SM）；32 个 bank 的共享内存子系统；以及 Tensor Core 带来的算力量级跃变（GTX 980 的 4.6 TFLOPs → H100 的 1000 TFLOPs）。
  • 软件侧：__shared__（含 extern __shared__ 动态共享内存）、__syncthreads()、atomicAdd（可作用于 global 与 shared）、kernel decomposition（把一个算法拆成多次 kernel launch，用 kernel 边界充当全局同步）、协作装载（cooperative fetching）、__shfl_down_sync 类的 warp 级原语、以及”用 L 与 V 两个复用旋钮描述一个 kernel 的访存量”的分析方法。

• 在并行计算知识体系中的角色：本讲是从”会写并行程序”到”会做性能工程”的转折点。它把前几讲的 work-span / 算术强度 / Roofline / 带宽-延迟工具，第一次放到一台真实且极端的机器上做定量诊断：同一份矩阵乘，未优化版本是 0.25 FLOP/B 的访存灾难，调好 L 与 V 之后变成计算受限。后续的 Performance Optimization、异构与硬件专用化、并行深度学习（tensor core、算子融合、算子间复用）全部建立在本讲的三条判据和两个 case study 之上。

• 配套材料：
  • lectures/05-CUDA-programming.pdf（抽取文本 extracted/05-CUDA-programming.txt，共 83 页）：已公开，可在 https://www.cs.cmu.edu/~418/lectures/ 直接下载。讲义首页标题为 “Lecture 5 & 6: GPU Architecture & CUDA Programming”，即一份讲义覆盖 Sep 2 与 Sep 4 两次课（Fall 2026 日程表：https://www.cs.cmu.edu/~418/schedule.html，第 5 讲 Sep 2，第 6 讲 Sep 4 “GPU Architecture and CUDA Programming (continued)”）。本讲（Lecture 6）聚焦该讲义的后半部分：CUDA 到硬件的映射细节、Case Study 1（matmul）、Case Study 2（reduction）、以及 recap 中列出的五条 GPU 优化技术。讲义首页含历史学期的署名信息（如邀请讲者与历史页脚），属讲义沿用现象，不是错误。
  • doc_CUDA-recitation.pdf（抽取文本 extracted/doc_CUDA-recitation.txt，共 39 页）：CUDA recitation 讲义，已公开（页脚标注历史学期 Fall 2025，同为沿用）。内容为”怎么写 CUDA”与”怎么写快 CUDA”：三种函数限定符、dim3 与线程索引公式、host/device 内存管理的 7 个步骤、归约的四种渐进方案（kernel 当屏障 / __syncthreads() / 单 block / 多 block + 共享内存）、matmul 的寄存器分块与共享内存分块。
  • cs149_supp/gpuarch.txt：Stanford CS149（Fall 2025）Lecture 7 “GPU Architecture & CUDA Programming”（共 74 页），已公开的补充读物；本笔记中 V100 SM 的 sub-core 微架构、warp selector 与指令发射时序、work scheduler 的 7 步动态映射、”为什么必须为整个 block 预留执行上下文”、”persistent thread” 编程风格、histogram 原子操作合法性等细节来自该材料。
  • 讲课录像（Panopto / YouTube）：Fall 2026 日程表中被注释隐藏，属未发布。
  • Ed 讨论区、Autolab、Canvas：需登录，非公开。
  • Fall 2026 尚未在公开目录发布的讲义（Performance Analysis / Profiling、Transactional Memory、AI in System Design 等），其历史学期 PDF 位于 /afs/cs/academic/class/15418-*/public/ 之下，需要 CMU 登录，属未公开。
  • Fall 2026 授课教师为 Brian Railing 与 Dimitrios Skarlatos；课程由 Kayvon Fatahalian 创建。

2. 核心概念与硬件/软件架构图解

2.1 三条性能判据：本讲的全部实用内容

• 定义与目的：讲义的 recap 把 GPU 优化技术收敛成五条：coherent warps（warp 执行一致）、coalesced memory access（合并访存）、shared memory bank conflict（避免 bank 冲突）、warp level optimizations（warp 级优化，如 shuffle 代替共享内存）、Tensor Core（把矩阵乘交给专用单元）。它们共同回答一个问题：GPU 的峰值算力要靠”每个时钟都有满宽的 warp 在执行有用的指令”来兑现；任何让 lane 空转、让访存拆成多次事务、让数据反复穿越存储层次的行为，都是在这个兑现率上打折。

• 直观解释（”它是什么？”）：把 SM 想成一条 32 人并排作业的装配线（32 宽 SIMD）。三条判据分别对应三种浪费：
  1. 发散（divergence） = 32 个工位里只有 8 个在干活，其余 24 个被”遮罩（mask）”住站着不动 → 装配线照常开动，产出只有 1/4；
  2. 非合并访存（non-coalesced） = 32 个工人各自去仓库不同货架取一个小螺丝，仓库被迫发 4 趟车而不是 1 趟 → 同样的数据量，带宽利用率掉到 1/4；
  3. 共享内存 bank 冲突（bank conflict） = 车间的 32 个储物柜（bank）里，有 16 个人同时挤在第 0、2、4…号柜子前，另外 16 个柜子空着 → 存取要排队两轮。
  4. 再叠加一条数据复用：从仓库（HBM）取回的料要在车间（寄存器/共享内存）被尽可能多次使用，否则装配线的速度完全由仓库门口决定。
• 架构/机制图解：

图 1：SM 内部结构 —— 4 个 sub-core、warp selector 与 SIMT 掩码

 +--------------------------------------------------------------------------+
 |            SM (Streaming Multiprocessor)  —— 一个"GPU 核"                 |
 |                                                                          |
 |  +---------------------------+   +---------------------------+           |
 |  | Sub-core 0                |   | Sub-core 1                |  ... x4   |
 |  |   Warp Selector           |   |   Warp Selector           |           |
 |  |        |                  |   |        |                  |           |
 |  |   Fetch / Decode          |   |   Fetch / Decode          |           |
 |  |        |                  |   |        |                  |           |
 |  |   +----+----+----+----+   |   |   (每个 sub-core 每个时钟  |           |
 |  |   |fp32 |int |fp64|LSU |   |   |    选 1 个可运行 warp)    |           |
 |  |   |16   |16  |8   |    |   |   |                           |           |
 |  |   |lane |lane|lane|    |   |   |                           |           |
 |  |   +-----+----+----+----+   |   +---------------------------+           |
 |  |   [ Tensor Core unit ]    |                                           |
 |  |                           |                                           |
 |  |   寄存器堆 64 KB/sub-core |   一个 warp = 32 个线程的标量寄存器组      |
 |  |   （R0,R1,... 每线程一套）|   每 sub-core 最多交错 16 个 warp         |
 |  +---------------------------+                                           |
 |                                                                          |
 |  +------------------------------------------------------------------+    |
 |  | Shared memory + L1 cache 存储（V100: 128 KB，可配置划分）        |    |
 |  +------------------------------------------------------------------+    |
 +--------------------------------------------------------------------------+
        |                                             ^
        v                                             |
 +---------------------+   L2 cache (V100: 6 MB)  +------------------------+
 | Register / Shared   |<------------------------>|  HBM 全局内存          |
 | (片上, 高带宽低延迟)|                          |  V100: 16 GB, 900 GB/s |
 +---------------------+                          +------------------------+

关键性能特征（来自讲义与 CS149 补充材料给出的 V100 数据）：

• warp 发射节奏：一个 32 宽 fp32 SIMD 操作每 2 个时钟完成一次（16 lane × 2 时钟），fp64 是每 4 个时钟一次（8 lane）。因此”一个 warp 的一条 fp32 指令”占用功能单元 2 个时钟。
• SM 的每时钟动作（418 讲义描述）：从驻留的最多 64 个 warp（64 × 32 = 2048 个 CUDA 线程）中选出最多 4 个可运行 warp（对应 4 个 sub-core 的线程级并行），每个被选中的 warp 再选出最多 2 条可运行指令（指令级并行）。
• 寄存器预算：V100 每 SM 256 KB 寄存器 = 65536 个 32-bit 寄存器；若要把 64 个 warp × 32 线程 = 2048 个线程全部驻留，则每线程只有 32 个寄存器。这是”寄存器分块（register tiling）”的硬上限——V 越大，accumulator 越多，occupancy 越低。

图 2：SIMT —— 硬件动态检测 32 条标量指令流是否”恰好一致”

   32 个独立线程，各自一条标量指令流            硬件把它们合并成 1 条 32 宽 SIMD 指令
   --------------------------------------       ------------------------------------
   t0 : mul r0 r1 r2  ->  ...
   t1 : mul r0 r1 r2  ->  ...                        lane: 0  1  2  3 ... 31
   t2 : mul r0 r1 r2  ->  ...                        ALU : [A][A][A][A]...[A]   32 宽
   ...                                              数据: d0 d1 d2 d3 ... d31
   t31: mul r0 r1 r2  ->  ...                        => 1 条指令 / 全宽 / coherent

   一旦某个 lane 走了别的分支：
   t6 : (分支到其他代码)                              lane: 0  1  2  3  4  5  6 ... 31
                                                     ALU : [A][A][A][A][A][A][-][A]..[A]
                                                                          ^
                                                          lane 6 被 mask（空转） —— divergent
   之后两条分支路径必须**串行**执行：先跑 if 体（mask 掉 else 的 lane），再跑 else 体。
   warp 的执行时间 ≈ 各分支路径用时之和，效率 ≈ 有用 lane 数 / (32 × 路径数)

• 性能后果：coherence（一致执行）是 GPU 高效的前提，divergence 应当被最小化。注意 CUDA 的 warp 是实现细节而非编程模型的一部分（例外是 __shfl_*、vote 等 warp 级内建操作）：ISPC 在编译期把程序编译成 SIMD 指令，而 GPU 是运行期由硬件动态判断这 32 个独立线程是否恰好走在同一条指令上。

2.2 线程块调度：work scheduler 如何在时间上填满 GPU

• 定义与目的：CUDA 的核心假设是 “线程块可以按任意顺序执行，块与块之间没有依赖”；GPU 用一个硬件 work scheduler 把 block 动态映射到 SM，映射时必须尊重编译产物中记录的资源需求（每块多少线程、每块多少共享内存、每线程多少本地数据）。这样同一份 CUDA 程序无需改动就能跑在 16 SM 的 GTX 980 和 132 SM 的 H100 上——程序里没有任何 “num_cores” 概念。

• 直观解释（”它是什么？”）：这是课程里反复出现的 “工人池 + 任务分解”设计模式的又一个实例：GPU 的 SM 是固定数量的工人（数目由硬件决定，与请求数无关），线程块是待办任务，work scheduler 是派工头。同样的模式出现在 ISPC 的 task 运行时（为每个 CPU 超线程创建一个 pthread 并长期保活）和 Web 服务器的线程池（线程数由核数决定，不由并发请求数决定）中。

图 3：一个 kernel 的 1000 个 block 在”假想双核 GPU”上的动态映射（Gantt 视图）

 假想双核 GPU（讲义的简化模型）：
   每个 core: 执行上下文可容纳 384 个 CUDA 线程（12 warps）、共享内存 1.5 KB
 convolve 的 block 需求: 128 线程 + 130 floats = 520 B 共享内存
 => 每 core 只能同时驻留 2 个 block: 2 x 520 B = 1040 B <= 1536 B
                                    3 x 520 B = 1560 B >  1536 B  (放不下第三个)

 时间 -->
 Core 0: [ Blk0 ][ Blk0 ][ Blk4 ][ Blk4 ][ Blk8 ][ Blk8 ] ...   (NEXT 指针: 0,2,4,6,8,...)
         [ Blk2 ][ Blk2 ][ Blk6 ][ Blk6 ][ Blk10] ...           (交错映射: 0,1,2,3,... 轮流)
 Core 1: [ Blk1 ][ Blk1 ][ Blk5 ][ Blk5 ][ Blk9 ] ...
         [ Blk3 ][ Blk3 ][ Blk7 ][ Blk7 ][ Blk11] ...
                ^
                └─ Blk0 完成后（Step 4），其 128 个执行上下文与 520 B 共享内存立刻被释放，
                   scheduler 在 Step 5 把 Blk4 映射到 contexts 0-127

 EXECUTE: convolve | ARGS: N, input, output | NUM_BLOCKS: 1000 | NEXT=4 | TOTAL=1000

• 关键机制与性能含义：
  1. block 的粒度是”全有或全无”：一个 block 的所有线程必须同时获得执行上下文（寄存器状态），因为块内线程之间可能有依赖（最简单的例子是 __syncthreads()）。CS149 的追问很关键：为什么不能先把线程 0-127 跑完再跑 128-255？因为块内线程按语义就是并发的——”如果块内某线程可运行，它最终一定会被运行（不会死锁）”。
  2. 资源决定 occupancy：真正的 GPU 上，SM 的驻留 block 数由 warp 执行上下文（64 个 warp）、共享内存（V100 96/128 KB 级，H100 256 KB）、寄存器总量三者取最小决定。
  3. 动态调度天然做负载均衡：块大小相同、块数远多于 SM 数时，靠”先完成先派新块”即可平滑负载不均；但块数太少（< SM 数 × 驻留数）会造成尾部空转。

2.3 精确的执行语义边界：块内并发、块间任意序

• 定义与目的：CUDA 的执行语义是两种模型的混合，混淆它会导致两类典型错误。
  • block 之间：任何顺序、逻辑并发、无依赖 → 对应数据并行模型的 forall（机器无关，由系统调度到任意数量的核上）。
  • block 内部：线程真正并发运行、共享地址空间、可协作 → 对应 SPMD 共享地址空间模型（很像一个 ISPC gang）。

• 直观解释（”它是什么？”）：把 GPU 想成一栋楼。块 = 一个项目组，组内成员在同一间办公室（SM + 共享内存）里实时协作、可以开会同步（__syncthreads()）；块之间 = 不同项目组，公司只保证”每个组最终都会被安排办公室”，但不保证谁先谁后、也不保证同时存在。所以：组内可以互相等，跨组不能互相等——一等的就是死锁。

• 重要推论：
  • 合法：所有 CUDA 线程用 atomicAdd 更新同一个 global 数组（histogram 例子）。原子操作用作互斥，不影响实现的调度自由——因为无论块以什么顺序执行，结果都对（加法可交换），并且不存在”等别的块”的行为。
  • 不合法（可能死锁）：block 0 设置 myFlag，block 1 自旋等待 myFlag != 0。若 GPU 只有 1 个 SM 且只能驻留 1 个 block，且实现先跑了 block 1，则 block 1 永久自旋；而 block 0 永远得不到执行资源。CUDA 没有全局同步（no global synchronization）：一方面为高 SM 数的 GPU 建造全局屏障硬件代价高，另一方面当 #blocks > #SM × #resident blocks 时全局屏障必然死锁。
  • 正确替代：kernel 分解（kernel decomposition）——用 kernel launch 边界充当全局同步；讲义强调 kernel launch 的硬件/软件开销很低，而且”所有层次的代码长得一样”，只是把 for 循环搬到 host 侧。

图 4：用 kernel 边界代替全局屏障（归约的两阶段分解）

  单个 kernel 内做不到:
  ---------------------------------------------------------------
  Block0  Block1  Block2  ...  BlockN      <- 块间无法 __syncthreads()
  (x)     (x)     (x)          (x)         跨块自旋 = 可能死锁
  ---------------------------------------------------------------

  分解为两次 launch:
  阶段 1 (grid = 4096 blocks)              阶段 2 (grid = 1 block)
  ---------------------------------------------------------------
  Block0  Block1  ...  Block4095           Block0
   |        |            |                   |
  部分和    部分和       部分和   ────────>  把 4096 个部分和再归约
   |        |            |                   |
  partial[0..4095]                           out[0]
  ---------------------------------------------------------------
       ^ kernel launch #1 结束 = 一次隐式全局屏障（所有块都完成）
       ^ kernel launch #2 结束 = 一次隐式全局屏障

  反例（recitation 中的 "Idea 1"）：把树的每一层都做成一次 launch
     d=1,2,4,...,2^23  =>  24 次 launch。正确，但每次 launch 数微秒级固定开销，
     24 次 x ~5 us ≈ 120 us，远超 67 MB 数据本身 75 us 的搬运时间 => 被 launch 开销吃掉

2.4 访存效率之一：global memory 的 coalescing

• 定义与目的：coalescing（合并访存） 指”一个 warp 的多个线程访问连续的内存地址”，从而把 32 个 lane 的请求合并成尽可能少的内存事务（memory transaction），最大化 HBM 带宽利用率。一个 warp 的 32 个线程各取 4 字节且地址连续时，正好覆盖 128 字节——这是硬件最喜欢的一次全宽事务。

• 直观解释（”它是什么？”）：快递取件。32 个人要取 32 个相邻柜子里的包裹——仓库管理员开一次门、走一趟就全部取回（coalesced）；如果这 32 个人要取的包裹分散在 32 个不同货架，管理员就得跑 32 趟（non-coalesced），实际带宽利用率可能掉到 1/4 甚至更低，而且每一趟都会把整条 cache line 搬回来，浪费的是带宽而不是容量。

图 5：合并访存 vs 非合并访存（以及共享内存的 bank 映射）

 [A] Coalesced（最优）          地址: 0 1 2 3 | 4 5 6 7 | 8 9 10 11  (以 4B 字为单位)
     warp 的 32 个 lane 连续访问
     t0 t1 t2 t3   t4 t5 t6 t7   t8 t9 t10 t11 ...
     |  |  |  |    |  |  |  |    |  |   |   |
     +--+--+--+    +--+--+--+    +--+---+---+     128 B = 1 次全宽事务
     第一次 load    第二次 load    第三次 load

 [B] Non-coalesced（次优）     warp 内 stride = 4 个字（16 B）
     t0 -> word 0     t1 -> word 4     t2 -> word 8    t3 -> word 12
     |                |                |               |
     +----------------+----------------+---------------+   每 lane 落在不同 128 B 段
     一次 warp 请求被拆成 4 段 => 有效带宽 1/4；但搬回的字节数不变 => 浪费 4 倍带宽

 [C] 共享内存 bank（V100/H100：32 个 bank，每 bank 4 B，每时钟可服务 1 个字）
     无冲突:  lanes t0..t7 -> words 0,1,2,3,4,5,6,7
        B0  B1  B2  B3  B4  B5  B6  B7 ... B31
        t0  t1  t2  t3  t4  t5  t6  t7 ...  -
        => 每个 bank 服务 1 个 lane，1 个周期完成，满带宽

     冲突:    lanes t0..t7 -> words 0,2,4,6,8,10,12,14   (stride 2)
        B0     B1   B2     B3   B4     B5   B6     B7
        t0,t4  -    t1,t5  -    t2,t6  -    t3,t7  -
        => 4 个 bank 各服务 2 个 lane，其余 4 个 bank 空闲
           => 2 路 bank 冲突，访问串行化为 2 个周期，带宽减半

• 数值化含义：如果一次 warp 访存被拆成 4 个事务，那么”每字访问的有效成本”是原来的 4 倍；在 V100（900 GB/s、12.7 TFLOPs、机器平衡点 14.2 FLOP/B）上，一个本该 0.15 ms 的 kernel 会退化到 0.6 ms——正好抵消掉 4 倍的算力优化。

2.5 三级存储层次与延迟/带宽/容量

• 定义与目的：CUDA 对 kernel 可见三个地址空间：per-thread private（寄存器/本地内存）、per-block shared memory、device global memory。它们不是”三块内存”，而是程序中三种不同 locality 的体现；把数据放到正确的层次，是 GPU 性能工程的核心动作。

• 两个关键事实：
  1. host 与 device 地址空间彼此不可直接访问。要移动数据必须显式调用 cudaMemcpy（参数中的方向常量 cudaMemcpyHostToDevice / cudaMemcpyDeviceToHost）。这一条正是”消息传递模型”的味道——所以讲师会反复问”CUDA 是共享地址空间模型还是消息传递模型？”答案取决于你站在哪一层看：host↔device 之间像消息传递，block 内部是共享地址空间。
  2. __shared__ 存在的唯一理由是让 block 内的线程协作：它把”每个线程各自从 global 取 3 个元素（3×128 次 load）”变成”全体协作取 130 个元素（130 次 load），再从片上反复读”。1D 卷积版本 1 → 版本 2 的整个收益就是这个 3N → (N + 2·blocks) 的访存缩减。

2.6 Case Study 1 的几何：矩阵乘的两级复用

• 定义与目的：C = A × B（M×K 乘 K×N 得 M×N）。朴素实现是访存灾难：每个线程算 1 个输出元素、每个元素需要读 A 的一整行与 B 的一整列（2N 次 global 访问），N² 个线程 → 总访存 2N³，而总计算只有 2N³ FLOP → 算术强度 0.25 FLOP/B，比机器平衡点（V100：12.7 TFLOPs / 900 GB/s = 14.2 FLOP/B）低 57 倍。两个”复用旋钮”解决它：
  • V（thread-level register tiling，寄存器分块）：一个线程算 V×V 个输出，把 2V 次加载的数据用于 V² 次乘加 → 总访存降到 2N³/V。
  • L（block-level shared memory tiling，共享内存分块）：一个 block 算 L×L 个输出，A/B 的 L×K 与 K×L 条带被整个 block 复用 → 总访存降到 2N³/L。
• 直观解释（”它是什么？”）：做菜类比。朴素做法是每做一道菜就跑一趟菜市场买齐所有原料（每个输出元素都重新读 A 行、B 列）。L 旋钮 = 让一个厨师班组统一采购一批原料放在班组自己的冰箱（shared memory）里，全班一起用；V 旋钮 = 每个厨师一次拿一小篮（寄存器 a[V]、b[V]）原料，同时做 V² 道菜（一次拿料，多次下锅）。L 减少”跑市场的总趟数”，V 减少”从冰箱里取料的总次数”。

图 6：矩阵乘的两级分块（block 级 L×L + 线程级 V×V）

                 K                                  N
        A  (M x K)                          B  (K x N)
   +-------------------+               +-------------------------+
   |                   |               |        ^                |
   |   +-------+       |               |   +----+----+           |
 M |   | L x S |<--------------------->|   | S x L   |           |
   |   | 瓦片  |       |               |   |  瓦片   |           |
   |   +-------+       |               |   +---------+           |
   |   ^ block 协同装载|               |   ^ block 协同装载      |
   +-------------------+               +-------------------------+

   共享内存: sA[S][L], sB[S][L]   (S = K 方向的分块步长)

                 N
        C  (M x N)
   +-------------------------+
   |      |      |           |
   |  +---+---+  |           |     一个 block 负责 L x L 的输出瓦片
   |  | L x L |  |           |     瓦片内每个线程负责一个 V x V 小方块
 M |  |  瓦片 |  |           |     线程内 4 个 (V=2) accumulator 提供 4 条独立
   |  +---+---+  |           |     的 FMA 依赖链 -> 指令级并行 (ILP)
   |      |      |           |
   +-------------------------+
              ^
              +-- 每个线程: 2V 次 sA/sB 加载  ->  V^2 次 FMA  ->  访存:计算 = 2/V

• 定量结论（讲义给出的公式）：

注解：表中”算术强度”按”2 FLOP / 4 B”的比率推出（每个 4 字节浮点访问贡献 2 FLOP）；共享内存分块后 L 负责砍 global 流量，V 负责砍 shared 流量——两个旋钮各管一层，这就是讲义那句 “L and V are the two reuse knobs!”。

• 协作装载（cooperative fetching）：把 sA[:,:] = A[...] 这种”整块赋值”翻译成每个线程搬一片：

   nthreads = blockDim.x * blockDim.y
   tid      = threadIdx.y * blockDim.x + threadIdx.x
   for (j = 0; j < L*S / nthreads; ++j) {
       y = (j*nthreads + tid) / L
       x = (j*nthreads + tid) % L
       s[y][x] = A[k + y][yblock*L + x]        // 连续 tid -> 连续 x => 合并访存
   }

关键点：让 tid 映射到连续的 x（列）上，global 读取就是连续的 → coalesced；随之而来的 __syncthreads() 必须成对出现：装载前一个屏障（防止上一个 k 迭代还在读的时候被覆盖）、装载后一个屏障（保证所有数据就位后才能开始计算）。

2.7 Case Study 2 的几何：树形归约与三种寻址

• 定义与目的：归约（reduction）是 normalization、softmax 等 ML 算子的底层原语，任务是把 n 个元素求和。串行版是 n 次依赖加法（Span = n）；树形归约把依赖链缩短到 log n：for (d=1; d<n; d<<=1) for (i=0; i<n; i+=2d) A[i] += A[i+d]。

• 直观解释（”它是什么？”）：淘汰赛。串行求和像”一个裁判逐个听 n 个人报数”，要 n 轮；树形归约像”两两配对比赛，胜者（和）进入下一轮”，只要 log₂n 轮。但比赛必须一轮一轮同步——CUDA 不保证块内线程 lockstep 执行，所以每轮之间必须插 __syncthreads()（块内屏障），或用 kernel 边界（跨块屏障）。这解释了 recitation 里”为什么需要两个 __syncthreads()“以及”为什么只加一个 __syncthreads() 的版本是错的”。

图 7：归约的三种寻址方式（同一棵树的三种”排队方式”）

 [1] Interleaved addressing（交错寻址, reduce0）   if (tid % (2s) == 0) sdata[tid] += sdata[tid+s]
     s=1: 线程  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
          活跃  T  F  T  F  T  F  T  F  T  F  T  F  T  F  T  F      <- 半数 lane 空转
     s=2: 活跃  T  F  F  F  T  F  F  F  T  F  F  F  T  F  F  F      <- 1/4 有效
     s=4: 活跃  T  F  F  F  F  F  F  F  T  F  F  F  F  F  F  F      <- 1/8 有效
     问题: warp 高度发散; 且奇偶 lane 访问 sdata[0],[2],[4]... => 2 路 bank 冲突

 [2] Strided index + 非发散分支        index = 2*s*threadIdx.x; if (index < blockDim.x) ...
     s=1: 活跃  T  T  T  T  T  T  T  T  F  F  F  F  F  F  F  F      <- 前 8 个 lane 连续活跃
     s=2: 活跃  T  T  T  T  F  F  F  F  F  F  F  F  F  F  F  F
     s=4: 活跃  T  T  F  F  F  F  F  F  F  F  F  F  F  F  F  F
     改善: 分支结果在 warp 内"前缀连续", 后期线程统一为 false; 但访问 sdata[0],[2],[4]... 仍 stride 2

 [3] Sequential (reversed) addressing（顺序寻址, 最终版）
     for (s = blockDim.x/2; s > 0; s /= 2) { if (threadIdx.x < s) sdata[tid] += sdata[tid+s]; }
     s=256: 活跃 lane 0..255  (8 个 warp 满负荷!)   访问 sdata[tid] 与 sdata[tid+256]
     s=128: 活跃 lane 0..127  (4 个 warp)
     s=64 : 活跃 lane 0..63   (2 个 warp)
     s=32 : 活跃 lane 0..31   (1 个 warp, 满)
     s=16,8,4,2,1: 单 warp 内前缀活跃 (不可避免的尾部), 但访问 tid, tid+1, ... 连续
     结果: 完全合并访存 + 无 bank 冲突 (相邻 lane 落在相邻 bank)

 归约树示意 (n=16, 顺序寻址):
   第1轮: [a0+a8][a1+a9][a2+a10][a3+a11][a4+a12][a5+a13][a6+a14][a7+a15]
   第2轮: [b0+b4][b1+b5][b2+b6][b3+b7]
   第3轮: [c0+c2][c1+c3]
   第4轮: [d0+d1]                      总轮数 = log2(16) = 4

• 性能含义：三种写法工作量相同（Work ≈ n 次加法），差别全在”这 n 次加法有多少个 lane 是有效的、访存有没有被拆成多个事务”。这正是”同样的 Work/Span，性能差 5–10 倍”的教科书案例。

2.8 算力演进与 Tensor Core

• 定义与目的：Tensor Core 是 SM 内的矩阵乘加专用单元（matrix multiplication unit）。它把”小尺寸矩阵乘加”固化成一条指令，用远高于通用 SIMD lane 的密度完成乘加，是 H100 峰值算力从 GTX 980 的 4.6 TFLOPs 跃升到 1000 TFLOPs 的主要原因。

• 直观解释（”它是什么？”）：通用 fp32 lane 像一群会做任何算术的通用工人；Tensor Core 像一台专门的”批量乘法打孔机”——它只会做矩阵乘加这一件事，但一次冲压就是一个 4×4（乃至更大）的乘加块。代价是你必须把数据摆成它要求的形状（分块、对齐、数据类型），”喂料”的组织工作落到了程序员/编译器身上。

• 注意三次”没变”：从 GTX 980 到 H100，SM 的基本形态没变、每 SM 的 warp 数上限（64）没变、warp 宽度（32）没变。变的是 SM 数量（16 → 132）、共享内存容量（96 KB → 256 KB）与专用单元（Tensor Core）。这解释了为什么本讲的三条判据在十年后依然成立。

3. 代码示例与性能分析

3.1 示例 A：两级分块的 CUDA 矩阵乘（Case Study 1，完整可运行）

// mm_tiled.cu  ——  C = A * B  (M x K 乘 K x N), 两级分块: block 级 L x L + thread 级 V x V
// 编译(发布优化): nvcc -O3 -arch=sm_80 -lineinfo mm_tiled.cu -o mm_tiled
// 运行:          ./mm_tiled
// 说明: 报告 GFLOP/s 与"有效全球访存带宽"，并与 CPU 朴素实现校验正确性。

#include <cuda_runtime.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>

// ---------- 分块参数（两个复用旋钮） ----------
#define BM 64          // block 计算的输出瓦片行数  L = 64
#define BN 64          // block 计算的输出瓦片列数  L = 64
#define BK 16          // K 方向分块步长 S = 16
#define V  4           // 每个线程算 V x V 个输出    V = 4

#define THREADS_X (BN / V)          // 16
#define THREADS_Y (BM / V)          // 16   => block = 16 x 16 = 256 线程

#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(1);                                                          \
        }                                                                     \
    } while (0)

__global__ void mm_tiled(const float* __restrict__ A,
                         const float* __restrict__ B,
                         float* __restrict__ C,
                         int M, int N, int K)
{
    __shared__ float sA[BM][BK];    // A 的 BM x BK 瓦片: 64 x 16 x 4 B = 4 KB
    __shared__ float sB[BK][BN];    // B 的 BK x BN 瓦片: 16 x 64 x 4 B = 4 KB

    const int ty = threadIdx.y;
    const int tx = threadIdx.x;
    const int tid = ty * THREADS_X + tx;
    const int nthreads = THREADS_X * THREADS_Y;   // 256

    // 该线程负责的 C 子块左上角（V x V）
    const int crow = blockIdx.y * BM + ty * V;
    const int ccol = blockIdx.x * BN + tx * V;

    float acc[V][V];
#pragma unroll
    for (int i = 0; i < V; ++i)
#pragma unroll
        for (int j = 0; j < V; ++j) acc[i][j] = 0.0f;

    float a[V], b[V];

    for (int ko = 0; ko < K; ko += BK) {
        // ---- 协作装载 A 的 BM x BK 瓦片: 每个线程搬 4 个元素, tid -> 连续列 => 合并访存 ----
#pragma unroll
        for (int t = tid; t < BM * BK; t += nthreads) {
            const int r = t / BK;
            const int c = t % BK;
            sA[r][c] = A[(size_t)(blockIdx.y * BM + r) * K + (ko + c)];
        }
        // ---- 协作装载 B 的 BK x BN 瓦片 ----
#pragma unroll
        for (int t = tid; t < BK * BN; t += nthreads) {
            const int r = t / BN;
            const int c = t % BN;
            sB[r][c] = B[(size_t)(ko + r) * N + (blockIdx.x * BN + c)];
        }
        __syncthreads();                    // 数据就位前不许算

#pragma unroll
        for (int k = 0; k < BK; ++k) {
#pragma unroll
            for (int i = 0; i < V; ++i) a[i] = sA[ty * V + i][k];
#pragma unroll
            for (int j = 0; j < V; ++j) b[j] = sB[k][tx * V + j];
#pragma unroll
            for (int i = 0; i < V; ++i)
#pragma unroll
                for (int j = 0; j < V; ++j)
                    acc[i][j] += a[i] * b[j];    // V^2 条独立 FMA 链 => ILP
        }
        __syncthreads();                    // 下一轮装载前，保证所有线程都读完本瓦片
    }

#pragma unroll
    for (int i = 0; i < V; ++i)
#pragma unroll
        for (int j = 0; j < V; ++j)
            C[(size_t)(crow + i) * N + (ccol + j)] = acc[i][j];
}

// ---------------- CPU 参考实现（仅用于校验小规模结果） ----------------
static void mm_cpu(const float* A, const float* B, float* C, int M, int N, int K)
{
    for (int i = 0; i < M; ++i)
        for (int j = 0; j < N; ++j) {
            double s = 0.0;
            for (int k = 0; k < K; ++k) s += (double)A[(size_t)i * K + k] * B[(size_t)k * N + j];
            C[(size_t)i * N + j] = (float)s;
        }
}

int main()
{
    // ---- 1) 小规模正确性校验: 256 x 256 x 256 ----
    {
        const int M = 256, N = 256, K = 256;
        const size_t bA = (size_t)M * K * sizeof(float);
        const size_t bB = (size_t)K * N * sizeof(float);
        const size_t bC = (size_t)M * N * sizeof(float);
        float *hA = (float*)malloc(bA), *hB = (float*)malloc(bB);
        float *hC = (float*)malloc(bC), *ref = (float*)malloc(bC);
        for (size_t i = 0; i < (size_t)M * K; ++i) hA[i] = (float)((i * 7) % 13) * 0.125f - 0.5f;
        for (size_t i = 0; i < (size_t)K * N; ++i) hB[i] = (float)((i * 5) % 11) * 0.125f - 0.5f;

        float *dA, *dB, *dC;
        CUDA_CHECK(cudaMalloc(&dA, bA));
        CUDA_CHECK(cudaMalloc(&dB, bB));
        CUDA_CHECK(cudaMalloc(&dC, bC));
        CUDA_CHECK(cudaMemcpy(dA, hA, bA, cudaMemcpyHostToDevice));
        CUDA_CHECK(cudaMemcpy(dB, hB, bB, cudaMemcpyHostToDevice));

        dim3 block(THREADS_X, THREADS_Y);              // (16, 16) = 256 线程
        dim3 grid(N / BN, M / BM);                     // (4, 4) = 16 个 block
        mm_tiled<<<grid, block>>>(dA, dB, dC, M, N, K);
        CUDA_CHECK(cudaGetLastError());
        CUDA_CHECK(cudaDeviceSynchronize());
        CUDA_CHECK(cudaMemcpy(hC, dC, bC, cudaMemcpyDeviceToHost));
        mm_cpu(hA, hB, ref, M, N, K);

        double maxerr = 0.0;
        for (size_t i = 0; i < (size_t)M * N; ++i)
            maxerr = fmax(maxerr, fabs((double)hC[i] - ref[i]));
        printf("[verify] 256x256x256  最大绝对误差 = %.6g  (K=256 时 fp32 累加误差量级 ~1e-3)\n", maxerr);

        cudaFree(dA); cudaFree(dB); cudaFree(dC);
        free(hA); free(hB); free(hC); free(ref);
    }

    // ---- 2) 性能测量: 1024 x 1024 x 1024 ----
    {
        const int M = 1024, N = 1024, K = 1024;
        const size_t bA = (size_t)M * K * sizeof(float);
        const size_t bB = (size_t)K * N * sizeof(float);
        const size_t bC = (size_t)M * N * sizeof(float);
        float *hA = (float*)malloc(bA), *hB = (float*)malloc(bB), *hC = (float*)malloc(bC);
        for (size_t i = 0; i < (size_t)M * K; ++i) hA[i] = (float)((i * 7) % 13) * 0.01f;
        for (size_t i = 0; i < (size_t)K * N; ++i) hB[i] = (float)((i * 5) % 11) * 0.01f;

        float *dA, *dB, *dC;
        CUDA_CHECK(cudaMalloc(&dA, bA));
        CUDA_CHECK(cudaMalloc(&dB, bB));
        CUDA_CHECK(cudaMalloc(&dC, bC));
        CUDA_CHECK(cudaMemcpy(dA, hA, bA, cudaMemcpyHostToDevice));
        CUDA_CHECK(cudaMemcpy(dB, hB, bB, cudaMemcpyHostToDevice));

        dim3 block(THREADS_X, THREADS_Y);              // (16, 16) = 256 线程
        dim3 grid(N / BN, M / BM);                     // (16, 16) = 256 个 block

        mm_tiled<<<grid, block>>>(dA, dB, dC, M, N, K);  // 预热
        CUDA_CHECK(cudaDeviceSynchronize());

        cudaEvent_t e0, e1;
        CUDA_CHECK(cudaEventCreate(&e0));
        CUDA_CHECK(cudaEventCreate(&e1));
        const int iters = 20;
        CUDA_CHECK(cudaEventRecord(e0));
        for (int it = 0; it < iters; ++it)
            mm_tiled<<<grid, block>>>(dA, dB, dC, M, N, K);
        CUDA_CHECK(cudaEventRecord(e1));
        CUDA_CHECK(cudaEventSynchronize(e1));
        float ms = 0.f;
        CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
        ms /= iters;

        const double flops = 2.0 * (double)M * N * K;              // 2N^3
        const double gflops = flops / (ms * 1e-3) / 1e9;
        // 该 kernel 在 "至 L2" 层面搬运的 global 字节数: 2N^3/L * 4 B
        const double globalBytes = 2.0 * (double)M * N * K / BM * 4.0;
        printf("[perf  ] N=%d  L=%d  V=%d  平均耗时 = %.3f ms\n", M, BM, V, ms);
        printf("[perf  ] %.1f GFLOP/s ; global(L2 级)流量 %.2f MB/次 -> %.1f GB/s\n",
               gflops, globalBytes / 1e6, globalBytes / (ms * 1e-3) / 1e9);

        CUDA_CHECK(cudaMemcpy(hC, dC, bC, cudaMemcpyDeviceToHost));
        printf("[check ] C[0]=%.4f  C[last]=%.4f\n", hC[0], hC[M * N - 1]);

        cudaFree(dA); cudaFree(dB); cudaFree(dC);
        free(hA); free(hB); free(hC);
    }
    return 0;
}

【代码做什么？】

1. host 侧：分配并初始化 A、B（device 端），分配 C；设置 block = (BN/V, BM/V) = (16,16) = 256 线程、grid = (N/BN, M/BM) 个 block；预热一次后连跑 20 次取平均；用 cudaEvent 计时；最后把 C 拷回 host 做小规模校验。
2. device 侧外层循环（K 方向分块）：每次迭代把 A 的 BM×BK 瓦片与 B 的 BK×BN 瓦片由 256 个线程协作装载进共享内存（每个线程搬 BM*BK/256 = 4 个元素），__syncthreads() 后进入内层。
3. device 侧内层循环（K 瓦片内）：对 BK 个 k 值，每个线程从 sA 取 V=4 个 A 元素、从 sB 取 V=4 个 B 元素，做 V²=16 次 FMA 累加到 acc[4][4]；结束时再 __syncthreads() 保护瓦片不被下一轮提前覆盖。
4. 索引到线程的映射：crow = blockIdx.y*BM + ty*V，ccol = blockIdx.x*BN + tx*V——同一个 warp 里的连续 tx 落在连续的 C 列上，因此最后写回 C 时也是合并访存；协作装载时 tid 也对齐到连续列，保证读 A/B 合并。

【并行机制与性能解说】

• 线程/block 的创建与工作分配：一次 launch 产生 16×16 = 256 个 block、每块 256 线程，共 65536 个 CUDA 线程，映射为 256 / 32 = 8 个 warp/block。硬件 work scheduler 把 block 动态派到 132（H100）个 SM 上，程序本身不关心 SM 数量。
• 共享数据的处理：A/B 瓦片是所有线程的只读共享数据，放在 __shared__（每 block 8 KB）里被 block 内 256 个线程各读 4 次；accumulator acc[4][4]、a[4]、b[4] 是线程私有数据，放在寄存器里——它们承载了全部复用收益，也决定了寄存器压力（16 + 4 + 4 = 24 个浮点寄存器，加索引寻址约 40 个）。
• Work / Span / 并行度：
  • Work = M·N·K 次 FMA = 2N³ FLOP（N=1024 时为 2.147 GFLOP）。分块不改变 Work，只改变访存量。
  • Span（关键路径）：每个输出元素的累加链在 k 上是串行的，长度 = K 次 FMA（分散在 K/BK 次外层迭代里，每次迭代间还夹两次 __syncthreads()）→ Span ≈ K·(FMA 延迟) + (K/BK)·2·(屏障延迟)。取 FMA 延迟 4 周期、屏障 ~30 周期、K=1024、BK=16：1024×4 + 64×2×30 = 4096 + 3840 ≈ 7936 周期 ≈ 6.4 µs @1.245 GHz。
  • 并行度 = Work / Span ≈ 2.147e9 FLOP / (7936 周期 × 32 lane-FMA/周期…)。用更直接的算法：并行度 = 可并行的工作量 / 关键路径工作量 = 2N³ / (2K) = N² = 1.05e6 个独立 FMA 链（每个输出元素一条链），远大于硬件并发线程数（V100：163840）→ 并行度充足，瓶颈不在并行度。
• 瓶颈诊断（N=1024，V100 参数：12.7 TFLOPs、900 GB/s、共享带宽 ≈ 12.75 TB/s）：
  • 计算下界：2.147e9 / 12.75e12 = 0.168 ms
  • global（至 L2）流量：2N³/L × 4 B = 2.147e9/64 × 4 = 134 MB → 0.149 ms
  • 共享内存流量：2N³/V × 4 B = 2.147e9/4 × 4 = 2.147 GB → 0.168 ms
  • 三者几乎相等（0.149 / 0.168 / 0.168 ms），说明 L=64, V=4 正好把 kernel 推到计算受限的平衡点；再增大 L 收益递减（global 已经不再是限制），再增大 V 会把寄存器吃光、occupancy 掉下来。
  • 若把 V 降为 2：共享流量翻倍到 0.337 ms，kernel 立刻变成共享内存带宽受限；若把 L 降为 32 而 V=4：global 流量翻倍到 0.298 ms，变成 L2/HBM 受限。这就是”两个旋钮各管一层”的实测含义。
  • 另需注意：N=1024 时 A、B 各 4 MB，合起来接近 V100 的 6 MB L2，因此真实 DRAM 流量远小于上面按”每次 block 都回 HBM”计算的 134 MB；换句话说 L2 抹掉了一部分 global 流量，但L2 带宽本身成为新的上限，结论（继续增大 L 收益有限）不变。

3.2 示例 B：两阶段并行归约（Case Study 2，完整可运行）

// reduce.cu ——  大数组求和: 阶段1 (多 block + 共享内存树 + warp shuffle) + 阶段2 (单 block 收尾)
// 编译(发布优化): nvcc -O3 -arch=sm_80 reduce.cu -o reduce
// 运行:          ./reduce
// 关键点: CUDA 无全局屏障 => 用 kernel 分解代替; 顺序寻址(sequential addressing)保证无 bank 冲突

#include <cuda_runtime.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#define BLK        512
#define FULL_MASK  0xffffffffu

#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(1);                                                          \
        }                                                                     \
    } while (0)

// ---- warp 级归约: 用 shuffle 代替共享内存, 不消耗 shared 带宽也不需要 __syncthreads ----
__device__ __forceinline__ float warp_reduce_sum(float v)
{
#pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        v += __shfl_down_sync(FULL_MASK, v, offset);
    return v;
}

// ---- 阶段 1: 每个 block 用 grid-stride 读取一块数据, 块内两级树形归约 ----
__global__ void reduce_stage1(const float* __restrict__ in,
                              float* __restrict__ partial,
                              long long n)
{
    __shared__ float sdata[BLK / 32];        // 512/32 = 16 个 warp 部分和, 64 B

    const int  tid    = threadIdx.x;
    const long long stride = (long long)gridDim.x * blockDim.x;

    float sum = 0.0f;
    for (long long i = (long long)blockIdx.x * blockDim.x + tid; i < n; i += stride)
        sum += in[i];                        // 连续 lane 读连续地址 => 合并访存

    sum = warp_reduce_sum(sum);              // 第一级: warp 内 shuffle 归约
    if ((tid & 31) == 0) sdata[tid >> 5] = sum;
    __syncthreads();                         // 第二级前必须同步: 块内无 lockstep 保证

    if (tid < 32) {                          // 只让 0 号 warp 参与(整个 warp 都在, mask 合法)
        float v = (tid < BLK / 32) ? sdata[tid] : 0.0f;
        v = warp_reduce_sum(v);
        if (tid == 0) partial[blockIdx.x] = v;
    }
}

// ---- 阶段 2: 单 block 把 nblocks 个部分和收尾 ----
__global__ void reduce_stage2(const float* __restrict__ partial,
                              float* __restrict__ out,
                              int nblocks)
{
    __shared__ float sdata[BLK];

    const int tid = threadIdx.x;
    float sum = 0.0f;
    for (int i = tid; i < nblocks; i += blockDim.x)
        sum += partial[i];
    sdata[tid] = sum;
    __syncthreads();

    // 顺序寻址(sequential addressing): 相邻 lane 访问相邻元素 -> 无 bank 冲突 + 完全合并
    for (int s = blockDim.x >> 1; s > 0; s >>= 1) {
        if (tid < s) sdata[tid] += sdata[tid + s];
        __syncthreads();
    }
    if (tid == 0) out[0] = sdata[0];
}

int main()
{
    const long long n = 1LL << 24;                       // 16,777,216 个 float = 64 MiB
    const int gridsz  = 4096;                            // 4096 个 block x 512 线程
    const size_t bytes = (size_t)n * sizeof(float);

    float* h = (float*)malloc(bytes);
    double exact = 0.0;
    for (long long i = 0; i < n; ++i) {
        h[i] = (float)((i % 1024) * 1e-3);               // 确定性数据, 便于校验
        exact += (double)h[i];
    }

    float *d_in, *d_partial, *d_out;
    CUDA_CHECK(cudaMalloc(&d_in, bytes));
    CUDA_CHECK(cudaMalloc(&d_partial, (size_t)gridsz * sizeof(float)));
    CUDA_CHECK(cudaMalloc(&d_out, sizeof(float)));
    CUDA_CHECK(cudaMemcpy(d_in, h, bytes, cudaMemcpyHostToDevice));

    // 预热
    reduce_stage1<<<gridsz, BLK>>>(d_in, d_partial, n);
    reduce_stage2<<<1, BLK>>>(d_partial, d_out, gridsz);
    CUDA_CHECK(cudaDeviceSynchronize());

    cudaEvent_t e0, e1;
    CUDA_CHECK(cudaEventCreate(&e0));
    CUDA_CHECK(cudaEventCreate(&e1));
    const int iters = 50;
    CUDA_CHECK(cudaEventRecord(e0));
    for (int it = 0; it < iters; ++it) {
        reduce_stage1<<<gridsz, BLK>>>(d_in, d_partial, n);
        reduce_stage2<<<1, BLK>>>(d_partial, d_out, gridsz);
    }
    CUDA_CHECK(cudaEventRecord(e1));
    CUDA_CHECK(cudaEventSynchronize(e1));

    float ms = 0.f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
    ms /= iters;

    float gpu_sum = 0.f;
    CUDA_CHECK(cudaMemcpy(&gpu_sum, d_out, sizeof(float), cudaMemcpyDeviceToHost));

    const double gib = (double)bytes / 1e9;
    printf("[verify] GPU=%.4f  CPU(double)=%.4f  相对误差=%.3e\n",
           gpu_sum, exact, fabs((double)gpu_sum - exact) / exact);
    printf("[perf  ] 耗时 = %.4f ms ; 读 %.2f MB -> 有效带宽 = %.1f GB/s\n",
           ms, bytes / 1e6, bytes / (ms * 1e-3) / 1e9);
    printf("[model ] 12.7 TFLOPs 下纯计算下界 = %.4f ms ; 900 GB/s 下读带宽下界 = %.4f ms\n",
           (double)(n - 1) / 12.75e12 * 1e3, (double)bytes / 900e9 * 1e3);
    (void)gib;

    cudaFree(d_in); cudaFree(d_partial); cudaFree(d_out);
    free(h);
    return 0;
}

【代码做什么？】

1. host 侧：造一个有确定性的 n = 2^24 浮点数组（64 MiB），同时在 CPU 上用 double 累加出精确值作参考；分配 d_in、d_partial[4096]、d_out[1]；预热后连跑 50 次两阶段流程并计时。
2. 阶段 1 kernel：4096 个 block × 512 线程 = 2,097,152 个 CUDA 线程（映射为 65536 个 warp）。每个线程用 grid-stride 循环读 n / (grid×block) = 8 个元素并累加到私有寄存器 sum；然后两级树形归约：先在 warp 内用 __shfl_down_sync 把 32 个 lane 的值折成 1 个（无需共享内存、无需屏障），再把 16 个 warp 的部分和放进 64 B 的 sdata，__syncthreads() 后由 0 号 warp 收尾，写 partial[blockIdx.x]。
3. 阶段 2 kernel：只有一个 block（这就是”用 kernel 分解替代全局同步”的落点：阶段 1 结束是一次隐式全局屏障），512 个线程把 4096 个部分和读进来放进共享内存，用顺序寻址做 log2(512) = 9 轮树形归约，最后由 tid == 0 写出总和的唯一结果。

【并行机制与性能解说】

• 线程/向量/warp 的创建与工作分配：launch 语义是”批量创建”——一次 launch 声明 gridsz × BLK 个线程；硬件把它们按 32 个一组打包成 warp（每 block 16 个 warp），work scheduler 按资源把 block 派到 SM；每个 SM 每时钟最多选 4 个 warp 发射、每 warp 最多选 2 条指令（线程级 + 指令级并行）。
• 共享数据的处理：sdata 是块内共享的中间结果（16 个 float 或 512 个 float）；块间通信只能经 global memory 的 partial 数组，并用 kernel 边界保证可见性/顺序。warp 内归约刻意不走共享内存——shuffle 直接在寄存器间搬数据，省掉 bank 访问与 5 次屏障。
• Work / Span / 并行度：
  • Work = n - 1 = 16,777,215 次加法 ≈ 1.68e7 FLOP（阶段 2 的 4096+511 次可忽略）。这是work-optimal 的：树形归约总加法数 = n-1，与串行相同。
  • Span = 阶段 1 内：(每线程 8 次串行累加) + (warp 内 5 步 shuffle) + (1 次 __syncthreads) + (warp 0 的 5 步 shuffle) ≈ 8·4 + 5·(shuffle 延迟 ~10) + 30 + 50 ≈ 162 周期，再加阶段 2 的 grid-stride 8 次 + 9 轮 (共享读+加+屏障) ≈ 8·20 + 9·(30+30) ≈ 700 周期，以及一次 kernel 边界（约 1–5 µs 量级）。以周期计，Span ≈ 900 周期 ≈ 0.72 µs（不含 launch 开销）。
  • 并行度 = Work / Span ≈ 1.68e7 / 900 ≈ 1.9e4（按”加法个数 / 关键路径加法个数”的等价估算则是 n / log n ≈ 1.0e6）。无论哪种算法，并行度都远大于硬件能同时运行的线程数（V100：163,840），所以”并行度不足”不是问题。
• 瓶颈诊断（数值见程序输出）：
  • 计算下界 = 1.68e7 / 12.75e12 = 1.3 µs；
  • 访存下界 = 64 MiB / 900 GB/s ≈ 74.6 µs；
  • 两者相差 57 倍 → 这个 kernel 是纯带宽受限，优化目标是”把 HBM 流量压到 1 遍读 + 1 遍写”，任何额外的中间数组往返都是纯亏损。
  • 反面教材（recitation 的 Idea 1）：如果把树的每一层都做成一次 launch（24 层 = 24 次 launch），即使每层 kernel 本身很快，24 × 约 5 µs ≈ 120 µs 的 launch 开销就已超过 74.6 µs 的数据搬运下界；再加上每层都要把中间数组写回 global 再读回，实际时间会到 200 µs 以上。这就是”正确但很慢”的典型。
  • 共享内存与屏障成本：顺序寻址让每轮共享访问都落在不同 bank（无冲突，1 周期），但 __syncthreads() 本身有 ~几十周期的延迟，9 轮 × 2 次同步（阶段 2）在只有 1 个 block 时无法被其他 block 掩盖——好在阶段 2 的数据量只有 4096 个 float，占比极小。
  • 负载不均：n = 2^24 除以 4096×512 = 2^21 正好整除（每线程 8 个元素），无尾部；若 n 不是整除，grid-stride 会让前几个 block 多干一轮，此时应让grid 大小略小于”SM × 每 SM 驻留 block 数”的整数倍并依赖 grid-stride 循环平衡负载。

3.3 示例 C：共享内存原子操作实现直方图（原子性与调度自由的边界）

// histogram.cu ——  用 shared memory 原子操作建直方图, 再原子汇总到 global
// 编译(发布优化): nvcc -O3 -arch=sm_80 histogram.cu -o histogram
// 运行:          ./histogram
// 要点: (1) 块内原子化到共享内存 -> 大幅减少 global 原子竞争
//       (2) 这个用法合法: 原子操作只用于互斥, 不约束 block 的调度顺序

#include <cuda_runtime.h>
#include <cstdio>
#include <cstdlib>

#define NBINS 256
#define BLK   256

#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(1);                                                          \
        }                                                                     \
    } while (0)

__global__ void histogram_shared(const unsigned char* __restrict__ data,
                                 long long n,
                                 unsigned int* __restrict__ hist)
{
    __shared__ unsigned int sHist[NBINS];        // 256 x 4 B = 1 KB

    for (int i = threadIdx.x; i < NBINS; i += blockDim.x)
        sHist[i] = 0u;
    __syncthreads();                             // 清零对所有线程可见后才能开始计数

    const long long stride = (long long)gridDim.x * blockDim.x;
    for (long long i = (long long)blockIdx.x * blockDim.x + threadIdx.x; i < n; i += stride) {
        atomicAdd(&sHist[data[i]], 1u);          // 块内竞争发生在片上, 延迟远低于 global
    }
    __syncthreads();                             // 保证所有计数完成后再汇总

    for (int i = threadIdx.x; i < NBINS; i += blockDim.x)
        atomicAdd(&hist[i], sHist[i]);           // 每个 block 对每个 bin 只做一次 global 原子加
}

int main()
{
    const long long n = 1LL << 26;                       // 64 Mi 个字节值
    unsigned char* h = (unsigned char*)malloc((size_t)n);
    unsigned int ref[NBINS] = {0};
    for (long long i = 0; i < n; ++i) {
        h[i] = (unsigned char)((i * 31 + 7) % NBINS);    // 确定性数据
        ref[h[i]]++;
    }

    unsigned char* d_data;
    unsigned int*  d_hist;
    CUDA_CHECK(cudaMalloc(&d_data, (size_t)n));
    CUDA_CHECK(cudaMalloc(&d_hist, NBINS * sizeof(unsigned int)));
    CUDA_CHECK(cudaMemcpy(d_data, h, (size_t)n, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemset(d_hist, 0, NBINS * sizeof(unsigned int)));

    const int gridsz = 1024;
    histogram_shared<<<gridsz, BLK>>>(d_data, n, d_hist);
    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaDeviceSynchronize());

    unsigned int out[NBINS];
    CUDA_CHECK(cudaMemcpy(out, d_hist, sizeof(out), cudaMemcpyDeviceToHost));

    long long bad = 0;
    for (int i = 0; i < NBINS; ++i) if (out[i] != ref[i]) bad++;
    printf("[verify] 不匹配的 bin 数 = %lld / %d ;  bin0: GPU=%u CPU=%u\n", bad, NBINS, out[0], ref[0]);

    cudaEvent_t e0, e1;
    CUDA_CHECK(cudaEventCreate(&e0));
    CUDA_CHECK(cudaEventCreate(&e1));
    const int iters = 20;
    CUDA_CHECK(cudaEventRecord(e0));
    for (int it = 0; it < iters; ++it)
        histogram_shared<<<gridsz, BLK>>>(d_data, n, d_hist);
    CUDA_CHECK(cudaEventRecord(e1));
    CUDA_CHECK(cudaEventSynchronize(e1));
    float ms = 0.f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
    ms /= iters;

    printf("[perf  ] %.4f ms ; 读 %.2f MB -> %.1f GB/s\n",
           ms, (double)n / 1e6, (double)n / (ms * 1e-3) / 1e9);

    cudaFree(d_data); cudaFree(d_hist); free(h);
    return 0;
}

【代码做什么？】 每个 block 先在共享内存里放 256 个计数器（1 KB），清零并同步；然后用 grid-stride 循环读数据，对 sHist[value] 做块内原子加；最后每个 block 只对每个 bin 做一次 global atomicAdd。结果是：global 原子操作次数从 n 次（64 Mi 次）降到 blocks × NBINS = 1024 × 256 = 262144 次。

【并行机制与性能解说】

• 原子操作与调度自由：这是一个合法的 CUDA 用法——原子操作只充当互斥，无论 block 以何种顺序、是否并发执行，结果都相同，因此不约束 work scheduler 的调度自由。反例：让 block 1 自旋等待 block 0 写的标志位（while(atomicAdd(&flag,0)==0){}），在”只能驻留 1 个 block 的 SM”上若先跑 block 1 就永久死锁——CUDA 只承诺”块可以任意顺序执行”，从不承诺”块之间没有依赖”能被满足。
• Work / Span / 并行度：
  • Work = n 次原子加（外加每个 block 的清零与汇总，各 O(NBINS)）→ 6.7e7 次原子操作。
  • Span = 最坏情况是”所有元素都落在同一个 bin”→ 该 bin 上的 n 次原子加被串行化，Span = Θ(n)；此时并行度 = Work/Span ≈ 1，整个 kernel 退化为串行。实测数据接近均匀分布，因此真实的 Span 由”每个 bin 上的冲突次数”决定 ≈ n / NBINS（每 bin 约 26 万次串行原子）。这就是基数/冲突决定的批判路径：并行度 = 元素数 / 最大 bin 的冲突数。
  • 并行度 = Work / Span ≈ NBINS 量级（均匀分布时），远小于硬件线程数 → 这个 kernel 的性能上限由”最热 bin 的原子串行化”决定，而不是由线程数决定。缓解手段：让每个 warp/block 拥有私有副本（本示例的共享内存版本）、或把 bin 数减少到能被 warp 内 shuffle 处理。
• 瓶颈与实测预期：读 64 MiB 的带宽下界 = 67.1e6/900e9 = 74.6 µs；一旦共享内存上的原子冲突严重（随机 bin → 多次重试、bank 争用），实测会明显高于此值。注意 sHist[NBINS] 是 256 个 4 B 字，正好每 bank 8 个字，随机 bin 会让同一 warp 的多个 lane 撞到同一 bank，这是”共享内存原子争用”与”bank 冲突”叠加的典型案例。

3.4 示例 D：CPU 侧 OpenMP 参考实现（Work/Span 对照与 roofline 定位）

// reduce_omp.cpp ——  同一归约任务的 CPU 多核版本, 用来做"机器平衡点"对照
// 编译(发布优化): g++ -O3 -fopenmp -march=native reduce_omp.cpp -o reduce_omp
// 运行:          OMP_NUM_THREADS=16 ./reduce_omp

#include <omp.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>

int main()
{
    const long long n = 1LL << 24;                 // 与 CUDA 版同为 16,777,216 个 float (64 MiB)
    float* a = (float*)aligned_alloc(64, (size_t)n * sizeof(float));
    double exact = 0.0;
    for (long long i = 0; i < n; ++i) { a[i] = (float)((i % 1024) * 1e-3); exact += (double)a[i]; }

    double sum = 0.0, t0 = 0.0, t1 = 0.0;
    const int iters = 20;

    // 预热
    sum = 0.0;
#pragma omp parallel for reduction(+ : sum) schedule(static)
    for (long long i = 0; i < n; ++i) sum += (double)a[i];

    sum = 0.0;
    t0 = omp_get_wtime();
    for (int it = 0; it < iters; ++it) {
        double s = 0.0;
#pragma omp parallel for reduction(+ : s) schedule(static)
        for (long long i = 0; i < n; ++i) s += (double)a[i];
        sum = s;                                   // 避免被优化掉
    }
    t1 = omp_get_wtime();

    const double ms = (t1 - t0) / iters * 1e3;
    printf("[verify] OMP=%.4f  CPU(double)=%.4f  相对误差=%.3e\n",
           sum, exact, fabs(sum - exact) / exact);
    printf("[perf  ] 线程数=%d  耗时=%.4f ms  读 %.2f MB -> 有效带宽=%.2f GB/s\n",
           omp_get_max_threads(), ms, (double)n * 4 / 1e6,
           (double)n * 4 / (ms * 1e-3) / 1e9);
    printf("[expect] 64 MiB / 20 GB/s(单路 DDR) = %.2f ms ;  / 100 GB/s(双路服务器) = %.2f ms\n",
           (double)n * 4 / 20e9 * 1e3, (double)n * 4 / 100e9 * 1e3);
    free(a);
    return 0;
}

【代码做什么？】 用 #pragma omp parallel for reduction(+:s) schedule(static) 把同一归约任务分给 16 个线程，每个线程把自己那块连续区间（schedule(static) 意味着每个线程拿到一大段连续地址）累加进私有 s，最后由 OpenMP 运行时的归约把 16 个 s 合并（内部就是一棵小树）。计时用 omp_get_wtime()，跑 20 次取平均。

【并行机制与性能解说】

• 并行机制：OpenMP 在进入并行区时唤醒一个线程池（线程数默认 = 核数/超线程数），用 fork-join 模式把循环的迭代块静态切给各线程；每个线程的 s 是线程私有变量（不在同一 cache line 上，避免伪共享——OpenMP reduction 的实现会做 padding）。这与 GPU 的”grid-stride + 树形归约”在结构上是同构的：私有部分和 → 树形合并。
• Work / Span / 并行度：Work = n−1 次加法；Span = n/P（每线程串行段） + log₂P（归约树）；并行度 = (n-1) / (n/P + log₂P) ≈ P = 16。这正是 CPU 与 GPU 的分野：CPU 版并行度被核数（16）限制，GPU 版并行度可达 10⁶ 量级——但因为任务本身带宽受限，两者都会撞在同一堵墙（DRAM 带宽）上。
• 瓶颈诊断：单路 DDR5 约 20 GB/s 时，64 MiB / 20 GB/s = 3.35 ms；双路服务器约 100 GB/s 时是 0.67 ms。16 个核再多也压不穿内存墙——这直接解释了为什么 GPU 要用 HBM（900 GB/s ~ 1 TB/s）：只有当带宽提高 10–50 倍，算力才有意义。

4. 性能模型与复杂度分析

4.1 Work-Span 在 GPU 上的映射

• Work（W） = 程序的总操作数；Span（D） = 最长依赖链；并行度 = W/D。
• 在 GPU 上，”并行度”要与两个硬件常量比较：
  • 可驻留线程总数 = SM 数 × 每 SM warp 数 × 32。例如 V100：80 × 64 × 32 = 163,840；H100：132 × 64 × 32 = 270,336。
  • 每时钟可发射的 warp 指令数 = SM 数 × 4 sub-core（V100：80 × 4 = 320 条 warp 指令/时钟）。
• 判据：若 并行度 ≫ 可驻留线程数，则”并行度”不是瓶颈，应转而分析带宽 / 算术强度 / 同步开销；若 并行度 < 可驻留线程数，则 GPU 会大量空转（这正是直方图原子化、稀疏不规则算法的病根）。
• Roofline（屋顶线）的三条线：
  1. 计算屋顶 = 峰值算力 P_peak；
  2. 带宽屋顶 = BW × I（I = 算术强度 FLOP/Byte）；
  3. 机器平衡点 I* = P_peak / BW——低于它 = 带宽受限，高于它 = 计算受限。

4.2 算例一：机器平衡点与算术强度阶梯（matmul，N=1024）

假设参数（取其自讲义与 CS149 补充材料的 V100 数据）：

P_peak(fp32) = 80 SM × 4 sub-core × 16 lane × 1.245 GHz × 2 FLOP/FMA = 12.75 TFLOP/s
BW_HBM                                                        = 900 GB/s
BW_shared = 80 SM × 32 bank × 4 B × 1.245 GHz                 = 12.75 TB/s
机器平衡点 I* = 12.75e12 / 900e9                              = 14.2 FLOP/Byte

算例：N = 1024 的稠密矩阵乘，总计算量 2N³ = 2.147 GFLOP，计算下界：

T_compute = 2.147e9 / 12.75e12 = 0.168 ms

各版本的访存量与时间下界：

同一张表换到共享内存维度看（共享带宽 12.75 TB/s）：

结论：T ≥ max(0.168, 9.543/L [ms], 0.674/V [ms])。取 L=64, V=4 时三项分别是 0.168 / 0.149 / 0.168 ms，三堵墙同时被顶到——这就是”两个复用旋钮调到位”的定量标志。继续调参只会让另一堵墙变得更矮：性能工程就是找到这几堵墙的交点。

4.3 算例二：1D 卷积的带宽下界（共享内存暂存的收益）

假设：N = 2²⁰ = 1,048,576 个输出，output[i] = (input[i]+input[i+1]+input[i+2])/3，THREADS_PER_BLK = 128，每个 block 需要 128+2 = 130 个输入。

计算量对照：3N ≈ 3.15 MFLOP（2 次加 + 1 次除，除算 1 FLOP）→ 3.15e6/12.75e12 = 0.25 µs。计算时间只占访存时间的 1.3%，是极端带宽受限的例子；同时要注意版本 2 引入了一次 __syncthreads()（约几十周期）与 520 B/block 的共享内存占用，若每线程工作量太小（本算例每线程只算 1 个输出），屏障与 block 启动开销会开始显著。

再看另一个更朴素的对照（把 CPU 侧数字放进同一个模型）：

CPU: 16 核 × 3.0 GHz × 8 宽 SIMD × 2 (FMA)              = 768 GFLOP/s
     单路 DDR5 带宽                                     ≈   20 GB/s
     机器平衡点 = 768e9 / 20e9                           = 38.4 FLOP/Byte

     一个 256 MB 的数组, 单次遍历(读+写各 256 MB)在 20 GB/s 下至少需
     512 MB / 20 GB/s = 25.6 ms      （只读一次则 256 MB / 20 GB/s = 12.8 ms）
     同期 GPU: 512 MB / 900 GB/s     = 0.57 ms          -> 快 45 倍

这解释了为什么”带宽”在 GPU 编程里比”算力”更常成为瓶颈：GPU 的机器平衡点（14.2）比 CPU 的（38.4）更低，因为它的带宽相对算力更充裕，但也更容易被”多读几遍数据”的写法挥霍掉。

4.4 算例三：归约的 Work/Span 与”kernel 分解层数”的代价

假设：n = 2²⁴ = 16,777,216 个 float = 64 MiB，V100 参数，BLK = 512，grid = 4096。

Work  = n - 1        = 16,777,215 次加法  = 1.68e7 FLOP
Span  ≈ (n/线程数 的串行段) + (log2(32) shuffle) + (log2(512/BLK_warp) 树) + kernel 边界
      = 8 次串行加 + 5 步 + 5 步 + 1 次 launch 边界   ≈ 900 周期 ≈ 0.72 µs
可驻留线程数 = 80 SM × 64 warp × 32 = 163,840
同时可执行的独立部分和 = 4096 block × 512 线程 = 2,097,152

一个常被误用的比较：若直接用 Work/Span = 1.68e7 / 900 ≈ 1.9e4 与 163,840 相比，会得到”并行度不够”的错误结论——因为这里的 Span 以周期为单位、Work 以操作数为单位，两者不同量纲，不能直接相比。正确的比较是同量纲的”同时可执行的独立工作单元数”：4096 × 512 = 2,097,152 个独立部分和，远大于 163,840 个可驻留线程，因此并行度充足。真正决定性能的是带宽：

T_bandwidth = 64 MiB / 900 GB/s = 74.6 µs      <- 实际能达到的最好成绩
T_compute   = 1.68e7 / 12.75e12 = 1.3 µs
比值        = 57×                              <- 纯带宽受限

kernel 分解层数的代价（recitation 的四种归约方案对照）：

算一笔”层数 vs 带宽”的账：Idea 1 在 24 层里，第 i 层读写 n/2^i 个元素，总流量 ≈ 2n(1+1/2+1/4+…) ≈ 4n 个元素 = 4 × 64 MiB = 256 MiB → 单是流量就 284 µs，再加 120 µs 的 launch 开销 → 约 400 µs，是理想值的 5 倍以上。而示例 B 的流量只有 n 个元素（读） + 4096 个部分和（写）≈ 64 MiB → 74.6 µs。同一个算法、同一个 Work/Span，两种分解方式差 5 倍。

4.5 算例四：延迟隐藏的定量要求（Little’s Law）

内存延迟不会因为”有很多线程”自动消失，它必须被并发请求数掩盖：

Little 定律:  并发请求数(in flight) = 带宽 × 延迟
假设: HBM 往返延迟 ≈ 500 周期, V100 时钟 1.245 GHz, 带宽 900 GB/s
  延迟(秒) = 500 / 1.245e9 = 4.016e-7 s
  需要 in-flight 字节 = 900e9 × 4.016e-7 = 361 KB
  摊到 80 个 SM: 约 4.5 KB / SM 必须同时在路上

每个 warp 的一次合并访存 = 32 lane × 4 B = 128 B
  => 每个 SM 需要约 4.5 KB / 128 B ≈ 36 条 warp 访存同时在途
  => 该 SM 至少要能驻留并保持发射 ~36 个 warp（V100 上限 64 个 warp）

结论：occupancy（驻留 warp 数）是”延迟隐藏能力”的代理指标。这也解释了寄存器分块的代价：V100 每 SM 只有 65536 个寄存器；若一个 kernel 用 40 个寄存器/线程，则最多驻留 65536/40 = 1638 个线程 = 51 个 warp（而不是 64），occupancy ≈ 80%；若 V=8 的 matmul 用到 80+ 个寄存器，驻留线程降到 819 = 25 个 warp，可能低于隐藏 500 周期延迟所需的 36 个 warp，于是”减少了访存却跑得更慢”。

一个具体的自检公式：

某 kernel 每线程寄存器的编译结果 R  ->  可驻留 warp 数 W = 65536 / (R × 32)
若 W < 需要掩盖的并发访存数(≈ 带宽×延迟 / (SM数 × 128 B))  ->  带宽用不满

4.6 Amdahl 与”kernel 边界”的隐含税

设串行部分（host 侧准备 + kernel 启动 + 结果回收）占比为 s，则加速比 S = 1/(s + (1-s)/P)。在 GPU 上，s 的常见来源是每次 launch 的固定开销与host↔device 拷贝：

假设 launch 开销 5 µs, 数据搬运 8.45 MB / 900 GB/s = 9.4 µs (卷积例)
若把算法拆成 24 次 launch: 固定开销 = 120 µs
   相对"数据搬运 9.4 µs" 的串行占比 s = 120/(120+9.4) = 92.7%
   => 无论 GPU 有多少个 SM, 加速比上限 = 1/s ≈ 1.08 !

这就是”kernel 分解要少而大“的定量理由：能用 2 次 launch 解决，就绝不用 24 次。

5. 关键要点

1. 三条判据决定 CUDA kernel 的成败：coherent warps（warp 内不发散）、coalesced memory access（合并访存）、高数据复用（register tiling + shared memory tiling）。它们不是三个独立技巧，而是同一个目标——让每个时钟的每个 lane 都在做有用的事，并且这些 lane 需要的字节已经被搬到片上。
2. 两个复用旋钮分管两层存储：L（block 级共享内存瓦片）把 global 流量降到 2N³/L，V（线程级寄存器瓦片）把 shared 流量降到 2N³/V。性能上界 T ≥ max(T_compute, 2N³·4/(L·BW_HBM), 2N³·4/(V·BW_shared))；调到三堵墙等高时就是该硬件上的最优分块。
3. CUDA 没有全局屏障，跨块同步只能靠 kernel 分解；而每次 launch 都有固定开销，所以分解要”少而大”。任何”块间自旋等待”的代码在”驻留 block 数 < block 总数”的机器上都可能死锁——因为 block 可以被按任意顺序执行。
4. 优化的第一步永远是定位瓶颈类型：先算算术强度 I 与机器平衡点 I* = P_peak/BW；I < I* 就去减少字节搬动（分块、合并访存、缩短数据类型），I > I* 才去优化指令数（向量化、减少发散、用 Tensor Core）。在 I < I* 时优化计算是纯粹的浪费。
5. 访存模式比访存次数同样重要：32 lane × 4 B 的连续访问 = 1 次全宽事务（最优）；stride 访问会被拆成多次事务（有效带宽除以拆分数）；共享内存 stride-2 访问产生 2 路 bank 冲突（带宽减半）。顺序寻址（sequential addressing）同时解决合并访存与 bank 冲突，是归约类 kernel 的标准收尾。
6. 延迟必须被并发掩盖：由 in-flight 字节 = 带宽 × 延迟 可算出所需并发访存数（V100 约 36 条 warp 访存/SM），这直接约束了寄存器分块的规模上限——occupancy 与复用率是一对必须权衡的矛盾。

6. 常见陷阱与注意事项

• 把 __syncthreads() 当作跨 block 的屏障：它只对同一 block 内的线程有效。recitation 里”在 kernel 里直接加 __syncthreads() 就以为能跨块同步”的写法是错误的（版本 2 的错误原因）；跨块只能靠 kernel 边界，或换用 cooperative groups 的 grid sync（需要特殊的 launch 方式与驻留保证）。
• 在 block 之间做自旋等待（spin-wait）：while (atomicAdd(&flag,0)==0) {} 这类代码在”SM 只能驻留 1 个 block 且先跑了等待方”时永久死锁。CUDA 只保证”块可以任意顺序执行”，从不保证”并发执行”或”依赖能被满足”。用原子操作做计数/直方图是合法的（只用于互斥），用原子操作做块间握手是不合法的。
• 忽略 warp 内的发散代价：if (tid % (2*s) == 0) 这种写法让一半 lane 在第一步就空转、后续步骤活跃率按 1/4、1/8 递减（示例中对 16 个 lane 的三个步骤有效利用率只有 (8+4+2)/(16×3) = 29%）。改成”strided index + 非发散分支”（index = 2*s*threadIdx.x; if (index < blockDim.x)）或顺序寻址，可以让前面的整 warp 满负荷工作，只把不可避免的部分留在尾部。
• 把 shared memory 当成”免费”的：共享内存有 32 个 bank、每 bank 每时钟只能服务 1 个字。stride 访问（如 sdata[2*tid]）会产生 bank 冲突使带宽减半；float 数组按 sdata[tid*stride] 访问更会出现严重冲突。解决办法：连续访问、或者把共享内存数组的列数 padding 成奇数（例如 [BLK][BLK+1]）来打散 bank 映射。
• 为了减少 global 访问而无节制地增大 V，结果 occupancy 崩掉：V×V 个 accumulator 全部占用寄存器；V=8 时光 accumulator 就要 64 个寄存器/线程，加上寻址会突破 80，使每 SM 驻留 warp 数从 64 掉到 20 多，延迟无法隐藏，反而变慢。判断依据是 4.5 节的 W = 65536/(R×32) 与 Little 定律给出的并发需求。
• 忽视”写回”也是一次访存，以及忽视 L2 的存在：写 C 矩阵（N² × 4 B）在 N 大时是实打实的 DRAM 流量；而当 N 小到 A、B 能装进 L2（V100 6 MB）时，按”每次 block 都回 HBM”估算的流量会严重高估时间——真实瓶颈会从 HBM 带宽变成 L2 带宽。分析时必须分清”至 L2 的流量”与”至 DRAM 的流量”。
• 把 CUDA 的线程当成 pthread：CUDA 线程没有”创建/销毁”的运行时语义，硬件只保存寄存器状态；但一个 block 的所有线程必须同时获得执行上下文（因为它们按语义就是并发的，且可能有 __syncthreads() 依赖），所以”先跑线程 0-127 再跑 128-255”是不允许的。同样，块内线程与 warp 的关系（连续 32 个 tid 属于同一 warp）、以及”warp 不是编程模型的一部分而是实现细节”这一点，都必须在写代码时记在心里。

7. 思考题（带答案）

问题 1（分块参数的定量选择）：某 GPU 的 fp32 峰值算力 P_peak = 20 TFLOP/s，HBM 带宽 BW = 1 TB/s，共享内存总带宽 BW_shared = 20 TB/s。要算 N = 2048 的稠密矩阵乘，请：(a) 求机器平衡点 I*；(b) 写出以 L（block 瓦片边长）和 V（线程瓦片边长）表示的三个时间下界；(c) 选一组 (L, V)，使得三项大致平衡，并说明寄存器约束；(d) 若有人把 V 从 4 提到 8，会出现什么问题？

【答案】 (a) I* = P_peak / BW = 20e12 / 1e12 = 20 FLOP/Byte。低于 20 FLOP/B 的写法就是带宽受限。

(b) 总计算量 2N³ = 2×2048³ = 1.718e10 FLOP。

• 计算下界：T_comp = 1.718e10 / 20e12 = 0.859 ms
• global 下界：总访问次数 2N³/L，字节数 8N³/L = 6.87e10 / L B → T_global = 6.87e10/(L × 1e12) s = 68.7/L ms
• shared 下界：总访问次数 2N³/V，字节数 8N³/V → T_shared = 6.87e10/(V × 20e12) s = 3.43/V ms

(c) 令三者相等：

• 68.7/L = 0.859 → L ≈ 80（取 64 或 128，实践中取 64 或 128 这样的 2 的幂）
• 3.43/V = 0.859 → V ≈ 4 取 L = 64, V = 4：T_comp = 0.859 ms，T_global = 68.7/64 = 1.073 ms，T_shared = 3.43/4 = 0.858 ms → 上界 max = 1.073 ms，与计算下界同量级（差 25%），已经接近最优。若想再压 global，取 L = 128, V = 4：T_global = 0.537 ms < 0.859 ms，此时 T = max(0.859, 0.537, 0.858) = 0.859 ms，正好卡在计算屋顶上——这就是”三堵墙等高”的解。 寄存器约束：acc[V][V] = 16 个浮点寄存器 + a[4] + b[4] = 8 个 + 索引/指针约 12–16 个 ≈ 40 个寄存器/线程。共享内存 / block：(BM×BK + BK×BN) × 4 B，若 BK=16、L=128 则为 (128×16 + 16×128)×4 = 16 KB——需与 SM 的共享内存容量（H100 为 256 KB）和寄存器文件（若 64K 寄存器/SM，40 reg/thread → 每 SM 最多 1638 线程 ≈ 51 warp，occupancy ≈ 80%）一起权衡。

(d) V 从 4 提到 8：acc 需要 64 个寄存器，加 a[8]+b[8] 与寻址将超过 80–90 个寄存器/线程。后果有两层：(1) 寄存器溢出（register spilling），多出来的 accumulator 被放到 local memory（本质是 global memory），访存量反而暴涨；(2) 即使不溢出，W = 65536/(R×32) 会让可驻留 warp 数从 ~51 掉到 ~22，低于隐藏内存延迟所需的并发 warp 数，于是带宽用不满、kernel 变慢。而 T_shared 从 0.858 ms 降到 0.429 ms 的收益此时完全用不上，因为 T_comp = 0.859 ms 才是屋顶。结论：V 不是越大越好，它要停在”共享内存不再是最慢的那堵墙”的位置。

问题 2（归约的正确性与性能）：下面这段 kernel 想把长度为 n 的数组求和（g_odata[blockIdx.x] 存每个 block 的部分和），但它既可能算错、又慢。请指出两处正确性问题与两处性能问题，并给出修改方案。

__global__ void reduce_bad(const float* g_idata, float* g_odata, int n) {
    extern __shared__ float sdata[];
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    sdata[tid] = g_idata[i];
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0)
            sdata[tid] += sdata[tid + s];
    }
    g_odata[blockIdx.x] = sdata[0];
}

【答案】

正确性问题 1：缺少 __syncthreads()。 CUDA 不保证 block 内线程 lockstep 执行，而在同一轮里，sdata[tid+s] 可能正被另一个线程写入（第 s 轮里 tid 读取的 tid+s 正是另一个线程的位置），下一轮又要读本轮的结果。没有屏障，就会出现读写竞争，结果不确定。修改：在进入循环前加一次 __syncthreads()（保证装载完成），并在循环体末尾每轮加一次 __syncthreads()。注意：循环体内的屏障是必需的，且必须放在”所有线程都会执行到”的位置——如果把它放进 if (tid % (2*s) == 0) 里，就会因为分支内屏障而使 warp 死锁。

正确性问题 2：越界访问（数组边界）。 当 blockDim.x 不是 2 的幂、或 n 不是 blockDim.x 的整数倍时，最后一个 block 的线程会读 g_idata[i] 越界（sdata[tid+s] 在 tid+s ≥ blockDim.x 时也越界）。修改：装载时加 sdata[tid] = (i < n) ? g_idata[i] : 0.f;，并保证 blockDim.x 为 2 的幂（或把归约循环条件收紧）。

性能问题 1：warp 高度发散。 if (tid % (2*s) == 0) 的第一个 warp（lane 0–31）在 s=1 时只有偶数 lane 活跃（有效利用率 1/2），s=2 时 1/4，s=4 时 1/8……整个前几轮的 lane 利用率是 (16+8+4+2+1)/(32×5) = 31/160 ≈ 19%。修改：改成顺序寻址（reversed loop）：

for (unsigned int s = blockDim.x / 2; s > 0; s /= 2) {
    if (tid < s) sdata[tid] += sdata[tid + s];
    __syncthreads();
}

这样第一轮有 blockDim.x/2 个线程活跃（多个完整 warp 满负荷），只有最后 5 轮落在单个 warp 内，且那 5 轮的分支是”前缀活跃”，代价最小。

性能问题 2：共享内存 bank 冲突 + 非合并访存。 sdata[tid] 与 sdata[tid+s] 在交错寻址下，活跃 lane 的地址是 0, 2s, 4s, ...，相邻活跃 lane 相隔 2s 个字 → 当 s=1 时就是 stride-2 访问，落在偶数 bank 上，产生 2 路 bank 冲突（例如 lane 0–7 访问 word 0,2,4,6,8,10,12,14 → B0 服务 t0,t4，B2 服务 t1,t5 …，4 个 bank 各服务 2 个 lane，另外 4 个 bank 空闲，访问串行化为 2 个周期，带宽减半）。顺序寻址下 sdata[tid]（tid=0..s-1 连续）与 sdata[tid+s]（连续区间）都是连续访问 → 无 bank 冲突。同时 g_idata[i] 是连续 tid 读连续地址，本来就是合并访存——这一点在两版中都成立，要保留。

补充（可选的进一步优化）：把 warp_reduce_sum 用 __shfl_down_sync 做最后 32 个元素的归约（省掉 5 轮共享内存访问与 5 次屏障），并在 blockDim.x >= 64 时先做”每个线程读多个元素”的 grid-stride 循环，把 Work 与访存比提高；最后用一个 block 的第二个 kernel（或 atomicAdd）汇总各 block 的部分和——因为CUDA 没有全局屏障。

问题 3（Roofline 应用与决策）：某 CUDA kernel 处理一个 n = 2^24 的 float 数组（64 MiB），做以下三件事之一：

• (A) 求和（1 次读，1 个标量输出）；
• (B) y[i] = a*x[i] + y[i]（读 x 与 y，写 y）；
• (C) 10 阶多项式求值 y[i] = c0 + x[i]*(c1 + x[i]*(... ))（读 x，写 y，每次 10 个 FMA）。

GPU 参数：P_peak = 12.75 TFLOP/s（fp32 FMA），BW = 900 GB/s，I* = 14.2 FLOP/B。请分别算出算术强度、判断瓶颈类型、给出理论最短时间，并说明应该优先优化什么。

【答案】

统一取 n = 16,777,216 个元素。

(A) 求和

• 流量：读 n × 4 B = 67.1 MB（输出 4 B 可忽略）。
• 计算：n − 1 ≈ 1.68e7 FLOP（加法）。
• 算术强度 I = 1.68e7 / 67.1e6 B = 0.25 FLOP/B ≪ I* = 14.2 → 极端带宽受限。
• T = max(67.1e6/900e9, 1.68e7/12.75e12) = max(74.6 µs, 1.3 µs) = 74.6 µs。
• 优先优化：任何”减少字节”的手段——用多阶段 kernel 避免中间数组往返、用 float4 向量化访存提高事务效率、确保合并访存、必要时改用更窄的数据类型（bf16/fp16 会把流量砍半，但要注意精度）。优化计算毫无意义。

(B) y[i] = a*x[i] + y[i]

• 流量：读 x（67.1 MB）+ 读 y（67.1 MB）+ 写 y（67.1 MB）= 201.3 MB。
• 计算：2n = 3.36e7 FLOP（1 次 FMA）。
• I = 3.36e7 / 201.3e6 = 0.167 FLOP/B → 带宽受限，且比 (A) 更严重（要搬 3 份数据）。
• T = max(201.3e6/900e9, 3.36e7/12.75e12) = max(223.7 µs, 2.6 µs) = 223.7 µs。
• 优先优化：把 y 留在片上（寄存器分块：每个线程一次读入 k 个 y、算完再写回，减少 y 的重复读写）、提高每次访存的元素数（float4）、把 x 与 y 的访问合并成连续的大事务。若整个数组能驻留 L2 也可以显著降低 DRAM 流量。

(C) 10 阶多项式

• 流量：读 x（67.1 MB）+ 写 y（67.1 MB）= 134.2 MB（系数在常量内存/立即数里，可忽略）。
• 计算：10 FMA × n = 1.68e8 FMA = 3.36e8 FLOP。
• I = 3.36e8 / 134.2e6 = 2.5 FLOP/B，仍 < 14.2 → 仍偏带宽受限，但已接近临界（差 5.7 倍）。
• T = max(134.2e6/900e9, 3.36e8/12.75e12) = max(149.1 µs, 26.4 µs) = 149.1 µs。
• 优先优化：仍然是减少字节（合并访存、向量化、避免 y 的额外读），但当 K 继续增大（例如 50 阶）时 I = 12.5，会逼近 I*，此时 ILP（Horner 链是串行依赖！需要拆成 2–4 条独立链）与取指/发射效率 就变成新的瓶颈——这正是”从带宽受限转向计算受限”的分界线。

总结：三者都落在带宽受限区，但受限程度依次递减（0.25 → 0.167 → 2.5 FLOP/B）。这直接给出优化优先级：先问”能不能少搬字节”，再问”能不能算得更快”；只有当 I > I* 时（本例中需要 K ≥ 57 阶多项式这类高复用计算）才轮到 FMA 吞吐、ILP 与 Tensor Core 的优化。

问题 4（执行语义判断）：判断下列 CUDA 代码片段的合法性与正确性，并说明理由。

// (1) 直方图: 所有 block 用原子操作更新同一 global 数组
__global__ void k1(int* A, int* counts, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) atomicAdd(&counts[A[i]], 1);
}

// (2) 块间握手: block 1 等待 block 0 完成
__global__ void k2(int* myFlag) {
    if (blockIdx.x == 0) { /* 干活 */ atomicAdd(myFlag, 1); }
    else { while (atomicAdd(myFlag, 0) == 0) { } /* 干活 */ }
}

// (3) 块内屏障放在分支里
__global__ void k3(float* A, int n) {
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    if (i < n) {
        A[i] *= 2.0f;
        __syncthreads();     // 只有部分线程会执行到
    }
}

【答案】

(1) 合法且正确。 原子操作用于互斥（保证对 counts[A[i]] 的读-改-写不可分割），它不约束 work scheduler 的调度自由：无论 block 以任何顺序、任何并发度执行，最终计数都相同（因为加法可交换并结合）。这是 CUDA 中”跨 block 通信”的少数合法形式之一。注意性能问题仍然存在：同一个 bin 上的原子加会串行化，热点 bin 决定 Span，因此需要 block 私有的共享内存副本或分片计数来降低争用。

(2) 不合法/可能死锁。 它依赖”block 0 与 block 1 并发执行”这一CUDA 从未承诺的语义。CUDA 只保证 block 可以按任意顺序执行（无依赖假设）。在”只有 1 个 SM 且每 SM 只能驻留 1 个 block”的 GPU 上，如果调度器先运行 block 1，block 1 会永远自旋（因为 block 0 得不到资源），整个 kernel 挂死。正确做法是用 kernel 分解：把”干活”和”等待后继续干活”放到两次 launch 中，利用 kernel 边界作为全局同步点。

(3) 不合法（未定义行为，可能挂死）。 __syncthreads() 是块内全部线程的屏障：它要求 block 中的每一个线程都到达该点。这里屏障被放在 if (i < n) 分支内，当 n 不是 blockDim.x 的整数倍时，最后一个 block 中部分线程不会执行到屏障（并且同一个 warp 内的 lane 会分叉），导致：(a) 在不支持独立线程调度的旧架构上直接永久挂起；(b) 在 Volta 及以后虽然支持 per-thread 调度，但 __syncthreads() 语义仍要求所有非退出线程到达，且 warp 内分叉执行屏障是明确定义的未定义行为。修改：把 __syncthreads() 移到分支外面（所有线程都会经过的公共路径上），分支只包住真正的计算；或者根本不给这个 kernel 加屏障（它本来不需要——每个线程只写自己的 A[i]）。

总结：判断一段 CUDA 代码是否合法，只看两个问题——(i) 它是否假设了 block 之间的执行顺序或并发性？(ii) 它是否让块的全体线程无法一致地到达每一个 __syncthreads()？ 任何一条答”是”，就必须重写。