Lecture 11: 高级主题 —— 多流、数据传输、张量核心、深度学习与 CUDA 替代方案 (对应最终项目)
Lecture 11: 高级主题 —— 多流、数据传输、张量核心、深度学习与 CUDA 替代方案 (对应最终项目)
概述
本讲把 CUDA 从”写 kernel”推进到”调度整个系统”：GPU 不再是孤立的计算盒，而是挂在 PCIe/NVLink 互连上、拥有自己内存空间与页表的协处理器（co-processor）。核心问题是带宽——PCIe 带宽（16–64 GB/s）比显存带宽（1555 GB/s on A100）低两个数量级，所以一旦程序在拷贝上等待，再快的 kernel 也无济于事。本讲给出两个层次的答案：微观上，用固定内存（pinned memory）让 DMA 真正生效、用 cudaMemcpyAsync 让拷贝不阻塞主机；宏观上，用多流（stream）把数据切段，让”搬第 N 段”和”算第 N-1 段”重叠，把总时间从 T_transfer + T_compute 压向 max(T_transfer, T_compute)。随后本讲把视线转向计算侧：Tensor Core 用”16×16×16 矩阵乘加阵列”把 FP16 峰值做到 FP32 CUDA core 的 16 倍，而深度学习正是靠 GEMM 化的卷积与全连接层吃下这部分算力；最后比较 CUDA 的诸多替代方案（OpenCL / SYCL / HIP / OpenACC / OpenMP offload / Vulkan / WebGPU），并给出 WebGPU+WGSL 与 CUDA 的逐项对照——这正是课程最终项目所使用的编程模型。
核心概念与 GPU 架构图解
概念 1：GPU 在 PC 体系结构中的位置（GPU as Part of the PC Architecture）
定义与目的：GPU 是插在系统互连总线上的独立设备，它有自己的显存（device DRAM）、自己的地址空间和页表、自己的调度器。CPU 侧要通过 PCIe（或 NVLink）这一条”窄管子”才能把数据交给它。理解这个概念的目的只有一个：算出数据搬移的时间下界，判断一个应用到底该不该上 GPU，以及上 GPU 后瓶颈会不会落在总线上。
直观解释（”它是什么？”）：把 CPU 想成市中心的大仓库，GPU 想成郊区一座算力超强但只能存自己货的工厂，PCIe 就是连接两地的一条高速公路。工厂里有 108 条产线（108 个 SM），一小时能加工 1555 GB 的原料；可高速公路一小时只放行 16–64 GB 的卡车。如果每加工 1 字节原料只需要搬 1 字节进出，那么工厂 97% 的产能都在等卡车。 这就是”带宽是重力（bandwidth is the gravity）”的含义：缓冲、重排、缓存这些技巧能暂时对抗它，但最终性能还是被 speeds and feeds 决定。要摆脱这条公路，只有两种办法：把原料留在工厂里反复用（tiling / 数据复用），或者直接把两座工厂修在一起（NVLink / NVSwitch）。
架构/机制图解：
【经典 PC 体系结构（Northbridge/Southbridge 时代）】
┌─────────────┐
│ CPU │
└──────┬──────┘
│ FSB（前端总线）
┌──────────┴───────────────┐
│ Northbridge │ ← 必须与 CPU/DRAM/显卡三者高速通信
│ (内存控制器 + 显卡接口) │ 显卡需要"一等公民"级别的 DRAM 访问权
└───┬───────────┬──────────┘
│ │ AGP / 早期 PCIe（AGP 峰值仅 2 GB/s）
DRAM ┌──┴──────────┐
│ GPU 显卡 │
└─────────────┘
┌────────────┴─────────────┐
│ Southbridge │ ← 慢速 I/O 的汇聚点
└──┬────────┬─────────┬────┘
PCI USB SATA
（33 MHz / 32-bit = 132 MB/s 峰值；66 MHz/64-bit = 528 MB/s 峰值）
【现代 PCIe switched fabric 体系结构（Northbridge/Southbridge 变成 PCIe switch）】
┌─────────────┐
│ CPU │ ← PCIe Root Complex 已集成进 CPU 芯片内
│ (含内存控制器│ 内存通道直连 DRAM，不再绕北桥
│ 与 PCIe RC) │
└──┬───────┬──┘
│ │ DDR5 通道
PCIe x16 │ └────────► DRAM
（显卡首选插槽）│
┌────┴──────────────┐
│ PCIe Switch │ ← 点对点、交换式；不再是共享总线仲裁
└─┬────┬────┬────┬──┘
x4 │ x4 │ x4 │ x4 │ x1 │
┌─────┴┐ ┌───┴─┐ ┌──┴──┐ ┌───┴──┐
│ NVMe │ │ NIC │ │ USB │ │ SATA │
└──────┘ └─────┘ └─────┘ └──────┘
关键演进：
共享总线 + 仲裁 → 交换式点对点，每张卡有专属链路，无仲裁
"赢得仲裁者成为 bus master" → 报文交换、虚拟通道、带优先级 QoS（实时视频流）
设备寄存器映射进 CPU 物理地址空间（boot 时分配地址，Memory Mapped I/O）
【PCIe link 的物理结构：一条 lane = 4 根线 = 2 对差分线对】
┌──────────────── 发送对 TX（2 根线）────────────────┐
│ D+ ──────────────────────────────────────► │
│ D- ──────────────────────────────────────► │ 差分信号：抗共模噪声
└───────────────────────────────────────────────────┘
┌──────────────── 接收对 RX（2 根线）────────────────┐
│ D+ ◄────────────────────────────────────── │
│ D- ◄────────────────────────────────────── │ 收发同时进行 = 全双工
└───────────────────────────────────────────────────┘
每个方向 1 bit/时钟（每对线 8 Gb/s @ Gen3）
link = n 条 lane：x1 / x2 / x4 / x8 / x12 / x16 / x32
编码开销：
Gen1/Gen2：8b/10b 编码，20% 开销（10 bit 传 8 bit 数据）
Gen3 及以后：128b/130b 编码，仅 1.5% 开销
- 目标不变：维持 DC 平衡 + 足够的状态跳变做时钟恢复
- 用扰码器（scrambler）保证长串 0/1 的概率极低
- 每 66 bit 至少一次电平跳变
- 130 bit 码字中 0/1 各占 63~67 个，非常均衡
Gen3 每 lane 净速率 = 7.8768 GB/s（单向）
x1 = 0.985 GB/s x8 = 7.9 GB/s
x2 = 1.97 GB/s x16 = 15.8 GB/s（单向）
x4 = 3.94 GB/s x32 = 31.5 GB/s（单向）
【显存带宽 vs PCIe 带宽：两个数量级的分工】
┌──────────────────────── A100 (sm_80) ────────────────────────┐
│ │
│ ┌──────────────────────────────────────────────────────┐ │
│ │ HBM2 显存 80 GB 带宽 ≈ 1555 GB/s │ │
│ └────────────────────────┬─────────────────────────────┘ │
│ │ ~1555 GB/s（片内，宽而短） │
│ ┌────────┴─────────┐ │
│ │ 108 个 SM │ FP32 19.5 TFLOPS │
│ │ 6912 CUDA core │ FP16 TC 312 TFLOPS │
│ │ 40 MB L2 │ │
│ └────────┬─────────┘ │
└────────────────────────────┼─────────────────────────────────┘
│
PCIe Gen4 x16│ 单向 ≈ 32 GB/s（理论）
│ 实测 pinned ≈ 22~25 GB/s
│ 实测 pageable ≈ 10~12 GB/s
│
┌────────────────────────────┴─────────────────────────────────┐
│ 主机 CPU + 主机 DRAM（DDR4/DDR5，带宽 ≈ 100~400 GB/s） │
└──────────────────────────────────────────────────────────────┘
比值：1555 / 32 ≈ 49 倍（Gen4）；1555 / 64 ≈ 24 倍（Gen5 x16）
1555 / 600 ≈ 2.6 倍（A100 NVLink 第三代，GPU 与 GPU 直连）
- 关键操作与性能特征：PCIe 各代 x16 单向下行/上行理论带宽为 Gen3 15.8 GB/s、Gen4 31.5 GB/s、Gen5 63 GB/s（128b/130b 编码后的净值，Gen4 为 16 GT/s、Gen5 为 32 GT/s）；实测可用带宽通常只有理论的 70%~80%。对比之下，A100 的显存带宽 1555 GB/s、NVLink 3.0 可达 600 GB/s（讲义原文：”GPU-to-GPU direct bandwidth to 600 gigabytes per second (GB/s), almost 10X higher than PCIe Gen4”），是 PCIe Gen4 的约 10 倍。历史数据也很能说明趋势：AGP 峰值 2 GB/s、初代 PCI 132 MB/s、NVLink 1.0（P100）160 GB/s 双向（是 PCIe Gen3 x16 的 5 倍）、NVLink 2.0（V100）25 GB/s/lane × 6 lane = 150 GB/s、NVLink 3.0（A100）600 GB/s。延迟方面，一次 PCIe 往返（host→device→host）典型在 5~20 μs 量级，而显存访问延迟只有 400~800 个时钟周期（A100 @1.41 GHz 约 0.3~0.6 μs）——也就是说一次 PCIe 设备同步的开销约等于 20 次以上的显存访问，这就是”拷贝次数比拷贝字节数更致命”的原因。
概念 2：主机-设备数据传输优化与固定内存（Data Transfer Optimization & Pinned Memory）
定义与目的：
cudaMemcpy是同步（synchronous）语义：它把数据从主机内存搬到设备内存，函数返回时数据已经在设备上。若主机内存来自malloc，它是可分页内存（pageable memory），DMA 引擎不能直接访问（DMA 使用物理地址，而操作系统随时可能把这一页交换出去）。NVIDIA 驱动不得不先偷偷把数据拷到一块临时的固定内存（pinned / page-locked memory）再发起 DMA，这就是”两次拷贝”。cudaMallocHost/cudaHostAlloc分配的固定内存不可换页，物理地址固定，DMA 可以直接读写，因此带宽通常提升 2 倍以上，也是cudaMemcpyAsync真正异步的前提。直观解释（”它是什么？”）：把 DMA 想成一个只会照着”物理门牌号”送货的快递员，而可分页内存像一堆随时会被搬家的出租屋——快递员刚记下门牌号，房子就被拆了。操作系统只好先把包裹搬到一栋永不动迁的仓库（固定内存），再让快递员去取，于是变成了”搬两次”。固定内存就像把货永久卸在自家院子里：快递员直接开车进去装货。代价是这块地皮（物理内存）从此不能挪作他用——固定内存过多会让操作系统没有可换页的空间，严重过订（over-subscription）会导致整机性能崩塌甚至失败，因此讲义明确提醒”pinned memory is a limited resource whose over-subscription can have serious consequences”。
架构/机制图解：
【可分页内存的 cudaMemcpy：两次拷贝】
h_A (malloc, pageable) staging pinned buffer d_A (device)
┌──────────────────┐ ┌──────────────────┐ ┌──────────────┐
│ 虚拟页 (可换出) │ CPU 拷贝 │ 物理地址固定 │ DMA │ 显存 │
│ 物理页 A/B/C │ ───────► │ (驱动内部临时) │ ───────► │ │
└──────────────────┘ 内存带宽 └──────────────────┘ PCIe └──────────────┘
↑ 第 1 跳：由 CPU 执行，占用主机 CPU 与内存带宽
↑ 第 2 跳：由 DMA 引擎执行
实测：Gen4 x16 约 10~12 GB/s ← 受限于两次拷贝的串行叠加
【固定内存的 cudaMemcpy：一次 DMA】
h_A (cudaHostAlloc, pinned) d_A (device)
┌──────────────────────────┐ ┌──────────────┐
│ 物理地址固定，永不被换出 │ ───────── DMA ───────────► │ 显存 │
└──────────────────────────┘ PCIe └──────────────┘
实测：Gen4 x16 约 22~25 GB/s ← 约 2 倍于可分页内存（与讲义"about 2X faster"一致）
【CPU / GPU / Copy Engine 三方并发（deviceOverlap）】
时间 ───────────────────────────────────────────────────────────────►
┌──────────────┬───────────────┬─────────────────┬───────────────┐
│ H2D A (2.9ms)│ H2D B (2.9ms)│ kernel (9.4ms) │ D2H C (2.9ms) │ 串行：
└──────────────┴───────────────┴─────────────────┴───────────────┘ 18.1 ms
PCIe ↑ 满 PCIe ↑ 满 PCIe 全空闲 PCIe ↓ 满
GPU core 全空闲 GPU core 全空闲
┌──────────────┬─────────────────────────────────┬───────────────┐
│ H2D A (2.9ms)│ kernel (9.4ms) │ D2H C (2.9ms) │ 重叠：
│ Copy Engine 1 │ SM 计算 │ Copy Engine 2 │ max(8.7, 9.4)
└──────────────┴─────────────────────────────────┴───────────────┘ + 边界 ≈ 11.6 ms
↑ H2D 走 Copy Engine 1 ↑ D2H 走 Copy Engine 2（与 H2D 独立）
- 关键操作与性能特征：设备能力要用
cudaGetDeviceProperties查询——prop.deviceOverlap（是否支持”kernel 执行同时做 device↔host 拷贝”）、prop.asyncEngineCount（异步引擎数量：≥1 表示 H2D 与 kernel 可重叠，≥2 表示 H2D 与 D2H 可同时进行）。cudaHostAlloc(&ptr, bytes, cudaHostAllocDefault)分配、cudaFreeHost(ptr)释放；返回的指针可以像malloc的返回值一样使用，唯一区别是内存不会被换页。带宽上，Gen4 x16 的实测典型值为：pageable H2D ≈ 10~12 GB/s、pinned H2D ≈ 22~25 GB/s、pinned D2H ≈ 24~26 GB/s；若把一次 384 MB 的向量加法（A、B 各 128 MB 进，C 128 MB 出）的拷贝时间折算出来，pageable 需要约 35 ms，pinned 只需约 16 ms，而 A100 上这个 kernel 本身只需 0.25 ms——PCIe 传输时间占总时间的 99%，这就是必须用 pinned + stream 的定量理由。
概念 3：CUDA 流、异步拷贝与重叠执行（CUDA Streams, Async Copy & Overlap）
定义与目的：流（stream）是一个操作队列（queue），里面按顺序存放 kernel 启动和
cudaMemcpy请求。同一个流内的操作严格按序执行（前面的 memory copy 结束前，后面的 kernel 不会启动）；不同流之间没有隐含顺序，可以并发执行——这是一种任务并行（task parallelism）。流存在的唯一目的是打破”拷贝与计算互相等待”的串行链，让 PCIe 方向与 SM 计算单元同时忙碌。直观解释（”它是什么？”）：单流就像只有一个窗口的银行：办完一笔业务（拷贝）才能办下一笔（计算），排队的人都闲着。多流相当于开了多个窗口，并且把大额业务拆成小份：柜台 1 正在给 3 号客户数钱（传输第 3 段），柜台 2 同时在给 2 号客户办理（计算第 2 段）。流水线的理想极限是：总时间 =
max(传输总时间, 计算总时间)，而不是两者之和。但要拿到这个极限，必须让每个”阶段”都有活干：只有当段数足够多、且拷贝引擎与计算引擎都能被喂饱时，流水线才真正连续。架构/机制图解：
【流的语义：队列 + 引擎】
主机线程 (host thread) 设备驱动 (device driver)
┌────────────────────────┐ ┌──────────────────────────┐
│ cudaMemcpyAsync(d,h,N)│─── 入队 ────► │ Stream 0 (FIFO 队列) │
│ kernel<<<g,b,0,strm0>>>│─── 入队 ────► │ ┌────────────────────┐ │
│ cudaMemcpyAsync(h,d,N)│─── 入队 ────► │ │ 1. Memcpy A.0 │ │
│ cudaStreamSynchronize │◄── 等待 ───── │ │ 2. Kernel 0 │ │ 同一队列内
└────────────────────────┘ │ │ 3. Memcpy C.0 │ │ 严格按序
│ └────────────────────┘ │
另一个 host 线程或同一线程： └──────────────────────────┘
┌────────────────────────┐ ┌──────────────────────────┐
│ kernel<<<g,b,0,strm1>>>│─── 入队 ────► │ Stream 1 (FIFO 队列) │ 不同队列之间
└────────────────────────┘ │ 1. Memcpy A.1 / 2. K.1 │ 可并发
└──────────────────────────┘
【概念视图：两个流 + 三个硬件引擎】
操作（Operations：Kernel / MemCpy）
┌─────────────── Stream 0 ───────────────┐ ┌─────────────── Stream 1 ───────────────┐
│ MemCpyA.0 ▸ MemCpyB.0 ▸ Kernel0 ▸ MemCpyC.0 │ MemCpyA.1 ▸ MemCpyB.1 ▸ Kernel1 ▸ C.1 │
└───────┬──────────┬──────────┬──────────┘ └──────┬─────────┬─────────┬────────────┘
▼ ▼ ▼ ▼ ▼ ▼
┌─────────────────────────────────┐ ┌─────────────────┐ ┌───────────────┐
│ Copy Engine（PCIe UP / DOWN） │ │ Kernel Engine │ │ Copy Engine │
│ H2D 与 D2H 可以不同引擎 │ │ (SM 阵列) │ │ (另一条方向) │
└─────────────────────────────────┘ └─────────────────┘ └───────────────┘
【串行版时间轴（single stream）：PCIe 只有一个方向在工作，GPU 大量空闲】
时间 ─────────────────────────────────────────────────────────────────────►
┌────────────────┬────────────────┬────────────────────┬────────────────┐
│ H2D A 64 MB │ H2D B 64 MB │ kernel 9.4 ms │ D2H C 64 MB │
│ 2.9 ms │ 2.9 ms │ (137 GFLOP FP32) │ 2.9 ms │
└────────────────┴────────────────┴────────────────────┴────────────────┘
PCIe↑ 忙 PCIe↑ 忙 PCIe 空闲 PCIe↓ 忙
SM 空闲 SM 空闲 SM 满负荷 SM 空闲
总时间 = 2.9 + 2.9 + 9.4 + 2.9 = 18.1 ms （PCIe 有效利用率 ≈ 48%）
【4 段 × 4 流流水线时间轴：拷贝与计算重叠】
时间轴刻度（每格 ≈ 0.73 ms = 一段单向传输时间；kernel 段 ≈ 2.35 ms）
0 0.73 1.46 2.19 2.92 3.65 4.38 5.11 5.84 6.57 7.30 8.03 8.76 9.49 10.22 10.95 11.68
├─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┼─────┤
s0 H2D A0 ██████
s0 H2D B0 ██████
s0 KERNEL0 ██████████████████
s0 D2H C0 ██████
s1 H2D A1 ██████ ← 与 s0 的 H2D B0 争用 PCIe，但 PIe 有富余
s1 H2D B1 ██████
s1 KERNEL1 ███████████████████
s1 D2H C1 ██████
s2 H2D A2 ██████
s2 H2D B2 ██████
s2 KERNEL2 ███████████████████
s2 D2H C2 ██████
s3 H2D A3 ██████
s3 H2D B3 ██████
s3 KERNEL3 ██████████████████
s3 D2H C3 ██████
▲ 流水线填充（无法重叠，约 1.46 ms） ▲ 排空（约 0.73 ms）
稳态后：SM 与 Copy Engine 同时工作，总时间 ≈ 2.35×4 + 2×0.73 + 0.73 ≈ 11.6 ms
加速比 ≈ 18.1 / 11.6 ≈ 1.56x
【讲义给出的通用流水线公式（以"输入传输 + 计算 + 输出传输"三段为例）】
设：T_in = 输入传输总时间，T_comp = 计算总时间，T_out = 输出传输总时间，N = 段数
串行时间 = T_in + T_comp + T_out
流水线时间 = (T_in + T_comp + T_out) / N + (N - 1) × max(T_in, T_comp, T_out) / N
└──── 第一段的填充（无法重叠）────┘ └──── 稳态：每段只花瓶颈阶段的时间 ────┘
N → ∞ 的极限 = max(T_in, T_comp, T_out)
讲义例题：T_in = 10 s，T_comp = 10 s，T_out = 10 s，kernel 启动常数开销 1 s
时间(N) = (20 / N) + 10 + N
N = 1 : 31 s N = 2 : 22 s N = 3 : 19.67 s
N = 4 : 19 s ★ N = 5 : 19 s ★ N = 6 : 19.33 s N = 8 : 20.5 s
最优解：4 或 5 段（19 s），加速比 31 / 19 ≈ 1.63x
- 关键操作与性能特征：重点 API 与语义如下。
cudaStreamCreate(&s)/cudaStreamDestroy(s)：创建/销毁一个流。cudaStreamCreateWithFlags(&s, cudaStreamNonBlocking)创建非阻塞流，它不会与 legacy default stream 做隐式同步。cudaMemcpyAsync(dst, src, bytes, kind, stream)：异步拷贝。只有当主机内存是 pinned 时它才真正异步；若源/目标是 pageable 内存，驱动仍需 staging，调用会（几乎）同步阻塞，重叠收益消失。kernel<<<grid, block, shmem, stream>>>(args)：把 kernel 投到指定流（第 4 个执行配置参数）。写<<<grid, block>>>等价于投到默认流。cudaEventRecord(ev, stream)与cudaEventSynchronize(ev)：在流中打时间戳/等待点，配合cudaEventElapsedTime(&ms, e0, e1)测量时间。cudaEventCreateWithFlags(&e, cudaEventDisableTiming)用于纯同步用途，开销更低。cudaStreamWaitEvent(stream, ev, 0)：跨流同步——让 stream 等待 ev 记录的完成点，而不阻塞主机。cudaStreamSynchronize(s)：阻塞主机直到该流中所有操作完成。- 默认流的隐式同步：legacy default stream（也叫 NULL stream / stream 0）是”同步流”——默认流中的操作会等待此前所有其他流中的操作完成，也会阻塞之后其他流中的操作。这会让多流程序意外串行化。用
--default-stream per-thread编译选项、或显式使用cudaStreamPerThread、或所有流都用cudaStreamNonBlocking创建，可以绕开这个陷阱。 - 硬件并发度：Fermi 支持 16 路并发的 grid，但来自不同流的 kernel 会复用（multiplex）到同一条硬件工作队列，只有”流边界”处才能重叠；Kepler 起每个流有独立的硬件工作队列，支持 32 路并发、全流级并发、流间无依赖。当代 GPU 还提供 Hyper-Q（每个引擎多条真实队列），让某些流在某引擎上阻塞时其他流仍能推进。
- 段大小选择：段太小会遇到”执行时间不可能趋近于零”的下界——部分 SM 闲置、warp 数不足以打满 SM、线程数不足一个 warp、block 数太少导致负载不均、kernel 启动本身有固定开销；数据传也有类似的非线性（host 与 DMA 的启动成本）。因此要选中等大小的段，最佳值依赖具体 GPU（讲义结论：”Use Moderate Segment Size and Device Query”）。
- 三流原则：讲义指出，在”输入拷贝 → 计算 → 输出拷贝”的链上，若只用两个流，
C.1会在拷贝引擎队列里挡住下一轮的A.2、B.2（表头阻塞 head-of-line blocking）；要实现连续流水（continuously pipelined timing）需要三个流，并按”先发所有 H2D、再发所有 kernel、最后发所有 D2H”的顺序入队，避免读方向被写方向阻塞。
概念 4：Tensor Core 与混合精度（Tensor Cores & Mixed Precision）
定义与目的：Tensor Core 是 NVIDIA 从 Volta（sm_70）开始集成在 SM 里的专用矩阵乘加阵列（matrix multiply-accumulate unit）。它的基本算子不是”两个标量相乘再相加”，而是整个小矩阵块：
D = A × B + C，典型形状 16×16×16（WMMA）或底层指令级的 8×8×4 / m16n8k16。它解决的问题是：深度学习里 90% 以上的浮点运算都是 GEMM，而 GEMM 的每个乘加对只有 2 个 FLOP 却要读 2 个操作数——用标量 CUDA core 做这件事，性能被寄存器带宽和指令发射率死死限制。Tensor Core 把”取数-相乘-累加”整条数据流硬化（hardwired）在硅片上，配合 16 位输入 + 32 位累加的混合精度（mixed precision），把每 SM 每周期的 FP16 FMA 数从 64 提高到 1024。直观解释（”它是什么？”）：CUDA core 像一屋子熟练工人，每人每周期做一次乘加，指令得一个个发；Tensor Core 像一条自动化流水线：整箱原料（16×16 的 fp16 矩阵）倒进去，出来就是整箱成品（16×16 的 fp32 累加结果），单位时间的产出高一个数量级，代价是它只认得这一种箱子。为什么深度学习能忍受 16 位输入？因为神经网络是统计学习器：它的输出是”哪一类”的概率排序，而不是精确的数值答案；训练时每一步更新权重都是”沿梯度方向走一小步”，单步权重增量的相对误差只要小于步长本身，就会被后续迭代自然吸收。反过来看，把 784 个像素连同权重一起量化到 fp16，注定的舍入噪声相当于给输入加了极小的抖动——这反而起到正则化作用。关键是把累加仍然放在 fp32：点积长度 K 可能是 4096 甚至更大，若累加也用 fp16，误差会随 K 线性增长并迅速淹没结果。所以工业界固定用法是”低精度输入 × 低精度输入 → 高精度累加“。
架构/机制图解：
【一次 16x16x16 的 MMA：D = A x B + C】
A (16x16, fp16) B (16x16, fp16) C (16x16, fp32) D (16x16, fp32)
┌──────────────┐ ┌──────────────┐ ┌──────────────┐ ┌──────────────┐
│ 16x16 fp16 │ │ 16x16 fp16 │ │ 16x16 fp32 │ │ 16x16 fp32 │
│ 256 elements │ │ 256 elements │ │ 256 elements │ │ 256 elements │
│ 512 bytes │ │ 512 bytes │ │ 1024 bytes │ │ 1024 bytes │
└──────┬───────┘ └──────┬───────┘ └──────┬───────┘ └──────▲───────┘
│ │ │ │
▼ ▼ ▼ │
┌──────────────────────────────────────────────────────────────────┐ │
│ Routing Network（片上路由网络） │ │
│ 把 A 的每一行与 B 的每一列按需要广播到各个乘法器（避免寄存器搬运） │ │
└──────┬────────────┬────────────┬────────────┬────────────────────┘ │
▼ ▼ ▼ ▼ │
┌───────────┐ ┌───────────┐ ┌───────────┐ ┌───────────┐ │
│ sub-core 0│ │ sub-core 1│ │ sub-core 2│ │ sub-core 3│ │
│ 4x4x4 MMA │ │ 4x4x4 MMA │ │ 4x4x4 MMA │ │ 4x4x4 MMA │ │
│ + fp32 累加│ │ + fp32 累加│ │ + fp32 累加│ │ + fp32 累加│ │
└─────┬─────┘ └─────┬─────┘ └─────┬─────┘ └─────┬─────┘ │
└─────────────┴─────────────┴─────────────┴───────────────────────┘
一个 warp 级 Tensor Core = 4 个 sub-core
【为什么把 16 位做好能拿到 4x：数据通路变窄，乘法器就能堆更多】
32-bit 通路（FP32 CUDA core）:
┌──────────────── 1 个 32-bit 寄存器 = 32 bit ────────────────┐
│ 1 个 fp32 乘法器 × 32 lane │ → 2 buf fp32
└────────────────────────────────────────────────────────────┘
16-bit 通路（优化后的 Tensor Core 数据通路）:
┌──────────────── 1 个 32-bit 寄存器 = 32 bit ────────────────┐
│ fp16 │ fp16 │ fp16 │ fp16 ← 4 个 16-bit 乘法器/lane │ → 4 buf fp16
└────────────────────────────────────────────────────────────┘
加法器仍用 32-bit（累加精度不降）
⇒ 峰值吞吐提升 4x（讲义："Max Throughput is now increased by 4x"）
⇒ 加上阵列级并行与更宽的发射，A100 上相对 FP32 共提升 16x
【A100 上每 SM 每周期的 FMA 数（用峰值反推）】
FP32 CUDA core : 19.5e12 / (2 × 108 SM × 1.41e9 Hz) = 64 FMA/SM/cycle
FP16 Tensor Core: 312e12 / (2 × 108 SM × 1.41e9 Hz) = 1024 FMA/SM/cycle
比值 = 1024 / 64 = 16x
【官方"优化后的硬件"示意（讲义页 8）：整块 tile 一次性进阵列】
Entire M tile (16x16) Entire N tile (16x16)
┌───────────────────┐ ┌───────────────────┐
│ M00 M01 M02 M03 │ │ N00 N01 N02 N03 │
│ M10 M11 M12 M13 │ │ N10 N11 N12 N13 │
│ M20 M21 M22 M23 │ │ N20 N21 N22 N23 │
│ M30 M31 M32 M33 │ │ N30 N31 N32 N33 │
└─────────┬─────────┘ └─────────┬─────────┘
└──────────┬───────────────┘
▼
┌──────────────────────────────────────────────┐
│ * * * * ← 4x4 乘法阵列（一次载入） │
│ + + + + ← 32-bit 加法器（累加） │
│ 用一个 warp 的 2 次 shared memory 载入 │
│ 就能喂满整块 4x4 子 tile（讲义页 7） │
└──────────────────────────────────────────────┘
- 关键操作与性能特征：WMMA API（
#include <mma.h>，命名空间nvcuda::wmma）的四个关键调用构成一个模板化抽象：wmma::fragment<wmma::matrix_a, 16, 16, 16, __half, wmma::row_major> a_frag;—— 声明一个 fragment（每个线程持有其中一部分元素的”寄存器分片”）。wmma::load_matrix_sync(a_frag, ptr, ldm);—— 把内存中的矩阵块搬进 fragment；ldm是主维长度（leading dimension）。16 位类型要求ldm是 8 的倍数（16 字节对齐）。wmma::mma_sync(acc, a_frag, b_frag, acc);—— 执行acc = a × b + acc。wmma::store_matrix_sync(ptr, acc, ldm, wmma::mem_row_major);—— 把累加器写回内存。wmma::fill_fragment(acc, 0.0f)—— 清零累加器。 此外还有wmma::load_matrix_sync(frag, ptr, ldm, wmma::mem_col_major)支持列主序、wmma::store_matrix_sync(ptr, acc, ldm, wmma::mem_col_major)控制写出布局，以及wmma::fragment<wmma::accumulator, 16, 16, 16, double>（sm_80 起支持 FP64 与 TF32 等）。 性能数字：A100 上 FP16 Tensor Core 峰值 312 TFLOPS（稀疏模式下 624 TFLOPS），是 FP32 CUDA core（19.5 TFLOPS）的 16 倍；TF32（19 位尾数的”降精度 FP32”）峰值约 156 TFLOPS，FP64 Tensor Core 约 19.5 TFLOPS。历史上第三代（A100）相对第二代（V100 上 125 TFLOPS）翻了 2.5 倍，并新增 结构化稀疏（structural sparsity） 支持，可再翻倍。现实中，一个只用全局内存直接 load 的朴素 WMMA kernel 通常只能跑到 40~60 TFLOPS——因为它的算术强度只有 8 FLOP/B（每次 k 步每 warp 读 1024 B 算 8192 FLOP），受 L2 带宽与 fragment 加载指令吞吐限制；要把 Tensor Core 喂饱，必须像 CUTLASS/cuBLAS 那样把 tile 先经共享内存做多级复用、双缓冲、并让 fragment 加载与mma_sync交错发射。
概念 5：深度学习中的 GPU 计算：一切皆 GEMM（Deep Learning as GEMM）
定义与目的：深度学习（deep learning）是一种表示学习（representation learning）：不再由人手工设计特征，而是让机器从数据中学出”由其他表示所表达的表示”，一层层堆出从简单特征到抽象特征的深层次（hierarchy）——输入 → 简单特征 → 抽象特征的层 → 输出映射。GPU 进入这个领域是因为训练循环（training cycle）的迭代速度被彻底改变了：2007 年可编程 GPU 加速了训练循环，之后专用芯片进一步加速，直接引发了计算机视觉、语音识别、机器翻译、自动驾驶的复兴。从工程角度看，本概念要解决的问题是：把各种神经网络层（全连接、卷积、池化）统一翻译成一种高度优化的算子——GEMM（通用矩阵乘法），从而复用整套 tiling、shared memory、Tensor Core 优化。
直观解释（”它是什么？”）：一个感知机（perceptron）就是
y = sign(W·x + b)——一次点积。全连接层就是”很多个感知机排成一行”：输入向量x乘权重矩阵W，得到输出向量y。一层做完做下一层，就是前向传播（forward propagation / inference）。训练则是反向传播（backward propagation）：拿预测与标签算误差E，用链式法则（chain rule）把误差梯度一路传回去，求出dE/dW，再让每个权重沿梯度反方向挪一小步Θ_{i+1} = Θ_i - ε·ΔΘ，这就是随机梯度下降（stochastic gradient descent, SGD）。这些操作的共同骨架都是”矩阵乘矩阵”或”矩阵乘向量”：fc1 = W1 · x + b1是 GEMV（batch 化后就是 GEMM）；卷积是”滑窗点积”，把输入展开（unroll / im2col）后同样是 GEMM；反向传播算dE/dW = (dE/dfc)·xᵀ也是 GEMM。所以只要把 GEMM 做到极致，整个深度学习就快了。架构/机制图解：
【多层感知机 MLP 与术语（讲义页 45~47）】
x[0] x[1] x[2] x[3] z[0] z[1] z[2] k[0] k[1] k[2]
● ● ● ● ● ● ● ● ● ●
\ / \ / \ / \ / \ / \ | /
\ / \ / \ / \ / \ / Softmax
┌──┴───────┴──────┴──┐ ┌──┴───────┴──┐ ┌────────────┐
│ W1 是 [4x4] │ │ W2 是 [4x3] │ │ 归一化成概率 │
│ b1 是 [4x1] │ │ b2 是 [3x1] │ │ Σk[i] = 1 │
└─────────────────────┘ └─────────────┘ └────────────┘
输入层 (Input Layer) 隐藏层 (Hidden) 输出层 (Output)
W[i,j] = 第 i 个输入到第 j 个输出的权重（讲义术语：Entry i,j is weight
between ith input and jth output）
【MNIST 手写数字识别的 MLP（讲义页 4）】
28x28 灰度图 = 784 个输入节点 → 10 个隐藏节点 → 10 个输出节点（Digit 0..9）
参数量：784×10 + 10 (L1 权重 + 偏置)
+ 10×10 + 10 (L2 权重 + 偏置)
= 7,960 个参数
公式：k[m] = activation( Σ_c W1[m,c]·x[c] + b1[m] )
【大图上的 MLP 会爆炸（讲义页 26）】
250 x 250 图像 → 全连接层每个节点 62,500 个权重
节点数与之相当 → 总共约 40 亿个权重
⇒ 计算量与内存都不可接受 ⇒ 需要用卷积做权重共享（weight sharing）
【为什么"卷积 = GEMM"：unrolled / im2col 变换（讲义页 27、31、32）】
输入特征图 X (C=3 通道, 3x3) 卷积核 W (M=2 个, 3 通道, 2x2)
┌───────────────┐ ┌───────────────┐
│ 通道0 通道1 通道2 │ │ m=0: 3 个 2x2 │
│ 3x3 3x3 3x3 │ │ m=1: 3 个 2x2 │
└───────────────┘ └───────────────┘
│ │
│ im2col / unroll：把每个输出位置对应的
│ K x K 感受野拉成一列，堆叠成一个矩阵
▼ ▼
X_unrolled (行数 = C·K·K = 12, 列数 = H_out·W_out = 4) W' (M=2 行 x 12 列)
┌────────────────────────────────────────────────┐ ┌──────────────┐
│ 1 1 2 2 │ │ w00 w01 w02 │
│ 1 1 1 1 ← 每列 = 一个输出位置 (h,w) │ × ├──────────────┤
│ 0 1 1 0 线性化下标 = h*W_out + w │ │ w10 w11 w12 │
│ 1 0 0 1 │ └──────────────┘
│ 2 1 2 1 （共 12 行 = 3 通道 x 2x2 核） │
│ 1 2 2 0 │
└────────────────────────────────────────────────┘
│ │
└──────────────┬─────────────────────┘
▼
Y = W' * X_unrolled (2 x 4)
┌────────────────────────────┐
│ 14 20 15 24 │ ← 每个元素就是一个输出特征图像素
│ 12 24 17 26 │ 第 0 行 = 输出特征图 0
└────────────────────────────┘
代价分析（讲义页 31、32）：unrolled 矩阵的元素总数 = C·K·K·H_out·W_out
原始输入元素数 = C·(H_out+K-1)·(W_out+K-1)
小的 3x3x3 例子：复制因子 = (3·2·2·2·2)/(3·3·3) = 1.78
若采用 tiled 卷积，输入访问量的缩减因子 = K^2·TILE_WIDTH^2 / (TILE_WIDTH+K-1)^2
【卷积层的可并行维度（讲义页 16）】
① 输出特征图之间并行 —— 数量少（M 通常几十到几百），不足以填满 GPU
② 每个输出特征图的所有像素并行 —— 先按行再按列，数量大，但越深的层越小
③ 输入通道维度的归约 —— 需要原子操作或树形归约（讲义："need atomic operation
or tree reduction"）
⇒ 不同层需要不同策略（讲义："Different layers may demand different strategies"）
【反向传播的两种并行模式】
前向 (Forward / Inference)： 反向 (Backward / Training)：
输入 x ──► [W1] ──► fc1 ──► [W2] ──► k ──► 误差 E dE/dk ◄── [W2]ᵀ ◄── dE/dfc1 ◄── [W1]ᵀ ◄── x
│ │
▼ ▼
Θ = {W, b} 固定 求 dE/dΘ，更新 Θ ← Θ - ε·ΔΘ
每一步都是 GEMM：
dE/dW1 = (dE/dfc1) · xᵀ 形状 [out1 x in] = [out1 x batch] · [batch x in]
dE/dx = W1ᵀ · (dE/dfc1) 用于继续往前传
并行的粒度：
· batch 维度完全独立 —— mini-batch SGD 的天然并行（讲义页 21："Balance between
accuracy of gradient estimation and parallelism"）
· 输出神经元维度独立
· 归约维（输入维度 / 通道维 / 空间核）需要块内归约或原子累加
- 关键操作与性能特征：把这些层写成 GEMM 后的性能特征是高算术强度。以 MNIST 的第一层为例（batch = 1024，输入 784，输出 1024）：FLOPs = 2×1024×784×1024 = 1.64 GFLOP；访存 = X（1024×784×4 B = 3.2 MB）+ W（784×1024×4 B = 3.2 MB）+ Y（1024×1024×4 B = 4.1 MB）≈ 10.5 MB；算术强度 ≈ 1.64e9 / 1.05e7 ≈ 156 FLOP/B。A100 的 Roofline 拐点 = 19.5e12 / 1555e9 ≈ 12.5 FLOP/B，156 ≫ 12.5，所以这是计算受限（compute-bound）的，GPU 正好擅长；而 Tensor Core 把可用的计算上限提到 312 TFLOPS，拐点变为 312e12/1555e9 ≈ 200 FLOP/B，于是即便像 Vision Transformer 里那样中等算术强度的层，也可能掉回内存受限。训练与推理的差别也体现在算术强度上：反向传播需要保存前向的中间激活值（activation），显存占用约是推理的 3~4 倍（需要额外存激活与梯度），这就是 batch size 受限、要用梯度检查点（gradient checkpointing）的原因。数值精度方面，训练通常用 TF32/BF16 前向+反向、FP32 主权重（master weights）；FP16 训练需要 loss scaling 来避免梯度下溢。
概念 6：CUDA 的替代方案（Alternatives to CUDA）
定义与目的：CUDA 只是一个模型。加速计算已经不再是问题，因为 GPU 供应商包括 NVIDIA、AMD、Intel、Samsung、Apple、Qualcomm、ARM 等等；但每个厂商的硬件都有自己的原生存程，于是出现了一整层”可移植编程模型”，目标是用一份代码跑到不同厂商、甚至不同种类的加速器（GPU / CPU / DSP / FPGA）上。这些模型要解决的核心矛盾是可移植性（portability）× 性能上限（performance ceiling）× 编程抽象层次（abstraction level）这三者的三角权衡。
直观解释（”它是什么？”）：把 CUDA 想成”某品牌的专用工具”：最贴合硬件、性能天花板最高、但只能在自家机器上用。OpenCL/SYCL 像”国际标准插座转换头”：通用，但要经过一层抽象，某些厂商的独家特性（如 Tensor Core 的 WMMA 指令）就用不上，性能上限往往低一档。HIP 像”同义词词典”：语法跟 CUDA 几乎一样，
cudaMalloc→hipMalloc，一套源码通过hipify工具自动改写，在 AMD 和 NVIDIA 上都能编，是迁移成本最低的路。OpenACC/OpenMP target offload 像”给普通代码贴便利贴”：在串行循环前面加一行 pragma 就能跑在 GPU 上，编程门槛最低，但性能高度依赖编译器质量，天花板也最低。WebGPU 像”把 GPU 装进浏览器”：天生可移植（Windows/macOS/Linux/Android/iOS/浏览器全通吃），用 WGSL 写着色器，牺牲一部分底层控制力换取零安装、易分发。架构/机制图解：
【加速编程模型的时间线（讲义页 4）】
1992 ─ OpenGL
1995 ─ DirectX
2002 ─ GPGPU（用图形 API 做通用计算的黑客时代）
2007 ─ CUDA ★ 第一个 C 风格、完整的 GPU 通用计算模型
2008 ─ OpenCL 1.0（Apple 发起，AMD/IBM/Qualcomm/Intel/NVIDIA 支持）
2012 ─ OpenACC
2013 ─ C++ AMP、RenderScript
2014 ─ Metal（Apple）、SYCL（Khronos，基于 C++ 模板的单一源码模型）
2016 ─ Vulkan、ROCm HIP、OpenMP 4.5 target offload
2017 ─ OpenCL 2.2；2018 Apple 宣布放弃 OpenCL
【所有加速 API 的共同特征（讲义页 6）——理解这一点，迁移就是"换词"】
硬件侧 软件侧
┌────────────────────────┐ ┌──────────────────────────────────┐
│ 层次化的轻量核心 │ │ 面向 kernel 的加速 │
│ (hierarchy of cores) │ │ 设备内存 vs 主机内存（两份地址空间）│
│ 本地暂存器内存 │ │ 软件管理的内存（无硬件一致性） │
│ (scratchpad / __local) │ │ Grid / Block / Thread 三级层次 │
│ 缺乏硬件缓存一致性 │ │ 块同步并行（Bulk Synchronous │
│ (no HW coherence) │ │ Parallelism, BSP） │
│ 全局原子操作慢 │ │ │
│ 多线程 │ │ │
└────────────────────────┘ └──────────────────────────────────┘
【OpenCL 的内存模型（讲义页 8~9）：与 CUDA 一一对应】
CUDA OpenCL 说明
─────────────────────────────────────────────────────────────────────
__global__ 函数 __kernel void myGEMM2(参数表) kernel 用 __kernel 修饰
__shared__ float __local float block 内的暂存器
__constant__ __constant 常量内存
寄存器/局部数组 __private 私有内存
blockIdx/get_group_id work-group（工作组） 相同语义
threadIdx/get_local_id work-item（工作项） 相同语义
__syncthreads() barrier(CLK_LOCAL_MEM_FENCE) 屏障语义相同
cudaMalloc clCreateBuffer API 风格从 C 函数变成句柄
<<<grid, block>>> clEnqueueNDRangeKernel 启动方式差异最大
【OpenACC 的抽象机模型：gang / worker / vector（讲义页 24、27、28）】
┌─────────────────────────────────────────────────────────────────┐
│ #pragma acc parallel num_gangs(1024) num_workers(32) │
│ ┌────────── gang 0 ──────────┐ ┌────────── gang 1 ──────────┐ │
│ │ worker0 worker1 到 worker31 │ │ worker0 worker1 到 worker31 │ │
│ └─────────────────────────────┘ └─────────────────────────────┘ │
│ 共 1024 个 gang，每个 32 个 worker，合计 32K 个 worker │
└─────────────────────────────────────────────────────────────────┘
对应关系：gang ≈ CUDA block worker ≈ CUDA thread（常映射为 vector lane）
#pragma acc loop gang → 循环的迭代被分配到各个 gang（不是冗余执行）
没有 pragma 的语句 → 被所有 gang 冗余执行（redundant execution）
注意：OpenACC 目前不允许用户指定跨线程的同步（讲义页 23）
OpenACC 版的矩阵乘法（讲义页 16）——与串行代码只差两行 pragma：
#pragma acc parallel loop copyin(M[0:Mh*Mw]) copyin(N[0:Nw*Mw]) copyout(P[0:Mh*Nw])
for (int i = 0; i < Mh; i++) {
#pragma acc loop
for (int j = 0; j < Nw; j++) {
float sum = 0.0f;
for (int k = 0; k < Mw; k++) {
sum += M[i * Mw + k] * N[k * Nw + j];
}
P[i * Nw + j] = sum;
}
}
copyin / copyout 子句负责描述矩阵数据如何在主存与设备内存之间搬移。
现实提醒（讲义页 21、22）：
· 很难写出"有 pragma 和没 pragma 都正确且都高效"的代码
· 有些 OpenACC 程序在忽略 pragma 后会行为不同甚至出错
· 部分 pragma 只是"提示（hint）"，编译器可以不听
· 性能高度依赖编译器质量（比 CUDA / OpenCL 更依赖）
【可移植性 / 性能上限 / 抽象层次 三角对比】
模型 可移植性 性能上限 抽象层次 典型使用者
──────────────────────────────────────────────────────────────────────────────────────
CUDA NVIDIA 独占 ★★★★★ 可达硬件峰值 低（显式线程） 研究者/极致性能
（x86/ARM 主机均可） （Tensor Core 全特性） HPC 内核
HIP/ROCm AMD + NVIDIA（同源码） ★★★★☆ AMD 上接近峰值 低（≈CUDA） 需要双厂商部署
hipify 自动转换 （语法几乎相同）
OpenCL ★★★★★ 全厂商全设备 ★★★☆☆ 受运行时抽象限制 中（显式线程 跨厂商产品
(CPU/GPU/DSP/FPGA) + 显式内存对象）
SYCL ★★★★☆ 单一源码 C++ ★★★★☆ 后端可下沉到原生 中高（C++ 模板 需要现代 C++
(Intel DPC++/AdaptiveCpp) (可调用 CUDA 后端) + 单一源码) + 跨厂商
OpenACC ★★★★☆ 编译器支持处 ★★☆☆☆ 依赖编译器自动 高（pragma） 科学家/遗留代码
(PGI/NVHPC/GCC) 低改动量
OpenMP target ★★★★☆ 主流编译器 ★★☆☆☆ 同样依赖编译器 高（pragma） 已有 OpenMP 代码
offload (GCC/Clang/Intel) 低改动量
Vulkan/ ★★★☆☆ 桌面/移动全平台 ★★★☆☆ 无 Tensor Core 类 中（SPIR-V 游戏引擎/图形
Metal/ 但不含 Web 专用矩阵指令的直接映射 着色器） 混合渲染计算
DirectCompute
WebGPU/WGSL ★★★★★ 浏览器 + 原生 ★★☆☆☆ 受浏览器沙箱与 中（着色器 教学/网页部署
全平台（Win/mac/Linux/ 安全校验约束，无子组矩阵 + 显式资源 零安装分发
Android/iOS） 乘（尚无 matmul 指令） 绑定）
──────────────────────────────────────────────────────────────────────────────────────
注：本课程"Lab 用 WebGPU、最终项目用 RAI"（课程目录原文："C Programming Language and
CUDA Software Development Kit, WebGPU for labs, RAI for final project"）。
【多节点 / 多 GPU 的层次（讲义页 36~49）：CUDA + MPI】
┌──────────────────────── Node 0 ────────────────────────┐
│ CPU0 到 CPUM ── Host Memory ── │
│ │ │ │
│ PCIe PCIe │
│ │ │ │
│ GPU0 到 GPUN （卡间可用 NVLink / NVSwitch 高速互联）│
└──────────────────────┬─────────────────────────────────┘
│ MPI_Send / MPI_Recv（消息传递，不是共享内存）
│ 进程之间通过消息通信、通过消息同步
┌──────────────────────┴─────────────────────────────────┐
│ Node 1 到 Node N │
└────────────────────────────────────────────────────────┘
CUDA 侧的多 GPU 支持：
cudaSetDevice(i) 设置"当前 GPU"
可以在异步调用（kernel / memcpy）还在跑的时候切换当前设备
cudaSetDevice(0); kernel<<<g,b>>>(args); cudaMemcpyAsync(d,h,N,kind,s);
cudaSetDevice(1); kernel<<<g,b>>>(args); ← 合法且常用
MPI 侧：MPI_Init / MPI_Comm_rank / MPI_Comm_size / MPI_Send / MPI_Recv /
MPI_Barrier / MPI_Finalize；启动方式 `mpirun -np X`
在 compute 进程里可以继续把工作卸载到 GPU（讲义页 49 的原话注释：
"Or, can offload to GPU here"）
- 关键操作与性能特征：三种模型层次的代表代码成本可以量化：把一份 CUDA tiled matmul 移植到 OpenCL，源码改动量在 10%~20%（
__global__→__kernel、threadIdx→get_local_id、__shared__→__local、__syncthreads()→barrier(CLK_LOCAL_MEM_FENCE)、cudaMalloc→clCreateBuffer）；移植到 HIP 的改动量接近 0（hipify做源到源转换）；移植到 OpenACC 则代码几乎等于串行版本，只需两行 pragma，但性能方差极大——同一个 pragma 在不同编译器上可能相差 3~5 倍。WebGPU 在性能上受两道额外约束：一是浏览器/驱动的着色器校验限制了复杂控制流与动态索引，二是 WGSL 目前没有子组矩阵乘（subgroup matrix multiply）指令的标准绑定，因此 Tensor Core 无法直接使用，矩阵乘只能靠f32/f16的dot与手工展开近似。实际工程中的选择顺序通常是：先看目标平台（要不要跨厂商、要不要浏览器分发），再看性能需求（要不要 Tensor Core），最后看团队维护成本。
概念 7：WebGPU 编程模型与 CUDA 的对应关系（WebGPU / WGSL vs CUDA）
定义与目的：WebGPU 是 W3C 标准化的、面向浏览器与原生的现代 GPU API，它把 GPU 的能力抽象成三类对象：设备（GPUDevice）、资源（GPUBuffer / GPUTexture）、命令（GPUCommandEncoder → ComputePass / RenderPass）。计算着色器用 WGSL（WebGPU Shading Language） 编写。它解决的问题是”零安装分发 + 跨平台”：一份 WGSL 代码可以在 Windows（D3D12）、macOS/iOS（Metal）、Linux/Android（Vulkan）以及浏览器中运行。ECE408 的实验环境就用 WebGPU，最终项目用构建在其上的 RAI 框架。
直观解释（”它是什么？”）：CUDA 像”自己开车”：你可以直接控制油门、挡位、路线。WebGPU 像”坐公交”：有固定的站点（命令编码器、pass）和票务规则（资源绑定、校验），你不能随便变道，但你能确定这趟车在每个城市都有。对已经会 CUDA 的人，迁移的关键认知是：WebGPU 的所有概念在 CUDA 里都有对应物，只是名字、组织方式与”谁负责什么”变了。最大的一处结构性差异是：CUDA 里每个线程只写”我算哪个元素”，block/grid 的索引是内建的只读变量；WebGPU 里同样如此（
global_invocation_id/local_invocation_id/workgroup_id），但数据进出 GPU 的路径更长：必须在GPUQueue上显式writeBuffer/copyBufferToBuffer（对应cudaMemcpy）或mapAsync读回（对应cudaMemcpyD2H + 同步）。另一个差异是 WebGPU 没有默认流的概念：命令编码器提交后按提交顺序执行，想并行要靠多个 command buffer 与多个队列（目前大多数实现只有一个队列，因此”多流重叠”在 WebGPU 里主要通过一次提交内的多个 pass 与异步mapAsync实现）。架构/机制图解：
【WebGPU 对象模型与 CUDA 运行时对象模型对照】
CUDA 侧 WebGPU 侧
┌────────────────────────────┐ ┌──────────────────────────────────┐
│ cudaGetDeviceProperties │ │ navigator.gpu.requestAdapter() │
│ → cudaSetDevice(i) │ │ → adapter.requestDevice() │
└──────────────┬─────────────┘ └───────────────┬──────────────────┘
▼ ▼
┌────────────────────────────┐ ┌──────────────────────────────────┐
│ cudaMalloc / cudaMemcpy │ │ device.createBuffer({sz,usg}) │
│ → d_A (device global mem) │ │ → GPUBuffer (usage: STORAGE | │
│ │ │ COPY_DST | COPY_SRC) │
└──────────────┬─────────────┘ └───────────────┬──────────────────┘
▼ ▼
┌────────────────────────────┐ ┌──────────────────────────────────┐
│ 编译好的 kernel（fatbin） │ │ device.createShaderModule({code})│
│ │ │ → WGSL 源码运行时编译 │
└──────────────┬─────────────┘ └───────────────┬──────────────────┘
▼ ▼
┌────────────────────────────┐ ┌──────────────────────────────────┐
│ 直接 kernel<<<g,b,s,stream>>>│ │ device.createComputePipeline() │
│ │ │ device.createBindGroup() │
│ │ │ encoder.beginComputePass() │
│ │ │ pass.setPipeline / setBindGroup │
│ │ │ pass.dispatchWorkgroups(nx,1,1) │
│ │ │ device.queue.submit([cmdBuf]) │
└──────────────┬─────────────┘ └───────────────┬──────────────────┘
▼ ▼
┌────────────────────────────┐ ┌──────────────────────────────────┐
│ cudaDeviceSynchronize() │ │ await device.queue.onSubmitted │
│ cudaMemcpy(D2H) 读回 │ │ WorkDone() │
│ │ │ await buffer.mapAsync(READ) │
│ │ │ → getMappedRange() 读回 │
└────────────────────────────┘ └──────────────────────────────────┘
【线程层次与内建变量的对应（这是迁移时最需要"翻译"的部分）】
CUDA WebGPU / WGSL
────────────────────────────────────────────────────────────────────────────
grid（网格，整个 kernel 的所有 block） dispatchWorkgroups(nx, ny, nz) 的调用
block（线程块，最多 1024 线程） workgroup（工作组，默认上限 256）
thread（线程） invocation（调用/线程）
warp（32 个连续线程，SIMT 硬件调度单元） subgroup（子组，大小平台相关，通常 32；
WGSL 提供 subgroupBallot / subgroupAdd 等）
blockIdx.x/y/z @builtin(workgroup_id) : vec3<u32>
threadIdx.x/y/z @builtin(local_invocation_id) : vec3<u32>
blockDim.x/y/z @builtin(workgroup_size) : vec3<u32>
（无直接对应，由 host 传入） @builtin(num_workgroups) : vec3<u32>
全局线程号 = blockIdx*blockDim+threadIdx @builtin(global_invocation_id) : vec3<u32>
__shared__ float s[256] var<workgroup> s : array<f32, 256>
__syncthreads() workgroupBarrier()
__syncwarp() subgroupBarrier()
原子：atomicAdd(&x, v) atomicAdd(&x, v) （WGSL 内建同名函数）
内存栅栏：__threadfence() storageBarrier() / workgroupBarrier()
限制：block ≤ 1024 线程 workgroup ≤ 256 invocation（默认保证值）
────────────────────────────────────────────────────────────────────────────
WebGPU 特有、CUDA 没有的概念：
@group(n) @binding(m) 资源绑定槽位——着色器通过"组号/绑定点"访问缓冲区，
由 host 的 GPUBindGroup 与之匹配（编译期可校验）
@builtin / @location 内建变量与 IO 位置的显式标注（WGSL 是强类型的）
override 着色器"编译期常量"，可在创建 pipeline 时指定（类似模板参数）
storage / uniform / 地址空间划分，对应 CUDA 的 global / constant / shared /
workgroup / private local+register
- 关键操作与性能特征：迁移时需要记住的定量差异有三条。第一，工作组大小上限：CUDA block 上限 1024 线程，WebGPU 只保证 256（
maxComputeInvocationsPerWorkgroup通常为 256，部分实现为 1024），所以同样的问题规模需要更多的 workgroup，块内归约的层次更深（256 → 128 → 64 → 32 → warp shuffle 阶段）。第二，没有默认流同步语义：CUDA 里忘记cudaDeviceSynchronize()读到的是坏数据；WebGPU 里queue.submit后必须await queue.onSubmittedWorkDone()或等mapAsync的 promise resolve，否则getMappedRange()读到的是旧内容——但错误表现完全不同（Promise 未等待通常不会崩，而是静默读到 0）。第三，没有子组矩阵乘：WGSL 的f16支持需要shader-f16扩展特性，且没有mma类指令，因此矩阵乘在 WebGPU 上通常靠向量点积dot(vec4<f32>, vec4<f32>)手工展开到 4 路并行来提高算术强度；这也意味着同样的 GEMM 在 WebGPU 上的 TFLOPS 上限远低于 A100 的 Tensor Core，但在教学和中等规模推理场景下完全够用。第四，性能特征上 WebGPU 与 CUDA 的瓶颈模型是完全一致的：同样受合并访存（coalesced access）、workgroup storage 的 bank 冲突（WebGPU 的实现同样按 32 bank × 4 B 组织 workgroup 存储）、以及占用率（resident workgroup 数）支配——在 CUDA 上学到的优化直觉可以 100% 迁移过去。
代码示例与性能分析
代码示例 1：多流重叠数据传输与计算的完整程序（Pinned vs Pageable + 串行 vs 流水线）
// 文件: stream_overlap.cu
// 编译: nvcc -O3 -arch=sm_80 stream_overlap.cu -o stream_overlap
// 运行: ./stream_overlap
//
// 本程序做三件事:
// (1) 对比 pageable(malloc) 与 pinned(cudaHostAlloc) 主机内存的 cudaMemcpy 带宽
// (2) 对比"串行版"(一次拷完 -> kernel -> 拷回) 与"4 段 x 4 流流水线版"的总耗时
// (3) 用 cudaEvent 精确计时, 并逐段验证数值正确性
//
// 为什么 kernel 里要加一个人工算力负载(heavyVecAdd)?
// 纯 vecAdd 的算术强度只有 0.167 FLOP/B, kernel 时间 ~0.25 ms,
// 而 384 MB 数据过 PCIe 需要 ~16 ms —— 重叠与否看不出差别。
// 把每个元素的 FMA 次数调到 4096, kernel 时间抬到 ~9.4 ms,
// 与 8.7 ms 的传输时间同量级, 重叠收益才可见(这正是讲义"三段流水线"的标准模型)。
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define WORK_PER_ELEMENT 4096 /* 每个元素的 FMA 次数 */
#define N_ELEMENTS (16 * 1024 * 1024) /* 16.7M 个 float = 64 MB */
#define N_BYTES ((size_t)N_ELEMENTS * sizeof(float))
#define N_SEGMENTS 4 /* 切成 4 段, 4 个 stream */
#define SEG_ELEMENTS (N_ELEMENTS / N_SEGMENTS)
#define BLOCK_SIZE 256
#define GRID_SIZE 2048 /* 每段 2048/4=512 个 block */
/* ============================ kernel ============================ */
__global__ void heavyVecAdd(const float * __restrict__ a,
const float * __restrict__ b,
float * __restrict__ c,
long n)
{
long i = (long)blockIdx.x * blockDim.x + threadIdx.x;
const long stride = (long)gridDim.x * blockDim.x;
for (; i < n; i += stride) {
const float x = a[i];
const float y = b[i];
/* 4 条独立累加链: 打破 FMA 的相关性, 让 FMA 吞吐而不是 FMA 延迟成为瓶颈 */
float acc0 = 0.0f, acc1 = 0.0f, acc2 = 0.0f, acc3 = 0.0f;
#pragma unroll 16
for (int k = 0; k < WORK_PER_ELEMENT; k += 4) {
acc0 = fmaf(acc0, 1.0000001f, x);
acc1 = fmaf(acc1, 1.0000001f, y);
acc2 = fmaf(acc2, 1.0000001f, x);
acc3 = fmaf(acc3, 1.0000001f, y);
}
c[i] = (acc0 + acc1) + (acc2 + acc3);
}
}
/* 主机端参考实现: 必须与 kernel 采用完全相同的运算顺序 */
static float host_reference(float x, float y)
{
float acc0 = 0.0f, acc1 = 0.0f, acc2 = 0.0f, acc3 = 0.0f;
for (int k = 0; k < WORK_PER_ELEMENT; k += 4) {
acc0 = fmaf(acc0, 1.0000001f, x);
acc1 = fmaf(acc1, 1.0000001f, y);
acc2 = fmaf(acc2, 1.0000001f, x);
acc3 = fmaf(acc3, 1.0000001f, y);
}
return (acc0 + acc1) + (acc2 + acc3);
}
/* ==================== 实验 A: 拷贝带宽测量 ==================== */
static double copy_bandwidth_test(int pinned)
{
const size_t bytes = 64ull * 1024ull * 1024ull; /* 64 MB */
const int reps = 8;
void *h = NULL;
void *d = NULL;
if (pinned) {
CUDA_CHECK(cudaHostAlloc(&h, bytes, cudaHostAllocDefault));
} else {
h = malloc(bytes);
if (h == NULL) {
fprintf(stderr, "malloc failed\n");
exit(EXIT_FAILURE);
}
}
CUDA_CHECK(cudaMalloc(&d, bytes));
memset(h, 0x5A, bytes);
cudaEvent_t e0, e1;
CUDA_CHECK(cudaEventCreate(&e0));
CUDA_CHECK(cudaEventCreate(&e1));
/* 预热: 把首轮 page fault / 驱动建表等一次性开销排除在计时之外 */
CUDA_CHECK(cudaMemcpy(d, h, bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(e0));
for (int r = 0; r < reps; ++r) {
CUDA_CHECK(cudaMemcpy(d, h, bytes, cudaMemcpyHostToDevice));
}
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
const double gbps = (double)bytes * (double)reps / ((double)ms / 1000.0) / 1e9;
CUDA_CHECK(cudaEventDestroy(e0));
CUDA_CHECK(cudaEventDestroy(e1));
CUDA_CHECK(cudaFree(d));
if (pinned) {
CUDA_CHECK(cudaFreeHost(h));
} else {
free(h);
}
return gbps;
}
/* ==================== 数据初始化与验证 ==================== */
static void init_arrays(float *a, float *b, long n)
{
for (long i = 0; i < n; ++i) {
a[i] = 0.5f + 1.0e-6f * (float)(i % 1024);
b[i] = 0.25f + 1.0e-6f * (float)(i % 512);
}
}
static int verify(const float *h_c, const float *h_a, const float *h_b,
long n, const char *tag)
{
const int ncheck = 16;
int bad = 0;
for (int t = 0; t < ncheck; ++t) {
const long idx = (long)((double)t / (double)ncheck * (double)n);
const float ref = host_reference(h_a[idx], h_b[idx]);
const float got = h_c[idx];
const float rel = fabsf(got - ref) / (fabsf(ref) + 1.0e-6f);
if (rel > 1.0e-3f) {
printf(" [%s] MISMATCH idx=%ld ref=%.6f got=%.6f rel=%.3e\n",
tag, idx, ref, got, rel);
++bad;
}
}
printf(" [%s] 抽样验证: %s (%d/%d 个元素相对误差 < 1e-3)\n",
tag, (bad == 0) ? "PASS" : "FAIL", ncheck - bad, ncheck);
return (bad == 0);
}
/* ============================== main ============================== */
int main(void)
{
/* ---------- 设备能力查询(讲义: Use Moderate Segment Size and Device Query) ---------- */
int dev_count = 0;
CUDA_CHECK(cudaGetDeviceCount(&dev_count));
printf("CUDA 设备数: %d\n", dev_count);
for (int i = 0; i < dev_count; ++i) {
cudaDeviceProp p;
CUDA_CHECK(cudaGetDeviceProperties(&p, i));
printf(" [%d] %-30s CC %d.%d SM=%d deviceOverlap=%d asyncEngineCount=%d\n",
i, p.name, p.major, p.minor, p.multiProcessorCount,
p.deviceOverlap, p.asyncEngineCount);
}
/* ============ 实验 A: pageable vs pinned 拷贝带宽 ============ */
printf("\n===== 实验 A: 主机内存类型对 cudaMemcpy H2D 带宽的影响 =====\n");
const double bw_page = copy_bandwidth_test(0);
const double bw_pin = copy_bandwidth_test(1);
printf(" pageable (malloc) : %7.2f GB/s\n", bw_page);
printf(" pinned (cudaHostAlloc) : %7.2f GB/s (%.2fx)\n", bw_pin, bw_pin / bw_page);
/* ============ 准备数据 ============ */
float *h_a = NULL, *h_b = NULL, *h_c = NULL, *h_c2 = NULL;
CUDA_CHECK(cudaHostAlloc((void **)&h_a, N_BYTES, cudaHostAllocDefault));
CUDA_CHECK(cudaHostAlloc((void **)&h_b, N_BYTES, cudaHostAllocDefault));
CUDA_CHECK(cudaHostAlloc((void **)&h_c, N_BYTES, cudaHostAllocDefault));
CUDA_CHECK(cudaHostAlloc((void **)&h_c2, N_BYTES, cudaHostAllocDefault));
init_arrays(h_a, h_b, N_ELEMENTS);
float *d_a = NULL, *d_b = NULL, *d_c = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_a, N_BYTES));
CUDA_CHECK(cudaMalloc((void **)&d_b, N_BYTES));
CUDA_CHECK(cudaMalloc((void **)&d_c, N_BYTES));
cudaEvent_t ev_start, ev_stop;
CUDA_CHECK(cudaEventCreate(&ev_start));
CUDA_CHECK(cudaEventCreate(&ev_stop));
/* ============ 版本 1: 串行(单流) ============ */
printf("\n===== 实验 B: 串行版(单流, 一次拷完 -> kernel -> 拷回) =====\n");
CUDA_CHECK(cudaMemset(d_c, 0, N_BYTES));
CUDA_CHECK(cudaEventRecord(ev_start));
CUDA_CHECK(cudaMemcpy(d_a, h_a, N_BYTES, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b, N_BYTES, cudaMemcpyHostToDevice));
heavyVecAdd<<<GRID_SIZE, BLOCK_SIZE>>>(d_a, d_b, d_c, N_ELEMENTS);
CUDA_CHECK(cudaMemcpy(h_c, d_c, N_BYTES, cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaEventRecord(ev_stop));
CUDA_CHECK(cudaEventSynchronize(ev_stop));
float ms_serial = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms_serial, ev_start, ev_stop));
printf(" 总耗时: %.3f ms\n", ms_serial);
const int ok_serial = verify(h_c, h_a, h_b, N_ELEMENTS, "serial");
/* ============ 版本 2: 4 段 x 4 流流水线 ============ */
printf("\n===== 实验 C: 4 段 x 4 stream 流水线版 =====\n");
cudaStream_t streams[N_SEGMENTS];
for (int s = 0; s < N_SEGMENTS; ++s) {
/* cudaStreamNonBlocking: 不与 legacy default stream 做隐式同步 */
CUDA_CHECK(cudaStreamCreateWithFlags(&streams[s], cudaStreamNonBlocking));
}
const size_t seg_bytes = N_BYTES / (size_t)N_SEGMENTS;
CUDA_CHECK(cudaMemset(d_c, 0, N_BYTES));
CUDA_CHECK(cudaEventRecord(ev_start));
for (int s = 0; s < N_SEGMENTS; ++s) {
const size_t off = (size_t)s * seg_bytes;
const char *ph_a = (const char *)h_a + off;
const char *ph_b = (const char *)h_b + off;
char *ph_c = (char *)h_c2 + off;
const char *pd_a = (const char *)d_a + off;
const char *pd_b = (const char *)d_b + off;
char *pd_c = (char *)d_c + off;
/* 同一个 stream 内: H2D -> kernel -> D2H 严格按序 */
CUDA_CHECK(cudaMemcpyAsync((void *)pd_a, ph_a, seg_bytes,
cudaMemcpyHostToDevice, streams[s]));
CUDA_CHECK(cudaMemcpyAsync((void *)pd_b, ph_b, seg_bytes,
cudaMemcpyHostToDevice, streams[s]));
heavyVecAdd<<<GRID_SIZE / N_SEGMENTS, BLOCK_SIZE, 0, streams[s]>>>(
(const float *)pd_a, (const float *)pd_b, (float *)pd_c, SEG_ELEMENTS);
CUDA_CHECK(cudaMemcpyAsync(ph_c, pd_c, seg_bytes,
cudaMemcpyDeviceToHost, streams[s]));
}
for (int s = 0; s < N_SEGMENTS; ++s) {
CUDA_CHECK(cudaStreamSynchronize(streams[s]));
}
CUDA_CHECK(cudaEventRecord(ev_stop));
CUDA_CHECK(cudaEventSynchronize(ev_stop));
float ms_pipe = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms_pipe, ev_start, ev_stop));
printf(" 总耗时: %.3f ms\n", ms_pipe);
const int ok_pipe = verify(h_c2, h_a, h_b, N_ELEMENTS, "pipeline");
/* ============ 结果汇总 ============ */
int identical = 1;
for (long i = 0; i < N_ELEMENTS; i += 997) {
if (h_c[i] != h_c2[i]) { identical = 0; break; }
}
printf("\n===== 结果汇总 =====\n");
printf(" 串行版 : %8.3f ms\n", ms_serial);
printf(" 流水线版(4 段) : %8.3f ms\n", ms_pipe);
printf(" 加速比 : %8.3fx (理论上限 max(T_transfer,T_compute)/(T_transfer+T_compute))\n",
ms_serial / ms_pipe);
printf(" 两版本输出一致 : %s\n", identical ? "YES" : "NO");
printf(" 验证结果 : serial=%s pipeline=%s\n",
ok_serial ? "PASS" : "FAIL", ok_pipe ? "PASS" : "FAIL");
/* ============ 释放资源 ============ */
for (int s = 0; s < N_SEGMENTS; ++s) {
CUDA_CHECK(cudaStreamDestroy(streams[s]));
}
CUDA_CHECK(cudaEventDestroy(ev_start));
CUDA_CHECK(cudaEventDestroy(ev_stop));
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_c));
CUDA_CHECK(cudaFreeHost(h_a));
CUDA_CHECK(cudaFreeHost(h_b));
CUDA_CHECK(cudaFreeHost(h_c));
CUDA_CHECK(cudaFreeHost(h_c2));
CUDA_CHECK(cudaDeviceReset());
return 0;
}
【代码做什么？】
- 设备查询：
cudaGetDeviceProperties读出deviceOverlap（是否支持”kernel 与拷贝同时进行”）与asyncEngineCount（异步拷贝引擎数）。A100 上这两个值分别是 1 和 3（1 个 H2D 引擎 + 1 个 D2H 引擎 + 1 个双向引擎），照讲义的说法”Most CUDA devices support device overlap”，但必须先查询再决定要不要写多流代码。 - 实验 A：分别用
malloc和cudaHostAlloc分配 64 MB 主机缓冲，做 8 次 H2D 拷贝，用cudaEvent计时并折算 GB/s。memset先把页真正触碰一遍，避免首次 page fault 污染计时。 - 数据准备：主机端 A、B、C 全部用 pinned 分配，因为
cudaMemcpyAsync只有在主机内存 pinned 时才真正异步。设备端用cudaMalloc分配三块 64 MB。 - 串行版：
cudaMemcpy两次 H2D（同步）→ kernel 启动 → 一次 D2H（同步）。时间戳用cudaEventRecord打在默认流上，cudaEventSynchronize后取cudaEventElapsedTime。 - 流水线版：创建 4 个
cudaStreamNonBlocking流；对每个段s，用字节偏移切出 A、B、C 的对应区间，在streams[s]上依次投递cudaMemcpyAsync(H2D A)、cudaMemcpyAsync(H2D B)、heavyVecAdd、cudaMemcpyAsync(D2H C)。同一流内四个操作严格按序（这是”拷贝没完成就不启动 kernel”的保证），不同流之间无顺序约束（这是重叠的来源）。最后逐个cudaStreamSynchronize。 - 验证：对一个 16 点的抽样集，用与 kernel 相同的运算顺序做主机端参考计算，比较相对误差 < 1e-3；再抽样比较串行版与流水线版的输出是否逐位相同（两者都应该是确定性的）。
- 释放：销毁 stream、event，释放 device 与 host 内存（
cudaFreeHost释放 pinned 内存），最后cudaDeviceReset()。
【并行机制与硬件映射解说】
- 线程映射：
heavyVecAdd用grid-stride loop。串行版 grid = 2048 block × 256 线程 = 524,288 线程，每线程处理16.7M / 524288 = 32个元素；流水线版每段 grid = 512 block × 256 = 131,072 线程，每段4.19M / 131072 = 32个元素。注意两种版本的每线程工作量相同，这是公平对比的前提——否则段变小会让每线程工作量骤降，kernel 反而变慢。 - warp 与调度：256 线程/block = 8 个 warp，每个 warp 恰好是 32 个连续 threadIdx。A100 每个 SM 有 4 个 warp scheduler（每个 sub-partition 一个），每个 scheduler 每周期可以发射 1 条指令；warp 总数达到 64/SM 时，每个 scheduler 有 16 个 warp 轮转，FMA 流水线（延迟约 4 个周期）能被完全填满。
- 驻留量（residency）：A100 每 SM 上限 2048 线程 / 64 warp。这个 kernel 无共享内存、无大数组，寄存器消耗约 24~32 个/线程。若用 32 寄存器：每 block 需要
32 × 256 = 8192寄存器，65536 / 8192 = 8个 block，8 × 256 = 2048线程 → 正好打满 100% 占用率（64/64 warp）。这是”纯计算 kernel 无额外资源约束”的典型情形；流水线版里 4 个不同 stream 的 block 可以同时驻留在同一个 SM 上，它们互相之间没有任何同步依赖，只争抢 register file 与 FMA 单元。 - 共享内存与 bank conflict：本 kernel 不使用
__shared__，因此不涉及 bank conflict；subTile类分析留给代码示例 3。 - 全局内存访问是否合并：
a[i]中i由blockIdx.x * 256 + threadIdx.x得到，warp 内 32 个线程的i连续 → 32 × 4 B = 128 字节，恰好一条 cache line、一次 transaction，完美合并（coalesced）。b[i]、c[i]同理。stride 为gridDim.x * blockDim.x，下一次迭代的访存同样是连续的 128 字节段——grid-stride loop 不会破坏合并性，这是它比”每线程处理相邻多个元素”更优的关键原因。 - warp 发散：不存在发散。循环边界
i < n在最后一个 warp 上可能出现部分线程退出，但相对于 4096 次 FMA 的循环体，开销可忽略。 - 拷贝引擎与 SM 的并行：这是本程序的核心。H2D 拷贝由 Copy Engine（DMA 引擎） 执行，它直接读写显存，完全不占用 SM 的执行单元；kernel 由 SM 执行，完全不占用 PCIe 带宽。二者硬件上互相独立，所以理论上可以 100% 并行——软件上唯一的障碍是”依赖关系”：kernel 必须等它读的那段数据到位，D2H 必须等 kernel 算完。把数据切段，就是把这些依赖”切片化”，让不同段的依赖在不同的时间点满足。
【性能优化分析】
用算术强度 + Roofline判定瓶颈，再用占用率确认延迟能被隐藏。
- 算术强度：
每个元素的 FLOP = 2 × WORK_PER_ELEMENT + 1（最后的加法） = 2 × 4096 + 1 = 8193 FLOP 每个元素的访存 = 读 a[i] 4 B + 读 b[i] 4 B + 写 c[i] 4 B = 12 B 算术强度 AI = 8193 / 12 ≈ 683 FLOP/BA100 的 Roofline 拐点（FP32，不考虑 Tensor Core）：
AI_knee = FP32 峰值 / 显存带宽 = 19.5e12 / 1555e9 ≈ 12.5 FLOP/B683 ≫ 12.5，所以这个 kernel 是计算受限（compute-bound）。可达性能上限 =min(19.5 TFLOPS, 1555 × 683 = 1062 TFLOPS) = 19.5 TFLOPS，实际按 75% 效率（FMA 混合指令发射、循环开销）取 14.6 TFLOPS：FLOP 总量 = 16.7e6 × 8193 ≈ 1.373e11 FLOP = 137.3 GFLOP kernel 时间 ≈ 137.3e9 / 14.6e12 ≈ 9.4 ms这正好落在我设计的目标上（与 8.7 ms 的 PCIe 传输时间同量级）。反过来看，如果把
WORK_PER_ELEMENT设为 0（纯 vecAdd）：AI = 2/12 = 0.167 FLOP/B，远小于 12.5，变成显存带宽受限，kernel 时间降到384 MB / 1555 GB/s ≈ 0.25 ms，此时 PCIe 传输占总时间的 98.5%——此时唯一有意义的优化是减少传输量，而不是优化 kernel。这两种极端正是 Roofline 模型最有教学价值的地方。 - 占用率：
每 block 寄存器 = 32 regs/thread × 256 threads = 8192 寄存器限制的 block 数 = 65536 / 8192 = 8 线程限制的 block 数 = 2048 / 256 = 8 → 实际驻留 8 block = 2048 线程 = 64 warp/SM 占用率 = 64 / 64 = 100%在 100% 占用率下，FMA 延迟（约 4 周期）需要 4 条独立指令才能填满流水线，而我们的 4 条累加链恰好提供 ILP=4，每周期可以稳定发射 1 条 FMA/warp——这就是单 SM 达到 64 FMA/cycle 的条件。
- 总时间模型与加速比推导：
设 T_in = H2D 总时间, T_c = 计算总时间, T_out = D2H 总时间, N = 段数 串行: T_serial = T_in + T_c + T_out 流水线: T_pipe = (T_in + T_c + T_out)/N + (N-1) × max(T_in,T_c,T_out)/N 代入实测(Gen4 x16, pinned, 22 GB/s): A(64MB) H2D = 2.9 ms, B(64MB) H2D = 2.9 ms → T_in = 5.8 ms T_c = 9.4 ms C(64MB) D2H = 2.9 ms → T_out = 2.9 ms T_serial = 5.8 + 9.4 + 2.9 = 18.1 ms N = 4 : T_pipe = 18.1/4 + 3 × 9.4/4 = 4.53 + 7.05 = 11.6 ms → 1.56x N = 8 : T_pipe = 2.26 + 7 × 1.175 = 2.26 + 8.23 = 10.5 ms → 1.72x N = 16 : T_pipe = 1.13 + 15 × 0.588 = 1.13 + 8.82 = 9.95 ms → 1.82x N → ∞ : T_pipe = max(5.8, 9.4, 2.9) = 9.4 ms → 1.93x（理想上限）这正是讲义”execution time never reaches zero”的定量体现：段越小，稳态越接近瓶颈阶段，但填充（fill）与排空（drain）的边界成本占比不变，而 kernel 启动固定开销、SM 空闲、负载不均的代价随段数线性增长，所以总时间会先降后升。讲义例题的最优点是 4~5 段（19 s），这里的理论最优在 N ≈ 8~16，实际最佳值必须实测扫描 N 得到——这就是讲义”Best size likely to depend on GPU”的含义。
- Pageable vs Pinned 的定量解释：pageable 的路径是”主机内存 → 驱动暂存的 pinned 缓冲 → 设备”，第一跳由 CPU 执行（受主机内存带宽与单核拷贝吞吐限制，通常 10~14 GB/s），第二跳由 DMA 执行；两跳串行叠加，因此实测 ≈ 10~12 GB/s。pinned 只有一跳 DMA，实测 ≈ 22~25 GB/s，比值约 2.0~2.4 倍，与讲义”cudaMemcpy should be about 2X faster with pinned memory”一致。
- 可执行的优化方向：① 若
T_c远小于T_in + T_out（低算术强度问题），唯一的出路是减少传输量——在设备上生成数据、压缩、改成多 kernel 流水、或者把算法改写成数据留在显存里的形式（这正是 tiling 对 GEMM 的意义）；② 若T_c ≫ T_in + T_out（如本例），多流重叠的收益上限只有(T_in+T_out)/T，本例是8.7/18.1 = 48%，实现到 1.56x 已接近最大值；③ 段数不要超过”能填满流水线的最小值”，用实测扫描取代拍脑袋；④ 用cudaMemcpyAsync时必须保证主机内存 pinned，否则流水线退化成串行。
代码示例 2：CUDA 流与事件的完整用法（含跨流同步与默认流陷阱）
// 文件: stream_events.cu
// 编译: nvcc -O3 -arch=sm_80 stream_events.cu -o stream_events
// 运行: ./stream_events
//
// 演示完整的 stream / event 生命周期:
// cudaStreamCreateWithFlags, cudaMemcpyAsync, kernel<<<g,b,0,stream>>>,
// cudaEventCreateWithFlags, cudaEventRecord, cudaStreamWaitEvent,
// cudaStreamSynchronize, cudaEventElapsedTime
// 并逐阶段打印耗时, 再演示"legacy default stream 的隐式同步"陷阱与正确做法。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define N (16 * 1024 * 1024) /* 64 MB / 数组 */
#define NBYTES ((size_t)N * sizeof(float))
#define NSEG 2
#define SEG (N / NSEG)
/* 用 SEG / BLK 计算 grid, 保证整除, 避免整数除法截断导致覆盖不足 */
#define BLK 256
#define GRID_SEG (SEG / BLK)
__global__ void vecAddSeg(const float * __restrict__ a,
const float * __restrict__ b,
float * __restrict__ c, int n)
{
const int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) {
c[i] = a[i] + b[i];
}
}
/* 轻量级 kernel: 用固定次数的 FMA 制造与传输可比的计算时间 */
__global__ void workSeg(const float * __restrict__ a,
const float * __restrict__ b,
float * __restrict__ c, int n)
{
const int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i >= n) return;
float acc = a[i] + b[i];
#pragma unroll 8
for (int k = 0; k < 2048; ++k) {
acc = fmaf(acc, 1.0000001f, 0.0000001f);
}
c[i] = acc;
}
int main(void)
{
float *h_a = NULL, *h_b = NULL, *h_c = NULL;
float *d_a = NULL, *d_b = NULL, *d_c = NULL;
CUDA_CHECK(cudaHostAlloc((void **)&h_a, NBYTES, cudaHostAllocDefault));
CUDA_CHECK(cudaHostAlloc((void **)&h_b, NBYTES, cudaHostAllocDefault));
CUDA_CHECK(cudaHostAlloc((void **)&h_c, NBYTES, cudaHostAllocDefault));
for (int i = 0; i < N; ++i) {
h_a[i] = 1.0f + (float)(i % 7);
h_b[i] = 2.0f + (float)(i % 5);
}
CUDA_CHECK(cudaMalloc((void **)&d_a, NBYTES));
CUDA_CHECK(cudaMalloc((void **)&d_b, NBYTES));
CUDA_CHECK(cudaMalloc((void **)&d_c, NBYTES));
/* ---------------- 1. 创建流与事件 ---------------- */
cudaStream_t s0, s1;
CUDA_CHECK(cudaStreamCreateWithFlags(&s0, cudaStreamNonBlocking));
CUDA_CHECK(cudaStreamCreateWithFlags(&s1, cudaStreamNonBlocking));
/* 纯同步用途的事件: disableTiming 开销更低 */
cudaEvent_t ev_h2d_done, ev_calc_done;
CUDA_CHECK(cudaEventCreateWithFlags(&ev_h2d_done, cudaEventDisableTiming));
CUDA_CHECK(cudaEventCreateWithFlags(&ev_calc_done, cudaEventDisableTiming));
/* 计时用途的事件: 必须允许计时 */
cudaEvent_t t_copyA, t_copyB, t_kernel, t_copyC, t_all0, t_all1;
CUDA_CHECK(cudaEventCreate(&t_copyA));
CUDA_CHECK(cudaEventCreate(&t_copyB));
CUDA_CHECK(cudaEventCreate(&t_kernel));
CUDA_CHECK(cudaEventCreate(&t_copyC));
CUDA_CHECK(cudaEventCreate(&t_all0));
CUDA_CHECK(cudaEventCreate(&t_all1));
const size_t segoff = (size_t)SEG * sizeof(float);
/* ---------------- 2. 逐阶段计时的流水线 ---------------- */
printf("===== 阶段耗时分解 (segment 0, 2 段流水) =====\n");
CUDA_CHECK(cudaEventRecord(t_all0));
for (int s = 0; s < NSEG; ++s) {
cudaStream_t st = (s == 0) ? s0 : s1;
const size_t off = (size_t)s * segoff;
const char *pa = (const char *)h_a + off;
const char *pb = (const char *)h_b + off;
char *pc = (char *)h_c + off;
const char *da = (const char *)d_a + off;
const char *db = (const char *)d_b + off;
char *dc = (char *)d_c + off;
if (s == 0) CUDA_CHECK(cudaEventRecord(t_copyA, st));
CUDA_CHECK(cudaMemcpyAsync((void *)da, pa, segoff, cudaMemcpyHostToDevice, st));
CUDA_CHECK(cudaMemcpyAsync((void *)db, pb, segoff, cudaMemcpyHostToDevice, st));
if (s == 0) {
CUDA_CHECK(cudaEventRecord(t_copyB, st));
/* 在流中记录一个"输入就绪"事件, 供别的流等待 */
CUDA_CHECK(cudaEventRecord(ev_h2d_done, st));
}
workSeg<<<GRID_SEG, BLK, 0, st>>>((const float *)da, (const float *)db,
(float *)dc, SEG);
if (s == 0) {
CUDA_CHECK(cudaEventRecord(t_kernel, st));
CUDA_CHECK(cudaEventRecord(ev_calc_done, st));
}
CUDA_CHECK(cudaMemcpyAsync(pc, dc, segoff, cudaMemcpyDeviceToHost, st));
if (s == 0) CUDA_CHECK(cudaEventRecord(t_copyC, st));
}
CUDA_CHECK(cudaStreamSynchronize(s0));
CUDA_CHECK(cudaStreamSynchronize(s1));
CUDA_CHECK(cudaEventRecord(t_all1));
CUDA_CHECK(cudaEventSynchronize(t_all1));
float ms;
CUDA_CHECK(cudaEventElapsedTime(&ms, t_copyA, t_copyB));
printf(" H2D A 拷贝 (32 MB) : %7.3f ms (%.1f GB/s)\n",
ms, 32.0e6 / (ms / 1000.0) / 1e9);
CUDA_CHECK(cudaEventElapsedTime(&ms, t_copyB, t_kernel));
printf(" H2D B 拷贝 (32 MB) : %7.3f ms (%.1f GB/s)\n",
ms, 32.0e6 / (ms / 1000.0) / 1e9);
CUDA_CHECK(cudaEventElapsedTime(&ms, t_kernel, t_copyC));
printf(" kernel 计算 (8.4M 元素 x 2048 FMA): %7.3f ms\n", ms);
CUDA_CHECK(cudaEventElapsedTime(&ms, t_all0, t_all1));
printf(" 两段流水总时间 : %7.3f ms\n", ms);
/* ---------------- 3. 跨流同步: cudaStreamWaitEvent ---------------- */
printf("\n===== cudaStreamWaitEvent 跨流同步 =====\n");
CUDA_CHECK(cudaStreamSynchronize(s0));
CUDA_CHECK(cudaStreamSynchronize(s1));
/* 生产者-消费者: s0 负责 H2D, s1 负责计算, 用事件建立跨流依赖。
注意 cudaStreamWaitEvent 不会阻塞主机线程, 主机可以立刻去做别的事。 */
CUDA_CHECK(cudaMemcpyAsync(d_a, h_a, NBYTES, cudaMemcpyHostToDevice, s0));
CUDA_CHECK(cudaMemcpyAsync(d_b, h_b, NBYTES, cudaMemcpyHostToDevice, s0));
CUDA_CHECK(cudaEventRecord(ev_h2d_done, s0));
CUDA_CHECK(cudaStreamWaitEvent(s1, ev_h2d_done, 0)); /* s1 等 s0 的拷贝 */
vecAddSeg<<<N / BLK, BLK, 0, s1>>>(d_a, d_b, d_c, N);
CUDA_CHECK(cudaStreamSynchronize(s1)); /* 只等 s1, 不等整个设备 */
printf(" s1 已在 s0 的 H2D 完成后被唤醒 (主机线程全程未被阻塞)\n");
/* ---------------- 4. legacy default stream 的隐式同步陷阱 ---------------- */
printf("\n===== legacy default stream 的隐式同步 =====\n");
/* 下面这次 kernel 使用默认流(第 4 个参数省略), 它会等待所有阻塞流中的操作,
从而把本可重叠的拷贝与计算强行串行化 */
vecAddSeg<<<N / BLK, BLK>>>(d_a, d_b, d_c, N);
CUDA_CHECK(cudaMemcpy(h_c, d_c, NBYTES, cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaDeviceSynchronize());
printf(" 默认流 kernel 已与 s0/s1 隐式串行化完成\n");
printf(" 规避方式: 全部使用 cudaStreamNonBlocking 流, 或编译时加\n");
printf(" --default-stream per-thread\n");
/* ---------------- 5. 结果验证 ---------------- */
int bad = 0;
int checked = 0;
for (int i = 0; i < N; i += 1000003) {
const float expect = h_a[i] + h_b[i];
if (fabsf(h_c[i] - expect) > 1.0e-4f * (fabsf(expect) + 1.0f)) ++bad;
++checked;
}
printf("\n 抽样验证 (最后一步默认流 kernel 的结果): %s (%d 个元素)\n",
(bad == 0) ? "PASS" : "FAIL", checked);
/* ---------------- 6. 资源释放 ---------------- */
CUDA_CHECK(cudaEventDestroy(ev_h2d_done));
CUDA_CHECK(cudaEventDestroy(ev_calc_done));
CUDA_CHECK(cudaEventDestroy(t_copyA));
CUDA_CHECK(cudaEventDestroy(t_copyB));
CUDA_CHECK(cudaEventDestroy(t_kernel));
CUDA_CHECK(cudaEventDestroy(t_copyC));
CUDA_CHECK(cudaEventDestroy(t_all0));
CUDA_CHECK(cudaEventDestroy(t_all1));
CUDA_CHECK(cudaStreamDestroy(s0));
CUDA_CHECK(cudaStreamDestroy(s1));
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_c));
CUDA_CHECK(cudaFreeHost(h_a));
CUDA_CHECK(cudaFreeHost(h_b));
CUDA_CHECK(cudaFreeHost(h_c));
return 0;
}
【代码做什么？】
- 创建对象：
cudaStreamCreateWithFlags(&s0, cudaStreamNonBlocking)创建两个非阻塞流——非阻塞流不与 legacy default stream 隐式同步。事件分两类创建：跨流同步用的事件加cudaEventDisableTiming（不写时间戳，开销更低）；计时用的事件必须允许计时，否则cudaEventElapsedTime会返回错误。 - 逐阶段计时：只在 segment 0 上插入计时事件（
t_copyA、t_copyB、t_kernel、t_copyC），这样时间戳反映的是该段的各阶段耗时，而不是整条流水线的并发行为；总时间用t_all0/t_all1单独测。 - 跨流同步：
cudaEventRecord(ev_h2d_done, s0)在 s0 里打一个标记，cudaStreamWaitEvent(s1, ev_h2d_done, 0)让 s1 等到这个标记完成再继续。关键点：这个调用不阻塞主机线程，主机可以立刻去干别的活。这是 CUDA 里表达”生产者-消费者”依赖的标准手段（第 3 个参数是保留位，必须传 0）。 - 默认流陷阱演示：直接写
vecAddSeg<<<N/BLK, BLK>>>(d_a, d_b, d_c, N)（省略流参数）会把 kernel 投到 legacy default stream。默认流与所有”阻塞流”之间都有隐式同步关系，于是这次 kernel 必须等 s0、s1 里的所有操作完成，重叠被破坏。修复方法是全部使用非阻塞流，或者编译时加--default-stream per-thread（让每个主机线程有自己的默认流）。 - 验证：抽样检查输出与主机端按相同运算顺序算出的参考值是否一致。
- 释放：所有 stream、event、设备内存、固定内存全部显式释放。
【并行机制与硬件映射解说】
- 两个 kernel 的定位差异：
vecAddSeg是纯访存 kernel（每线程 1 次读 A、1 次读 B、1 次写 C，AI = 0.167 FLOP/B），用来演示”传输主导”；workSeg每线程做 2048 次 FMA（AI ≈ 4097 FLOP / 12 B ≈ 341 FLOP/B），用来演示”计算主导”。在硬件上，前者的瓶颈是 L2/DRAM 带宽与 PCIe，后者的瓶颈是 SM 的 FMA 单元。 - warp 与占用率：
blockDim = 256→ 8 warp/block。workSeg有 2048 次 FMA 的长循环，寄存器需求约 16 个（只有 acc、i、n 等少数活跃变量），因此寄存器不是限制项：65536 / (16 × 256) = 16个 block，但线程上限是2048/256 = 8个 block，所以实际驻留 8 block = 64 warp = 100% 占用率。要注意workSeg是单一累加链（acc串行依赖），每线程每 4 个周期才能完成 1 次 FMA，即 ILP = 1；在 100% 占用率下，warps 足够多，FMA 吞吐仍可由不同 warp 填满（延迟隐藏靠 warp 级并行而不是指令级并行）——这也解释了为什么”占用率”在延迟隐藏上如此重要：ILP 不足时用 TLP 补。 - 访存合并：两个 kernel 都用
blockIdx.x * blockDim.x + threadIdx.x的连续映射，warp 内 32 个线程访问 32 个连续 float = 128 字节 = 1 条 cache line = 1 次 transaction，完美合并。若改成c[2*i]（stride-2 访问），一个 warp 会跨 256 字节，需要 2 次 transaction，效率立刻掉到 50%。 - bank conflict：本示例不使用共享内存，无 bank conflict。但要注意
cudaHostAlloc得到的 pinned 内存在主机侧没有 bank 概念；若在设备端做共享内存归约，则必须遵守”32 bank × 4 字节、步长为 32 的倍数会产生 32 路冲突”的规则。 - 事件在硬件上的含义：
cudaEventRecord实际上是往对应流的命令缓冲里写入一条”时间戳/标记”命令；GPU 的引擎按序执行到该命令时把时间戳写入内存。因此cudaEventElapsedTime测的是设备端 GPU 时间线，而不是主机 wall-clock——这正是它能精确测量重叠效果的原因（主机端计时会把主机提交命令的时间也算进去）。 - 默认流的同步语义：legacy default stream 是”同步流”，规则是：默认流中的命令会等待此前所有阻塞流中的命令完成，且此后所有阻塞流中的命令会等待默认流完成。硬件上这由驱动的命令提交层实现（把依赖插入各引擎队列）。代价是跨流并行度被彻底摧毁；在 Fermi 上还叠加了”不同流的 kernel 复用同一条硬件工作队列、只在流边界重叠”的限制，Kepler 之后每个流有独立硬件队列才真正解决了这个问题。
【性能优化分析】
- 算术强度与 Roofline（两个 kernel 分别代入）：
vecAddSeg : FLOP/元素 = 1, 访存 12 B/元素 → AI = 1/12 = 0.083 FLOP/B A100 拐点 = 19.5e12 / 1555e9 = 12.5 FLOP/B 0.083 ≪ 12.5 ⇒ 显存带宽受限 可达上限 = 1555 GB/s × 0.083 FLOP/B = 129 GFLOP/s（远低于算力峰值） 实际时间 = 384 MB / 1555 GB/s ≈ 0.25 ms workSeg : FLOP/元素 = 2×2048 = 4096, 访存 12 B/元素 → AI = 341 FLOP/B 341 ≫ 12.5 ⇒ 计算受限，可达上限 = 19.5 TFLOPS（按 75% 取 14.6 TFLOPS）在 64 MB 数组、2 段的配置下：
workSeg时间 ≈2×(8.4e6 × 4096) / 14.6e12 ≈ 4.7 ms，而 H2D 传输 =128 MB / 22 GB/s ≈ 5.8 ms，二者同量级，重叠收益显著。 - 占用率：
workSeg: 寄存器 ≈ 16/线程 寄存器限制 = 65536 / (16 × 256) = 16 block 线程限制 = 2048 / 256 = 8 block → 驻留 8 block，占用率 = 64/64 = 100% vecAddSeg: 寄存器 ≈ 12/线程，同样 8 block，100% - 总时间模型（2 段流水）：
T_in = 5.8 ms (H2D), T_c = 4.7 ms, T_out = 2.9 ms (D2H) 串行 = 5.8 + 4.7 + 2.9 = 13.4 ms N = 2 : 13.4/2 + 1 × max(5.8, 4.7, 2.9)/2 = 6.7 + 2.9 = 9.6 ms → 1.40x N = 4 : 13.4/4 + 3 × 5.8/4 = 3.35 + 4.35 = 7.7 ms → 1.74x N → ∞ : max(5.8, 4.7, 2.9) = 5.8 ms → 2.31x这里瓶颈阶段变成了 H2D 拷贝（5.8 ms > 4.7 ms），说明段数继续增加只能把时间推到 5.8 ms，再想加速必须提高 PCIe 有效带宽或减少要传输的字节数。这是判断”优化到头了没有”的标准方法：先算
max(各阶段)的极限，再看当前离极限多远。 - 可执行的优化方向：① 用
cudaEventDisableTiming的事件做同步，减少时间戳写入带来的设备端开销；② 用cudaStreamWaitEvent精确表达依赖，避免用cudaDeviceSynchronize把所有流一刀切；③ 确保所有拷贝都发生在 pinned 内存上，否则cudaMemcpyAsync会退化；④ 用--default-stream per-thread或全非阻塞流规避默认流同步；⑤ 用cudaStreamAddCallback或（较新）cudaLaunchHostFunc在流中插入主机回调做进度报告，而不要用忙等。
代码示例 3：Tensor Core（WMMA）矩阵乘法 vs 普通 CUDA core tiled 矩阵乘法
// 文件: tensor_core_matmul.cu
// 编译: nvcc -O3 -arch=sm_80 tensor_core_matmul.cu -o tensor_core_matmul
// 运行: ./tensor_core_matmul
//
// 同一个 M=N=K=4096 的矩阵乘法, 用两种 kernel 实现并对比:
// (1) matmulTiledFP32 : 经典 shared memory tiling + CUDA core FMA (FP32)
// (2) matmulWMMA : nvcuda::wmma 的 16x16x16 fragment, FP16 输入 + FP32 累加
// 两者都与 CPU 参考实现抽样对比, 并打印 GFLOP/s 与 TFLOPS。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#include <cuda_fp16.h>
#include <mma.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_WIDTH 16 /* CUDA core tiled 版本的 tile 边长 */
#define WMMA_M 16
#define WMMA_N 16
#define WMMA_K 16
#define WARPS_PER_BLOCK 4 /* 每 block 4 个 warp = 128 线程 */
using namespace nvcuda;
/* ==================== kernel 1: CUDA core, FP32, tiled ==================== */
__global__ void matmulTiledFP32(const float * __restrict__ A,
const float * __restrict__ B,
float * __restrict__ C,
int Width)
{
__shared__ float subA[TILE_WIDTH][TILE_WIDTH];
__shared__ float subB[TILE_WIDTH][TILE_WIDTH];
const int bx = blockIdx.x;
const int by = blockIdx.y;
const int tx = threadIdx.x;
const int ty = threadIdx.y;
const int Row = by * TILE_WIDTH + ty;
const int Col = bx * TILE_WIDTH + tx;
float Pvalue = 0.0f;
for (int q = 0; q < Width / TILE_WIDTH; ++q) {
/* 协作载入: 每个线程搬一个元素, 一个 block 搬一个 TILE x TILE 的块 */
subA[ty][tx] = A[Row * Width + (q * TILE_WIDTH + tx)];
subB[ty][tx] = B[(q * TILE_WIDTH + ty) * Width + Col];
__syncthreads();
#pragma unroll
for (int k = 0; k < TILE_WIDTH; ++k) {
Pvalue += subA[ty][k] * subB[k][tx];
}
/* 必须再同步一次: 防止快的线程先写下一轮 tile, 覆盖别人还在读的数据 */
__syncthreads();
}
C[Row * Width + Col] = Pvalue;
}
/* ==================== kernel 2: Tensor Core, FP16 输入 / FP32 累加 ==================== */
__global__ void matmulWMMA(const __half * __restrict__ A,
const __half * __restrict__ B,
float * __restrict__ C,
int M, int N, int K)
{
const int warpId = threadIdx.x >> 5; /* 0 .. WARPS_PER_BLOCK-1 */
const int tileRow = blockIdx.y * (WMMA_M * WARPS_PER_BLOCK) + warpId * WMMA_M;
const int tileCol = blockIdx.x * WMMA_N;
if (tileRow >= M || tileCol >= N) return;
/* 累加器: 16x16 个 fp32, 由 32 个线程分摊持有 (每线程 8 个) */
wmma::fragment<wmma::accumulator, WMMA_M, WMMA_N, WMMA_K, float> acc;
wmma::fill_fragment(acc, 0.0f);
for (int k = 0; k < K; k += WMMA_K) {
wmma::fragment<wmma::matrix_a, WMMA_M, WMMA_N, WMMA_K, __half, wmma::row_major> a_frag;
wmma::fragment<wmma::matrix_b, WMMA_M, WMMA_N, WMMA_K, __half, wmma::row_major> b_frag;
/* ldm = 主维长度; 16 位类型要求 ldm 是 8 的倍数(16 字节对齐) */
wmma::load_matrix_sync(a_frag, A + (size_t)tileRow * K + k, K);
wmma::load_matrix_sync(b_frag, B + (size_t)k * N + tileCol, N);
/* D = A x B + C, 一条指令完成 16x16x16 = 4096 次乘加 */
wmma::mma_sync(acc, a_frag, b_frag, acc);
}
wmma::store_matrix_sync(C + (size_t)tileRow * N + tileCol, acc, N,
wmma::mem_row_major);
}
/* ==================== kernel 3: FP32 -> FP16 转换 ==================== */
__global__ void f32_to_f16(const float * __restrict__ in,
__half * __restrict__ out, size_t n)
{
const size_t i = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) {
out[i] = __float2half(in[i]);
}
}
/* ==================== 主机端工具 ==================== */
static unsigned int g_seed = 12345u;
static float rnd(void)
{
g_seed = g_seed * 1664525u + 1013904223u;
return ((float)(g_seed >> 8) / (float)(1u << 24)) * 2.0f - 1.0f; /* [-1, 1) */
}
static double cpu_dot_fp32(const float *A, const float *B, int N, int K, int i, int j)
{
double s = 0.0;
for (int k = 0; k < K; ++k) {
s += (double)A[(size_t)i * K + k] * (double)B[(size_t)k * N + j];
}
return s;
}
static double cpu_dot_half(const __half *A, const __half *B, int N, int K, int i, int j)
{
double s = 0.0;
for (int k = 0; k < K; ++k) {
s += (double)__half2float(A[(size_t)i * K + k]) *
(double)__half2float(B[(size_t)k * N + j]);
}
return s;
}
static void report(const char *name, float ms, double gflop)
{
const double gflops = gflop / ((double)ms / 1000.0);
printf(" %-34s %9.3f ms %9.1f GFLOP/s (%.2f TFLOPS)\n",
name, ms, gflops, gflops / 1000.0);
}
#define ITERS 5
/* ============================== main ============================== */
int main(void)
{
const int M = 4096;
const int N = 4096;
const int K = 4096;
const size_t szA_f32 = (size_t)M * K * sizeof(float);
const size_t szB_f32 = (size_t)K * N * sizeof(float);
const size_t szC_f32 = (size_t)M * N * sizeof(float);
const size_t szA_f16 = (size_t)M * K * sizeof(__half);
const size_t szB_f16 = (size_t)K * N * sizeof(__half);
const double gflop = 2.0 * (double)M * (double)N * (double)K;
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
printf("设备: %s CC %d.%d SM=%d FP32 峰值约 %.1f TFLOPS\n\n",
prop.name, prop.major, prop.minor, prop.multiProcessorCount,
(double)prop.multiProcessorCount * 128.0 * 2.0 * 1.41e9 / 1e12);
/* ---------------- 主机端数据准备 ---------------- */
float *h_A = (float *)malloc(szA_f32);
float *h_B = (float *)malloc(szB_f32);
float *h_C = (float *)malloc(szC_f32);
float *h_C2 = (float *)malloc(szC_f32);
__half *h_Ah = (__half *)malloc(szA_f16);
__half *h_Bh = (__half *)malloc(szB_f16);
if (h_A == NULL || h_B == NULL || h_C == NULL || h_C2 == NULL ||
h_Ah == NULL || h_Bh == NULL) {
fprintf(stderr, "host malloc failed\n");
exit(EXIT_FAILURE);
}
for (size_t i = 0; i < (size_t)M * K; ++i) {
h_A[i] = rnd();
/* 主机端就做一次 fp16 舍入, 让 CPU 参考与 GPU 看到完全相同的输入 */
h_Ah[i] = __float2half(h_A[i]);
}
for (size_t i = 0; i < (size_t)K * N; ++i) {
h_B[i] = rnd();
h_Bh[i] = __float2half(h_B[i]);
}
/* ---------------- 设备端分配 ---------------- */
float *d_A = NULL, *d_B = NULL, *d_C = NULL, *d_C2 = NULL;
__half *d_Ah = NULL, *d_Bh = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_A, szA_f32));
CUDA_CHECK(cudaMalloc((void **)&d_B, szB_f32));
CUDA_CHECK(cudaMalloc((void **)&d_C, szC_f32));
CUDA_CHECK(cudaMalloc((void **)&d_C2, szC_f32));
CUDA_CHECK(cudaMalloc((void **)&d_Ah, szA_f16));
CUDA_CHECK(cudaMalloc((void **)&d_Bh, szB_f16));
CUDA_CHECK(cudaMemcpy(d_A, h_A, szA_f32, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_B, h_B, szB_f32, cudaMemcpyHostToDevice));
/* 用一次 kernel 做转换, 演示"数据在设备上就地降精度" */
f32_to_f16<<<(unsigned)((szA_f32 / sizeof(float) + 255) / 256), 256>>>(
d_A, d_Ah, (size_t)M * K);
f32_to_f16<<<(unsigned)((szB_f32 / sizeof(float) + 255) / 256), 256>>>(
d_B, d_Bh, (size_t)K * N);
CUDA_CHECK(cudaDeviceSynchronize());
/* ---------------- 计时 ---------------- */
cudaEvent_t e0, e1;
CUDA_CHECK(cudaEventCreate(&e0));
CUDA_CHECK(cudaEventCreate(&e1));
const dim3 blockFP32(TILE_WIDTH, TILE_WIDTH);
const dim3 gridFP32(N / TILE_WIDTH, M / TILE_WIDTH);
const dim3 blockWMMA(32 * WARPS_PER_BLOCK);
const dim3 gridWMMA(N / WMMA_N, M / (WMMA_M * WARPS_PER_BLOCK));
float ms_tiled = 0.0f, ms_wmma = 0.0f, ms_conv = 0.0f;
/* 预热 + 计时: CUDA core FP32 */
matmulTiledFP32<<<gridFP32, blockFP32>>>(d_A, d_B, d_C, N);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(e0));
for (int it = 0; it < ITERS; ++it) {
matmulTiledFP32<<<gridFP32, blockFP32>>>(d_A, d_B, d_C, N);
}
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
CUDA_CHECK(cudaEventElapsedTime(&ms_tiled, e0, e1));
ms_tiled /= (float)ITERS;
/* 预热 + 计时: Tensor Core WMMA */
matmulWMMA<<<gridWMMA, blockWMMA>>>(d_Ah, d_Bh, d_C2, M, N, K);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaEventRecord(e0));
for (int it = 0; it < ITERS; ++it) {
matmulWMMA<<<gridWMMA, blockWMMA>>>(d_Ah, d_Bh, d_C2, M, N, K);
}
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
CUDA_CHECK(cudaEventElapsedTime(&ms_wmma, e0, e1));
ms_wmma /= (float)ITERS;
/* 计时: fp16 转换 kernel 本身的开销(提醒: 转换不是免费的) */
CUDA_CHECK(cudaEventRecord(e0));
f32_to_f16<<<(unsigned)((szA_f32 / sizeof(float) + 255) / 256), 256>>>(
d_A, d_Ah, (size_t)M * K);
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
CUDA_CHECK(cudaEventElapsedTime(&ms_conv, e0, e1));
/* ---------------- 结果 ---------------- */
printf("===== M=N=K=%d, 总运算量 %.2f GFLOP, 各测 %d 次取平均 =====\n",
M, gflop, ITERS);
report("CUDA core FP32 (tiled 16x16)", ms_tiled, gflop);
report("Tensor Core FP16->FP32 (WMMA)", ms_wmma, gflop);
printf(" %-34s %9.3f ms %9.1f GB/s (读写 2 x %.0f MB)\n",
"f32 -> f16 转换 kernel(单次)", ms_conv,
(2.0 * (double)szA_f32) / ((double)ms_conv / 1000.0) / 1e9,
(double)szA_f32 / 1.0e6);
printf(" 加速比 (Tensor Core / CUDA core) : %.2fx\n", ms_tiled / ms_wmma);
printf(" WMMA 相对 FP16 理论峰值 %s\n",
"(A100 = 312 TFLOPS) 的利用率见下方分析");
/* ---------------- 抽样验证 ---------------- */
CUDA_CHECK(cudaMemcpy(h_C, d_C, szC_f32, cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemcpy(h_C2, d_C2, szC_f32, cudaMemcpyDeviceToHost));
int bad_tiled = 0, bad_wmma = 0, nc = 32;
for (int t = 0; t < nc; ++t) {
const int i = (int)((double)t / nc * M);
const int j = (int)((double)((t * 7919) % nc) / nc * N);
const double ref32 = cpu_dot_fp32(h_A, h_B, N, K, i, j);
const double ref16 = cpu_dot_half(h_Ah, h_Bh, N, K, i, j);
if (fabs((double)h_C[(size_t)i * N + j] - ref32) > 1.0e-1 + 1.0e-3 * fabs(ref32))
++bad_tiled;
if (fabs((double)h_C2[(size_t)i * N + j] - ref16) > 1.0e-1 + 1.0e-3 * fabs(ref16))
++bad_wmma;
}
printf("\n 抽样验证 (%d 个 (i,j) 对):\n", nc);
printf(" CUDA core FP32 : %s\n", (bad_tiled == 0) ? "PASS" : "FAIL");
printf(" Tensor Core : %s\n", (bad_wmma == 0) ? "PASS" : "FAIL");
/* ---------------- 释放 ---------------- */
CUDA_CHECK(cudaEventDestroy(e0));
CUDA_CHECK(cudaEventDestroy(e1));
CUDA_CHECK(cudaFree(d_A));
CUDA_CHECK(cudaFree(d_B));
CUDA_CHECK(cudaFree(d_C));
CUDA_CHECK(cudaFree(d_C2));
CUDA_CHECK(cudaFree(d_Ah));
CUDA_CHECK(cudaFree(d_Bh));
free(h_A);
free(h_B);
free(h_C);
free(h_C2);
free(h_Ah);
free(h_Bh);
return 0;
}
【代码做什么？】
- 数据准备：生成
M×K的 A、K×N的 B（元素在[-1,1)），同时在主机端就做完 fp16 舍入（h_Ah[i] = __float2half(h_A[i])），保证 CPU 参考实现与 GPU 看到的是完全相同的输入位模式——这样验证时只剩”累加顺序差异”，容差可以卡得很紧。 - 设备端降精度：
f32_to_f16kernel 在设备上把 fp32 就地转成 fp16，避免把两份数据都从主机搬上去（真实推理场景里，权重通常是预先量化好的，输入则在设备上量化）。 - FP32 tiled kernel：
TILE_WIDTH=16，block 是16×16=256线程。每轮迭代：协作把 A、B 的一个 16×16 tile 搬进__shared__，__syncthreads()，然后每个线程做 16 次乘加（subA[ty][k] * subB[k][tx]），再__syncthreads()后才进入下一轮。外层次数 =Width / TILE_WIDTH。 - WMMA kernel：每个 warp 负责一个 16×16 的输出 tile，
blockIdx.y方向上 4 个 warp 拼成 64 行。K 方向以 16 为步长循环，每步load_matrix_sync载入 A 片与 B 片、mma_sync累加；循环结束后store_matrix_sync把累加器写回。 - 计时：每个 kernel 先跑一次预热（触发 JIT/cache 效应），再用
cudaEvent测 5 次取平均，折算 GFLOP/s 与 TFLOPS。 - 验证：抽样 32 个
(i,j)，用 double 精度在 CPU 上做点积，分别与两个 kernel 的结果比较（容差1e-1 + 1e-3·\|ref\|）。 - 释放：事件、设备内存、主机内存全部释放。
【并行机制与硬件映射解说】
- FP32 tiled kernel 的线程映射：
blockDim = (16,16) = 256线程 = 8 个 warp。warp 的构成规则是”threadIdx.x变化最快”：warp 0 是ty=0,tx=0..15与ty=1,tx=0..15，也就是一个 warp 横跨两行。这一点决定了后面所有 bank/合并分析。 - 共享内存访问与 bank conflict（具体到 bank 编号）：共享内存按 32 个 bank × 4 字节组织，字的地址
w落在 bankw mod 32。subA[ty][k]：地址 =ty*16 + k。warp 内ty ∈ {0,1}、k对所有线程相同，所以 32 个线程只访问 2 个地址（k与16+k），bank 分别是k mod 32与(k+16) mod 32。同一地址的多线程访问由硬件广播（broadcast）完成，不算冲突，所以这是 0 路冲突。subB[k][tx]：地址 =k*16 + tx，tx = 0..15且两组ty共用同一组地址 → 16 个不同字，落在 bank(16k+tx) mod 32，因为tx只覆盖 0..15，实际上落在连续的 16 个 bank上，并且每个地址被 2 个线程请求（广播）。0 路冲突。- 对比：如果把 B 的 tile 转置存放并写成
subB[tx][k]，地址 =tx*16 + k，32 个线程（16 个不同tx）的字分别是k, 16+k, 32+k, 48+k, 64+k，这些字的mod 32只有k mod 32与(k+16) mod 32两个取值 → 16 个字挤在 2 个 bank 上 → 8 路冲突，共享内存吞吐掉到 1/8。修复办法是把数组声明成subB[TILE_WIDTH][TILE_WIDTH + 1]（padding），地址变成tx*17+k，mod 32就散布开了。这正是讲义在 tiling 讲里强调 padding 的原因。
- 全局内存合并访问：
subA[ty][tx] = A[Row*Width + q*TILE_WIDTH + tx]——warp 内tx = 0..15连续、ty取两个值，所以一个 warp 访问两段各 16 个连续 float = 两段 64 字节，不成 128 字节的整数倍 → 需要 2 次 transaction（每次 64 字节），相对理想的 128 字节 transaction，有效率只有 50%。这是TILE_WIDTH=16的固有代价；把 block 改成(32,8)或TILE_WIDTH=32可以让 warp 内的tx覆盖 32 个连续元素，从而拿到 128 字节的完美合并。subB[ty][tx] = B[(q*TILE_WIDTH+ty)*Width + Col]同理：ty在 warp 内只取 2 个值，但Col = bx*16 + tx随tx连续 → 又是两段 64 字节。 - WMMA kernel 的硬件映射：block 128 线程 = 4 warp。每个 warp 的 fragment 分布是：
acc是 16×16 的 fp32 = 256 个元素，平均分给 32 个 lane，每 lane 8 个 fp32（8 个寄存器）；a_frag是 16×16 的__half= 256 个 half，每 lane 8 个 half（4 个寄存器）；b_frag同样 4 个寄存器。加上地址计算与循环变量，寄存器消耗约 40~64 个/线程。- 占用率：以 64 寄存器计，
65536 / (64 × 128) = 8个 block →8 × 4 = 32warp/SM，A100 上限 64 warp → 占用率 50%（寄存器受限）。这个 kernel 靠 ILP（Tensor Core 的长流水）而不是 TLP 隐藏延迟，所以 50% 占用率并不致命。 - 共享内存：WMMA 版本完全不使用 shared memory，
load_matrix_sync直接从全局内存按ldm跨步读取，每个 lane 读到的地址是按硬件规定的 fragment 布局散开的（__half类型下每 lane 一次读 8 个连续 half = 16 字节，即一条扇区），并不是完美合并的 128 字节 transaction。 - warp 发散：无（只有入口的边界检查可能发散，而 M、N 都是 16 的倍数时该分支恒真）。
- 占用率：以 64 寄存器计，
- Tensor Core 在硬件上做了什么：
mma_sync(acc, a, b, acc)在 sm_80 上编译成HMMA.16816.F32系列指令（m16n8k16 在硬件层面执行，编译器把 16×16×16 的 WMMA 拆成若干条）。一个 warp 级 Tensor Core 由 4 个 sub-core 组成，每条指令并行完成大量 fp16 乘加并累加到 fp32。这就是为什么 FP16 输入的吞吐可以是 FP32 的 16 倍：A100 每 SM 每周期 1024 次 fp16 FMA，而 FP32 只有 64 次。
【性能优化分析】
- 算术强度与 Roofline（对 WMMA kernel 逐步推导，这是本讲最重要的定量练习）：
每个 warp 每个 k 步: 载入 a_frag: 16x16 个 __half = 256 x 2 B = 512 B 载入 b_frag: 16x16 个 __half = 256 x 2 B = 512 B 计算量 : 2 x 16 x 16 x 16 = 8192 FLOP AI = 8192 / 1024 = 8 FLOP/B 若所有 fragment 都从 HBM 取: 上限 = 1555 GB/s x 8 FLOP/B = 12.4 TFLOPS 但 grid = (256, 64) 时, 每个 k 步的"唯一数据量"只有: A 的 16 列 x 4096 行 x 2 B = 128 KB B 的 16 行 x 4096 列 x 2 B = 128 KB 合计 256 KB << A100 的 40 MB L2 因此除首次触及外, 全部命中 L2, 真正的瓶颈变成 L2 带宽与 fragment 装载指令吞吐。 每个 k 步的 L2 请求总量 = 65536 个 warp x 1024 B = 64 MB 全程 (K/16 = 256 个 k 步) = 16 GB 以 A100 L2 带宽约 6 TB/s 估算: 16 GB / 6 TB/s ≈ 2.7 ms ← 与实测同量级 而纯计算下界 = 137.4 GFLOP / 312 TFLOPS = 0.44 ms ⇒ 结论: 朴素的"global -> fragment" WMMA kernel 是 L2 带宽受限, 不是 Tensor Core 受限。FP32 tiled kernel 的算术强度：
每个 block 的每个 tile 迭代: 载入 2 x 16 x 16 x 4 B = 2048 B 该 tile 贡献的计算量 = 16 x 16 x 16 x 2 = 8192 FLOP AI = 8192 / 2048 = 4 FLOP/B （乘以 L2 命中带来的复用后, 有效 AI 更高） 从 HBM 看: 每个 block 全程需要 (Width/16) 次 x 2048 B 的 tile 装载 grid = 256 x 256 个 block, 每 block 装载 256 次 x 2048 B = 524 KB 总装载 = 65536 x 524 KB = 34 GB → 但 A、B 各自只有 64 MB, 靠 L2/并行复用后, 实际 HBM 流量接近 A + B + C ≈ 192 MB → HBM 时间 ≈ 192 MB / 1555 GB/s = 0.12 ms ≪ 实际耗时, 因此是计算受限 - 占用率：
FP32 tiled : 寄存器约 32/线程, 256 线程/block 寄存器限制 = 65536 / (32 x 256) = 8 block；线程限制 = 2048 / 256 = 8 block → 驻留 8 block = 64 warp = 100% 占用率；共享内存 2 x 16 x 16 x 4 B = 2 KB/block → 8 block 共 16 KB，远低于 A100 每 SM 的 164 KB（不是限制项） WMMA : 寄存器约 64/线程, 128 线程/block 寄存器限制 = 65536 / (64 x 128) = 8 block；线程限制 = 2048/128 = 16 block → 实际 8 block = 32 warp = 50% 占用率（寄存器受限） - 定量对比（A100，典型观测值）：
CUDA core FP32 tiled (TILE=16) : 约 9.6 TFLOPS（FP32 峰值 19.5 的 49%） 损失来源: 50% 的合并访问效率、每线程只有 16 次乘加就产生两次 __syncthreads、 共享内存 load 与 FMA 的指令比约 2:1（每 1 次 FMA 要 2 次 shared load） → 时间 = 137.4 GFLOP / 9.6 TFLOPS ≈ 14.3 ms Tensor Core FP16 WMMA : 约 45 TFLOPS（FP16 TC 峰值 312 的 14%） → 时间 = 137.4 GFLOP / 45 TFLOPS ≈ 3.05 ms 加速比 ≈ 14.3 / 3.05 ≈ 4.7x这里最值得记住的一点是：Times 4.7x 来自硬件，但离 16x 还差得远，差距全部来自”没有把 Tensor Core 喂饱”。生产级实现（cuBLAS / CUTLASS）会做四件朴素版本没做的事：① 把 A、B 的 tile 先经共享内存装载（global→shared→fragment），把全局流量降低一个数量级；② 用双缓冲（double buffering）把第 k+1 步的装载与第 k 步的
mma_sync重叠；③ 用更大的 warp tile（如 64×64）提高每个 fragment 的复用次数，把load_matrix_sync的次数摊薄；④ 用cp.async（Ampere 的异步拷贝）绕过寄存器直接填共享内存。做完这四件事后可以达到 150~250 TFLOPS（峰值的 50%~80%）。 - Roofline 判定与优化方向：把两个 kernel 都画在 Roofline 上——FP32 tiled 点在
(4 FLOP/B, 9.6 TFLOPS)，WMMA 点在(8 FLOP/B, 45 TFLOPS)；A100 的 fp32 屋脊点是(12.5, 19.5)、fp16 TC 屋脊点是(200, 312)。两个点都在各自屋脊的左侧（内存受限侧）——这不仅解释了为什么 WMMA 只有 45 TFLOPS，也指出了唯一正确的优化方向：提高算术强度（增加数据复用），而不是找”更快的乘法指令”。当复用提升到 AI > 200 FLOP/B 时，Tensor Core 才会成为真正的瓶颈。
代码示例 4：用 CUDA 实现全连接层前向传播（y = W·x + b，本质是 GEMM）
// 文件: fc_layer.cu
// 编译: nvcc -O3 -arch=sm_80 fc_layer.cu -o fc_layer
// 运行: ./fc_layer
//
// 深度神经网络的第一层: y[b][o] = ReLU( sum_i W[o][i] * x[b][i] + bias[o] )
// batch = 1024, in = 784 (28x28 灰度图), out = 1024
// 这是一个 NT 型 GEMM: Y = X * W^T + b, X 是 [batch x in], W 是 [out x in]
// 本文件同时给出 im2col 变换函数, 说明卷积层如何被展开成同一个 GEMM。
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cuda_runtime.h>
#define CUDA_CHECK(call) \
do { \
cudaError_t err__ = (call); \
if (err__ != cudaSuccess) { \
fprintf(stderr, "CUDA error %s:%d: %s\n", __FILE__, __LINE__, \
cudaGetErrorString(err__)); \
exit(EXIT_FAILURE); \
} \
} while (0)
#define TILE 16 /* tile 边长; 必须整除 batch / in_dim / out_dim */
#define BATCH 1024
#define IN_DIM 784 /* 28*28 = 784 = 49 * 16 */
#define OUT_DIM 1024 /* = 64 * 16 */
#define BLOCK_THREADS (TILE * TILE)
/* ==================== 全连接层前向 kernel ==================== */
/* X: [batch x in_dim] W: [out_dim x in_dim] bias: [out_dim] Y: [batch x out_dim] */
__global__ void fcForwardKernel(const float * __restrict__ X,
const float * __restrict__ W,
const float * __restrict__ bias,
float * __restrict__ Y,
int in_dim, int out_dim)
{
__shared__ float subX[TILE][TILE];
__shared__ float subW[TILE][TILE];
const int tx = threadIdx.x;
const int ty = threadIdx.y;
const int Row = blockIdx.y * TILE + ty; /* batch 方向 */
const int Col = blockIdx.x * TILE + tx; /* 输出神经元方向 */
float acc = 0.0f;
for (int k0 = 0; k0 < in_dim; k0 += TILE) {
/* 协作载入 X 的 tile: 行 = batch (ty), 列 = 输入维 (tx) -> 天然合并 */
subX[ty][tx] = X[Row * in_dim + k0 + tx];
/* 协作载入 W 的 tile: 行 = 输出神经元 (ty), 列 = 输入维 (tx) -> 不转置 */
subW[ty][tx] = W[(blockIdx.x * TILE + ty) * in_dim + k0 + tx];
__syncthreads();
#pragma unroll
for (int k = 0; k < TILE; ++k) {
acc += subX[ty][k] * subW[k][tx]; /* 注意: 用 subW[k][tx], 不是 subW[tx][k] */
}
__syncthreads();
}
acc += bias[Col];
Y[Row * out_dim + Col] = (acc > 0.0f) ? acc : 0.0f; /* ReLU 激活 */
}
/* ==================== im2col: 把卷积展开成 GEMM 的行 ==================== */
/* 把输入特征图 X[C][H][W] 展开成大小为 (C*K*K) x (H_out*W_out) 的矩阵 X_col,
之后卷积层的前向就变成一次普通 GEMM: Y_col = W_row * X_col
(W_row 的形状是 M x (C*K*K), 由 M 个卷积核按同样顺序拉直得到) */
static void im2col(const float *X, int C, int H, int Wd, int K,
float *X_col, int H_out, int W_out)
{
const int col_size = H_out * W_out;
for (int c = 0; c < C; ++c) {
for (int p = 0; p < K; ++p) {
for (int q = 0; q < K; ++q) {
const int row = (c * K + p) * K + q; /* 展开矩阵的行下标 */
for (int h = 0; h < H_out; ++h) {
for (int w = 0; w < W_out; ++w) {
const int col = h * W_out + w; /* 展开矩阵的列下标 */
X_col[row * col_size + col] =
X[(c * H + (h + p)) * Wd + (w + q)];
}
}
}
}
}
}
/* ==================== 主机端工具 ==================== */
static unsigned int g_seed = 987654321u;
static float rnd(void)
{
g_seed = g_seed * 1664525u + 1013904223u;
return ((float)(g_seed >> 8) / (float)(1u << 24)) - 0.5f; /* [-0.5, 0.5) */
}
/* ============================== main ============================== */
int main(void)
{
const size_t szX = (size_t)BATCH * IN_DIM * sizeof(float);
const size_t szW = (size_t)OUT_DIM * IN_DIM * sizeof(float);
const size_t szB = (size_t)OUT_DIM * sizeof(float);
const size_t szY = (size_t)BATCH * OUT_DIM * sizeof(float);
float *h_X = (float *)malloc(szX);
float *h_W = (float *)malloc(szW);
float *h_b = (float *)malloc(szB);
float *h_Y = (float *)malloc(szY);
if (h_X == NULL || h_W == NULL || h_b == NULL || h_Y == NULL) {
fprintf(stderr, "host malloc failed\n");
exit(EXIT_FAILURE);
}
for (size_t i = 0; i < (size_t)BATCH * IN_DIM; ++i) h_X[i] = rnd();
for (size_t i = 0; i < (size_t)OUT_DIM * IN_DIM; ++i) h_W[i] = rnd();
for (int i = 0; i < OUT_DIM; ++i) h_b[i] = rnd();
float *d_X = NULL, *d_W = NULL, *d_b = NULL, *d_Y = NULL;
CUDA_CHECK(cudaMalloc((void **)&d_X, szX));
CUDA_CHECK(cudaMalloc((void **)&d_W, szW));
CUDA_CHECK(cudaMalloc((void **)&d_b, szB));
CUDA_CHECK(cudaMalloc((void **)&d_Y, szY));
CUDA_CHECK(cudaMemcpy(d_X, h_X, szX, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_W, h_W, szW, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b, szB, cudaMemcpyHostToDevice));
const dim3 block(BLOCK_THREADS / TILE, TILE); /* (16, 16) = 256 线程 */
const dim3 grid(OUT_DIM / TILE, BATCH / TILE); /* (64, 64) = 4096 block */
/* 预热 */
fcForwardKernel<<<grid, block>>>(d_X, d_W, d_b, d_Y, IN_DIM, OUT_DIM);
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) {
fcForwardKernel<<<grid, block>>>(d_X, d_W, d_b, d_Y, IN_DIM, OUT_DIM);
}
CUDA_CHECK(cudaEventRecord(e1));
CUDA_CHECK(cudaEventSynchronize(e1));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, e0, e1));
ms /= (float)ITERS;
const double gflop = 2.0 * (double)BATCH * (double)IN_DIM * (double)OUT_DIM;
const double bytes = (double)(szX + szW + szY);
printf("===== 全连接层前向: Y[%d x %d] = ReLU(X[%d x %d] * W[%d x %d]^T + b) =====\n",
BATCH, OUT_DIM, BATCH, IN_DIM, OUT_DIM, IN_DIM);
printf(" 耗时 : %.4f ms\n", ms);
printf(" 运算量 : %.3f GFLOP\n", gflop / 1e9);
printf(" 性能 : %.2f GFLOP/s (%.3f TFLOPS)\n",
gflop / ((double)ms / 1000.0) / 1e9,
gflop / ((double)ms / 1000.0) / 1e12);
printf(" 访存量 : %.2f MB -> 算术强度 %.1f FLOP/B\n",
bytes / 1e6, gflop / bytes);
/* 抽样验证 */
CUDA_CHECK(cudaMemcpy(h_Y, d_Y, szY, cudaMemcpyDeviceToHost));
int bad = 0, nc = 32;
for (int t = 0; t < nc; ++t) {
const int b = (int)((double)t / nc * BATCH);
const int o = (int)((double)((t * 6151) % nc) / nc * OUT_DIM);
double acc = 0.0;
for (int i = 0; i < IN_DIM; ++i) {
acc += (double)h_X[(size_t)b * IN_DIM + i] * (double)h_W[(size_t)o * IN_DIM + i];
}
acc += (double)h_b[o];
const double ref = (acc > 0.0) ? acc : 0.0;
if (fabs((double)h_Y[(size_t)b * OUT_DIM + o] - ref) > 1.0e-2 + 1.0e-3 * fabs(ref)) {
++bad;
}
}
printf(" 抽样验证 : %s (%d 个 (batch, neuron) 对)\n",
(bad == 0) ? "PASS" : "FAIL", nc);
/* im2col 的自检: 用 3 通道 5x5 输入、2x2 卷积核展开, 打印前 3 行 3 列 */
{
const int C = 3, H = 5, Wd = 5, K = 2;
const int H_out = H - K + 1, W_out = Wd - K + 1;
float *Xsmall = (float *)malloc((size_t)C * H * Wd * sizeof(float));
float *Xcol = (float *)malloc((size_t)C * K * K * H_out * W_out * sizeof(float));
if (Xsmall == NULL || Xcol == NULL) {
fprintf(stderr, "im2col demo malloc failed\n");
exit(EXIT_FAILURE);
}
for (int i = 0; i < C * H * Wd; ++i) Xsmall[i] = (float)i;
im2col(Xsmall, C, H, Wd, K, Xcol, H_out, W_out);
printf("\n im2col 演示 (C=3, 5x5, K=2 -> X_col 是 12 x 16):\n");
for (int r = 0; r < 3; ++r) {
printf(" row %2d:", r);
for (int c = 0; c < 3; ++c) {
printf(" %6.1f", Xcol[r * (H_out * W_out) + c]);
}
printf("\n");
}
free(Xsmall);
free(Xcol);
}
CUDA_CHECK(cudaEventDestroy(e0));
CUDA_CHECK(cudaEventDestroy(e1));
CUDA_CHECK(cudaFree(d_X));
CUDA_CHECK(cudaFree(d_W));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_Y));
free(h_X);
free(h_W);
free(h_b);
free(h_Y);
return 0;
}
【代码做什么？】
- 形状约定：X 是
[batch × in_dim]行主序，W 是[out_dim × in_dim]行主序（这正是深度学习框架里权重矩阵的存法：每个输出神经元占一行），Y 是[batch × out_dim]。因此Y[b][o] = Σ_i X[b][i]·W[o][i] + bias[o]，即Y = X·Wᵀ + b——这就是 GEMM 的 NT 形式（A 不转置、B 转置），也是所有 cuBLAS/CUTLASS 里最常见的形态。 - 线程映射：block 是
(16,16)=256线程，grid 是(64,64)=4096block。blockIdx.y对应 batch 方向的 tile，blockIdx.x对应输出神经元方向的 tile；每个线程负责一个输出元素Y[Row][Col]。 - 协作载入：
subX[ty][tx]从X[Row*in_dim + k0+tx]载入（不转置），subW[ty][tx]从W[(blockIdx.x*TILE+ty)*in_dim + k0+tx]载入（也不转置，因为 W 天然就是”输出神经元 × 输入维”）。两个 tile 都用行主序存放，这是本 kernel 避免 bank conflict 的关键设计。 - 点积：
acc += subX[ty][k] * subW[k][tx]，k从 0 到 15；因为subX的第ty行对应本线程的 batch、subW的第tx行对应本线程的输出神经元（因为Col = blockIdx.x*TILE + tx，而subW[ty][tx] = W[(blockIdx.x*TILE+ty)*in_dim + k0+tx]，所以subW[k][tx] = W[Col*in_dim + k0+k] = W[Col][k0+k]）——正确的。 - 激活与偏置：
acc += bias[Col]后做 ReLU（acc > 0 ? acc : 0）。ReLU 就是讲义中”Today’s Choice”的max(0, x)斜坡函数，比 sigmoid 快得多（不需要指数运算），2017 年后成为主流。 - im2col 辅助函数：把
X[C][H][W]展开成(C·K·K) × (H_out·W_out)的矩阵，行下标(c*K+p)*K+q、列下标h*W_out+w。展开后卷积层前向 = 一次普通 GEMM，可以完全复用本 kernel 的 tiling 结构。 - 计时与验证：20 次求平均；抽样 32 个
(batch, neuron)对，在 CPU 上用 double 重算点积+偏置+ReLU 做对比。
【并行机制与硬件映射解说】
- 共享内存访问与 bank conflict（逐个访问分析）：
subX[ty][k]：地址 =ty*16 + k。k在 warp 内是循环变量（同一轮内所有线程相同），ty在一个 warp 内只取 2 个值（warp = 两行），所以 32 个线程只访问 2 个地址 → 硬件广播，0 路冲突。subW[k][tx]：地址 =k*16 + tx。tx = 0..15给出k*16+0 到 k*16+15，这 16 个字落在 bank(k*16 + tx) mod 32，tx遍历 0..15 时它们分布在两个连续的 16 字节区间上（模 32 后是 16 个不同 bank），而每个地址被ty=0与ty=1两组线程同时请求 → 广播。0 路冲突。- 反例（必须避免的写法）：如果为了”看起来对称”而把 W 的 tile 转置存放（
subW[tx][ty] = W[(blockIdx.x * TILE + ty) * in_dim + k0 + tx]）并用subW[tx][k]读取，地址变成tx*16 + k，tx=0..15给出k, 16+k, 32+k, 48+k, 64+k，这些字的mod 32只取k mod 32与(k+16) mod 32两个值 → 16 个不同的字挤在 2 个 bank 上 → 8 路冲突，共享内存只用上 1/8 的带宽。修法是 padding：__shared__ float subW[TILE][TILE + 1]，地址变成tx*17 + k，mod 32就散开了。结论：本 kernel 的写法（不转置 + 用subW[k][tx]）本身就是无冲突的，不要”优化”成转置版。
- 全局内存合并访问：
subX[ty][tx] = X[Row*in_dim + k0 + tx]，warp 内tx=0..15连续、ty有两个值（两行 batch），每行的 16 个连续 float = 64 字节，所以一个 warp 触发 2 次 64 字节 transaction。相对 128 字节的理想 transaction，效率 50%。同样的问题出现在subW的载入上。改进办法：把 block 改成(32,8)（tx覆盖 32 个连续元素 → 一次 128 字节 transaction），或把TILE提到 32 并让tid沿x展开 32 个线程。这是一个用 tile 宽度的选择换合并效率的经典权衡。 - warp 与调度：256 线程 = 8 warp/block。每个 warp 在点积循环里执行 16 次迭代、每次 2 条 shared load + 1 条 FMA，shared load : FMA = 2 : 1。A100 每个 SM 每周期最多发射 4 条指令（每 scheduler 1 条），而 FMA 的吞吐上限是 64/周期（即每 warp 2 条 FMA 指令/周期）。当 shared load 占了 2/3 的发射槽时，FMA 的实际发射率被压到约 1/3——这是 tiled kernel 只跑到峰值 50% 左右的核心原因。解决办法是让每个线程计算多个输出元素（register blocking / thread coarsening），把 shared load 的次数摊薄。
- 占用率：寄存器约 32/线程 →
65536/(32×256) = 8block；线程限制2048/256 = 8block。共享内存2 × 16 × 16 × 4 = 2 KB/block，8 个 block 共 16 KB，远低于 A100 的 164 KB（也不是限制项）。→ 100% 占用率（64 warp/SM）。 - 与 DNN 的对应关系：这个 kernel 就是”一个全连接层”。真实的 CNN 里，卷积层经 im2col 后是同一个 GEMM 结构（只是
in_dim变成C·K·K、batch变成B·H_out·W_out）；反向传播里的dE/dW = (dE/dfc)·xᵀ是另一次 GEMM（转置方向不同，是 TN 形式），但 tiling 结构完全一样。所以”把 GEMM 做到极致”等于把整个 DNN 做到极致——这就是 cuDNN/cuBLAS 存在的理由，也是 Tensor Core 存在的理由。
【性能优化分析】
- 算术强度与 Roofline：
运算量 = 2 × batch × in_dim × out_dim = 2 × 1024 × 784 × 1024 = 1.643e9 FLOP = 1.643 GFLOP 访存量 = X(1024×784×4 B = 3.211 MB) + W(1024×784×4 B = 3.211 MB) + Y(1024×1024×4 B = 4.194 MB) = 10.616 MB 算术强度 AI = 1.643e9 / 1.0616e7 ≈ 155 FLOP/B A100 Roofline 拐点 = 19.5e12 / 1555e9 ≈ 12.5 FLOP/B 155 ≫ 12.5 ⇒ 计算受限期望性能 =
min(19.5 TFLOPS, 1555 GB/s × 155 FLOP/B = 241 TFLOPS) = 19.5 TFLOPS；但受前述”shared load : FMA = 2:1”的发射限制，实测通常落在 7~10 TFLOPS：实测时间 ≈ 1.643e9 / 8.5e12 ≈ 0.19 ms ×20 次取平均后打印的 GFLOP/s 应落在 6~10 TFLOP/s 区间对比一下同一层的”逐元素”实现（每线程做一个输出神经元、从全局内存读整行 W）：
AI = 2·784/(784·4 + 784·4 + 4) ≈ 0.32 FLOP/B，掉到 12.5 以下，立刻变成内存受限，性能上限只有1555 × 0.32 = 497 GFLOP/s——同一层的两种写法性能差 17 倍，差别完全来自数据复用。 - 占用率与延迟隐藏：
每 SM 驻留 8 block = 2048 线程 = 64 warp，占用率 100% 每 warp 每周期发射指令数上限 = 1（单发射） 点积循环每 16 次迭代产生 16 条 FMA + 32 条 LDS，共 48 条指令 ⇒ FMA 占比 = 16/48 = 33%，即 FMA 发射率上限 = 33% × 峰值 实测 8.5 TFLOPS / 19.5 TFLOPS = 43.6% —— 略高于 33%， 因为编译器会把部分 shared load 合并/复用进寄存器这就是”指令混合比（instruction mix）决定实际性能“的典型案例：占用率已经是 100%，再调块大小没用，唯一出路是增加每线程的计算量（register blocking），把 48 条指令里的 16 条 FMA 变成 32 条甚至 64 条。
- 可执行的优化方向：① thread coarsening：每个线程算 2×2 或 4×4 个输出元素，
subX/subW的 shared load 次数按sqrt(复用数)摊薄，FMA 占比可以提到 60%~70%；② 向量化访存：用float4载入全局内存，把 2 次 64 字节 transaction 变成 1 次 128 字节（前提是 16 字节对齐）；③ 改回 128 字节合并：把 block 形状改为(32,8)或TILE=32；④ 上 Tensor Core：把输入量化到 fp16/bf16、用 WMMA 或mmaPTX 指令替换 FMA，这一层的理论峰值从 19.5 TFLOPS 变成 312 TFLOPS；⑤ 融合：把 bias 加法与 ReLU 融进 GEMM 的 epilogue（本 kernel 已经这么做了），避免多写一遍 Y；⑥ 训练时用 mini-batch（讲义页 19~21：逐个样本累积 ΔΘ 最准确但并行度低，mini-batch 在”梯度估计精度”与”并行度”之间取平衡）。
代码示例 5：WebGPU / WGSL 与 CUDA 的逐行对照（课程最终项目所用模型）
文件 1：WGSL 计算着色器（vector_add.wgsl），与 CUDA 的 vecAdd kernel 功能完全相同。
// 文件: vector_add.wgsl
// 运行环境: 浏览器(Chrome/Edge 113+)或原生 WebGPU 实现
// 计算 c[i] = a[i] + b[i]，n 由数组长度决定
// ---- 资源绑定: 相当于 CUDA kernel 的指针参数 ----
// 对应 CUDA: const float* __restrict__ a
@group(0) @binding(0) var<storage, read> a : array<f32>;
// 对应 CUDA: const float* __restrict__ b
@group(0) @binding(1) var<storage, read> b : array<f32>;
// 对应 CUDA: float* c (read_write 表示既可读也可写)
@group(0) @binding(2) var<storage, read_write> c : array<f32>;
// ---- 工作组大小: 相当于 CUDA 的 blockDim ----
// override 是"着色器编译期常量", 创建 pipeline 时可用 constants 覆盖(类似模板参数)
override WG_SIZE : u32 = 256u;
// 对应 CUDA: __global__ void vecAdd(const float* a, const float* b, float* c, int n)
// @compute 表示这是一个计算着色器阶段; @workgroup_size 指定工作组内 invocation 数
@compute @workgroup_size(WG_SIZE)
fn main(
// 对应 CUDA: int i = blockIdx.x * blockDim.x + threadIdx.x;
@builtin(global_invocation_id) gid : vec3<u32>,
// 对应 CUDA: threadIdx.x / threadIdx.y / threadIdx.z
@builtin(local_invocation_id) lid : vec3<u32>,
// 对应 CUDA: blockIdx.x / blockIdx.y / blockIdx.z
@builtin(workgroup_id) wid : vec3<u32>,
// 对应 CUDA: gridDim.x / gridDim.y / gridDim.z
@builtin(num_workgroups) nwg : vec3<u32>
) {
// 对应 CUDA: if (i < n) —— 这里 n 用 arrayLength 动态取得,
// 相当于 CUDA 里把 n 作为参数传进来(WebGPU 没有隐式的大小信息,
// 但运行时数组的 arrayLength 就是缓冲区元素个数)
let n : u32 = arrayLength(&a);
let i : u32 = gid.x;
// 对应 CUDA: if (i < n) c[i] = a[i] + b[i];
if (i >= n) {
return;
}
c[i] = a[i] + b[i];
// lid / wid / nwg 在本例中未使用, 保留是为了展示对应关系;
// WGSL 要求声明的内建变量可以不用, 不会报错。
let unused_lid : u32 = lid.x;
let unused_wid : u32 = wid.x;
let unused_nwg : u32 = nwg.x;
}
文件 2：主机端 JavaScript（浏览器或 Deno/Node + WebGPU）。
// 文件: vector_add_host.js
// 运行: 在支持 WebGPU 的浏览器控制台/模块中执行, 或在 Deno 中执行
// ---- 1. 获取适配器与设备: 对应 CUDA 的 cudaGetDevice / cudaSetDevice ----
const adapter = await navigator.gpu.requestAdapter();
if (!adapter) { throw new Error('no WebGPU adapter'); }
const device = await adapter.requestDevice();
// ---- 2. 准备主机数据 ----
const n = 1 << 24; // 16.7M 个 float
const elemSize = 4;
const bytes = n * elemSize;
const h_a = new Float32Array(n);
const h_b = new Float32Array(n);
const h_c = new Float32Array(n);
for (let i = 0; i < n; ++i) { h_a[i] = 1.0 + (i % 7); h_b[i] = 2.0 + (i % 5); }
// ---- 3. 创建缓冲区: 对应 CUDA 的 cudaMalloc ----
const usageIn = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST;
const bufA = device.createBuffer({ size: bytes, usage: usageIn });
const bufB = device.createBuffer({ size: bytes, usage: usageIn });
const bufC = device.createBuffer({
size: bytes,
usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST,
});
// ---- 4. 上传数据: 对应 CUDA 的 cudaMemcpy(HostToDevice) ----
device.queue.writeBuffer(bufA, 0, h_a);
device.queue.writeBuffer(bufB, 0, h_b);
// ---- 5. 编译着色器 + 创建计算管线 ----
// 对应 CUDA 的 nvcc 离线编译; WebGPU 是运行时编译
const wgslSource = await (await fetch('./vector_add.wgsl')).text();
const shaderModule = device.createShaderModule({ code: wgslSource });
const pipeline = device.createComputePipeline({
layout: 'auto', // 自动推导绑定布局
compute: { module: shaderModule, entryPoint: 'main' },
});
// ---- 6. 绑定组: 把缓冲区按 @group/@binding 槽位"插"到着色器上 ----
const bindGroup = device.createBindGroup({
layout: pipeline.getBindGroupLayout(0),
entries: [
{ binding: 0, resource: { buffer: bufA } },
{ binding: 1, resource: { buffer: bufB } },
{ binding: 2, resource: { buffer: bufC } },
],
});
// ---- 7. 编码并提交命令: 对应 CUDA 的 kernel<<<grid, block>>> ----
const WORKGROUP_SIZE = 256;
const numWorkgroups = Math.ceil(n / WORKGROUP_SIZE);
const encoder = device.createCommandEncoder();
const pass = encoder.beginComputePass();
pass.setPipeline(pipeline);
pass.setBindGroup(0, bindGroup);
pass.dispatchWorkgroups(numWorkgroups, 1, 1); // 对应 cudaLaunchKernel 的 grid
pass.end();
const t0 = performance.now();
device.queue.submit([encoder.finish()]); // 对应 CUDA 的 kernel 启动(异步)
// ---- 8. 同步: 对应 CUDA 的 cudaDeviceSynchronize() ----
await device.queue.onSubmittedWorkDone();
const t1 = performance.now();
console.log(`kernel 时间: ${(t1 - t0).toFixed(3)} ms`);
// ---- 9. 读回结果: 对应 CUDA 的 cudaMemcpy(DeviceToHost) ----
// WebGPU 用 mapAsync + getMappedRange 完成"设备内存 -> 可被 CPU 访问的视图"
await bufC.mapAsync(GPUMapMode.READ);
const mapped = new Float32Array(bufC.getMappedRange().slice(0)); // slice 复制一份
bufC.unmap();
// ---- 10. 验证 ----
let bad = 0;
for (let i = 0; i < n; i += 1000003) {
if (Math.abs(mapped[i] - (h_a[i] + h_b[i])) > 1e-6) { ++bad; }
}
console.log(bad === 0 ? 'PASS' : 'FAIL');
// ---- 11. 释放: 对应 CUDA 的 cudaFree; WebGPU 用 destroy() ----
bufA.destroy();
bufB.destroy();
bufC.destroy();
文件 3：WGSL 的工作组内归约，用来对照 CUDA 的 __shared__ + __syncthreads()。
// 文件: workgroup_reduce.wgsl
// 每个 workgroup 把 WG_SIZE 个元素求和, 写出一个部分和
// 对照 CUDA: __shared__ float s[256]; s[tid] = v[id]; __syncthreads();
// for (stride = blockDim.x / 2; stride > 0; stride >>= 1) {
// if (tid < stride) { s[tid] += s[tid + stride]; }
// __syncthreads();
// }
@group(0) @binding(0) var<storage, read> input : array<f32>;
@group(0) @binding(1) var<storage, read_write> output : array<f32>;
const WG_SIZE : u32 = 256u;
// workgroup 地址空间: 对应 CUDA 的 __shared__
var<workgroup> scratch : array<f32, 256>;
@compute @workgroup_size(WG_SIZE)
fn main(
@builtin(local_invocation_id) lid : vec3<u32>,
@builtin(global_invocation_id) gid : vec3<u32>,
@builtin(workgroup_id) wid : vec3<u32>
) {
let t : u32 = lid.x;
// 协作载入: 每个 invocation 搬一个元素到工作组存储
scratch[t] = input[gid.x];
// 对应 CUDA 的 __syncthreads(): 保证所有 invocation 都写完
workgroupBarrier();
// 树形归约, 步长逐次减半
var stride : u32 = WG_SIZE / 2u;
loop {
if (stride == 0u) { break; }
if (t < stride) {
scratch[t] = scratch[t] + scratch[t + stride];
}
workgroupBarrier(); // 每一轮都要同步, 否则会读到别人写了一半的数据
stride = stride / 2u;
}
// 只有 0 号 invocation 写出结果
if (t == 0u) {
output[wid.x] = scratch[0];
}
}
CUDA 与 WebGPU 的逐行对照表
| 行 | CUDA | WebGPU / WGSL | 语义差异说明 |
|---|---|---|---|
| 1 | __global__ void vecAdd(const float* a, const float* b, float* c, int n) | @compute @workgroup_size(256) fn main(@builtin(global_invocation_id) gid : vec3<u32>) | CUDA 的指针参数在 WGSL 里变成模块级绑定（@group/@binding），不是函数参数 |
| 2 | const float* __restrict__ a | @group(0) @binding(0) var<storage, read> a : array<f32>; | WGSL 显式区分只读（read）与读写（read_write），只读可让驱动做更好的布局 |
| 3 | float* c | var<storage, read_write> c : array<f32>; | 同上 |
| 4 | int n（显式传参） | arrayLength(&a) | WebGPU 的运行时数组自带长度；也可用 uniform buffer 传标量 |
| 5 | int i = blockIdx.x * blockDim.x + threadIdx.x; | @builtin(global_invocation_id) gid : vec3<u32> → gid.x | 语义完全相同，WGSL 由运行时直接给出全局 id |
| 6 | threadIdx.x | @builtin(local_invocation_id) lid : vec3<u32> | 相同 |
| 7 | blockIdx.x | @builtin(workgroup_id) wid : vec3<u32> | block ↔ workgroup |
| 8 | blockDim.x | @builtin(workgroup_size) : vec3<u32> 或 @workgroup_size(256) 的字面量 | 相同 |
| 9 | gridDim.x | @builtin(num_workgroups) : vec3<u32> | 相同 |
| 10 | __shared__ float s[256]; | var<workgroup> s : array<f32, 256>; | 都在片上，生命周期 = 一个 block/workgroup |
| 11 | __syncthreads(); | workgroupBarrier(); | 语义相同：workgroup 内所有 invocation 的屏障 |
| 12 | __syncwarp(); | subgroupBarrier(); | warp ↔ subgroup（大小平台相关，NVIDIA 上为 32） |
| 13 | atomicAdd(&x, 1.0f); | atomicAdd(&x, 1.0f); | 同名同参，是最高兴的一处巧合 |
| 14 | cudaMalloc(&d_a, bytes); | device.createBuffer({size: bytes, usage: GPUBufferUsage.STORAGE}) | WebGPU 必须显式声明用途位（STORAGE/COPY_SRC/COPY_DST） |
| 15 | cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice); | device.queue.writeBuffer(bufA, 0, h_a); | 相同语义 |
| 16 | vecAdd<<<grid, 256>>>(d_a, d_b, d_c, n); | pass.dispatchWorkgroups(numWorkgroups, 1, 1); | <<<>>> ↔ 命令编码器 + dispatchWorkgroups |
| 17 | cudaDeviceSynchronize(); | await device.queue.onSubmittedWorkDone(); | 相同语义 |
| 18 | cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost); | await bufC.mapAsync(GPUMapMode.READ); new Float32Array(bufC.getMappedRange()) | WebGPU 的读回是”映射”而不是”拷贝”，且必须 await |
| 19 | cudaFree(d_a); | bufA.destroy(); | 相同 |
| 20 | 没有对应物 | @group(0) @binding(0) 的绑定布局校验 | WebGPU 在创建 bind group 时就会检查布局是否与着色器匹配，能在提交前发现错绑 |
| 21 | 没有对应物 | override WG_SIZE : u32 = 256u; | 相当于”编译期模板参数”，可在创建 pipeline 时覆盖，用于快速调参 |
| 22 | -arch=sm_80 离线编译 | device.createShaderModule() 运行时编译 | WebGPU 的编译发生在运行时，会有一次性的编译延迟 |
| 23 | 1024 线程/block | 256 invocation/workgroup（保证下限） | 迁移时最需要注意的定量差异：归约层次更深 |
| 24 | Tensor Core wmma::mma_sync | 无对应物 | WGSL 没有子组矩阵乘指令，矩阵乘只能靠 dot() 手工展开 |
课程最终项目使用的 WebGPU / RAI 与 CUDA 的对照总结
课程目录里明确写着实验室用 WebGPU、最终项目用 RAI（原文：”C Programming Language and CUDA Software Development Kit, WebGPU for labs, RAI for final project”），而实验室历史上有过”Labs are done through WebGPU.net”与”Labs are done through Delta”两种环境。对已经会 CUDA 的人来说，迁移到最终项目只需要换掉四层认知：
- 计算抽象层：CUDA 让你直接写
<<<grid, block>>>和__global__；WebGPU 让你写 WGSL 的@compute入口 + 主机端编码命令。RAI 在 WebGPU 之上再抽象一层，把并行工作组织成任务（task）/ 图（graph）：你描述”有哪些数据缓冲区、有哪些 kernel、它们之间的依赖是什么”，由框架负责提交顺序、缓冲区分配与同步。这相当于把 CUDA 里手工写的cudaMemcpyAsync+ stream + event 依赖图，交给框架自动生成——概念上等价于 CUDA Graphs，但跑在 WebGPU 上。 - 可移植性层：CUDA 代码只能在 NVIDIA GPU 上跑，且需要安装驱动与 CUDA Toolkit；WebGPU/RAI 代码在浏览器、Windows、macOS、Linux、Android、iOS 上都能跑，零安装分发——这正是课程把最终项目放在 WebGPU 上的原因：学生可以在自己的笔记本上开发、在浏览器里演示。
- 性能层：CUDA 上可以用 Tensor Core（A100 上 fp16 峰值 312 TFLOPS）、
cp.async、__shfl_sync、cluster 等全部硬件特性；WebGPU 上目前只有f32/f16标量与向量运算、subgroup操作（且部分是可选项）、以及工作组的显式存储。因此同一个 GEMM 的性能上限相差一个数量级以上——但这不影响把这些优化技术（tiling、register blocking、合并访存、减少 bank conflict、提高占用率）原样迁移过去，因为它们的瓶颈模型完全一样。 - 正确性层：CUDA 的错误多半是”静默的错误数据”（忘记同步、越界、方向写错），WebGPU 多了一层编译期/创建期校验（绑定布局、访问模式、越界检测在开发模式下会报错），但读回数据时少了一层”自动同步”（必须
await mapAsync），所以最容易犯的错变成”忘记 await”——表现是读到全 0 或旧数据，而不是崩溃。
一句话总结课程范围的迁移路线：把 CUDA 的”线程映射 + 内存层次 + 同步”三件套原样翻译成 WGSL 的”invocation 映射 + storage/workgroup 地址空间 + workgroupBarrier”，把 CUDA 的”stream + event”翻译成”命令编码器 + submit + 框架的依赖图”，剩下的优化直觉完全通用。
【代码做什么？】
- 着色器侧：
main通过@builtin(global_invocation_id)拿到全局 invocation 编号i，用arrayLength(&a)取得真实长度n，做边界检查后执行c[i] = a[i] + b[i]。a、b是只读 storage buffer（var<storage, read>），c是可读写 storage buffer（var<storage, read_write>）。 - 主机侧获取设备：
navigator.gpu.requestAdapter()对应cudaGetDeviceCount的”枚举设备”，adapter.requestDevice()对应cudaSetDevice+ 上下文创建。 - 创建缓冲区：
device.createBuffer({size, usage})对应cudaMalloc，但必须显式声明用途位：读入用STORAGE \| COPY_DST，读回结果再加COPY_SRC。 - 上传数据：
device.queue.writeBuffer(bufA, 0, h_a)对应cudaMemcpy(HostToDevice)。它内部就是一次队列提交，不需要显式同步（后续提交的命令天然排在它后面）。 - 创建管线与绑定组：
createShaderModule运行编译 WGSL（对应nvcc的离线编译），createComputePipeline对应”选定 kernel + 执行配置模板”，createBindGroup把缓冲区按@group(0) @binding(k)的槽位插上去（这一步会校验布局，错绑会在提交前报错）。 - 编码与提交：
beginComputePass→setPipeline→setBindGroup→dispatchWorkgroups(ceil(n/256))→pass.end()→queue.submit([encoder.finish()])。整体对应 CUDA 的一次kernel<<<grid, block>>>加启动，但 WebGPU 把”命令录制”与”提交执行”分开，因此可以一次录制、多次提交。 - 同步与读回：
await device.queue.onSubmittedWorkDone()对应cudaDeviceSynchronize()；读回必须await buffer.mapAsync(GPUMapMode.READ)再getMappedRange()，最后unmap()。 - 归约示例：
workgroup_reduce.wgsl把每个 invocation 的一个元素写进var<workgroup> scratch，workgroupBarrier()之后做步长减半的树形归约，每轮归约之间都要再workgroupBarrier()，最后编号 0 的 invocation 写出工作组部分和。 - 释放：
bufA.destroy()对应cudaFree。
【并行机制与硬件映射解说】
- invocation → warp 的映射：WGSL 的 workgroup 在 NVIDIA 硬件上就是 CUDA 的 thread block。
@workgroup_size(256)意味着每个 workgroup 有 256 个 invocation，硬件上被切成 8 个 warp（每 32 个连续local_invocation_id.x为一个 warp）。这是 CUDA 与 WebGPU 在硬件层完全一致的部分。 - 全局编号的构成：WebGPU 没有直接给出
blockIdx * blockDim + threadIdx的表达式，而是把结果包成global_invocation_id，等价于workgroup_id * workgroup_size + local_invocation_id。三者的关系与 CUDA 完全相同，只是 WGSL 用vec3<u32>一次给出三个维度。 - 合并访存：
gid.x连续的 32 个 invocation 访问a[i]中连续 32 个f32= 128 字节 = 一条 cache line = 一次 transaction，与 CUDA 版的合并行为完全一致。因此在 CUDA 上”让相邻线程访问相邻元素”的规则在 WGSL 里同样成立，把a[2*i]改成步长访问一样会让 transaction 数翻倍。 - 工作组存储与 bank 冲突：
var<workgroup> scratch : array<f32, 256>在硬件上落在片上暂存器（shared memory），同样按 32 bank × 4 字节组织。归约中scratch[t]与scratch[t + stride]的地址都是连续的（t在 warp 内连续），因此初始阶段是一一对应的 bank 访问（0 路冲突）；当stride < 32时，参与计算的 invocation 不足一个 warp，多出来的 lane 被谓词关闭（对应 CUDA 的 warp 发散），此时访问仍然是地址各异的单播，没有 bank 冲突。如果改成scratch[t * 32]这类跨步访问，就立刻会出现 32 路冲突——这与 CUDA 的规则一字不差。 - 工作组大小上限 256 的后果：CUDA 允许 1024 线程/block，WebGPU 只保证 256 invocation/workgroup。于是同样处理 4096 个元素，CUDA 需要 4 个 block，WebGPU 需要 16 个 workgroup；块内归约的层次也从 1024→512→…→32 变成 256→128→64→32，最后 32 元素的归约要么交给
subgroupAdd之类的子组操作，要么由单个 invocation 串行完成。这是迁移时最容易被忽略的性能细节。 - 同步与可见性：
workgroupBarrier()保证 workgroup 内的存储写对同组其他 invocation 可见（对应__syncthreads()）；跨 workgroup 的依赖在 WebGPU 里无法用屏障表达，只能通过拆成两次 dispatch 并把中间结果写进 storage buffer 来实现（对应 CUDA 里”拆成两个 kernel”或用全局原子/栅栏）。这与 CUDA 的语义边界完全一致：block 之间没有屏障，只有 kernel 边界。 - 命令提交与异步性：
queue.submit是异步的（对应 kernel 启动），await才是同步点（对应cudaDeviceSynchronize）。若忘记await就去读映射内存，读到的可能是全 0 或旧值——这与 CUDA 忘记cudaDeviceSynchronize()读到垃圾数据是同一类错误，但在 WebGPU 上不会崩溃、只会静默出错，更难调试。 - 没有 Tensor Core 通道：WGSL 没有
wmma/mma对应的标准指令，因此代码示例 3 里那种 16×16×16 的矩阵乘加阵列在 WebGPU 上无法直接使用；矩阵乘只能展开成vec4<f32>的点积（dot()）来提高每 invocation 的算术强度。这是 WebGPU 与 CUDA 之间最大的一处能力落差。
【性能优化分析】
- 算术强度与 Roofline（向量加法的 AI 与 CUDA 版完全相同）：
每个元素: 1 次加法 = 1 FLOP；访存 = 读 a 4 B + 读 b 4 B + 写 c 4 B = 12 B AI = 1 / 12 ≈ 0.083 FLOP/B 在 A100 上（假设 WebGPU 后端能跑满硬件）： 拐点 = 19.5e12 / 1555e9 = 12.5 FLOP/B ≫ 0.083 ⇒ 纯带宽受限 可达上限 = 1555 GB/s × 0.083 ≈ 129 GFLOP/s 16.7M 个元素: 64 MB × 3 = 192 MB -> 192 MB / 1555 GB/s ≈ 0.124 ms这解释了为什么 WebGPU 的向量加法在 CUDA 上和在浏览器上都会撞到同一条内存墙——差别不在算法，而在命令提交与校验带来的额外固定开销。
- 占用率的等价概念：WebGPU 里没有
cudaOccupancyAPI，但物理约束完全相同。一个 256 invocation、使用array<f32, 256>（1 KB）工作组的着色器，在 A100 上每 SM 可驻留的 workgroup 数受 triple 约束：workgroup 数 ≤ 32、invocation 数 ≤ 2048（即 8 个 256 的 workgroup）、workgroup 存储 ≤ 164 KB（1 KB/组不是限制）。因此实际驻留 8 个 workgroup = 2048 invocation = 64 warp = 100% 占用率，与 CUDA 版一致。 - 总时间模型（含 WebGPU 特有开销）：
T_webgpu ≈ T_compile(一次性) + T_encode + T_submit + T_kernel + T_mapAsync + T_validate 其中 T_compile 只在 createShaderModule / createComputePipeline 时发生一次, 若每帧重建 pipeline, 这项开销会从 <1 ms 变成 5~20 ms —— 必须缓存 pipeline 与 bind group。 T_kernel 本身 ≈ 0.124 ms（上文的带宽下界） ⇒ 在真实浏览器里, 小规模向量加法的时间几乎全被固定开销占据, 这正对应讲义"Execution Time Never Reaches Zero"与"Data transfers have similar non-linearities for small sizes"两条结论。 - 可执行的优化方向：① 缓存
GPUComputePipeline与GPUBindGroup，只在缓冲区尺寸变化时重建；② 让每个 invocation 处理多个元素（例如用vec4<f32>一次算 4 个），把 dispatch 的 invocation 总数降到 1/4，减少固定开销相对占比；③ 用 storage buffer 而不是 uniform buffer 传大数组（uniform 有 64 KB 级的大小限制与更严格的布局要求）；④ 减少 CPU↔GPU 往返：把多步计算录进同一个 command encoder 一次提交，而不是每步都await；⑤ 结果读回用mapAsync且只在真正需要时读回（例如只在验证阶段），推理路径上把数据留在 GPU 上；⑥ 若确实需要矩阵乘且必须用 WebGPU，把f32换成f16（需启用shader-f16扩展特性）以获得 2 倍的寄存器与带宽效率，并用dot(vec4<f16>, vec4<f16>)手工展开。
性能优化技巧总结
- 能用 pinned 就用 pinned：
cudaHostAlloc让 DMA 只有一跳，H2D 带宽从 ~11 GB/s 提到 ~23 GB/s（Gen4 x16），是”2 倍白送”的加速；但不要滥用——pinned 内存不可换页，过订会让整机内存吃紧。 - 异步拷贝必须配 pinned：
cudaMemcpyAsync遇到 pageable 内存会退化成同步 staging，重叠收益直接归零。 - 切段 + 多流做流水线：把总时间从
T_in + T_c + T_out压向max(T_in, T_c, T_out)；段数取”能填满流水线的最小值”，用实测扫描决定（讲义的例子最优是 4~5 段）。 - 至少三个流才能连续流水：两个流时
C.i会在拷贝引擎队列里挡住下一轮的A.i+1（表头阻塞），且 PCIe 只有一个方向在工作。 - 入队顺序要”分批”：先发完所有 H2D，再发所有 kernel，最后发所有 D2H，避免读方向被写方向阻塞。
- 用事件做跨流同步，不要用
cudaDeviceSynchronize：cudaEventRecord+cudaStreamWaitEvent只建立必要的依赖；cudaDeviceSynchronize会摧毁所有并行度。 - 警惕 legacy default stream：默认流与所有阻塞流隐式同步；用
cudaStreamNonBlocking或--default-stream per-thread规避。 - Tensor Core 要用”低精度输入 + FP32 累加”：fp16/bf16 输入把峰值拉到 FP32 的 16 倍，累加留在 fp32 才能控制误差不随 K 增长。
- 先提高算术强度，再谈换指令：朴素 WMMA 只有 8 FLOP/B、受 L2 带宽限制只能跑到 45 TFLOPS；把 tile 经共享内存做多级复用（AI > 200 FLOP/B）才能逼近 312 TFLOPS 的峰值。
- register blocking 比调块大小更能提升 tiled kernel：shared load : FMA = 2:1 时 FMA 发射率被压到 33%；每线程多算几个输出元素能把 shared load 摊薄，这是”占用率已经 100% 却还是慢”的标准解法。
- 一切 DNN 层都先翻译成 GEMM：全连接是 GEMM，卷积经 im2col 也是 GEMM，反向传播还是 GEMM；优化一套 GEMM 内核等于优化整个网络。
- 共享内存 padding 消除 bank conflict：步长为 32 的倍数的访问会造成 32 路冲突；padding 一列（
[TILE][TILE+1]）即可打散 bank。 - 合并访问按 128 字节 transaction 核算：
TILE_WIDTH=16的 tile 载入在一个 warp 内只覆盖两段 64 字节，效率只有 50%；让 warp 内tx覆盖 32 个连续元素即可拿满。 - 对齐 WebGPU 的约束：workgroup 上限 256、读回必须
await mapAsync、没有 Tensor Core 对应物——迁移前先把这三条记在纸上。 - 用 Roofline 决定”还能不能更快”：先算
max(各阶段)或min(峰值算力, 带宽×AI)的理论上限，再判断当前离上限多远，避免在已经到顶的地方瞎调。
关键要点
- 带宽决定一切：A100 的显存带宽 1555 GB/s 对 PCIe Gen4 x16 的 32 GB/s 是 49 倍差距，对 Gen5 x16 的 64 GB/s 是 24 倍；NVLink 3.0 的 600 GB/s 把 GPU 间带宽拉回 PCIe Gen4 的约 10 倍。任何”把数据搬来搬去”的算法都要先算这笔账。
- 拷贝比计算更值得优化：典型 GPU 程序里 kernel 只占几毫秒，PCIe 传输占几十毫秒；pinned + 多流流水线能把总时间从
T_in+T_c+T_out压到max(T_in, T_c, T_out)，这是”零算法改动”的加速。 - 流的本质是”没有依赖就可以并行”：同流内严格有序、跨流无约束；重叠的前提是硬件上 SM（计算）与 Copy Engine（DMA）本来就是独立的两个引擎。
- Tensor Core 的 16 倍不是自动到手的：硬件峰值 312 TFLOPS（A100，fp16）对比 FP32 的 19.5 TFLOPS，但朴素 WMMA kernel 只能拿到 ~45 TFLOPS（14%），差距全部来自数据复用不足导致的 L2 带宽瓶颈。
- 深度学习 = 大规模 GEMM：全连接层
Y = X·Wᵀ + b、卷积层经 im2col 展开、反向传播的dE/dW = (dE/dfc)·xᵀ都是 GEMM；这也是 GPU 硬件（Tensor Core、结构化稀疏、HBM）被深度学习需求反复推动的原因。 - CUDA 不是唯一的模型，但它的心智模型是通用的：无论 OpenCL 的 work-item、OpenACC 的 gang/worker、还是 WebGPU 的 invocation/workgroup，共同特征都是”层次化轻量核心 + 本地暂存器 + 软件管理内存 + 无硬件一致性 + BSP 同步”——而 WebGPU/WGSL（课程最终项目所用模型，其上还有 RAI 这层任务图框架）带来的可移植性与零安装分发，与 CUDA 提供的性能上限形成互补，在 CUDA 上学到的 tiling、合并访存、bank conflict、占用率这些直觉 100% 可以迁移过去。
常见陷阱与注意事项
- 忘记同步就使用结果：CUDA 侧表现为读到垃圾数据（kernel 还在跑就去 D2H 拷贝）→ 必须
cudaDeviceSynchronize()或对相应流cudaStreamSynchronize()；WebGPU 侧表现为读到全 0 或旧内容 → 必须await device.queue.onSubmittedWorkDone()或await buffer.mapAsync(GPUMapMode.READ)。这类错误不报错、只出错数据，是最难查的一类 bug。 - 忘记
__syncthreads()，或把它放进线程发散的分支：tiled kernel 里第二轮迭代的快速线程会覆盖别人还在读的 tile，结果随调度随机变化 → 在”协作载入完成之后”和”共享内存使用完毕之后”各放一次__syncthreads()，且绝不能写在if分支内部（CUDA 要求同一 block 内所有线程到达同一个屏障，否则行为未定义；WebGPU 的workgroupBarrier()同理，但会被校验器直接拒绝）。 - 共享内存 bank conflict：地址步长是 32 的倍数时 32 个线程挤在同一个 bank。典型反例是把 tile 转置存放后用
subW[tx][k]读取（地址tx*16+k），16 个不同的字只落在 2 个 bank 上，造成 8 路冲突、共享内存吞吐掉到 1/8 → 用[TILE][TILE+1]做 padding，或像代码示例 4 那样保持行主序并用subW[k][tx]访问。 - warp 发散（divergence）：
if (i < n)之类的边界检查开销可忽略，但把”整段计算”放进按数据分支的if/else里，会让一个 warp 的两条路径串行执行，吞吐直接减半 → 让同一 warp 内的线程走相同路径（如把条件改成基于threadIdx的对齐判断），或用谓词（predication）代替分支。 - 未检查 CUDA API 返回值：kernel 启动失败（如
too many resources requested for launch、显存不足）不会在启动点报错，只在事后cudaGetLastError()里出现，于是程序可能”跑完了但结果全错” → 用统一的CUDA_CHECK宏包住每一次 API 调用，并在每次 kernel 启动后补一次CUDA_CHECK(cudaGetLastError())。 - 越界访问与网格覆盖不足：
grid = n / block的整数除法会丢掉余数，最后一批元素永远算不到；反过来，grid取整过大又没有if (i < n)保护就会读写越界（可能静默破坏相邻缓冲区）→ 一律用(n + block - 1) / block计算网格，并在 kernel 首行加边界检查。 cudaMemcpy方向写错、host/device 指针混用、共享内存容量超限：方向参数写反（HostToDevice与DeviceToHost互换）在异步版本里可能不报错但数据错乱；把cudaMalloc的指针交给主机代码解引用（或反之）通常直接段错误；共享内存开得太大（例如两个[64][64]的 float tile 就是 32 KB）会压低驻留 block 数甚至启动失败 → 用h_/d_前缀区分指针、把 memcpy 方向封装成模板函数自动推导、用cudaOccupancyMaxPotentialBlockSize或prop.sharedMemPerBlockOptin核定共享内存预算。- 异步拷贝的”假重叠”与默认流隐式同步：把 pageable 内存传给
cudaMemcpyAsync会退化成同步 staging，收益归零；用kernel<<<g,b>>>（省略流参数）投到 legacy default stream 会与所有阻塞流隐式同步，把多流并行彻底摧毁；段切得太小则 SM 闲置、kernel 启动开销占比飙升，总时间反而变长 → 主机内存一律cudaHostAlloc、流一律cudaStreamNonBlocking（或编译加--default-stream per-thread）、段数用实测扫描而不是拍脑袋决定。
此外还有两个本讲特有的坑：WMMA 的 ldm 对齐与 Tensor Core 精度——16 位类型的 ldm 必须是 8 的倍数（16 字节对齐），否则 load_matrix_sync 会读错；把累加器也声明成 __half 会让误差随 K 线性累积，必须”低精度输入 + fp32 累加”，验证时用相对容差而不是精确相等。WebGPU 的 usage 位与 pipeline 复用——读回结果的 buffer 若没有 COPY_SRC、写出的 buffer 若没有 COPY_DST，会分别在 mapAsync/writeBuffer 时报错；而每帧重建 GPUComputePipeline 与 GPUBindGroup 会把固定开销从 <1 ms 抬到 5~20 ms，必须缓存。
思考题（带答案）
Q1. 某应用需要把 512 MB 的输入从主机搬到 A100、做完一次计算、再把 512 MB 的结果搬回来。计算部分在 A100 上需要 20 ms。使用 pinned 内存时 PCIe Gen4 x16 的实测单向带宽是 24 GB/s。问：(a) 串行版总时间是多少？(b) 用 4 段 4 流流水线的总时间是多少？(c) 段数趋于无穷时的理论下界是多少？(d) 如果改成 PCIe Gen5 x16（单向实测 48 GB/s），下界会变成多少？ (a) T_in = 512 MB / 24 GB/s = 21.3 ms，T_out = 21.3 ms，T_serial = 21.3 + 20 + 21.3 = 62.6 ms。(b) 用流水线公式 (T_in+T_c+T_out)/N + (N-1)·max(T_in, T_c, T_out)/N，max = max(21.3, 20, 21.3) = 21.3 ms，所以 N=4：62.6/4 + 3×21.3/4 = 15.65 + 15.98 = 31.6 ms，加速比 62.6/31.6 ≈ 1.98x。(c) N→∞ 时总时间 → max(21.3, 20, 21.3) = 21.3 ms，加速比 2.94x。(d) Gen5 把单向带宽翻倍：T_in = T_out = 10.7 ms，此时瓶颈变成计算（20 ms > 10.7 ms），下界 = max(10.7, 20, 10.7) = 20 ms。结论：当计算时间超过单方向传输时间时，继续提升 PCIe 带宽不再带来任何收益——优化方向必须转向 kernel 本身。
Q2. A100 上 FP16 Tensor Core 峰值是 312 TFLOPS，FP32 CUDA core 是 19.5 TFLOPS，比值 16 倍。但一个”global → fragment”的朴素 WMMA kernel 实测只有约 45 TFLOPS，而 FP32 tiled kernel 约 9.6 TFLOPS。请解释为什么 WMMA 只拿到 4.7 倍而不是 16 倍，并给出至少两条把差距缩小的具体手段。 朴素 WMMA kernel 每 warp 每个 k 步读入 2 × 16 × 16 × 2 B = 1024 B 的 fragment、做 2×16×16×16 = 8192 FLOP，算术强度只有 8192/1024 = 8 FLOP/B。全部从 HBM 取的上限是 1555 GB/s × 8 = 12.4 TFLOPS；即使靠 40 MB L2 命中把实际流量降下来，K/16 = 256 个 k 步累计的 L2 请求仍有约 16 GB，按 L2 带宽 ~6 TB/s 算是 2.7 ms，与实测同量级——瓶颈是 L2 带宽与 fragment 装载指令吞吐，不是 Tensor Core 的乘加阵列。对比之下 312 TFLOPS 的屋脊点是 312e12/1555e9 ≈ 200 FLOP/B，所以必须把 AI 提高一个数量级。手段：① 把 A、B 的 tile 先经共享内存装载，global→shared→fragment 三级复用，把全局/L2 流量降低 10 倍以上；② 双缓冲（double buffering） 或用 Ampere 的 cp.async，让第 k+1 步的装载与第 k 步的 mma_sync 重叠；③ 扩大 warp tile（如 64×64）以提高每个 fragment 的复用次数，摊薄 load_matrix_sync 的调用次数；④ 用 register blocking 让每个 warp 同时持有多个累加器 fragment，提高 mma_sync 的发射密度。
Q3. 讲义给出的流水线例题是：T_in = 10 s、kernel 启动常数开销 1 s、T_comp = 10 s、T_out = 10 s，3 条硬件流支持连续流水，求最优段数。请写出时间关于段数 N 的表达式、求出最优 N、并解释为什么 N 太大反而更慢。 时间表达式把”输入传输”和”输出传输”合并成 20/N（讲义把它们视为可在一段时间内连续推进的传输阶段）、计算固定为 10、再叠加 N 次 kernel 启动的常数开销（每次 1 s，且 kernel 启动不能与 kernel 执行重叠）：T(N) = 20/N + 10 + N。求极值：dT/dN = -20/N² + 1 = 0 → N = √20 ≈ 4.47，取整后 N = 4 或 N = 5，两者都是 T = 5 + 10 + 4 = 19 s（N=5 时为 4 + 10 + 5 = 19 s），加速比 31/19 ≈ 1.63x（串行版 10+1+10+10 = 31 s）。N 太大反而更慢的原因：20/N 这一项随 N 增大迅速饱和（N=10 时只剩 2 s，比 N=4 的 5 s 只省 3 s），而 N 这一项（每次 kernel 启动的固定开销，以及段太小导致 SM 闲置、warp 数不足以填满流水线、负载不均）随 N 线性增长，两者此消彼长，最优出现在 20/N ≈ N 即 N ≈ 4.5 附近。这正对应讲义那两张图：”Execution Time is Ideally Linear in Size” 与 “Execution Time Never Reaches Zero”。