Lecture 5: GPU Architecture and CUDA Programming

目录 · ← l4 · l6 →

Lecture 5: GPU Architecture and CUDA Programming

1. 章节标题与概述

Lecture 5: GPU Architecture and CUDA Programming(GPU 体系结构与 CUDA 编程)
  • 本讲核心问题:GPU 为什么能把”同一条计算(kernel)作用在海量数据上”做得比 CPU 快十倍以上?要回答这个问题,必须同时说清两件事:(1) CUDA 编程抽象是什么(grid / thread block / thread 三级线程层次、分散的 host/device 地址空间、块内共享内存与 __syncthreads()),以及 (2) 这些抽象在现代 NVIDIA GPU 上如何被实现(线程块被 work scheduler 动态映射到 SM、block 内的线程按 32 个一组组成 warp 做 SIMT 执行、资源(寄存器/共享内存/执行上下文)决定一块 SM 上能同时驻留多少 block)。讲义把整讲收束为一句性能判据:要写好 CUDA,程序必须是”coherent execution + coalesced memory access + 高数据复用”——即 warp 内不发散、访存连续、每个字节从全局内存搬到片上后被尽量多地重复使用。

  • 涉及的主要硬件/软件机制
    • 软件侧(CUDA 抽象)__global__ / __device__ / __host__ 三种函数限定符与 kernel launch 语法 kernel<<<gridDim, blockDim, shmemBytes>>>(args)dim3、内置变量 threadIdx/blockIdx/blockDim/gridDimcudaMalloc/cudaMemcpy/cudaFree__shared__(静态与 extern 动态共享内存)、__syncthreads()、原子操作(atomicAdd 等,作用于 global 与 shared)、以及”用多次 kernel launch 代替全局同步”的 kernel decomposition 手法。
    • 硬件侧(GPU 微架构):HBM(High Bandwidth Memory,高带宽显存,讲义标注 ~1 TB/s 量级)、SM(Streaming Multiprocessor,流多处理器,即”GPU 核”)、SM 内部的 sub-core / warp selector / 寄存器堆 / SIMD 功能单元 / LSU(load-store unit)/ Tensor Core、片上的 shared memory + L1 存储、L2 与全局内存、以及 GPU work scheduler(线程块调度器)。讲义还给出两代硬件的对照:GTX 980(2014,16 SM)→ H100(2022,132 SM、共享内存 256 KB/SM、峰值算力从 4.6 TFLOPs 提升到 1000 TFLOPs,主要来自 Tensor Core)。
    • 抽象与实现的接口:编译产物(CUDA device binary 里除了指令文本,还带有”资源需求”清单:每块多少线程、每线程多少本地数据、每块多少共享内存)——这正好是 work scheduler 做资源匹配的依据。
  • 在并行计算知识体系中的角色:本讲是全课程从”CPU 侧并行”跳到”加速器侧并行”的转折点。它承接第 3、4 讲的 SIMD / 多线程 / SPMD / 数据并行模型,把它们放进一个真实而极端的机器里:GPU 用”很多个弱核心 + 硬件多线程(latency hiding)+ 32 宽 SIMT”换取吞吐量。往后看,本讲的 warp/coalescing/共享内存分块三条线索会直接支撑后续的性能优化(Performance Optimization)、异构与硬件专用化、以及并行深度学习(tensor core、算子融合、数据复用)等讲座;实验课里 CUDA 的 matmul/reduction 优化也几乎全部是本讲两个 case study 的延伸。

  • 配套材料
    • lectures/05-CUDA-programming.pdf(对应抽取文本 extracted/05-CUDA-programming.txt,共 83 页幻灯片):已公开,可在 https://www.cs.cmu.edu/~418/lectures/ 下直接下载。Fall 2026 日程表(https://www.cs.cmu.edu/~418/schedule.html)中,第 5 讲位于 Sep 2,主题写作 “GPU Architecture and CUDA Programming”,第 6 讲(Sep 4)是它的续讲。讲义首页标题写作 “Lecture 5 & 6: GPU Architecture & CUDA Programming”,即一份讲义覆盖两次课,属正常现象。
    • doc_CUDA-recitation.pdf(对应抽取文本 extracted/doc_CUDA-recitation.txt,共 39 页):CUDA recitation(习题/辅导课)讲义,已公开。它把 CUDA 从”怎么写”讲到”怎么写快”:三种函数限定符、dim3 与线程索引、host/device 内存管理的 7 个步骤、归约的四种渐进实现(kernel-as-barrier / __syncthreads() / 单 block / 多 block + 共享内存),以及 matmul 的两级分块。该讲义页脚标注的是历史学期(Fall 2025),属讲义沿用现象,不是错误。
    • cs149_supp/gpuarch.txt:Stanford CS149(Fall 2025)Lecture 7 “GPU Architecture & CUDA Programming”(共 74 页)——已公开的补充读物,也是本笔记中 V100 SM 微架构、warp selector、work scheduler 逐步调度、persistent thread 等细节的来源。
    • 讲课录像(Panopto / YouTube):Fall 2026 日程表中被注释隐藏,属 未发布
    • Ed 讨论区、Autolab、Canvas:需登录,非公开。
    • 历史学期 PDF(如 Performance Analysis/Profiling、Transactional Memory 等 Fall 2026 尚未在公开目录发布的讲义)位于 /afs/cs/academic/class/15418-*/public/ 之下,需要 CMU 登录,属未公开
    • Fall 2026 授课教师为 Brian Railing 与 Dimitrios Skarlatos;课程由 Kayvon Fatahalian 创建。

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

2.1 GPU 的来路:从 3D 渲染管线到 GPGPU

  • 定义与目的:GPU 原本只做一件事——3D 渲染。渲染任务的定义可以压成一句话:给定场景描述(三角形网格、材质、光源、相机),计算网格中每个三角形对图像中每个像素外观的贡献。这条定义天然是数据并行的:同一个”着色程序”被套用在海量顶点流、片元流、像素流上。理解这段历史的目的,是理解 CUDA 为什么长成现在这个样子——它的抽象是从图形管线里”长出来”的,而不是凭空设计的。
  • 直观解释(”它是什么?”):把渲染管线想成一条印刷流水线:先送来一堆”点坐标”(vertex),拼成”图章”(primitive,三角形),把图章盖到纸上决定盖住哪些格子(rasterization,光栅化 → fragment),对每个格子算颜色(fragment processing / shader),最后写进成品(pixel operations)。程序员能插手的地方是”算颜色”这一站:他写一个小程序(shader),管线把这个小程序自动套用到整条片元流上。类比:流水线上的喷漆工序——工厂(管线)决定工件怎么运,喷漆配方(shader)由你提供,配方会被自动应用到每一个经过的工件上。
  • 架构/机制图解:早期 GPU 只有一条通路——图形管线;2007 年 NVIDIA Tesla 架构才第一次给 GPU 装上”compute mode“接口:应用可以分配显存缓冲区、上传一个 kernel 二进制、然后说”以 SPMD 方式跑 N 份这个 kernel”(launch(myKernel, N))。这个操作比 drawPrimitives(vertex_buffer) 简单得多,从此 GPU 有了两条入口。

图 1:从图形管线接口到 compute 接口(GPGPU 的诞生)

  [A] 图形接口:2007 年之前通往 GPU 硬件的唯一通路
      应用(经图形驱动)提交 "draw" 命令 + shader 二进制 + 顶点缓冲

      vertex buffer → vertex proc → primitive gen → rasterization → [fragment shader] → pixel op → framebuffer
                                                          │                ▲
                                                          │                └─ 硬件把程序员写的 shader
                                                          │                   "map" 到整条片元流上
                                                          └─ 关键:一次 shader = 一个片元/像素

  [B] compute 接口:2007 年 NVIDIA Tesla 架构打开的第二个入口(GPGPU / CUDA)

      cudaMalloc / cudaMemcpy  →  launch(myKernel, N)  →  以 SPMD 方式跑 N 份 kernel
      对比:launch(myKernel, N) 这个操作比 drawPrimitives(vertexBuffer) 简单得多
            GPU 从此不再只能执行"图形管线计算"
  • 性能特征与历史坐标:2001–2003 年的人发现”GPU 是极快的处理器,它在做同一种计算(shader),作用于大量数据(顶点/片元/像素)“——这正是 90 年代超级计算机上的数据并行。于是出现了一个著名的 hack:把 OpenGL 的输出图像尺寸设成数组尺寸(如 512×512),画两个刚好铺满屏幕的三角形,于是”一个像素 = 一个数组元素 = 一次 shader 计算“,GPU 就被当成数据并行机用了。这就是 GPGPU(General-Purpose computing on GPU)的起点;随后 Brook(2004,Stanford 图形实验室)把”抽象 GPU 为数据并行处理器”变成一门语言:kernel void scale(...) + scale(amount, input_stream, output_stream),编译器再把流程序翻译成 drawTriangles 之类的图形命令。CUDA(2007)是把这套思路正式固化下来的产物:“C-like 语言 + 低抽象距离(low abstraction distance)”

2.2 SIMD 与 SIMT:warp、coherence 与 divergence

  • 定义与目的:SIMD(Single Instruction, Multiple Data,单指令多数据)指”同一条指令被广播到多个 ALU 上并行执行“;用增加 ALU 数量的办法提升算力,代价是控制逻辑被多个 lane 共享,因此一旦 lane 之间指令不同,就有人必须空转。CUDA 的执行载体是 warp:一个 block 里连续的 32 个 CUDA thread 被绑定为一个 warp,在硬件上以 32 宽 SIMD 方式推进。NVIDIA 把这套”32 条独立标量指令流被硬件动态检测出’恰好相同’后合并执行”的机制叫 SIMT(Single Instruction, Multiple Thread,单指令多线程)
  • 直观解释(”它是什么?”)龙舟类比。32 名桨手(warp 里的 32 个线程)共用一面鼓(warp selector + fetch/decode):鼓点一响,32 支桨同时入水(一条指令同时作用在 32 个数据上),这就是 coherent execution(一致执行)。但如果有人想向左划、有人想向右划(if (x > 0) ... else ...),船不可能同时向两边走——只能先左划一遍(把想右划的人按住不出力,mask 掉),再右划一遍,这就是 divergent execution(发散执行),代价是效率按分支数成比例下降。SIMT 与 ISPC 的 gang 很像,但有一个关键区别:ISPC 在编译期就把程序编译成 SIMD 指令,而 GPU 是在运行期由硬件动态判断 32 个独立线程是否恰好走在同一条指令上;换言之,warp 不是 CUDA 编程模型的一部分,而是实现细节(例外是 __shfl_*、vote 这类 intra-warp 原语)。

图 2:SIMT 执行流与分支发散(divergence)

  warp 的 32 条线程共享一条指令流              if (x > 0) { x = 2*x; } else { x = exp(x,5); }
  ┌────────────────────────────────────┐
  │  warp 指令流(逐条取出、广播)      │        lane:  0    1    2    3  ...  31
  │  00: LDG  r0, [r5]                 │        x>0:  T    F    T    F   ...  T
  │  01: FSETP.GT p0, r0, 0            │              │    │    │    │        │
  │  02: @p0 FMUL r1, r0, 2.0          │◀───────┐     └────┴──┬─┴────┴────────┘
  │  03: @!p0 ...exp...                │        │            │
  │  04: STG  [r6], r1                 │        │      ┌─────┴─────┐
  └──────────────┬─────────────────────┘        │      ▼           ▼
                 │  同一条指令 → 32 个 ALU       │  ┌────────┐  ┌────────┐
                 ▼  同时执行(coherent)         │  │ pass A │  │ pass B │
      ┌───────────────────────────┐              │  │ 掩码   │  │ 掩码   │
      │ ALU ALU ALU ... ALU (×32) │──────────────┘  │ T F T F│  │ F T F T│
      └───────────────────────────┘                 │ 有效/无效 lane 交替 → 有效算力减半
                                                    └────────┘  └────────┘
                    发散代价 ≈ 分支路径数 × 指令数(此处 2 条路径 → 最坏 2 倍)
  • 性能特征:warp 内 32 个 lane 走同一路径时,一条 32 宽 SIMD 指令就能服务全部 32 个 CUDA 线程;一旦发散,硬件按掩码(mask)分多次执行,有效吞吐量下降为”活跃 lane 数 / 32”。讲义明确把”coherence execution”列为”高效使用 GPU 的必要条件”,把”divergent execution”列为”应在 CUDA 程序中最小化的东西”。

2.3 CUDA 执行模型:grid / thread block / thread 三级层次

  • 定义与目的:CUDA 用一个两级层次 + 三级索引来把一个数据并行问题切成”可独立调度的工作单元”:线程(thread)→ 线程块(block)→ 网格(grid)。这样做的目的是让同一个程序在 6 核 GPU 和 16 核 GPU 上都正确、且尽量高效地运行,而不需要程序里出现 num_cores——把”有多少个核”这件事交给硬件调度器决定。
  • 直观解释(”它是什么?”)建筑工地类比grid 是”整个工程项目”,block 是一个”施工班组”,thread 是”班组里的一名工人”。项目经理(host 代码)只说”这个活总共要切成 8000 个班组,每班 128 人,各自负责自己编号那一小段”,而从不说”派到几号楼”——现场调度员(GPU work scheduler)看哪栋楼(SM)还有空位、资源够不够,就把班组塞进去。班组之间互不通信、无先后依赖(”thread blocks can be executed in any order”),所以怎么排都不会错;而同一个班组内部的 128 名工人是真正同时在场的(他们要靠共享内存和 __syncthreads() 协作),所以硬件必须为整班人同时准备工位(寄存器 + 执行上下文)。

图 3:CUDA 线程层次 → 硬件映射

   Host(CPU,串行 C/C++ 程序)
   ┌───────────────────────────────────────────────────────────────────────────────┐
   │  dim3 threadsPerBlock(4,3,1);  dim3 numBlocks(Nx/4, Ny/3, 1);                 │
   │  matrixAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);   ← "launch a grid"       │
   │  调用在所有线程结束前不返回(implicit barrier across all threads)              │
   └───────────────────────────────┬───────────────────────────────────────────────┘
                                   │ ① 下发 kernel 命令(含 NUM_BLOCKS、资源需求)
                                   ▼
                        ┌───────────────────────┐
                        │   GPU Work Scheduler  │  ② 按"资源是否放得下"动态分配 block
                        └───────────┬───────────┘
        ┌───────────────────────────┼───────────────────────────┐
        ▼                           ▼                           ▼
  ┌─────────── SM 0 ──────────┐ ┌────── SM 1 ──────┐      ┌──── SM k ────┐
  │ Registers / Shared / L1   │ │                  │      │              │
  │ ┌─ Block(0,0) ──────────┐ │ │ ┌─ Block(1,0) ─┐ │ ...  │              │
  │ │ 128 threads = 4 warps │ │ │ │              │ │      │              │
  │ │ w0: t0..t31  w1: t32..│ │ │ │              │ │      │              │
  │ │ w2: t64..t95 w3: t96..│ │ │ │              │ │      │              │
  │ │ shared mem: 520 B     │ │ │ │              │ │      │              │
  │ └───────────────────────┘ │ │ └──────────────┘ │      │              │
  └───────────────────────────┘ └──────────────────┘      └──────────────┘
     ▲ gridDim / blockIdx 决定"我负责数据的哪一块";blockDim / threadIdx 决定"我在块里的编号"
     │ 全局线程号:i = blockIdx.x*blockDim.x + threadIdx.x(二维时 j 同理)
  • 语法与语义要点(讲义给的例子):__global__ 表示 kernel 函数在 GPU 上运行、由 host 调用;__device__ 是只能被 device 代码调用的函数;__host__ 是普通 CPU 函数。kernel 的调用数量由程序显式写出,而不是由数据集合的大小自动决定(这与图形着色器的 map(shader, stream) 不同)——因此当 Nx=11blockDim.x=4 不能整除时,块数要向上取整,kernel 内部必须写 if (i < Nx && j < Ny) 做越界保护。host 与 device 代码的分离是程序员在源码里静态划定的

2.4 GPU 硬件:SM 微架构与 warp selector

  • 定义与目的:SM(Streaming Multiprocessor)是 GPU 的”核”。它内部由若干个 sub-core 组成,每个 sub-core 有自己的 warp selector、取指/译码单元、矢量功能单元、寄存器堆分片。SM 的设计目的是用大量硬件线程(warp)填满功能单元的空隙:功能单元的延迟以时钟计,但访存/依赖造成的停顿以数十上百时钟计,只有同时驻留几十个 warp、每个时钟挑一个”可发射”的,才能把吞吐量顶到峰值。
  • 直观解释(”它是什么?”)机场安检类比。sub-core 是安检通道,warp selector 是”挑下一位旅客上前”的安检员,寄存器堆是每个旅客的随身行李格。如果当前这位旅客行李还没过机(load 未返回),安检员不会傻等,而是立刻换下一位(换个 warp)——这就是硬件多线程(也称 interleaved multithreading)隐藏延迟的方式。另外,sub-core 上的功能单元天生是 SIMD 的:以 CS149 补充讲义里的 V100 sub-core 为例,每个 sub-core 有 16 个 fp32 ALU(”32 宽 SIMD 操作每 2 个时钟一次”)、16 个 int ALU、8 个 fp64 单元(”32 宽操作每 4 个时钟一次”)、一组 LSU 和 Tensor Core。

图 4:一个 SM 的内部结构(以 V100 量级的 SM 为参照)

 ┌────────────────────── Streaming Multiprocessor (SM) ───────────────────────┐
 │  ┌── sub-core 0 ──────────┐ ┌── sub-core 1 ──────────┐                     │
 │  │ Warp Selector          │ │ Warp Selector          │      ...×4 sub-core │
 │  │  ├─ 每时钟选 1 个 warp  │ │                        │                     │
 │  │  └─ 16 个 warp 的上下文 │ │                        │                     │
 │  │ ┌────────────────────┐ │ │                        │                     │
 │  │ │ fp32 SIMD unit ×16 │ │ │                        │  64 KB 寄存器/sub-core │
 │  │ │ int  SIMD unit ×16 │ │ │                        │  256 KB 寄存器/SM      │
 │  │ │ fp64 SIMD unit ×8  │ │ │                        │  寄存器按 (最多)64 个   │
 │  │ │ Load/Store unit    │ │ │                        │  warp 切分             │
 │  │ │ Tensor Core unit   │ │ │                        │                     │
 │  │ └────────────────────┘ │ │                        │                     │
 │  │ Fetch / Decode         │ │                        │                     │
 │  └────────────────────────┘ └────────────────────────┘                     │
 │  ┌───────────────────────────────────────────────────────────────────────┐  │
 │  │  片上存储:Shared Memory + L1 cache(V100 约 128 KB;H100 共享内存256KB)│  │
 │  └───────────────────────────────────────────────────────────────────────┘  │
 └──────────────────────────────────┬──────────────────────────────────────────┘
                                    │ L2 cache(V100: 6 MB)
                                    ▼
                     GPU memory / HBM(V100: 16 GB,≈900 GB/s;教学讲义标注 ~1 TB/s 量级)
  • warp 的调度语义与状态:block 一旦开始执行,”block 内所有线程就都已存在且已分配寄存器状态“(这是 CUDA 的语义约束,不是实现自由)。原因是 block 内线程之间可以有依赖——最简单的依赖就是 __syncthreads()。因此不能”先跑完线程 0–127,再跑线程 128–255”来”实现”一个 256 线程的 block:如果用 128 个线程先去等 256 线程的 barrier,就会死锁。CUDA 的语义保证是:block 内线程并发运行,只要一个线程是 runnable 的,它最终一定会被运行(不会饿死)

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

  • 定义与目的:GPU 的 work scheduler 是硬件里的一块专用逻辑,负责把 grid 里的线程块按资源需求映射到各个 SM 上。它存在的意义是:让同一份 CUDA 程序在任意核数的 GPU 上都能跑,且尽量跑满。
  • 直观解释(”它是什么?”)“班车 + 座位”类比。kernel 编译产物里写着”我这个班有多少人、需要多少共享内存、每人占多少寄存器”;调度器就像调度员,看见哪辆车(SM)上还有连续的座位(执行上下文)和放行李的地方(共享内存),就安排一个班组上去。放不下的班组就在门口等,绝不硬塞。“问题 → 子问题(任务)→ 工作线程池”这个分解-分配模式在本课程已经出现过多次:ISPC 的 task 启动、Web 服务器的线程池(线程数由核数决定,而不是由请求数决定)都是同一模式,这里是它在硬件里的版本。

图 5:线程块调度的时间线(讲义用”两个 SM 的虚构 GPU、1000 个 block、每块 128 线程 + 520 B 共享内存”逐步演示)

   时间 →
   SM0 上下文: [0..127] [128..255]        (每 SM 仅 384 线程上下文 / 1.5 KB 共享内存)
   ┌──────────────────────────────────────────────────────────────────────────┐
   │ t0  host 下发: EXECUTE=convolve, NUM_BLOCKS=1000, 需求 128 thr / 520 B     │
   │ t1  scheduler: Block0 → SM0(contexts 0-127,  shared 0x000)   NEXT=1      │
   │ t2            : Block1 → SM1(contexts 0-127,  shared 0x000)   NEXT=2      │
   │ t3            : Block2 → SM0(contexts 128-255, shared 0x520)  NEXT=3      │
   │ t4            : Block3 → SM1(contexts 128-255, shared 0x520)  NEXT=4      │
   │ t5  第三个 block 放不下(3×520 B = 1560 B > 1.5 KB)→ 必须等占位者完成     │
   │ t6  Block0 完成 → 资源回收 → Block4 → SM0(contexts 0-127)     NEXT=5      │
   │ t7  Block2 完成 → 资源回收 → Block5 → SM0(contexts 128-255)   NEXT=6      │
   └──────────────────────────────────────────────────────────────────────────┘
   关键点:block 的"起飞顺序"与编号无关(可以任意序);决定并行度的是**资源**而非线程数
  • 性能特征:因为 block 会被”分波(wave)”调度,每块占用的资源越多,同一 SM 上能同时驻留的 block 就越少,隐藏延迟的能力(occupancy)就越差。讲义给出的资源账本非常具体:一个 128 线程、130 个 float 共享内存的卷积 block,需要的资源是 128 个线程 + 520 B 共享内存 + 每线程若干字节本地内存。把 block 改成 256 线程就会让每块共享内存翻倍(258 floats),在同样的 SM 上可驻留的 block 数减半——“共享内存”和”寄存器”是限制并行度的两种硬约束

2.6 CUDA 存储模型:三种地址空间与两条拷贝路径

  • 定义与目的:CUDA 显式暴露不同的局部性层级给程序员:per-thread 私有(寄存器/本地内存)、per-block 共享(__shared__)、per-program 全局(global)。暴露它们的目的是让程序员能把”复用距离”控制在片上:全局内存慢而大,共享内存快而小,寄存器最快但容量由寄存器文件决定(这一设计哲学与第 3 讲”存储层次反映程序的局部性”完全一致)。
  • 直观解释(”它是什么?”)车间类比。全局内存是”中央仓库”(巨大、远、单次往返几百个时钟);共享内存是”班组共用的小白板”(就挂在车间墙上,几十个时钟,但只有 96–256 KB,且只有同一个班组的成员能看);寄存器是”工人手上的工具”(1 个时钟,但一人只有几十到二百多个格子)。block 之外的人看不到你的白板,所以跨 block 通信只能走仓库(全局内存),而且没有全局屏障可用。

图 6:CUDA 内存层次与数据复用(以 matmul 分块为例)

   Host (CPU, 串行)                    Device 全局内存 (GPU 所有线程可见, HBM, 慢而大)
   ┌──────────────────────────┐        ┌────────────────────────────────┐
   │  float* A = new float[N] │        │  deviceA   latency: 数百时钟    │
   │  cudaMalloc(&deviceA, ..)│◀──────▶│  deviceB   V100: 16 GB @900GB/s │
   └──────────────────────────┘ memcpy └────────────────────────────────┘
        Host↔Device 拷贝                      ▲                │
        (PCIe / NVLink, 最贵的一环)   coalesced│                │结果写回
                                              │                │
              ┌───────────────────────────────┴────────────────┴─────────────┐
              │  __shared__ sA[S][L] / sB[S][L]   (per-block, 片上 ~20-30 时钟) │
              │  一个 block 只把 LxL 瓦片从全局搬进来一次, 之后复用 L 次         │
              └───────────────────────────────┬──────────────────────────────┘
                                              │ 每线程取 V 个元素
                                 ┌────────────▼─────────────┐
                                 │ 寄存器 a[V], b[V], c[V][V]│ (per-thread, 1 时钟)
                                 │ V*V 次 FMA / 2V 次取数     │
                                 └──────────────────────────┘
   复用阶梯(越小越快越贵):寄存器 > 共享内存 > L2 > 全局内存/HBM
   设计规则:把"复用距离"压到最内层;把"全局流量"降到与分块尺寸成反比

表 1:CUDA 各地址空间对照(量级为公开规格的常见近似值)

地址空间声明方式作用域 / 生命周期典型延迟量级容量约束谁可以读写
寄存器 register自动变量(kernel 内局部标量)单线程 / 线程生命期~1 个时钟每线程有限(常见上限 255 个寄存器)仅本线程
本地内存 local溢出的自动变量/数组单线程 / 线程生命期类似全局(实际在 DRAM,走 L1/L2 缓存)由寄存器压力决定仅本线程
共享内存 shared__shared__ float s[L][L];extern __shared__block 内全部线程 / block 生命期~20–30 个时钟每 SM 96 KB(教学用 GTX 980)→ 168 KB(A100)→ 256 KB(H100)同 block 所有线程
全局内存 globalcudaMalloc + kernel 指针参数全程序 / 直到 cudaFree~400–800 个时钟(HBM)显存容量(V100: 16 GB)所有线程,host 需经 cudaMemcpy
常量/只读空间*__constant__ / __ldg全程序 / 只读命中常量缓存时很快64 KB 常量缓存所有线程只读

* 常量/只读空间属于 CUDA 提供的另一类片上缓存路径,讲义正文只强调了”per-thread / per-block / global”三类,这里作为补充。

  • 主从地址空间与拷贝路径:host 与 device 是两个互不可见的地址空间——device 代码不能解引用 host 指针,host 代码也不能直接读写 device 指针指向的内容(deviceA[i] 在 host 里是非法操作)。数据必须通过 cudaMemcpy(dst, src, bytes, kind) 显式搬运(cudaMemcpyHostToDevice / DeviceToHost / DeviceToDevice)。recitation 讲义把它写成一个 7 步流程:① host 分配 → ② 初始化输入 → ③ cudaMalloc 分配 device 内存 → ④ cudaMemcpy 上传 → ⑤ launch kernel → ⑥ cudaMemcpy 取回结果 → ⑦ cudaFree + free。CUDA 后来也提供 Unified Virtual Memory(统一虚拟内存)来简化这件事,但性能模型没有变:跨 PCIe/NVLink 的搬运仍然是全程序里最贵的环节之一(见 §4 的数值算例)。

2.7 访存效率:coalescing 与共享内存 bank conflict

  • 定义与目的:GPU 的内存系统为”多个线程访问连续地址“做了优化:一个 warp 的 32 个 lane 如果访问的是连续的内存字,硬件把它们合并成尽量少的、宽度接近 128 B 的事务,这叫 coalesced access(合并访存);如果访问是散乱的(stride 很大),同一个 warp 的请求会被拆成很多个事务,带宽利用率骤降。共享内存则把存储切成 bank(通常 32 个 bank,4 字节宽):一个 warp 的 32 次访问如果都落在不同 bank,一个周期就能完成;如果两个 lane 撞在同一个 bank 的不同字上,就叫 bank conflict(体冲突),访问被串行化。
  • 直观解释(”它是什么?”):coalescing 像快递合单:32 个人要寄东西到相邻门牌号,快递员一次拉走一整箱(1 个事务);如果他们要寄到城市各处,快递员就得跑 32 趟(32 个事务)。bank conflict 像银行 8 个柜台的叫号:如果 32 个人被均匀分到 8 个柜台,各办各的(无冲突,全速);如果一半人全挤在 2 个柜台,其他人干等(冲突,串行化,带宽减半)。

图 7:coalesced 与 non-coalesced 访存;无冲突与有冲突的共享内存访问

  ① 全局内存:合并访存(optimal)
     地址:   0    4    8   12   16   20   24   28
           ┌────┬────┬────┬────┬────┬────┬────┬────┐
     第1次 │ t0 │ t1 │ t2 │ t3 │ t4 │ t5 │ t6 │ t7 │  → 1 个事务服务 8 个 lane
           └────┴────┴────┴────┴────┴────┴────┴────┘
  ② 全局内存:非合并访存(suboptimal,地址间隔 = 4 个字)
           ┌────┬────┬────┬────┬────┬────┬────┬────┐
     第1次 │ t0 │    │    │    │ t1 │    │    │    │  → 每个 lane 都要一个事务
           └────┴────┴────┴────┴────┴────┴────┴────┘

  ③ 共享内存无冲突:lane t0..t7 访问 word 0,1,2,3,4,5,6,7
        B0  B1  B2  B3  B4  B5  B6  B7
        t0  t1  t2  t3  t4  t5  t6  t7     每个 bank 服务 1 个 lane → 1 个周期,满带宽
  ④ 共享内存 2-way 冲突:lane t0..t7 访问 word 0,2,4,6,8,10,12,14
        B0     B1   B2     B3   B4     B5   B6     B7
        t0,t4   —   t1,t5   —   t2,t6   —   t3,t7   —    4 个 bank 各服务 2 个 lane
                                                          → 串行化 2 个周期,带宽减半
  • 性能特征:这两条规则解释了 CUDA 优化的绝大部分收益来源。归约的”交错寻址(interleaved addressing)”版本之所以慢,正是因为它同时踩了两个坑:s 很小时 warp 内 tid % (2s) == 0 的判定让活跃 lane 稀疏分布(发散),而且 sdata[tid] + sdata[tid+s] 的第二次访问步长是 2(2-way bank conflict)。改成”顺序寻址”(s < blockDim.x/2 递减、判据 tid < s)后,活跃 lane 集中在前缀(warp 级别不发散,整 warp 可以整体跳过),且两次访问都是连续地址(无 conflict)。

2.8 同步:__syncthreads() 与”CUDA 没有全局屏障”

  • 定义与目的:block 内部需要协作(例如”大家先把瓦片搬进共享内存,再开始算”),所以 CUDA 提供了 block 级屏障 __syncthreads()、原子操作(atomicAdd 等,可作用于 global 与 shared),以及”kernel 返回时对所有线程的隐式屏障”。但没有 block 之间的同步——这不是疏漏,而是设计选择。
  • 直观解释(”它是什么?”)__syncthreads()班组内部的”点名”:所有人都到齐了才能继续(仅限同一车间/同一 block)。而跨车间没有点名机制:车间之间只能通过”下一班次”(下一个 kernel)来对表——kernel 的结束/启动本身就是一个全局同步点,代价极低(讲义强调”kernel launch has very low hardware/software overhead”)。
  • 为什么没有全局同步:讲义给了两条理由:(1) 在 SM 数量极多的 GPU 上做硬件全局屏障代价高昂;(2) 当 block 数 > SM 数 × 每 SM 可驻留 block 数时,等待其它 block 的 block 会永远等不到(死锁)。补救办法是把计算分解成多个 kernel(kernel decomposition):所有层次用同一段代码逻辑,靠多次 launch 串起”全局同步”。
同步手段作用范围谁来用代价典型用途
__syncthreads()同一 thread block 内所有线程block 内协作数十个时钟(warp 到齐即可放行)共享内存”装载 → 使用”之间
原子操作 atomicAdd/Inc/Exch/...global 或 shared 地址任意线程冲突时串行化,吞吐随冲突度下降直方图计数、全局部分和累加
kernel 边界(隐式屏障)整个 gridhost 通过多次 launch每次 launch 数微秒级开销取代不存在的全局屏障(多级归约)
自旋 + 全局标志(实践中危险跨 block手工实现可能死锁/饿死只在”所有 block 保证同时驻留”时才勉强成立

2.9 Tensor Core:为什么 H100 的峰值能到 1000 TFLOPs

  • 定义与目的:Tensor Core 是 SM 里的矩阵乘加专用单元(讲义原话:”Matrix multiplication unit in SMM”)。它把”小矩阵乘(如 4×4 或更大)的乘加阵列”做成硬件,使单位时钟完成的 MAC 数量级远高于通用 fp32 SIMD 单元——这正是从 GTX 980 的 4.6 TFLOPs 到 H100 的 1000 TFLOPs 的主要原因(讲义明确说明”mainly because of tensor cores”)。
  • 直观解释(”它是什么?”):通用 SIMD 单元像通用机床,什么都能做但一次只能加工一个工件;Tensor Core 像专用冲压模具:只做”矩阵乘”这一种活,但一次冲出一整排。代价是你必须有足够多的矩阵乘数据喂给它,否则模具空转——这正是后续”算术强度”分析要量化的东西:算力越强,达到峰值所需的算术强度(FLOP/Byte)就越高。
  • 硬件演进对照(讲义 slide “GTX 980 (2014) → H100 (2022)”)
参数GTX 980 (2014)H100 (2022)说明
时钟频率1.1 GHz1.11 GHz几乎没变
每 SM 可驻留 warp 数6464没变
warp 宽度(线程/warp)3232没变
每 SM 共享内存96 KB256 KB(A100 为 168 KB)片上容量随工艺增长
SM 数量161328.25×
峰值算力4.6 TFLOPs1000 TFLOPs(主要靠 tensor core)≈217×

3. 代码示例与性能分析

3.1 示例 A:1D 卷积的两个版本(直接全局访存 vs 共享内存暂存)

代码convolve_1d.cu)——对应讲义 “1D convolution version 1 / version 2”:

// 编译: nvcc -O3 -arch=sm_70 convolve_1d.cu -o convolve_1d
//       (release 优化用 -O3;需要行号信息可加 -lineinfo;只做调试用 -g -G,会关闭优化)
// 运行: ./convolve_1d 1048576
#include <cuda_runtime.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#define THREADS_PER_BLK 128

#define CUDA_CHECK(call)                                                     \
  do {                                                                       \
    cudaError_t err__ = (call);                                              \
    if (err__ != cudaSuccess) {                                              \
      fprintf(stderr, "CUDA error %s at %s:%d\n",                            \
              cudaGetErrorString(err__), __FILE__, __LINE__);                \
      exit(EXIT_FAILURE);                                                    \
    }                                                                        \
  } while (0)

// ---------- Version 1: 每个线程直接读 3 个全局元素 ----------
__global__ void convolve_v1(int N, const float* __restrict__ input,
                            float* __restrict__ output) {
  int index = blockIdx.x * blockDim.x + threadIdx.x;  // 全局线程号 = 输出下标
  if (index >= N) return;                             // 越界保护
  float result = 0.0f;
#pragma unroll
  for (int i = 0; i < 3; i++) result += input[index + i];
  output[index] = result / 3.0f;
}

// ---------- Version 2: 一个 block 协作把 support 区搬进共享内存 ----------
__global__ void convolve_v2(int N, const float* __restrict__ input,
                            float* __restrict__ output) {
  __shared__ float support[THREADS_PER_BLK + 2];       // 每块分配 130 个 float
  int index = blockIdx.x * blockDim.x + threadIdx.x;
  int tid = threadIdx.x;
  support[tid] = input[index];                         // 130 次装载取代 3*128 次
  if (tid < 2) support[THREADS_PER_BLK + tid] = input[index + THREADS_PER_BLK];
  __syncthreads();                                     // 全 block 屏障
  float result = 0.0f;
#pragma unroll
  for (int i = 0; i < 3; i++) result += support[tid + i];
  if (index < N) output[index] = result / 3.0f;
}

int main(int argc, char** argv) {
  const int N = (argc > 1) ? atoi(argv[1]) : (1 << 20);   // 默认 1M 个输出元素
  // 输入需要 N+2 个元素;再额外留 2*THREADS_PER_BLK 的尾部填充,
  // 这样即使 N 不是 block 大小的整数倍,块内"偷看"后续 128 个元素的读也不会越界
  const int inLen = N + 2 + 2 * THREADS_PER_BLK;
  const size_t inBytes = sizeof(float) * inLen;
  const size_t outBytes = sizeof(float) * N;

  float* h_in = (float*)malloc(inBytes);
  float* h_out = (float*)malloc(outBytes);
  float* h_ref = (float*)malloc(outBytes);
  for (int i = 0; i < inLen; i++) h_in[i] = sinf(0.001f * i);
  for (int i = 0; i < N; i++)     h_ref[i] = (h_in[i] + h_in[i + 1] + h_in[i + 2]) / 3.0f;

  float *d_in = nullptr, *d_out = nullptr;
  CUDA_CHECK(cudaMalloc(&d_in, inBytes));
  CUDA_CHECK(cudaMalloc(&d_out, outBytes));
  CUDA_CHECK(cudaMemcpy(d_in, h_in, inBytes, cudaMemcpyHostToDevice));

  dim3 block(THREADS_PER_BLK);
  dim3 grid((N + THREADS_PER_BLK - 1) / THREADS_PER_BLK);  // 向上取整
  cudaEvent_t t0, t1;
  CUDA_CHECK(cudaEventCreate(&t0));
  CUDA_CHECK(cudaEventCreate(&t1));

  for (int which = 1; which <= 2; ++which) {
    CUDA_CHECK(cudaMemset(d_out, 0, outBytes));
    CUDA_CHECK(cudaEventRecord(t0));
    for (int rep = 0; rep < 20; ++rep) {                  // 多次重复取平均
      if (which == 1) convolve_v1<<<grid, block>>>(N, d_in, d_out);
      else            convolve_v2<<<grid, block>>>(N, d_in, d_out);
    }
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    float ms = 0.f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, t0, t1));
    CUDA_CHECK(cudaMemcpy(h_out, d_out, outBytes, cudaMemcpyDeviceToHost));

    double maxErr = 0.0;
    for (int i = 0; i < N; i++) maxErr = fmax(maxErr, fabs(h_out[i] - h_ref[i]));
    // 有效带宽 = (逻辑上的读 + 写) 字节数 / 时间(不计尾部填充)
    double gb = (double)(sizeof(float) * (N + 2) + outBytes) / 1e9;
    printf("version %d: %.3f ms/iter, 有效带宽 %.1f GB/s, maxErr=%.2e\n",
           which, ms / 20.0, gb / (ms / 20.0 / 1e3), maxErr);
  }

  CUDA_CHECK(cudaFree(d_in));
  CUDA_CHECK(cudaFree(d_out));
  free(h_in); free(h_out); free(h_ref);
  return 0;
}
  • 【代码做什么?】 host 端按 recitation 讲义给出的 7 步流程走:分配 host 内存 → 初始化 → cudaMalloc device 内存 → cudaMemcpy 上传 → launch kernel(convolve<<<N/THREADS_PER_BLK, THREADS_PER_BLK>>>)→ 拷回 → 释放。kernel 里 index = blockIdx.x * blockDim.x + threadIdx.x 是全局线程号,一个线程算一个输出元素 output[index] = (input[index]+input[index+1]+input[index+2])/3。版本 2 的差别在于:block 的 128 个线程先协作把 130 个输入值(本块需要的 support 区)搬进共享内存__syncthreads() 之后各自从共享内存里取自己需要的 3 个值。
  • 【并行机制与性能解说】:kernel 启动后,每个 block(4 个 warp)被 work scheduler 分配到某个 SM 上;SM 的每个 sub-core 每时钟从驻留 warp 中挑一个可发射的 warp,warp 里 32 个 lane 同时处理 32 个相邻输出。访存方面,input[index]完全连续的,所以 v1 的三次装载都是 coalesced;v2 的 support[tid] = input[index] 同样连续,但装载指令数从每线程 3 次降到每 block 130 次(3×128 = 384 → 130,约 2.95×)。
    • Work(总工作量):Θ(N)。v1 每线程 3 次全局 load + 2 次 add + 1 次 store;v2 每 block 130 次全局 load + 各自 3 次共享 load。Work 的阶不变,改变的只是”从哪里取数”。
    • Span(关键路径):Θ(1)。每个线程内部的依赖链长度是常数(3 次加法),block 内唯一的同步点是 1 个 __syncthreads()。若把 span 记作”最长依赖链 + 同步点数”,本 kernel 的 span = O(1) 不随 N 增长。
    • 并行度 = Work/Span = Θ(N):N = 2²⁰ ≈ 1.05×10⁶ 个可并行线程。
    • 可扩展性上限:并行度不是瓶颈——V100 可同时驻留 80 SM × 64 warp × 32 = 163,840 个 CUDA 线程,本问题需要约 6.4 “波”。真正的瓶颈是内存带宽:v1 与 v2 的 DRAM 流量几乎相同(每个输入元素仍要从 HBM 取一次、每个输出写一次),所以 v2 不会降低 DRAM 流量,它降低的是 L1/共享内存通路上的请求数与指令数。当访存模式本身不合并、或同一数据被线程反复从 L1 请求导致 L1 带宽成为瓶颈时,v2 才真正救场。

3.2 示例 B:分块矩阵乘(Case Study 1:两层复用)

代码mm_tiled.cu)——对应讲义 “Optimization 1: thread-level register tiling” + “Optimization 2: block-level shared memory tiling”:

// 编译: nvcc -O3 -arch=sm_70 mm_tiled.cu -o mm_tiled
// 运行: ./mm_tiled 1024
#include <cuda_runtime.h>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#define L   32   // block 级瓦片尺寸: 一个 block 计算 C 的 L x L 子矩阵
#define BLK 16   // block 维度: BLK x BLK = 256 个线程
#define V   2    // thread 级寄存器瓦片: 每个线程算 V x V 个元素

#define CUDA_CHECK(call)                                                     \
  do {                                                                       \
    cudaError_t e__ = (call);                                                \
    if (e__ != cudaSuccess) { fprintf(stderr, "%s\n", cudaGetErrorString(e__)); exit(1);} \
  } while (0)

// A, B, C 都是 N x N 的行主序矩阵
__global__ void mm_tiled(const float* __restrict__ A,
                         const float* __restrict__ B,
                         float* __restrict__ C, int N) {
  __shared__ float sA[L][L];
  __shared__ float sB[L][L];

  const int tx = threadIdx.x, ty = threadIdx.y;
  const int tid = ty * BLK + tx;
  const int rowBase = blockIdx.y * L;       // 本 block 负责 C 的 [rowBase, colBase)
  const int colBase = blockIdx.x * L;
  const int nThreads = BLK * BLK;

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

  for (int ko = 0; ko < N; ko += L) {
#pragma unroll
    for (int j = 0; j < L * L / nThreads; ++j) {     // 协作装载: 每线程装 4 个元素
      int idx = j * nThreads + tid;
      int r = idx / L, cc = idx % L;
      sA[r][cc] = A[(rowBase + r) * N + (ko + cc)];  // A 的行瓦片
      sB[r][cc] = B[(ko + r) * N + (colBase + cc)];  // B 的列瓦片
    }
    __syncthreads();                                 // 装载完成才能被使用

#pragma unroll
    for (int k = 0; k < L; ++k) {                    // 在瓦片上做 L 次外积
#pragma unroll
      for (int i = 0; i < V; ++i) a[i] = sA[ty * V + i][k];
#pragma unroll
      for (int j = 0; j < V; ++j) b[j] = sB[k][tx * V + j];
#pragma unroll
      for (int i = 0; i < V; ++i)
#pragma unroll
        for (int j = 0; j < V; ++j) c[i][j] += a[i] * b[j];
    }
    __syncthreads();                                 // 下一轮装载前必须用完
  }

#pragma unroll
  for (int i = 0; i < V; ++i)
#pragma unroll
    for (int j = 0; j < V; ++j)
      C[(rowBase + ty * V + i) * N + (colBase + tx * V + j)] = c[i][j];
}

int main(int argc, char** argv) {
  const int N = (argc > 1) ? atoi(argv[1]) : 1024;
  const size_t bytes = sizeof(float) * (size_t)N * N;
  float *hA = (float*)malloc(bytes), *hB = (float*)malloc(bytes), *hC = (float*)malloc(bytes);
  for (int i = 0; i < N * N; i++) { hA[i] = 1.0f; hB[i] = 2.0f; }   // 结果应为 2N

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

  dim3 block(BLK, BLK);                    // 256 线程/块
  dim3 grid(N / L, N / L);                 // (N/L)^2 个块
  cudaEvent_t t0, t1; cudaEventCreate(&t0); cudaEventCreate(&t1);
  CUDA_CHECK(cudaEventRecord(t0));
  for (int rep = 0; rep < 5; ++rep) mm_tiled<<<grid, block>>>(dA, dB, dC, N);
  CUDA_CHECK(cudaEventRecord(t1));
  CUDA_CHECK(cudaEventSynchronize(t1));
  float ms = 0.f; cudaEventElapsedTime(&ms, t0, t1);
  CUDA_CHECK(cudaMemcpy(hC, dC, bytes, cudaMemcpyDeviceToHost));

  double flops = 2.0 * N * N * N;          // 每个 C 元素 N 次 MAC = 2N 次浮点运算
  printf("N=%d: %.3f ms, %.1f GFLOP/s, C[0]=%f (应为 %d)\n",
         N, ms / 5.0, flops / (ms / 5.0 / 1e3) / 1e9, hC[0], 2 * N);
  return 0;
}
  • 【代码做什么?】 每个 block 负责 C 的一个 L×L 子矩阵(L=32),块内 256 个线程各负责一个 V×V(V=2)的小块,即 16×16 个线程 × 2×2 = 32×32。主循环以 S=L 为步长沿 K 维推进:每个 block 协作把 A 的 L×S 行瓦片与 B 的 S×L 行瓦片搬进共享内存L*L/nThreads = 4 个元素/线程,访存完全连续 → coalesced),屏障后每个线程从共享内存里反复取数,做 L 次”外积”更新自己的 V×V 累加器,最后一次性写回 C。
  • 【并行机制与性能解说】:grid 有 (N/L)² 个 block,全部可乱序调度;块内 256 线程 = 8 个 warp,每个 warp 的 32 个 lane 协同装载 128 B 连续数据(完美合并)。真正的”复用”发生在两级:共享内存层复用 L 倍(同一个 A/B 瓦片被 block 内所有线程用),寄存器层复用 V 倍(同一个 a[i] 被 V 个累加器用、同一个 b[j] 被 V 个累加器用)。
    • 讲义给的流量公式(本讲的定量核心):
      • Strawman(每线程算一个元素):全局访存 = 2N³ 次;算术强度极低。
      • 寄存器分块后:每线程 2NV 次 → 总全局访存 2N³/V
      • 再叠加 block 级共享内存分块后:每 block 2NL 次 × (N²/L²) 个 block → 总全局访存 2N³/L;而共享内存访问总数是 2N³/V
      • 结论:L 与 V 是”两个复用旋钮”(讲义原话:”L and V are the two reuse knobs!”)。
    • Work:Θ(N³) 次乘加(2N³ FLOPs)。N = 1024 时为 2.15×10⁹ FLOPs。分块不改变 Work,只改变访存
    • Span:Θ(N)。关键路径是单个累加器 c[i][j] 沿 k 的 N 次串行累加(寄存器里是依赖链),这是本 kernel 无法回避的深度。
    • 并行度 = Work/Span = 2N³/N = Θ(N²) ≈ 2.1×10⁶,远大于硬件可同时驻留的线程数,因此并行度不是瓶颈,带宽/共享内存带宽才是(见 §4 的 Roofline 算例)。
    • 共享内存是否会成为新瓶颈:按讲义公式,共享内存访问次数 = 2N³/V,即”每次 FMA 对应 2/V 次共享取数“。V=2 时相当于”1 次 FMA 配 1 次共享取数”:在 V100 的 sub-core 里,一条 32 宽共享 load 若完全无冲突需要 1 个 LSU 周期,而一条 32 宽 fp32 FMA 在 16 宽 ALU 上要 2 个周期,于是 LSU 利用率约 50%——还没到瓶颈,但已经不低;V 增大(4 或 8)会把该比率降到 1/4、1/8,同时把寄存器压力抬高(c[V][V] 需要 V² 个寄存器)。这就是”分块尺寸是资源与带宽之间的取舍”的具体含义。
    • 讲义伪代码的排版问题:讲义第 64/65 页的伪代码里 b[:] = sA[ki, ...] 应为 sB[ki, ...](B 的瓦片来自 sB)。上面给出的实现按语义修正了这一点。

3.3 示例 C:并行归约(Case Study 2:树形归约 + 顺序寻址 + warp shuffle)

代码reduce.cu)——综合讲义 “reduce0 交错寻址 → strided index → sequential addressing” 三代改进与 CS149 的 two-pass 分解:

// 编译: nvcc -O3 -arch=sm_70 reduce.cu -o reduce
// 运行: ./reduce 16777216
#include <cuda_runtime.h>
#include <cstdio>
#include <cstdlib>

#define BLK 256

#define CUDA_CHECK(call)                                                     \
  do {                                                                       \
    cudaError_t e__ = (call);                                                \
    if (e__ != cudaSuccess) { fprintf(stderr, "%s\n", cudaGetErrorString(e__)); exit(1);} \
  } while (0)

// warp 内归约:用 __shfl_down_sync 在 32 个 lane 之间交换寄存器,无需共享内存与屏障
__inline__ __device__ float warpReduceSum(float v) {
#pragma unroll
  for (int offset = 16; offset > 0; offset >>= 1)
    v += __shfl_down_sync(0xffffffffu, v, offset);
  return v;
}

// 第一/第二阶段共用的 kernel:每个线程读 2 个元素,块内树形归约,写 1 个部分和
__global__ void reduce_kernel(const float* __restrict__ in,
                              float* __restrict__ out, int n) {
  extern __shared__ float sdata[];                       // 动态共享内存,长度 = blockDim.x
  const unsigned tid = threadIdx.x;
  const unsigned i = blockIdx.x * (blockDim.x * 2) + tid;

  float v = 0.0f;
  if (i < n) v = in[i];                                   // 连续访存 → coalesced
  if (i + blockDim.x < n) v += in[i + blockDim.x];        // 第二个元素,同样连续
  sdata[tid] = v;
  __syncthreads();

  // 阶段 1:顺序寻址(sequential addressing),活跃线程集中在低位前缀
  for (unsigned s = blockDim.x >> 1; s > 32; s >>= 1) {
    if (tid < s) sdata[tid] += sdata[tid + s];            // 无 bank conflict(连续地址)
    __syncthreads();                                      // 每轮都需要屏障
  }
  // 阶段 2:剩下 64 个部分和,用 1 个 warp + shuffle 收尾(不需要再屏障)
  if (tid < 32) {
    float w = sdata[tid] + ((tid + 32 < blockDim.x) ? sdata[tid + 32] : 0.0f);
    w = warpReduceSum(w);
    if (tid == 0) out[blockIdx.x] = w;
  }
}

int main(int argc, char** argv) {
  const int n = (argc > 1) ? atoi(argv[1]) : (1 << 24);
  const size_t bytes = sizeof(float) * n;
  float* h_in = (float*)malloc(bytes);
  double ref = 0.0;
  for (int i = 0; i < n; i++) { h_in[i] = 1.0f; ref += 1.0f; }

  float* d_src; float* d_dst;
  CUDA_CHECK(cudaMalloc(&d_src, bytes));
  const int firstBlocks = (n + 2 * BLK - 1) / (2 * BLK);
  CUDA_CHECK(cudaMalloc(&d_dst, sizeof(float) * (firstBlocks + 1)));
  CUDA_CHECK(cudaMemcpy(d_src, h_in, bytes, cudaMemcpyHostToDevice));

  cudaEvent_t t0, t1; cudaEventCreate(&t0); cudaEventCreate(&t1);
  CUDA_CHECK(cudaEventRecord(t0));

  // 多级 kernel 分解:每一级都是"全局同步点",直到只剩 1 个元素
  int cur = n, launched = 0;
  float* src = d_src; float* dst = d_dst;
  while (cur > 1) {
    int blocks = (cur + 2 * BLK - 1) / (2 * BLK);
    reduce_kernel<<<blocks, BLK, BLK * sizeof(float)>>>(src, dst, cur);
    CUDA_CHECK(cudaGetLastError());
    ++launched;
    cur = blocks;
    float* tmp = src; src = dst; dst = tmp;               // swap 双缓冲
  }
  CUDA_CHECK(cudaDeviceSynchronize());
  CUDA_CHECK(cudaEventRecord(t1));
  CUDA_CHECK(cudaEventSynchronize(t1));
  float ms = 0.f; cudaEventElapsedTime(&ms, t0, t1);

  float total = 0.f;
  CUDA_CHECK(cudaMemcpy(&total, src, sizeof(float), cudaMemcpyDeviceToHost));
  printf("n=%d: %.3f ms(%d 次 kernel launch,%.1f GB/s),sum=%f 期望=%f\n",
         n, ms, launched, (double)bytes / (ms * 1e-3) / 1e9, total, ref);
  return 0;
}
  • 【代码做什么?】 归约是归一化、softmax 等算子的基础原语(讲义原话)。这里用两级分解:第一级把 n 个元素切成 2*BLK 一组的块,每块树形归约成 1 个部分和写回全局内存;第二级(以及后续级)继续对部分和做同样的事,直到只剩 1 个数。因为 CUDA 没有跨 block 同步,多级之间的”全局同步”由 kernel 边界提供。recitation 讲义列出了四种做归约的思路:把每轮循环做成一个 kernel(正确但慢,launch 开销大)、一个 kernel 内用 __syncthreads()错误——它只保证块内)、只启动 1 个 block(正确但没法用满 GPU,适合小数据)、以及”块内 __syncthreads() + 跨块 kernel”(本示例采用的折中)。
  • 【并行机制与性能解说】:第一级有 n/(2·256) = 32768 个 block(n=16M),每个块 8 个 warp;sdata[tid] += sdata[tid+s]顺序寻址,两次访问都落在连续地址上,无 bank conflict;等到 s < 32 时,只有 warp 0 需要继续工作,其余 7 个 warp 因判据整体为假而直接退出,不再发射指令(这正是”stride 版 vs 顺序版”的关键差别),最后用 __shfl_down_sync 在 warp 内做无共享内存的收尾。
    • Work:Θ(n) 次加法(n−1 次)。n=2²⁴ 时 ≈ 1.68×10⁷ 次加法。
    • Span:Θ(log n)(树高)+ O(1) 次 kernel 边界。n=2²⁴ 时块内树高 log₂(256)=8 层,跨级约 24/8 ≈ 3 层 kernel,总同步深度约 11–12 个”阶段”。
    • 并行度 = Θ(n/log n) ≈ 7×10⁵(n=16M;span 取树高 log₂n = 24 层)。
    • 瓶颈:这个 kernel 是纯带宽受限。n=16M 个 float = 64 MB,V100 上 HBM 带宽 900 GB/s → 至少 64 MB / 900 GB/s ≈ 71 µs;而算力需求只有 1.68×10⁷ 次加法(算术强度 1.68e7 FLOP / 6.71e7 B ≈ 0.25 FLOP/Byte),在 12.7 TFLOPs 机器上只要约 1.3 µs——两者相差 50 倍以上。两次以上的 kernel launch(每次数微秒)在总时间中占 5–15%,所以”减少 kernel 数、用 grid-stride loop 或 persistent thread 把多级合并“是这类 kernel 的常见优化方向。
    • 版本对比(讲义三代实现)
版本归约判据 / 索引warp 内是否发散共享内存 bank每轮 warp 活跃情况结论
V1 交错寻址if (tid % (2*s) == 0) sdata[tid] += sdata[tid+s](活跃 lane 稀疏:T F T F…)2-way conflict(第二次访问步长 2)每个 warp 都要为 T/F 两条路径发指令最慢,两个坑都踩
V2 strided indexindex = 2*s*tid; if (index < blockDim.x) ...否(活跃 lane 连续在前缀)仍有步长 s 的访问模式部分 warp 整体不活跃快于 V1
V3 顺序寻址for (s = blockDim.x/2; s > 0; s /= 2) if (tid < s) sdata[tid] += sdata[tid+s]无冲突(两次都是连续地址)高位 warp 整块退出,不再发射讲义给出的推荐形态
V4 warp shuffle 收尾前 31 → 32 用 __shfl_down_sync完全不使用共享内存只剩 1 个 warp 干活减少屏障次数与共享内存流量

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

4.1 Work-Span 模型在 GPU 上的映射

  • Work(W):整个计算的总操作数(与处理器数量无关)。Span(S,关键路径长度):必须串行执行的最长依赖链。并行度(parallelism)= W / S,表示”用无限多执行单元时最多能利用的并发度”。
  • 在 GPU 上,硬件能提供的”处理器数量”不是 CPU 的核数 p,而是可同时驻留的 CUDA 线程数P_HW = SM 数 × 每 SM 可驻留 warp 数 × 32。V100 上是 80 × 64 × 32 = 163,840
  • 因此 GPU 的加速上限有两道门:min(W/S, P_HW) 决定”并行度够不够”,而带宽/共享内存/指令发射带宽决定”喂不喂得饱”。本讲两个 case study 都撞在第二道门上
程序WorkSpan并行度 W/S实际瓶颈说明
1D 卷积(N 个输出)Θ(N)Θ(1)Θ(N) ≈ 10⁶HBM 带宽每字节只做约 1 次运算
分块 matmul(N×N)Θ(N³) = 2.15×10⁹ FLOPsΘ(N) = 1024Θ(N²) ≈ 2.1×10⁶全局访存 + 共享内存带宽分块提高算术强度
树形归约(n 个元素)Θ(n) ≈ 1.68×10⁷Θ(log n) ≈ 24 层Θ(n/log n) ≈ 7×10⁵HBM 带宽 + kernel 启动开销每字节只做 0.25 次运算

4.2 Amdahl 定律与”The Free Lunch Is Over”的 GPU 版本

若程序中不可并行的比例为 f,则 p 个执行单元的加速比上限是

S(p) = 1 / ( f + (1 - f) / p )         p → ∞ 时 S_max = 1 / f
  • CPU 上的 f:串行段、I/O、锁竞争。
  • GPU 上的 f 至少多了三项:(1) host↔device 的 cudaMemcpy(慢通路,且常常无法与计算重叠);(2) 每次 kernel launch 的固定开销(数微秒);(3) 无法并行的 host 侧逻辑(分配、校验、启动)。
  • 数值算例(N=1024 的 fp32 matmul,V100 量级)
    • 计算量 2.15×10⁹ FLOPs。kernel 理想时间 = 2.15e9 / 12.7e12 ≈ 169 µs(100% 峰值;按 §4.3 的实际可达约 300 µs)。
    • 数据量:A、B、C 各 1024²×4 B = 4 MB,共 12 MB。若通过 PCIe Gen3 x16(≈16 GB/s 实测)搬运:12 MB / 16 GB/s ≈ 750 µs
    • 于是 f = 750/(750+300) ≈ 0.71,加速比上限 1/0.71 ≈ 1.4——GPU 的 12.7 TFLOPs 被数据搬运吃掉了。这正是 CUDA 优化的第一条工程准则:让数据留在显存里、把多次 kernel 连成流水(避免中间结果来回搬运),也正是后续讲座里”算子融合 / 减少 host-device 往返”的动机。
    • 若改成”数据已在显存、只做纯计算”(f→0),加速比就由带宽与算力决定,此时 GPU 相对单核 CPU(假设 20 GFLOP/s 有效)可达数百倍。

4.3 算术强度与 Roofline 模型(本讲最重要的定量工具)

算术强度(arithmetic intensity) AI = 浮点运算次数 / 从全局内存搬运的字节数(FLOP/Byte)。Roofline 模型给出可达性能:

可达性能 = min( 峰值算力 P_peak ,  AI × 内存带宽 BW )
"拐点"算术强度 = P_peak / BW

数值算例 1:matmul 的三种实现(N = 1024,V100:P_peak = 12.7 TFLOPs,BW = 900 GB/s)

  • 拐点 AI = 12.7e12 / 900e9 ≈ 14.1 FLOP/Byte
  • Strawman(每线程算一个元素):全局访存 2N³ 次 × 4 B = 8.6×10⁹ B,FLOPs 2.15×10⁹ → AI = 0.25 FLOP/Byte(= 每字节 1/4 次运算)。可达性能 = 0.25 × 900 GB/s = 225 GFLOP/s,只有峰值的 1.8%
  • 寄存器分块(V=2):总全局访存 2N³/V = 1.07×10⁹ 次 × 4 B = 4.3×10⁹ B → AI = 0.5,可达 ≈ 450 GFLOP/s。
  • 再叠加 block 共享内存分块(L=32):每 block 访存 2NL 次、共 (N/L)² 个 block → 总访存 2N³/L = 6.7×10⁷ 次 × 4 B = 268 MBAI = 2N³/(4·2N³/L) = L/4 = 8 FLOP/Byte。可达 = min(12.7 T, 8×900 G) = 7.2 TFLOP/s(仍是内存受限)。时间下界 = 268 MB / 900 GB/s ≈ 298 µs,而算力下界 169 µs → 268 MB 的搬运是主约束,有效效率约 57%。
  • 把 L 提到 64:AI = L/4 = 16 > 14.1 → 跨过拐点,变成算力受限:访存时间 134 MB/900 GB/s ≈ 149 µs < 169 µs,上限就是算力。共享内存需求 = 2 × 64 × 64 × 4 B = 32 KB/block,在 V100 的 128 KB 片上存储里可容纳。
  • 通用结论(简洁且好用):按讲义公式,float32 分块 matmul 的算术强度 ≈ AI = L/4 FLOP/Byte。因此”要撑满峰值算力需要多大的分块”可以由 L ≥ 4·P_peak/BW 反推。
  可达 GFLOP/s(示意,纵轴对数刻度)
   12.7 T ┤───────────────────────────────────────────────────┬───────────────  算力屋顶
          │                                                   │
          │                                          ┌────────┘
          │                                    ┌─────┘
    7.2 T ┤                              ┌─────┘   ● L=32: AI=8 → 7.2 TFLOP/s(内存受限)
          │                        ┌─────┘
          │                  ┌─────┘
          │            ┌─────┘   斜率 = 900 GB/s(带宽屋顶)
          │      ┌─────┘
   450 G  ┤  ┌───┘ ● V=2: AI=0.5
   225 G  ┤● │     strawman: AI=0.25
          └─┴────┴────┴────┴────┴────┴────┴────┴────┴────┴────► 算术强度 FLOP/Byte
            0.25  0.5   2    8   14.1  16   32   64
                                  ↑ 拐点 P_peak/BW = 14.1(V100 fp32)

数值算例 2:H100 为什么必须用 Tensor Core

  • 讲义给的 H100 峰值 1000 TFLOPs(tensor core),HBM 带宽标注 ~1 TB/s 量级。此时拐点 AI = 1000e12 / 1e12 = 1000 FLOP/Byte
  • 若按公开规格近似(132 SM、每 SM 128 个 fp32 FMA/时钟、1.11 GHz)估算 fp32 通用算力 ≈ 132×128×2×1.11e9 ≈ 37.5 TFLOP/s,拐点仍需 AI ≈ 37.5 FLOP/Byte,对应 L ≥ 150——远超共享内存与寄存器能维持的分块规模
  • 因此”算力越强,越需要靠降低数据搬运字节数(更深的复用、更低精度、更大分块、warp-level 的矩阵乘指令)来维持算术强度”。Tensor Core 把”每字节能换多少运算”整个抬高了一个量级:它做的是”在片上一小片寄存器/共享数据上完成更多 MAC”,等于在硬件里内置了更高阶的复用。

数值算例 3:延迟隐藏与 Little 定律

  • HBM 延迟量级 ≈ 500 ns,带宽 900 GB/s。要打满带宽,需要”在途(in-flight)”的字节数 = BW × latency = 900e9 × 500e-9 = 450 KB
  • 分摊到 80 个 SM:5.6 KB/SM,即约 1400 个 4 字节的未完成 load 每 SM。除以 4 个 sub-core 得 350 个 lane-load/sub-core ≈ 11 个满 warp 的 load
  • 而每个 sub-core 最多驻留 16 个 warp。结论:必须让每个驻留 warp 平均携带多个未完成 load(ILP),并且让 SM 上驻留足够多的 warp(occupancy),才可能打满 HBM。这从数量上解释了讲义里反复强调的三件事:访存要合并(1 个 warp 的 load 才等于 1 个 128 B 事务,而不是 32 个)不要用发散把 warp 切碎不要用过多共享内存/寄存器把可驻留 block 数压下去

4.4 一个把三项放在一起的定量对照表

指标公式1D 卷积(N=1M)分块 matmul(N=1024, L=32)归约(n=16M)
Work2N 次加法 + N 次除法 ≈ 3.1×10⁶ FLOPs2N³ FLOPs ≈ 2.15×10⁹n−1 次加法 ≈ 1.68×10⁷
SpanΘ(1)Θ(N) = 1024Θ(log n) ≈ 24
并行度W/S≈ 1.05×10⁶≈ 2.1×10⁶≈ 7×10⁵
全局内存流量≈ 8 MB(读 4 MB + 写 4 MB)268 MB(strawman 8.6 GB)≈ 64 MB(读一次)
算术强度FLOP/Byte≈ 0.378(strawman 0.25;L=64 时 16)≈ 0.25
带宽时间下界bytes / BW8 MB / 900 GB/s ≈ 8.9 µs268 MB / 900 GB/s ≈ 298 µs64 MB / 900 GB/s ≈ 71 µs
算力时间下界FLOPs / P_peak≈ 0.24 µs169 µs≈ 1.3 µs
实际受限方带宽(约 37× 于算力)带宽(约 1.8× 于算力)带宽(约 55× 于算力)
并行度是否足够vs 163,840约 6.4 波约 13 波约 4.3 波

5. 关键要点

  1. CUDA 的三个抽象支柱:线程层次 + 分层地址空间 + 块内屏障。 grid/block/thread 三级索引(blockIdx*blockDim+threadIdx)让”问题分解”与”机器规模”解耦(程序里没有 num_cores);host/device 分离的主从地址空间要求显式 cudaMemcpy__syncthreads() 只在 block 内有效,而 kernel 边界就是 CUDA 提供的唯一全局同步手段(多级 kernel 分解因此是最重要的结构性手法)。
  2. warp 是性能的基本单位(32 线程 SIMT)。 一个 block 的线程按 32 个一组绑定成 warp,block 数 = 线程数/32;warp 之内必须”同一条指令”,所以coherent execution 是高效使用 GPU 的必要条件,divergence 是主要性能损失源__syncthreads() 的存在意味着”block 内线程必须真的并发存在”,这解释了为什么硬件必须一次性为整个 block 分配寄存器与执行上下文。
  3. 两级分块 = 两级复用旋钮。 全局访存总量从 2N³(strawman)降到 2N³/V(寄存器分块)再降到 2N³/L(共享内存分块)。L 与 V 是唯二的复用杠杆,但它们同时受共享内存容量(每 SM 96 KB → 256 KB)与寄存器数量的双重限制,分块放大会牺牲 occupancy。
  4. 一切定量判断都可以用”算术强度 vs 拐点”来做。 Roofline:可达 = min(P_peak, AI×BW),拐点 = P_peak/BW(V100 fp32 ≈ 14.1 FLOP/Byte)。matmul 的 AI = L/4,L=32 时只有 8 → 内存受限;L=64 才跨过拐点。算力越强(tensor core),维持峰值所需的算术强度越高,这是 H100 时代编程方式的根本约束。
  5. 访存模式(coalescing / bank conflict)与延迟隐藏决定了”带宽能不能被吃满”。 一个 warp 访问连续地址才能合并成少量事务;共享内存 32 个 bank 只有”一个 lane 一个 bank”才无冲突(步长 2 的访问直接减半带宽)。同时按 Little 定律,要撑满 900 GB/s 需要约 450 KB 的在途数据,这必须靠”高 occupancy + 每 warp 多个未完成 load”来提供,因此不要为了分块把共享内存/寄存器用到没有余量。

6. 常见陷阱与注意事项

  • __syncthreads() 当成”全局屏障”,以及跨 block 自旋等待。 __syncthreads() 只同步同一个 thread block:recitation 讲义里的 reduce_2 就是典型错误——在多个 block 上跑、循环里放 __syncthreads(),跨 block 的依赖完全没有保证(”CUDA 不保证锁步执行”)。方向相反的错误同样致命:CS149 补充讲义里 block 0 自增 myFlag、block 1 自旋等 myFlag != 0 的写法,在”每 SM 只能驻留 1 个 block”的机器上若调度器先跑 block 1 就会永久死锁。正确做法是块内用 __syncthreads()、跨块用 kernel 边界(或只启动 1 个 block,仅在小数据时合理)。合法使用原子操作 ≠ 可以假设 block 之间有执行顺序;CUDA 只承诺”可以把 block 按任意顺序调度”。
  • 忽略发散(divergence)带来的浪费。 if (tid % (2*s) == 0) 这类判据让 warp 内活跃 lane 稀疏分布,硬件不得不为两条路径分别发射指令;正确做法是改成连续前缀判据(tid < s)或转成步长索引,让”整 warp 活跃/整 warp 不活跃”。注意:发散不一定错,但代价要算清楚(归约尾部只有 32 个线程干活是可接受的;如果一个 kernel 的每个 warp 长期只有 1/32 的 lane 有效,那就是设计问题)。
  • 把非合并访存当成小事,或以为”用了共享内存就一定会更快”。 一个 warp 的 32 个 lane 访问步长远大于 1 的地址(例如按列访问行主序矩阵)时,一个请求会被拆成 32 个事务,有效带宽掉到 1/32——矩阵乘里 B[k][y] 的列方向访问就是典型场景,分块到共享内存的一部分动机正是把非合并的全局访问换成合并的瓦片装载。但反过来,共享内存不是免费的:它有 32 个 bank,冲突会串行化访问,而且它是按 block 分配的稀缺资源。示例 A 的 v2 把全局 load 指令从 3N 降到 N+2,却没有减少 DRAM 流量(每个输入元素仍要从 HBM 取一次),所以纯带宽受限时几乎不会变快——它的收益在减少 L1 压力与指令数。
  • 忽视 host↔device 搬运与 kernel 启动的固定开销。 §4.2 的算例里,12 MB 的 PCIe 搬运(≈750 µs)可以轻松压过 300 µs 的计算,使加速比上限只有 1.4;此外每次 kernel 启动都有数微秒级开销,”为了用 kernel 边界做全局同步而拆成几十次 launch”会把开销累积成主要成本。优化的方向是让数据驻留显存、合并 kernel、用 grid-stride/persistent thread 减少 launch 次数。
  • 把 CUDA 线程当 pthread 来理解(抽象与实现的混淆)。 CUDA 线程与 pthread 在”逻辑控制流”这一层相似,但实现截然不同:pthread 的创建要分配栈与控制块、由 OS 调度;而 CUDA 的百万级线程不会真的分配百万份栈——block 只有在被调度到 SM 上时才”实化”寄存器状态与共享内存,且同一 block 内的所有线程必须同时实化。混淆这一点会得出错误结论(例如”线程越多越好”)。
  • 为分块把资源吃干净,导致 occupancy 崩塌或寄存器溢出(spill)。 c[V][V] 累加器随 V² 增长,共享内存随 L² 增长;一旦超过”每 SM 共享内存/寄存器除以每 block 占用”,可驻留 block 数下降,延迟隐藏能力变差,甚至寄存器溢出把数据挤到本地内存(实际走 DRAM),性能出现断崖。分块尺寸必须按目标 GPU 的资源上限反推。

7. 思考题(带答案)

问题 1(执行模型 / 正确性):recitation 讲义里有一个”错误”的归约实现:

__global__ void reduce_2(float* A, int length) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  for (int d = 1; d < length; d <<= 1) {
    if (i < length && i % (d << 1) == 0) A[i] += A[i + d];
    __syncthreads();                 // 这里加了块内屏障
  }
}

它比”不加任何屏障”的版本更”安全”了吗?为什么它仍然不正确?请给出两个正确的替代方案,并说明各自的适用场景。

【答案】 仍然不正确。__syncthreads() 只对同一个 thread block 内的线程有效,而这个 kernel 让多个 block 共同处理同一个数组 A:第 d 轮需要第 d−1 轮的全部结果,但第 d−1 轮的结果可能由另一个 block 写出的,块内屏障完全管不到它。所以不同 block 之间仍然存在未同步的读写依赖(数据竞争),结果不确定。两个正确方案:(a) 把每一轮做成一次 kernel launchfor (d=1; d<length; d<<=1) reduce_1<<<grid,block>>>(A, length, d);)——kernel 边界提供全局同步,语义正确;缺点是 launch 次数多、每次数微秒开销,适合 d 的轮数不多时(例如先把 length 降到较小值再逐级 kernel)。(b) 只用 1 个 block + __syncthreads()reduce_2_1):启动 <<<1, blockDim>>>,块内循环里用 for (; i < length; i += blockDim.x) 处理整个数组,屏障就是正确的全局屏障。它只在 length 较小(一个 block 能覆盖)时高效,因为一个 block 只能占一个 SM,GPU 其余部分闲置。(c) 更实用的第三种:混合方案——块内用 __syncthreads() 把每个 block 的一段数组归约成 1 个部分和,跨块再用 kernel launch(或让第二个 kernel 只有 1 个 block)合并,这正是本讲示例 C 采用的两级分解。

问题 2(性能模型 / 定量):某 GPU 的 fp32 峰值算力为 12.7 TFLOP/s,HBM 带宽 900 GB/s。你要实现 N=1024 的 fp32 矩阵乘 C=A×B,用共享内存分块,块尺寸 L×L。请回答:(1) 分块后的算术强度是多少?(2) 要使程序变成算力受限,L 至少要多大?(3) L=32 时的时间下界与可达 GFLOP/s 各是多少?(4) 若换成峰值 1000 TFLOP/s、带宽 1 TB/s 的 H100(假设仍用 fp32 通用单元),这个结论会怎样变化?

【答案】 (1) 全局访存总量是 2N³/L 次(每次 4 字节),浮点运算 2N³ 次,于是 AI = 2N³ / (2N³/L × 4 B) = L/4 FLOP/Byte每字节 1/4 次运算,即”分块把算术强度提升到 L/4”。 (2) 拐点算术强度 = P_peak / BW = 12.7e12 / 900e9 ≈ 14.1 FLOP/Byte。要跨过拐点需 L/4 ≥ 14.1,即 L ≥ 56.4,工程上取 L = 64(此时 AI = 16)。L=64 时共享内存需求为 2 × 64 × 64 × 4 B = 32 KB/block,在 96–256 KB 的片上空间内是可实现的。 (3) L=32 时:AI = 8 FLOP/Byte。可达性能 = min(12.7 T, 8 × 900 G) = 7.2 TFLOP/s(内存受限,只有峰值的 57%)。全局流量 = 2N³/L × 4 B = 2×1.074e9/32 × 4 B ≈ 268 MB,时间下界 = 268 MB / 900 GB/s ≈ 298 µs;此时算力下界为 2.15e9/12.7e12 ≈ 169 µs,所以带宽是主约束(298 µs > 169 µs)。 (4) 若峰值提高到约 37.5 TFLOP/s(132 SM × 128 fp32 FMA/时钟 × 2 × 1.11 GHz,公开规格近似)而带宽只有 ~1 TB/s,则拐点变成约 37.5 FLOP/Byte,需要 L ≥ 150——共享内存和寄存器都放不下这么大的瓦片。这说明:通用 fp32 SIMD 的算力提升会把程序推向”必然内存受限”的境地,因此硬件把复用作进了 Tensor Core(把”每字节的运算数”在单元内部放大),而软件必须走更深的复用、更低精度或 warp 级矩阵指令。这就是讲义”从 GTX 980 的 4.6 TFLOPs 到 H100 的 1000 TFLOPs,主要来自 tensor cores”这一句背后的性能模型含义。

问题 3(实现细节 / 资源与调度):假设卷积 kernel 的每个 block 需要 128 个线程和 130 个 float 的共享内存(520 B),而目标 GPU 的每个 SM 有 384 个线程的执行上下文和 1.5 KB 共享内存。请解释:(1) 一个 SM 上最多能同时驻留几个 block?(2) 如果把这个 kernel 的 block 改成 256 线程(共享内存相应变成 258 个 float),驻留情况会怎样变化?(3) 为什么硬件不允许”用 128 个线程分两批跑完一个 256 线程的 block”?

【答案】 (1) 两个约束同时起作用:线程上下文 384 / 128 = 3 个,共享内存 1536 B / 520 B = 2.95 → 2 个(第 3 个需要 1560 B > 1.5 KB,放不下)。因此取较小者:每 SM 最多 2 个 block(256 个 CUDA 线程、1040 B 共享内存),即共享内存是先撞上的那道墙。这也说明 occupancy 的瓶颈常常不是线程数而是片上存储。 (2) block 变成 256 线程后:线程上下文 384 / 256 = 1 个,共享内存 1536 / 1032 = 1 个 → 每 SM 只能驻留 1 个 block 的 256 个线程。相比 (1) 的 2×128 = 256 个线程,驻留的线程总数相同,但独立 block 数从 2 降到 1,warp 数从 8 降到 8(256/32),看上去一样,然而同步粒度变粗__syncthreads() 要让 8 个 warp 全部到齐,屏障等待的方差更大,且 block 结束时释放的是整块资源——调度灵活性下降。这解释了为什么”block 大小”是一个需要调参的实现细节。 (3) 不可以,因为 CUDA 语义规定 block 内的所有线程并发执行,且它们之间可以有依赖。最简单的依赖就是 __syncthreads():如果硬件只放得下 128 个线程,却”先跑 0–127 到 barrier、再跑 128–255”,那么先跑的那批会在 barrier 上等一批还没启动的线程,永久死锁。所以硬件必须在 block 开始时为整块预留寄存器与执行上下文;CUDA 的保证是”只要某线程是 runnable 的,它最终一定会被运行”。