Lecture 1: 课程概述、并行计算动机与 GPU 架构导论 (对应 Lab 0: Device Query)

目录 · ← l0 · l2 →

第二部分:分讲学习笔记

共 12 讲,严格按实验/讲座顺序组织。每讲结构统一为: 概述 → 核心概念与 GPU 架构图解 → 代码示例与性能分析 → 性能优化技巧总结 → 关键要点 → 常见陷阱与注意事项 → 思考题(带答案)。 所有 .cu 代码均为完整可编译程序,编译命令见各代码块头部注释,统一使用 nvcc -O3 -arch=sm_80


Lecture 1: 课程概述、并行计算动机与 GPU 架构导论 (对应 Lab 0: Device Query)

概述

本讲回答三个层层递进的问题:为什么要在 2004 年之后转向并行与异构计算(功耗墙、Dennard 缩放失效、多核普及、GPU 在算力与带宽上的双重优势);怎么衡量并行化的收益与极限(加速比、效率、可扩展性、Amdahl 定律、并行开销分类);以及GPU 长什么样(众核吞吐量处理器、SM(Streaming Multiprocessor,流多处理器)、SIMT 执行模型、延迟隐藏、内存带宽与算术强度、统一着色器/Tesla 架构如何让 CUDA 成为可能)。最后落到工程实践:CUDA 程序的五步模板(分配、拷贝入、启动、拷贝回、释放)、__global__ / __device__ / __host__ 限定符、<<<grid, block>>> 执行配置语法的含义,并用 Lab 0 的 Device Query 程序把上面所有抽象概念变成屏幕上可读的真实数字。

本讲统一使用的硬件档案(后续所有性能分析都以此为准,避免数字前后矛盾):

GPUSM 数每 SM warp 上限每 SM 线程上限每 SM 寄存器每 SM 共享内存显存带宽FP32 峰值
RTX 4090 (Ada, sm_89)12848153665536100 KB~1008 GB/s~82.6 TFLOPS
A100 (Ampere, sm_80)10864204865536164 KB~1555 GB/s~19.5 TFLOPS
H100 SXM (Hopper, sm_90)13264204865536228 KB~3350 GB/s~67 TFLOPS
RTX 2080 Ti (Turing, sm_75)683210246553664 KB~616 GB/s~13.4 TFLOPS

课程实验环境:UIUC 的 Lab 与 Project 通过 NCSA 的 Delta 超算提交,一个 Delta GPU 计算节点为 AMD Milan CPU(单路 64 核,~2.45 GHz,256 GB RAM)加 4 张 NVIDIA A40(48 GB GDDR6、696 GB/s 显存带宽、PCIe Gen4 64 GB/s、NVLink 112.5 GB/s 双向、300 W)。历史上课程也在 A100(sm_80)上运行;代码统一用 nvcc -arch=sm_80 编译,性能分析以 A100 为基准机,必要时对比 RTX 4090。

核心概念与 GPU 架构图解
功耗墙与并行计算的历史转向(The Power Wall and the 2004 Pivot)
  • 定义与目的:功耗墙(power wall)指 2000 年代中期芯片功耗密度增长到风冷散热无法承受的极限,使得”继续提高主频”这条自 1985 年以来最有效的性能提升路径突然失效。它解决的不是某个算法问题,而是整个产业必须换一条性能增长曲线的问题:从”把单个核做快”改为”把很多核做多”。

  • 直观解释(”它是什么?”):把 CPU 想象成一辆跑车。过去几十年工程师的做法是不断加大发动机排量(提高主频),每 18 个月车速翻倍。到了 2004 年,发动机已经热到会把自己熔掉——散热器(风冷)跟不上了,而继续加排量只会让整台车起火。于是厂商换了个思路:不再造一辆 300 km/h 的跑车,而是造 64 辆 80 km/h 的小车一起拉货。总运力(吞吐量)继续增长,但任何单件货物的送达速度(延迟)反而变慢了。这就是”多核转向”的本质代价。

  • 架构/机制图解:下面是 1985 至 2010 年主频与功耗的走势示意(讲义引用了 Karl Rupp 的 microprocessor-trend-data;对应讲义 “Frequency Scaled Too Fast 1993-2003” 与 “Total Processor Power Increased” 两页):

时钟频率 (MHz, 对数轴)                    功耗 / 芯片 (W, 对数轴)
 10000 |                        ......       1000 |
       |                  .....                   |              .....
  1000 |           .......                       100 |        ....
       |      .....                                 |    .....
   100 |  ....                                     10 | ...
       |..                                            |.
    10 |                                             1 |
       +----+----+----+----+----+----+----+----+      +----+----+----+----+
       85  87  89  91  93  95  97  99  01  03        85  87  89  91  93  95  97
                  ^                                     ^
                  |                                     |
       1993-2003 主频年增 ~40%/年              总功耗随之逼近风冷极限
                            \                         /
                             \                       /
                          ====  2004 年:功耗墙  ====
                          Intel 取消了一款产品(Tejas),
                          硬件转向并行:多核芯片成为标配。
                          今天已经很难买到单核的台式机/笔记本,
                          甚至很难买到单核的智能手机。
  • 关键操作与性能特征
    • 2004 年之前,CPU 性能提升几乎全部来自主频与指令级并行(ILP)。Pentium 4 “Prescott”(2004)达到 3.8 GHz、功耗约 115 W,继续往上做就要突破 150 W 风冷极限。
    • 1970–2004 年主频从 ~0.1 MHz 涨到 ~3.8 GHz(约 38000 倍),但 2004–2025 年主流桌面主频只从 3.8 GHz 涨到 ~5.5 GHz(约 1.4 倍)。主频红利在 2004 年一次性消失
    • 性能提升的重心由此转移到:多核更温和的主频大量使用向量执行同时使用延迟导向核与吞吐导向核3D 堆叠封装换取更大内存带宽。这五条正是讲义 “But Our Chips (Mostly) Still Don’t Burn — So what changed?” 一页给出的答案。
    • 量子化的后果:并行提供的是固定倍数的收益(受 Amdahl 定律限制),不会改变算法复杂度;没有任何代码或输入能扩展到无限资源(讲义 “Parallelism Does Not Affect Algorithmic Complexity”)。此外,一片芯片即使算力无穷,它仍然只有一片芯片的内存容量与 I/O 带宽(讲义 “Not Everything Can Fit on a Chip”)——这直接解释了后面”算术强度”为什么会成为 GPU 编程的核心约束。
    • 历史注脚:在 2000 年代后期之前,并行计算在学术界占比不到 5% 的研究群体,工业界关注度极低,因为”没有收益只有痛苦(No gain, just pain)”;老笑话是”只有研究生才写并行程序”。真正的转折发生在 2004 年(功耗墙)2007 年(CUDA 发布)
Dennard 缩放失效与延迟/吞吐量两种设计哲学(Dennard Scaling, Latency-Oriented vs Throughput-Oriented Design)
  • 定义与目的:Dennard 缩放(1974, JSSC)是 MOSFET 尺寸等比缩小(constant field scaling)的理想模型,它保证了”晶体管变小 → 更快 → 更省电”的良性循环。当它在 2004 年前后失效,芯片设计被迫一分为二:延迟导向(latency-oriented) 的 CPU 和 吞吐量导向(throughput-oriented) 的 GPU。理解这两者的取舍,是理解一切 CUDA 优化技巧的总纲。

  • 直观解释(”它是什么?”):CPU 像一个高级餐厅的服务员:一次只服务一桌,但每一步都极快——他记得客人要什么(分支预测),提前把菜端到桌边(数据前送),还会用小推车(大缓存)把常用调料备在手边。他处理”一道复杂工序”的单件延迟极低,但一次只能端一盘。GPU 像快餐店的一百个窗口:每个窗口的店员只会”装一个汉堡”这一招(简单控制、无分支预测),但他不介意排队(长延迟、深流水),只要队伍足够长,出餐的总速率就极高。点一份汉堡,快餐店比高级餐厅慢;点一万份汉堡,快餐店快十倍以上。

  • 架构/机制图解

        CPU:延迟导向设计 (Latency Oriented)          GPU:吞吐量导向设计 (Throughput Oriented)
   +-------------------------------------+      +-------------------------------------+
   |  Control (复杂控制)                 |      |  Control (极简控制)                 |
   |   分支预测 / 乱序 / 数据前送        |      |   无分支预测 / 无数据前送           |
   +-------------------------------------+      +-------------------------------------+
   |  Cache (大缓存)                     |      |  Cache (小缓存,只为提升吞吐)       |
   |   L1 32-48KB  L2 8-32MB  L3 32MB+   |      |   L1/Shared 合并  L2 40MB (A100)    |
   |   把长延迟访存变成短延迟 cache 命中 |      |   缓存不是为降低延迟,而是聚合请求  |
   +-------------------------------------+      +-------------------------------------+
   |  ALU  ALU  ALU  ALU                 |      |  ALU ALU ALU ALU ALU ALU ALU ALU    |
   |  (少量、超强、低延迟、高功耗)       |      |  (海量、简单、深流水、高能效)       |
   +-------------------------------------+      +-------------------------------------+
     主频高 (4-5.5 GHz)                           主频温和 (1.4-2.5 GHz)
     一次执行 4-8 条指令 (ILP)                    一次执行 32 条指令 (SIMT, 一个 warp)
     4-64 个硬件线程                              数万-数十万个硬件线程
  • 关键操作与性能特征
    • Dennard 缩放的数学:设 L → αL(典型 α ≈ 0.7),则有 VDD → αVDD、C → αC、I → αI。于是
      • 门延迟 Delay = C·VDD / I → α·CVDD/I,即延迟缩放为 α,频率 f → 1/α(每代约 ×1.43)
      • 单管功耗 P = C·V²·f → α·α²·(1/α) = α²(每代约 ×0.49)
      • 因此在芯片面积与总功耗不变的前提下,晶体管数量可以每代翻倍而功耗不涨——这就是摩尔定律能持续 30 年的物理基础。
    • 失效原因:阈值电压 Vth 不能随 VDD 等比下降(否则亚阈值漏电指数上升),约 2004 年 VDD 卡在 ~1.0–1.2 V,漏电功耗(静态功耗)在总功耗中占比快速上升。P 不再按 α² 下降,频率红利同时终结。
    • 数量级对比:一枚 10 代 Intel Core(10 核、14 nm)对比 NVIDIA GK110(2880 个 CUDA core、28 nm),后者在峰值吞吐量上远超前者,但单线程跑一段带分支的串行代码则慢 10 倍以上。讲义给出的结论是:”CPU 在串行部分比 GPU 快 10 倍以上,GPU 在并行部分比 CPU 快 10 倍以上;赢的策略是同时使用两者(Winning Strategies Use Both CPU and GPU)”。
    • 性能特征数值:CPU 访存延迟靠缓存压低到 ~4 cycle(L1 命中),但带宽只有 ~50–200 GB/s(双通道 DDR5-4800 约 76.8 GB/s,L2/L3 上百 GB/s);GPU 全局内存延迟高达 400–800 cycle,但 A100 显存带宽 1555 GB/s延迟高但对吞吐量不敏感——这正是 GPU 必须靠”海量线程”运转的原因。
加速比、效率与可扩展性(Speedup, Efficiency, Scalability)
  • 定义与目的:这三个量是并行计算界的”公共语言”。加速比回答”快了没有”,效率回答”快得值不值”,可扩展性回答”再堆机器还有用吗”。它们共同约束了所有后续优化的评价方式。

  • 直观解释(”它是什么?”):把并行化想象成请施工队盖房子
    • 加速比:本来你一个人盖要 100 天,请了 10 个人 12 天盖完,加速比 = 100/12 ≈ 8.3×。
    • 效率:你付了 10 个人的工钱,只拿到 8.3 个人的产出,效率 = 8.3/10 = 83%——有一部分工人在互相等待、传递工具、或者干着只有一个人能干的活(比如只有一个电闸)。
    • 可扩展性:请到 50 个人时,如果工地只有一条上料通道,大家开始互相挡路,加速比反而掉下来——这就是”不可扩展”。
  • 架构/机制图解(加速比 vs 处理器数的两条曲线,对应讲义 “Scalability Measures Effect of Parallel Overheads”):
Speedup
  140 |                                                        /
      |                                                      /   <-- 理想线性
  120 |                                                    /         (efficiency = 1)
      |                                                  /
  100 |                                                /
      |                                              /
   80 |                                            /
      |                                        ..--
   60 |                                   ..---'      <-- 实际曲线:P 增大后
      |                            ..----'                加速比开始"压平"
   40 |                    ..-----'
      |            ..------'
   20 |     ..----'
      |..--'
    0 +----+----+----+----+----+----+----+----+----+----+----+----> 处理器数 P
      0   10   20   30   40   50   60   70   80   90  100  110  120
                               ^
                    拐点:并行开销(通信、同步、负载不均)
                    开始与计算量可比。固定问题规模时,加速比必然压平。
  • 关键操作与性能特征
    • 定义Speedup(P) = T(sequential) / T(parallel, P),并且固定问题规模
    • T(sequential) 的选取有一条苛刻规则:必须是”在串行机器上最优的算法(未必是被并行化的那个算法)、针对串行机器优化过不含任何并行残留开销“。实践中很难做到,因此工程界的做法是——去找别人最好的实现当基线,而不是自己写一个。Hwu 教授的原话大意是:”没人会相信你为一个 baseline 付出了多少努力,即使你真的付出了。”(讲义 “Find (Don’t Write) a Competitive Baseline”)
    • 效率Efficiency(P) = Speedup(P) / P。付款方希望是 1;真实应用通常是”接近 1 但不等于 1”的某个不可忽略值。效率 > 1 称为超线性加速(superlinear speedup),极少见,可能来源是:CPU 上被容量限制的问题在 P 台机器上获得了额外资源(如总缓存容量随 P 增长,命中率提高),或者纯粹的运气(并行搜索恰好早找到答案)。
    • 单一 GPU 的 P 是什么? 是 1?是 SM 数?还是总 PE 数?讲义明确指出:对单个 GPU 而言效率的定义会含糊,实际做法是用资源占用率与该 GPU 的峰值资源对比来估计效率(这正好就是 Lab 0 Device Query 要做的事——查出峰值,然后看自己的配置用了多少)。
    • 可扩展性:问的是”在多少个处理器范围内加速比是线性的、效率是平的”。固定问题规模下,超过某个 P 之后加速比必然压平(除非让一部分处理器空转)。好的可扩展性 = 在你的机器上,直到可测的最大 P 都没有衰减
    • 面向不同场景的加速比变体(J. P. Singh, J. L. Hennessy, A. Gupta, IEEE Computer 26(7), 1993):
      • scaled speedup(规模扩展加速比):问题规模随 P 线性增长,好的结果是 1。
      • memory-constrained speedup(内存受限加速比):取”能装进(随 P 线性增长的)内存的最大问题”,只对 O(N) 算法成立。
      • time-constrained speedup(时间受限加速比):取”在我吃完午饭回来之前能算完的最大问题”。
    • 算法复杂度不变性:并行化把 T(N) 乘上一个常数因子(上限由 Amdahl 定律给出),不会改变 O(·)。因此”N 很大时并行一定赢”是错的——若算法本身是 O(N²),堆再多处理器也追不上一个 O(N log N) 的串行算法。讲义强调:并行编程的第一件事是选/造可扩展的算法,第二件事才是写 CUDA
Amdahl 定律与并行开销分类(Amdahl’s Law and Parallel Overhead)
  • 定义与目的:Amdahl 定律给出并行加速比的硬上界:加速比 ≤ 1 /(不可并行部分占比)。它解决的是”预期管理”问题——在上手写代码之前就告诉你值得投入多少。配套的”并行开销分类”则告诉你效率损失具体浪费在哪里。

  • 直观解释(”它是什么?”):想象 10 个人合作做一顿饭,但只有一个人能用炉子(烧菜占 75% 时间中的一部分必须串行,或者更贴切地说:切菜可以 10 人并行,但烧菜环节只能一个人守着一口锅)。即使切菜时间被压缩到 0,总时间也降不下来超过 1/(串行占比)。如果串行部分占 25%,那么无论请多少人,最多只能快 4 倍——第 5 个人的工资白付

  • 架构/机制图解

固定问题规模 N,串行占比 s = 0.25,可并行部分 1-s = 0.75

T(1)  = s + (1-s)              = 0.25 + 0.75 = 1.00
T(4)  = s + (1-s)/4            = 0.25 + 0.1875 = 0.4375   -> Speedup = 2.29
T(16) = s + (1-s)/16           = 0.25 + 0.0469 = 0.2969   -> Speedup = 3.37
T(∞)  = s + 0                  = 0.25          = 0.2500   -> Speedup = 4.00  <== 上界

时间条 (每格 = 0.05 单位):
P=1   [############ 串行 0.25 ############][ 可并行 0.75 ...........................................]
P=4   [############ 串行 0.25 ############][ 0.1875 ][ 空闲/等待 ...................................]
P=16  [############ 串行 0.25 ############][0.047][ 空闲 ...........................................]
P=∞   [############ 串行 0.25 ############][ 空 .................................................]
                                              ^
                          Amdahl 上界 = 1 / s = 1 / 0.25 = 4×
  • 关键操作与性能特征
    • 公式:设串行占比 s(0 ≤ s ≤ 1),处理器数 P,则 Speedup(P) = 1 / ( s + (1 - s) / P )lim(P→∞) = 1/s
    • 具体数字表(代入计算)
    sP=2P=4P=8P=16P=64P=1024上界 1/s
    0.25(并行化 75%)1.602.292.913.373.883.994.00
    0.10(并行化 90%)1.823.084.716.408.779.9110.00
    0.01(并行化 99%)1.983.887.4813.9139.2691.17100.00

    代入一例:s = 0.10、P = 64 时 1 / (0.10 + 0.90/64) = 1 / (0.10 + 0.0140625) = 1 / 0.1140625 = 8.77×,效率 = 8.77/64 = 13.7%——从 10% 串行代码到 40 核,效率就已经掉到 14%。这是 CUDA 优化中最需要记住的一条直觉:先消灭串行瓶颈,再谈并行技巧

    • 讲义给出的例子:如果只并行化了 75% 的代码,就永远拿不到超过 4× 的加速比
    • 并行开销分类(讲义 “Necessary/Good Sources” 与 “Bad Sources of Parallel Overhead”):
      • 必要且合理:(a) 搬数据(通信)——大多数并行代码的必要开销;(b) 做额外计算来省通信,例如 CNN 中在卷积之后立刻做池化(pooling)以减少”共享内存 → 全局内存”的流量,或 Kogge-Stone scan 中多做加法以减少 barrier 数量;(c) 争抢共享资源带来的优先级仲裁。
      • 纯粹浪费:(a) 干等长延迟事件(wait for long-latency events);(b) 看别人干活——GPU 上的分支发散(branch divergence)是最典型的例子;(c) 排成一列——不必要的串行化,例如过粗的同步、缺少私有化(privatization)、对共享硬件资源的时间相关访问、只用一个 CUDA stream
    • 负载均衡(Load Balance):讲义原话是”完成一个并行任务的总时间,由最慢的那个线程决定“。固定”每线程一份工作”最简单,但会导致负载不均(例如上/下三角矩阵的行长度不同)。常见解法是动态负载均衡:一个或多个工作队列,线程从队列取一块活、干完再取,先取大块后取小块,队列空了就从别的线程偷工作(work stealing)
    • 粒度(parallel grain size):即”每个线程做多少工作”。不同并行源的自然粒度不同:循环体、容器中的对象、矩阵的行/列/块/元素、图的节点/连通分量。粒度选择要比较各方案的负载方差与分支发散:矩阵元素通常工作量近似恒定(规整性好);而”循环体里带条件判断”、”每个对象调用复杂方法”、”上/下三角矩阵的行”、”图中节点的度、连通分量的大小”则方差很大。
    • 骨架式同步执行(bulk synchronous execution):HPC 与 CUDA 应用的主导风格。barrier 把代码切成时间区域(phase),每个 phase 通常只有 O(100) 行,交错与数据共享只发生在 phase 内部。好处是”调试一个区域比调试整个程序容易”;代价是会使用量在时间上相关(resource usage correlates),这是坏的——所有 warp 同时抢同一个执行单元。
  • 定义与目的:异构计算模型把”串行/控制密集部分”交给 CPU(host),把”数据并行/计算密集部分”交给 GPU(device),两者通过 PCIe/NVLink 互连。它解决的问题是:既然 CPU 与 GPU 各有所长,就不该二选一,而该让它们协作。代价是必须显式管理两套内存空间和它们之间的数据搬运。

  • 直观解释(”它是什么?”):CPU 是总部办公室(会开会、做决策、处理复杂流程),GPU 是远郊大工厂(只会重复做同一道工序,但产能惊人)。总部把原料用卡车(PCIe)运到工厂,工厂加工完再运回成品。卡车一趟只能运那么多,而且路上要花时间——如果工厂加工只要 1 分钟而运料要 40 分钟,那么优化工厂本身毫无意义。这就是 Lecture 19 的核心结论:”现代计算机系统的重力是带宽“(Bandwidth: The Gravity of Modern Computer Systems)。

  • 架构/机制图解

  CUDA 的规范五步结构与内存空间(对应讲义 "Canonical CUDA Program Structure")

  主机 (Host / CPU)                           设备 (Device / GPU)
  +----------------------+                    +----------------------------------+
  |  main() {            |                    |  __global__ void kernelOne(      |
  |                      |                    |      float* A, float* B, int N)  |
  |  (1) cudaMalloc  ----+--- 分配 global ----+--> [ Global Memory / DRAM ]      |
  |      (&d_A, bytes)   |      memory        |     d_A  d_B  d_C                |
  |                      |                    |                                  |
  |  (2) cudaMemcpy  ----+=== PCIe/NVLink ===>|     数据上行 (H2D)               |
  |      (d_A,h_A,H2D)   |                    |                                  |
  |                      |                    |                                  |
  |  (3) kernel<<<grid, -+--- 启动 grid ----->|   [SM0][SM1][SM2] .. [SM107]     |
  |      block>>>(d_A,   |     (异步!)         |     并行执行 kernelOne          |
  |      d_B, d_C, N);   |                    |     线程块被分发到各个 SM        |
  |  (4) cudaMemcpy  <---+=== PCIe/NVLink ====|     数据下行 (D2H)               |
  |      (h_C,d_C,D2H)   |                    |                                  |
  |                      |                    |                                  |
  |  (5) cudaFree(d_A) --+--- 释放 ---------->|     [ 空 ]                       |
  |  }                   |                    |                                  |
  +----------------------+                    +----------------------------------+
        主机内存                                    设备内存
   (可分页 pageable /                             (global memory,
    锁页 pinned)                                  容量 8-80 GB)
  经典 PC 架构 vs 现代 PCIe 架构(对应讲义 "Classic (Historical) PC Architecture" 与
  "Recent PCIe PC Architecture")

  【经典:Northbridge / Southbridge】
        CPU
         |
    +----+----+                 Northbridge(北桥):连接必须高速通信的三者
    | 北桥    |----------------- CPU / DRAM / 显卡
    +----+----+                 显卡需要"一等公民"级的 DRAM 访问权
         |                      早期 NVIDIA 卡走 AGP,峰值 2 GB/s
      DRAM
         |
    +----+----+
    | 南桥    |----------------- 慢速 I/O 集中器:PCI / USB / SATA / 音频
    +---------+
         |
    【原始 PCI 总线】33 MHz、32-bit、峰值 132 MB/s
                     后来 66 MHz、64-bit、峰值 528 MB/s
                     对设备而言上行带宽仍然只有约 256 MB/s 峰值
                     共享总线 + 仲裁:仲裁赢家成为 bus master
                     PCI 设备寄存器映射进 CPU 物理地址空间(memory-mapped I/O)

  【现代:PCIe 交换网络】
       CPU (内含 PCIe Root Complex / 控制器)
        |  PCIe x16 链路
    +---+-----------------------------+
    |  PCIe Switch (= 现代"北桥")     |
    +---+--------------+--------------+
        |              |
   [GPU 显卡]      [NVMe SSD]
   PCIe x16        PCIe x4
   点到点、交换式、无仲裁;报文组成虚拟通道;可对报文划分优先级做 QoS
  (例如实时视频流)。
  • 关键操作与性能特征
    • PCIe 编码与带宽:PCIe 用 128b/130b 编码(每 16 字节数据编成 130 bit,其中 1 与 0 个数相等,保证 DC 平衡与足够的跳变供时钟恢复),开销仅 1.5%;相比旧的 8b/10b 编码(20% 开销、要求任意 20-bit 流中 1/0 数量差 ≤ 2 且不允许连续 5 个以上相同位)大幅提升。128b/130b 用加扰器(scrambler)替代硬性游程限制,只需保证每 66 bit 至少一次跳变。若需要 2¹²⁸ 个码字(取自全部 2¹³⁰ 个 130-bit 模式),每个码字中 0 或 1 的个数必须在 63–67 之间——即码字相当平衡且富含 0-1 跳变。
    • 净数据率(PCIe 3.0,8 GT/s):每 lane 每方向 985 MB/s;链路可由 1/2/4/8/12/16 条 lane 组成(x1、x2、x4、x8、x16),每 lane 1 bit 宽(4 根线,两对差分,上下行同时且对称)。于是:

      链路宽度x1x2x4x8x16
      净带宽(每方向)985 MB/s1.97 GB/s3.94 GB/s7.9 GB/s15.8 GB/s
    • 代际:每代目标翻倍。PCIe 3.0(8 GT/s,极广泛部署)→ PCIe 4.0(16 GT/s,现代 AMD/Intel/IBM 系统,A40 即 PCIe Gen4,x16 双向 64 GB/s,单向约 32 GB/s)→ PCIe 5.0(32 GT/s,目前仅在极少数系统如 IBM Power10 上支持)。
    • NVLink 的数量级:NVLink 是 GPU 之间的直连互连,绕开 PCIe。Ampere 第三代 NVLink 的 GPU-GPU 单向带宽可达 600 GB/s,约为 PCIe Gen4 的近 10 倍;P100 上 8 卡混合立方网格(hybrid cube mesh)双向 160 GB/s,是 PCIe Gen3 x16 的 5 倍;V100 的 NVLink 2.0 为 25 GB/s/lane × 6 lane = 150 GB/s;Delta 上的 A40 为 112.5 GB/s 双向。IBM Power9 系统中 2 颗 Power9(各 10 核 80 线程 @4.02 GHz、256 GB DDR4)挂 4 张 V100(各 16 GB),GPU 间 NVLink 150 GB/s,每 CPU 到 GPU 120 GB/s,两 CPU 之间 X-bus 64 GB/s。
    • DMA 与锁页内存(pinned memory):PCIe 传输用 DMA(Direct Memory Access)才能吃满总线带宽。DMA 使用物理地址,因此需要一个不会被操作系统换出(page out)的缓冲区——否则 OS 可能在 DMA 读写期间把页换走,把别的虚拟页换到同一物理位置,导致数据损坏。这类”不能被换页”的内存叫锁页内存 / 页锁定内存(pinned / page-locked memory)
      • cudaHostAlloc(&ptr, bytes, cudaHostAllocDefault) 分配,用 cudaFreeHost(ptr) 释放;用法与 malloc 完全相同,唯一区别是 OS 不能对它换页。
      • cudaMemcpy() 的主机端缓冲区不是锁页的,运行时必须先在”分页内存 → 一块内部锁页的暂存区(staging buffer)”之间做一次额外的 CPU 拷贝,再启动 DMA。因此 cudaMemcpy 在源或目标为锁页内存时通常快约 2 倍
      • 锁页内存是有限资源,超额申请会有严重后果(可能拖垮整个系统),不要随手给大数组用。
    • 数据搬运的代价(本讲量化核心):以 A100(PCIe Gen4 x16,单向理论 32 GB/s,锁页实测约 25 GB/s)为例,向量加法 N = 2²⁴ 个 float:
      • 上行(A、B 两个数组):2 × 67.11 MB = 134.22 MB ÷ 25 GB/s ≈ 5.37 ms
      • kernel 计算:约 0.14 ms(见后文性能分析)
      • 下行(C 数组):67.11 MB ÷ 25 GB/s ≈ 2.68 ms
      • 传输总计 ≈ 8.05 ms,是 kernel 时间的 57 倍。若用可分页内存(实测约 6–8 GB/s),上行就要 ~19 ms,恶化到 130 倍以上。

      结论:当 kernel 的算术强度很低时,”把数据搬上搬下”就是整个应用的全部成本。工程对策是:让问题规模尽量大以摊薄固定开销;用锁页内存;用多 stream 把传输与计算重叠;更重要的是——尽量让数据一次上传后在 device 上被多个 kernel 反复使用(这正是后续 tiling、融合 kernel 等技术的动机)。

    • 讲义的产物与趋势:GPU 与 CPU 正在融合(同一封装、共享内存控制器),PC 世界正在”变平”(The PC world is becoming flatter),计算外包(outsourcing of computation)越来越容易。
众核吞吐量处理器:SM 与 SIMT 执行模型(Streaming Multiprocessor and SIMT)
  • 定义与目的:SM(Streaming Multiprocessor,流多处理器)是 GPU 内部可独立调度线程块的处理器核心;SIMT(Single Instruction, Multiple Threads,单指令多线程)是 CUDA 的执行模型,它把一个 warp 的 32 个线程组织成”共享同一条指令、但各自拥有独立寄存器和独立执行状态”的执行单元。SIMT 解决了”既要写标量风格的 C 代码,又要让硬件高效地一次驱动 32 个线程”这个核心矛盾。

  • 直观解释(”它是什么?”):把 SM 想象成一个 32 人组成的合唱团声部(warp):指挥(warp scheduler)举起一个手势(一条指令),32 个人同时唱同一个音。每个人有自己的嗓子(寄存器)和自己的乐谱页码(线程索引),所以”同一句歌词”唱出来可以是不同的音高——这就是 SIMT 与 SIMD 的区别:SIMD 是数据打包进向量寄存器,SIMT 是 32 个独立线程被硬件强行同步节奏。如果乐谱上写着”男高音唱这句,其他人闭嘴”,那么这个声部就得唱两遍(一次男高音、一次其他人)——这就是分支发散(divergence)

  • 架构/机制图解:自顶向下,从整卡到单个线程的层次结构:

  ============================ 一块 GPU(例如 A100,108 个 SM) ============================
   +-----------+  +-----------+  +-----------+            +------------+   +-------------+
   |   SM 0    |  |   SM 1    |  |   SM 2    |    ...     |  SM 106    |   |   SM 107    |
   +-----+-----+  +-----+-----+  +-----+-----+            +-----+------+   +------+------+
         |              |              |                        |                 |
         +--------------+--------------+-------- ... ----------+-----------------+
                                        |
                            +-----------+-----------+
                            |  L2 Cache (A100: 40MB) |  <-- 全芯片共享,所有
                            +-----------+-----------+       global 访问的汇聚点
                                        |
                            +-----------+-----------+
                            |  DRAM (HBM2, 1555 GB/s)|
                            +-----------------------+
                                         ^
                                         | PCIe Gen4 x16 / NVLink
                                    [ 主机 CPU + 主机内存 ]

  ============================ 单个 SM 内部(以 A100 sm_80 为例) ============================
   +-------------------------------------------------------------------------------------+
   |                           SM (Streaming Multiprocessor)                             |
   |                                                                                     |
   |  +----------------+  +----------------+  +----------------+  +----------------+     |
   |  | Warp Scheduler |  | Warp Scheduler |  | Warp Scheduler |  | Warp Scheduler |     |
   |  |      0         |  |      1         |  |      2         |  |      3         |     |
   |  | (每周期选 1 个 |  |                |  |                |  |                |     |
   |  |  就绪 warp 发射)|  |                |  |                |  |                |    |
   |  +-------+--------+  +-------+--------+  +-------+--------+  +-------+--------+     |
   |          |                   |                   |                   |              |
   |  +-------v-------------------v-------------------v-------------------v--------+     |
   |  |        寄存器堆 (Register File)  65536 × 32-bit = 256 KB / SM             |      |
   |  |        按 warp 静态划分;32 regs/thread 时正好支撑 2048 线程满占用        |      |
   |  +----------------------------------------------------------------------------+     |
   |                                                                                     |
   |  +-------------------+  +-------------------+  +-------------------+                |
   |  | 16 个 FP32 单元 ×4|  |  16 个 INT32 单元 |  |  LD/ST 单元 (访存) |               |
   |  | = 64 FP32 core/SM |  |                   |  |  发出 global/L1    |               |
   |  +-------------------+  +-------------------+  +---------+---------+                |
   |  +-------------------+  +-------------------+            |                          |
   |  | 4 个 SFU (超越函数)|  | Tensor Core (4 个) |           |                         |
   |  | sin/cos/rsqrt     |  | 4x4x4 / 16x8x16    |           |                          |
   |  +-------------------+  +-------------------+            |                          |
   |                                                           v                         |
   |  +--------------------------------------------------------------------------------+ |
   |  | L1 / Shared Memory  192 KB 统一体可配置(A100 最多 164 KB 作 shared 用)        ||
   |  |   - 共享内存:block 内显式共享,延迟 20-30 cycle,分 32 个 bank               |  |
   |  |   - L1 cache:隐式缓存 global/local,命中约 30 cycle                          |  |
   |  +--------------------------------------------------------------------------------+ |
   +-------------------------------------------------------------------------------------+
        A100 资源上限:2048 threads/SM = 64 warps/SM;32 blocks/SM;
                        65536 regs/SM;164 KB shared/SM;163 KB shared/block(需 opt-in)
  ============================ warp / block / grid 层次 ============================
   Grid(一次 kernel 启动的全部线程;最多 2^31-1 个 block 在 x 维)
    |
    +-- Block(0,0) --> 被调度到【某一个】SM 上,整个 block 生命周期不迁移
    |     |
    |     +-- Warp 0 = Thread(0..31)   <-- 硬件以 32 线程为单位调度
    |     +-- Warp 1 = Thread(32..63)
    |     +-- Warp 2 = Thread(64..95)   <-- 最后一个 warp 可能不满 32 线程,
    |     |                                  此时"空"的 lane 仍然占发射槽(浪费)
    |     +-- 共享内存 / __syncthreads() 只在这个 block 内有效
    |
    +-- Block(1,0) --> 可能在同一 SM,也可能在另一个 SM
    |
    +-- Block(N-1) --> block 之间的执行顺序【完全没有保证】

  索引计算(务必背下来):
     int i = blockIdx.x * blockDim.x + threadIdx.x;      // 一维全局线性索引
     gridDim.x  = 网格 x 维的 block 个数
     blockIdx.x = 本 block 在 x 维的编号(0 .. gridDim.x-1)
     blockDim.x = 本 block 在 x 维的线程数(= <<<..., block>>> 里的 block.x)
     threadIdx.x= 本线程在 block 内的 x 编号(0 .. blockDim.x-1)
  • 关键操作与性能特征
    • 线程是调度的基本单位,warp 是执行的原子单位:block 内的线程按 threadIdx 线性顺序每 32 个切成一个 warpthreadIdx.x = 0..31 是 warp 0,32..63 是 warp 1,以此类推。
    • SIMT 与发散:一个 warp 的 32 个 lane 共享同一个 PC(程序计数器)。遇到分支且条件在不同 lane 上不一致时,硬件串行化执行各分支路径,只让属于该路径的 lane 有效:
   if (tid % 2 == 0) { A; } else { B; }
   32 个 lane 的取值: 0 1 2 3 4 5 6 ... 31
   tid%2==0 的掩码:  1 0 1 0 1 0 1 ... 0

   第 1 遍:执行 A,掩码 = 101010...10  (16 个 lane 有效,16 个闲置)
   第 2 遍:执行 B,掩码 = 010101...01  (16 个 lane 有效,16 个闲置)
   --------------------------------------------------------------
   总时间 = t(A) + t(B),而非 max(t(A), t(B))
   利用率 = 16/32 = 50%  -->  这就是"看别人干活"式的纯浪费
  • block 之间无同步、无顺序保证:这是 CUDA 可扩展性的来源(block 数量可以远超 SM 数量、可以任意顺序调度),也是”跨 block 依赖必须靠多次 kernel 启动(或 2012 年起的动态并行 dynamic parallelism,允许 block 启动子 kernel 并等待其完成)来表达”的原因。
  • 调度与常驻量(用统一硬件档案代入)
GPUSM 数每 SM warp 上限每 SM 线程上限每 SM block 上限一次可常驻总线程
A10010864204832221 184
H100 SXM13264204832270 336
RTX 409012848153624196 608
RTX 2080 Ti683210241669 632
  • 寄存器预算:A100 每 SM 65536 个 32-bit 寄存器。要跑满 2048 线程,平均每线程只能用 65536 / 2048 = 32 个寄存器。若 nvcc 为你的 kernel 分配了 40 个寄存器,则每 SM 最多驻留 65536 / (40×32) = 51.2 → 51 个 warp,占用率降到 51/64 = 79.7%
  • 共享内存预算:A100 每 SM 164 KB 可作 shared 使用。若每个 block 用 48 KB,则每 SM 最多 3 个 block(164/48 = 3.4);即使线程数允许 8 个 256 线程的 block,shared 也把常驻量压到 3 个。
  • 历史题(讲义 Lecture 19 “Problem Solving” 页原文考点):某 GPU 的 compute capability(CC)1.3 的 SM 上限为:1024 threads/SM、8 blocks/SM、16384 registers/SM、16 KB shared memory/SM(即 32 warps/SM)。给定四个 kernel:
kernelthreads/blockregs/threadshared/block受线程限受寄存器限受共享内存限受 block 数限实际 block 数占用率
K1256208 KB416384/(20×256)=3.2→316/8=282512/1024 = 50%
K2128122 KB816384/(12×128)=10.6→1016/2=8881024/1024 = 100%
K3512324 KB216384/(32×512)=116/4=481512/1024 = 50%
K4646416 KB1616384/(64×64)=416/16=18164/1024 = 6.25%
**答案是 K2 达到最大占用率**。这道题把"**占用率 = 四类资源上限取最小值**"这条规则讲透了:`residentBlocks = min( byWarps, byThreads, byRegs, bySharedMem, byBlocksPerSM )`。
延迟隐藏与占用率(Latency Hiding and Occupancy)
  • 定义与目的:延迟隐藏(latency hiding)指当某个 warp 因为等待数据(访存、纹理、分支)而停顿时,硬件立刻切换到另一个就绪 warp 去发射指令,从而让执行单元始终有活干。占用率(occupancy)是”每 SM 实际常驻 warp 数 ÷ 每 SM 最大 warp 数”,它是延迟隐藏能力的上限指标。这套机制解决的是 GPU 最根本的工程矛盾:显存延迟长达 400–800 周期,而 GPU 只有很浅的缓存层次,靠什么不空转?

  • 直观解释(”它是什么?”):想象银行有 64 个窗口(warp 槽位),但每个客户办业务时都要等后台调档案(400–800 周期)。如果只开 1 个窗口,柜员每次都要干等档案送来,一天办不了几笔。如果 64 个窗口全开,柜员在每个窗口都留下”等待中”的客户,谁的材料到了就先办谁——没有人被加速,但整个大厅的吞吐量拉满了。这就是 GPU 的设计哲学:不减少延迟,而是用并发掩盖延迟。用讲义原话:”GPU 需要海量线程来容忍延迟(Require massive number of threads to tolerate latencies)”。

  • 架构/机制图解

  同一个 SM 上 4 个 warp scheduler 的时间线(每周期每个 scheduler 可发射 1 个 warp 的 1 条指令)

  周期:      0    1    2    3    4    5    6    7   ...  400  401  402
  ------------------------------------------------------------------------
  Sched 0:
    warp 0  [LDG].....................................................[用数据]
            ^ 发出全局加载,进入"未就绪"状态
    warp 32          [FADD]  <- 切换!warp 32 的就绪指令立刻被发射
    warp 1               [IMAD]
    warp 33                   [LDG]..................................
    warp 2                          [FADD]
    ...                              ...  (轮流发射,执行单元 100% 忙)

  ------------------------------------------------------------------------
  【关键】执行单元的利用率 = min(1, 可发射的 warp 数 / 需要填补的延迟槽位数)
          常驻 warp 越多 -> 越容易找到就绪 warp -> 延迟隐藏越充分

  【Little 定律:算清楚"要多少 warp 才够"】
      需要的在途字节数 = 内存延迟 × 内存带宽
      A100: 700 cycle / 1.41 GHz = 497 ns
            497 ns × 1555 GB/s = 772 KB 需要同时在途
            772 KB / 108 SM        = 7.15 KB / SM
            7.15 KB / 128 B 每 warp = 约 56 个 warp 的在途 load
            ==> A100 每 SM 最多 64 个 warp,用纯 load 打满带宽几乎需要满占用!

      H100: 700/1.98 GHz = 354 ns × 3350 GB/s = 1.184 MB
            1.184 MB / 132 SM / 128 B = 约 70 个 warp  <-- 超过 64!
            ==> H100 单靠占有率达到 100% 也不够,必须靠 ILP(每线程多条独立 load)

      RTX 4090: 700/2.52 GHz = 278 ns × 1008 GB/s = 280 KB
            280 KB / 128 SM / 128 B = 约 17 个 warp(上限 48)
            ==> 仅约 36% 占用率就足以打满带宽
  ------------------------------------------------------------------------
  延迟层级(写代码时必须记住的数字):
    寄存器            ~0 cycle        每线程私有,最快
    共享内存 / L1 命中 20-30 cycle     block 内共享,需 __syncthreads()
    L2 命中           ~200 cycle       全芯片共享(A100: 40 MB)
    全局内存 (DRAM)   400-800 cycle    "墙";必须靠并发掩盖
    PCIe H2D/D2H      微秒级          比 DRAM 再慢 1-2 个数量级
  • 关键操作与性能特征
    • 占用率的准确定义occupancy = 每 SM 实际常驻 warp 数 / 每 SM 最大 warp 数。A100 上 64/64 = 100%
    • 三级限制链每 SM 常驻 block 数 = min(线程上限, warp 上限, 寄存器上限, 共享内存上限, 硬件 block 上限)每 SM 常驻 warp 数 = 每 SM block 数 × 每 block 的 warp 数
    • 寄存器分配粒度会让理论值和实际值有出入:nvcc 按 warp 为单位、并以 256 个寄存器为粒度向上取整分配(不同 arch 略有差异),所以 65536 / (numRegs × 32) 的手算值可能与 cudaOccupancyMaxActiveBlocksPerMultiprocessor() 的返回值差 1 个 block。以运行时 API 为准,手算用于理解原理。
    • 占用率不是越高越好:占用率高只是”有足够多的 warp 可切换”。若每个线程可用寄存器变少导致寄存器溢出(register spilling,把变量挤到 local memory),性能反而暴跌。常见的工程折中是 50%–75% 占用率,配合每线程 2–4 条独立访存(ILP)。可以用 __launch_bounds__(maxThreadsPerBlock, minBlocksPerMultiprocessor) 告诉编译器”我要几个 block”,让它把寄存器压到预算内。
    • 实测建议:Lab 0 的 Device Query 应打印出”不同 blockDim 下的理论占用率表”,这样在后续 Lab 中调 blockDim 时能立刻看出瓶颈是寄存器、共享内存还是 block 数上限。
内存带宽与算术强度(Memory Bandwidth and Arithmetic Intensity)
  • 定义与目的:算术强度(arithmetic intensity, AI)指程序每搬运 1 字节数据所完成的浮点运算次数,单位 FLOP/Byte。它与机器的”平衡点”(peak FLOPS ÷ peak bandwidth)比较,就能在写代码之前判断程序是内存带宽受限还是计算受限。这是 GPU 编程中投入产出比最高的一次分析。

  • 直观解释(”它是什么?”):把 GPU 想象成一个超级厨师(算力)配一个很小的传菜窗口(带宽)。如果菜谱是”把一块豆腐切成丝”(每克原料要做很多刀 = 高算术强度),那么厨师的手速是瓶颈,窗口再小也无所谓。如果菜谱是”把一箱矿泉水从门口搬到仓库”(每克原料只做 0.08 次运算 = 极低算术强度),那么厨师再快也没用——瓶颈是门口那条通道。GPU 的厨师比 CPU 强 10 倍以上,但它的窗口(1555 GB/s)相对于它的手速(19.5 TFLOPS)其实很窄,所以绝大多数朴素 CUDA 程序都是”搬水工”,不是”切丝师傅”

  • 架构/机制图解

  Roofline 模型(A100: 峰值 FP32 = 19.5 TFLOPS,带宽 = 1555 GB/s)

  性能 (GFLOP/s, 对数)
   10^4 |========================================================= 计算屋顶 19 500 GFLOP/s
        |                                              /
        |                                            /  内存屋顶 (斜线, 斜率 = 带宽)
        |                                          /     每 1 Byte 换来 1555 GFLOP... 实为:
   10^3 |                                        /       y = AI × 1555 GB/s
        |                                      /
        |                                    /   <-- 折点 = 机器平衡点
    100 |                                  /          AI* = 19500/1555 = 12.5 FLOP/Byte
        |                                /
        |                              /
     10 |                            /
        |                          /
      1 |*  <-- 向量加法: AI = 0.083 FLOP/B, 只能跑到约 129 GFLOP/s
        +----+----+----+----+----+----+----+----+----+----+----+----> 算术强度 (FLOP/Byte)
           0.1  0.5   1    2    5   10  12.5  20   50  100  ...
                                ^
                    AI < 12.5 -> 内存带宽受限(Roofline 斜线段)
                    AI > 12.5 -> 计算受限(Roofline 水平段)

  顶点判断公式:  可达性能 = min( peakFLOPS, AI × peakBandwidth )
  机器平衡点:    AI* = peakFLOPS / peakBandwidth
                 A100     : 19500 GFLOP/s ÷ 1555 GB/s = 12.5 FLOP/Byte
                 H100 SXM : 67000 GFLOP/s ÷ 3350 GB/s = 20.0 FLOP/Byte
                 RTX 4090 : 82600 GFLOP/s ÷ 1008 GB/s = 82.0 FLOP/Byte
                 RTX 2080Ti: 13400 GFLOP/s ÷ 616 GB/s = 21.8 FLOP/Byte
  • 关键操作与性能特征
    • 合并访问(coalescing)决定了实际能拿到多少带宽:一个 warp 的 32 个线程若访问连续的 32 个 4 字节字,硬件会把它合成一次 128 字节的事务(transaction)= 一条 cache line;若访问步长是 32(strided),则 32 个线程散落在 32 条不同的 cache line 上,需要 32 次 128 字节事务 = 4096 字节才能满足 128 字节的有用数据,有效带宽掉到 1/32。这是 GPU 性能最常见的”隐形杀手”。
    • 典型 kernel 的算术强度(先算再写)
    kernel每元素访存每元素 FLOPAI (FLOP/B)A100 上判定
    向量加法 C=A+B12 B(读 A、B,写 C)10.083极度带宽受限(差 150 倍)
    归一化 A[i]*=s8 B(读+写)10.125带宽受限
    点积 / reduction4 B(只读)10.25带宽受限
    3×3 卷积 (每输出)36 B(9 次读)170.47带宽受限
    3×3 卷积 + tiling (每输出)约 4.5 B(摊薄后)17~3.8接近平衡点
    稠密矩阵乘 N×N×N每输出 8N B2N0.25朴素版带宽受限
    稠密矩阵乘 + tiling每输出约 8N/√tile B2N随 tile 增大可转为计算受限
    • 由 AI 直接算出预期耗时(A100,N = 2²⁴ 个 float)
      • 访存流量 = 3 × N × 4 B = 3 × 16 777 216 × 4 = 201.33 MB
      • 内存受限下界 = 201.33 MB ÷ 1555 GB/s = 129.5 µs
      • 计算量 = N × 1 FLOP = 16.78 MFLOP;计算下界 = 16.78e6 ÷ 19.5e12 = 0.86 µs
      • 结论:两者相差 150 倍,耗时完全由带宽决定。要达到 129.5 µs 需要 100% 带宽利用率——实践中典型只能跑到 85%–93%,即 140–152 µs。
    • “内存墙”的长期趋势:Amir Gholami 等人在 RiseLab 博客(2021)指出,过去 20 年 GPU/加速器的峰值算力增长远快于内存带宽增长(FLOPs/Byte 的平衡点持续上升,从 GPU 早期约 10 提升到 H100 的 ~20、RTX 4090 的 ~82)。这意味着”内存墙”越来越紧:未来的 CUDA 优化只会越来越偏向”减少字节搬运”,而非”减少指令数”。这条趋势是理解 tiling、kernel fusion、量化、稀疏化等所有现代技术的总纲。
统一着色器架构与 Tesla 架构:让 CUDA 成为可能的硬件基础(Unified Shader Architecture / Tesla)
  • 定义与目的:统一着色器架构(unified shader architecture)指把图形流水线中原本硬件独立、数量固定比例的顶点着色器(vertex shader)与像素/片段着色器(pixel/fragment shader)单元,合并成一组通用的可编程标量处理器(后来被称为 CUDA core),它们可以执行任意着色阶段、也可以执行通用计算。它在 2006 年(NVIDIA G80,Tesla 架构)落地,并在 2007 年通过 CUDA 暴露给程序员——没有统一着色器架构,就没有 CUDA

  • 直观解释(”它是什么?”):旧图形流水线像一条专用装配线:工位 A 只装轮子(顶点着色)、工位 B 只装车门(像素着色),而且必须按 1:3 的固定比例配人。如果今天的车型不需要那么多车门工序,工位 B 的人就只能闲着。统一着色器架构相当于把所有人都培训成多能工:今天需要 5 个人装轮子就派 5 个,明天需要 3 个装车门就派 3 个。更进一步——既然这些多能工什么都能干,那把非图形的通用计算(物理、线性代数、图像处理)交给他们做,不是理所当然吗? 这就是 GPGPU 与 CUDA 的诞生逻辑。

  • 架构/机制图解

  【2004-2005:分离的着色器(fixed-ratio shader)】
   +---------------------+     +---------------------+     +---------------------+
   | Vertex Shader 阵列  |     | Triangle Setup / Rasterizer |  | Pixel Shader 阵列 |
   | (VS 单元, 数量固定) |---->|                             |--->| (PS 单元, 固定) |
   +---------------------+     +---------------------+     +---------------------+
        VS : PS 的比例由硬件焊死(例如 1:3)。
        哪一段忙、哪一段闲,完全无法调配 -> 利用率低、且无法做通用计算。

  【2006:G80 / Tesla 统一着色器架构】
   +--------------------------------------------------------------------------+
   |  统一的标量处理器阵列(Unified Shader / 后来叫 CUDA core)               |
   |  同一个物理单元,既能跑 vertex shader、也能跑 pixel shader、             |
   |  还能跑程序员写的 __global__ kernel —— 三种用途共用同一批硬件。          |
   +--------------------------------------------------------------------------+
        G80 规模: 128 个 SP (Streaming Processor) 分布在 16 个 SM 上(每 SM 8 个)
        G80 频率: 1.35 GHz  |  每 SM 768 线程 (24 warp)  |  16 KB shared / SM
        显存带宽: 86.4 GB/s (GeForce 8800 GTX)

  【从 G80 到今天的谱系(CUDA compute capability)】
   2006 Tesla    G80      CC 1.0/1.1   共享内存首次对通用计算开放
   2008 Tesla    GT200    CC 1.3       <-- 讲义 Lecture 19 "Problem Solving" 用的就是这个
   2010 Fermi    GF100    CC 2.0       L1/L2 可配置、真正的缓存层次、ECC、并发 kernel
   2012 Kepler   GK110    CC 3.5       动态并行、Hyper-Q、更大寄存器堆 (65536/SM)
   2014 Maxwell  GM204    CC 5.2       每 SM 共享内存大幅提升、能效比飞跃
   2016 Pascal   GP100    CC 6.0       HBM2、NVLink、统一内存
   2017 Volta    GV100    CC 7.0       独立线程调度、Tensor Core 第一代
   2018 Turing   TU102    CC 7.5       RT Core、每 SM 64 FP32 lane
   2020 Ampere   GA100    CC 8.0       <-- ECE408 的 sm_80 目标:108 SM、164 KB shared/SM
   2020 Ampere   GA10x    CC 8.6       每 SM 128 FP32 lane(RTX 30 系)
   2022 Hopper   GH100    CC 9.0       132 SM、228 KB shared/SM、TMA、第四代 NVLink
   2022 Ada      AD102    CC 8.9       <-- RTX 4090:128 SM、1008 GB/s、82.6 TFLOPS
   2024 Blackwell GB100   CC 10.0      第五代 Tensor Core、FP4
  • 关键操作与性能特征
    • CUDA 发布的时间与动机:2007 年 NVIDIA 推出 CUDA,让”给这些硬件写通用程序”变得相当容易(讲义 “Meanwhile, in Another Part of the Industry”)。在此之前,要用 GPU 做通用计算必须把算法伪装成图形渲染(”用纹理做数组、用像素着色器做计算”),门槛极高。
    • G80 到 A100 的规模跃迁:SM 数 16 → 108(6.75 倍),每 SM 线程数 768 → 2048(2.7 倍),共享内存每 SM 16 KB → 164 KB(10 倍),显存带宽 86.4 GB/s → 1555 GB/s(18 倍),峰值 FP32 约 0.35 TFLOPS → 19.5 TFLOPS(56 倍)。算力增长快于带宽增长,这正是内存墙变紧的量化证据。
    • CC 1.3 的 SM 上限(历史基线):1024 threads/SM、8 blocks/SM、16384 registers/SM、16 KB shared/SM、warp = 32 线程(即 32 warp/SM)。对比 A100 的 2048 threads/SM、32 blocks/SM、65536 registers/SM、164 KB shared/SM(即 64 warp/SM)——每一维度都放大了 2 到 16 倍,但”占用率 = 各类资源上限取 min”这条规则 20 年未变。
    • 对 CUDA 程序员的实际含义:统一着色器架构把 GPU 变成了一台有很多简单核心、共享一个很宽的显存接口的 SIMT 机器。因此 CUDA 的性能模型始终是三条:(1) 制造足够多的独立工作(占用率高);(2) 让每个 warp 的访存合并且规整;(3) 把复用的数据搬到离计算更近的地方(寄存器 → 共享内存/L1 → L2 → DRAM 的层次)。这一讲之后的所有讲座,本质上都在细化这三条。
CUDA 五步模板、函数限定符与执行配置语法(The Five-Step Template, Qualifiers and Execution Configuration)
  • 定义与目的:五步模板是 CUDA 程序的骨架(分配 device 内存 → 拷贝入 → 启动 kernel → 拷贝回 → 释放),它在后续每一讲里都会重复出现。__global__ / __device__ / __host__ 三个限定符告诉 nvcc”这个函数编译给谁、谁能调用它”,<<<grid, block>>> 则是 CUDA 独有的执行配置(execution configuration)语法,用来描述”这次启动要用多少个 block、每个 block 多少个线程”。掌握这三件事,就能写出第一个可运行的 CUDA 程序。

  • 直观解释(”它是什么?”)
    • 五步模板 = 寄快递:先把仓库腾出来(cudaMalloc)、把货装上车(cudaMemcpy H2D)、工厂加工(kernel)、把成品运回(cudaMemcpy D2H)、把临时仓库退掉(cudaFree)。运费(PCIe)常常比加工费(kernel)贵得多。
    • 三个限定符 = 三种员工证__global__ 是”外派工程师”——由 CPU(host)招聘并派活,但在 GPU(device)上班,必须返回 void__device__ 是”工厂内部员工”——只能由 GPU 上的代码调用;__host__ 是”总部员工”——只能在 CPU 上跑。注意 __两个下划线字符,写成一个下划线会编译报错。
    • <<<grid, block>>> = 派工单grid 说明”开几条产线(多少个 block)”,block 说明”每条产线站几个人(每个 block 多少线程)”。24 小时运转也没关系——硬件会自动把 block 排队塞进 SM,block 数量可以远大于 SM 数量,这正是 CUDA “同一份代码在 10 个 SM 的旧卡和 132 个 SM 的新卡上都能跑满”的可扩展性来源。
  • 架构/机制图解
  ==================== CUDA 五步模板(后续每讲都以此为基础) ====================

   宿主程序 (main, 跑在 CPU)                       设备 (GPU)
   ------------------------------------------------------------------------------
   【准备】 host 侧 malloc + 初始化数据
            h_A, h_B, h_C

   【第 1 步】cudaMalloc(&d_A, bytes)  ------------>  在显存(global memory)分配
             cudaMalloc(&d_B, bytes)                 d_A, d_B, d_C
             cudaMalloc(&d_C, bytes)                 (返回的指针只能给 device 用!)

   【第 2 步】cudaMemcpy(d_A, h_A, bytes,  =========>  H2D(上行)
               cudaMemcpyHostToDevice)               走 PCIe/NVLink,用 DMA + 锁页内存
             cudaMemcpy(d_B, h_B, bytes,
               cudaMemcpyHostToDevice)               <<< 这一步常常是总耗时的大头 >>>

   【第 3 步】dim3 grid(ceil(N/256), 1, 1);   ------>  启动 grid
             dim3 block(256, 1, 1);                 硬件把 block 分发到各 SM
             vecAddKernel<<<grid, block>>>(          每个 SM 上多个 block 并发
                 d_A, d_B, d_C, N);                 block 内线程按 32 个一组切成 warp
             cudaGetLastError();   <-- 立刻查启动错误   ** 启动是异步的 **
             cudaDeviceSynchronize();                等待 device 全部完成

   【第 4 步】cudaMemcpy(h_C, d_C, bytes,  <=========  D2H(下行)
               cudaMemcpyDeviceToHost)

   【第 5 步】cudaFree(d_A); cudaFree(d_B);  -------->  释放显存
             cudaFree(d_C);
             free(h_A); free(h_B); free(h_C);        释放主机内存

   【第 6 步(可选但强烈推荐)】与 host 计算的 golden 结果对比,打印 PASS/FAIL
   ==============================================================================

   +-----------------------+---------------+----------------+------------------------------------------+
   | 限定符                | 执行位置      | 可调用者       | 备注                                     |
   +-----------------------+---------------+----------------+------------------------------------------+
   | __host__              | host (CPU)    | host           | 默认值;可省略不写                       |
   | __global__            | device (GPU)  | host 或 device | 必须返回 void;即 kernel;用 <<<>>> 启动 |
   | __device__            | device (GPU)  | device         | 设备端辅助函数                           |
   | __host__ + __device__ | host + device | 两者皆可       | 同一个函数编译两份                       |
   +-----------------------+---------------+----------------+------------------------------------------+
   补充限定符:
     __shared__   静态共享内存(block 内共享,生命周期 = block)
     __constant__ 常量内存(64 KB,全 device 只读,通过常量缓存广播)
     __restrict__ 告诉编译器指针不重叠,允许更激进的优化
     __launch_bounds__(maxThreads, minBlocks) 约束寄存器用量以控制占用率

  ==================== 执行配置语法 <<<grid, block>>> ============================
      vecAddKernel<<<Dg, Db, Ns, S>>>(d_A, d_B, d_C, N);
                     |   |   |   |
                     |   |   |   +-- S   : cudaStream_t(默认 0,即默认流)
                     |   |   +------ Ns  : 动态共享内存字节数(默认 0)
                     |   +---------- Db  : dim3 block —— 每个 block 的线程数
                     +-------------- Dg  : dim3 grid  —— block 的个数
     二者都可以直接写整数(等价于 dim3(n,1,1)):
         vecAddKernel<<<ceil(N/256.0), 256>>>(d_A, d_B, d_C, N);
     也可以显式检查边界:
         dim3 DimGrid(N/256, 1, 1);
         if (0 != (N % 256)) { DimGrid.x++; }   // 不足一个 block 也要补一个
         dim3 DimBlock(256, 1, 1);
         vecAddKernel<<<DimGrid, DimBlock>>>(d_A, d_B, d_C, N);
     上限(sm_80): blockDim.x ≤ 1024,blockDim.x*blockDim.y*blockDim.z ≤ 1024
                    gridDim.x ≤ 2^31-1,gridDim.y/z ≤ 65535

  ==================== nvcc 的编译流程 =========================================
      .cu 源文件
         |
         v
      [ nvcc ]  --->  分离 host 代码与 device 代码
         |                    |
         |                    +--> device 代码 -> PTX (虚拟 ISA)
         |                    |       |
         |                    |       v
         |                    |   [ ptxas ] -> SASS (真实机器码, 如 sm_80)
         |                    |       或 由驱动中的 JIT 编译器在运行时编译
         |                    |
         +--> host 代码 -----+--> [ 系统 C++ 编译器 / 链接器 ] -> 可执行文件
  • 关键操作与性能特征
    • kernel 启动是异步的(CUDA 1.0 起即如此):kernel<<<grid, block>>>(args) 语句返回时,kernel 很可能还没开始执行。主机若要使用 kernel 的结果,必须显式同步:cudaDeviceSynchronize()(等 device 全部完成)、cudaMemcpy(隐式同步该次拷贝)、或 cudaStreamSynchronize(stream)cudaEventSynchronize(event)
    • 异步是性能的来源:主机可以在 device 干活时继续准备下一批数据,或者启动另一个 stream 上的 kernel 与当前传输重叠。
    • cudaGetLastError() 检查启动是否成功:如果 grid/block 配置非法(例如 blockDim.x = 2048)或参数错误,错误只会记录在”最后一次错误”里,必须手动查询,否则会被后续某个成功的 API 调用悄悄清掉,导致你完全不知道 kernel 根本没有运行。
    • host 指针与 device 指针不可混用cudaMalloc 返回的指针只能用于 device 代码和 CUDA 运行时 API;把它传给 free() 或直接解引用会导致崩溃或不可预知的错误。反之,把 malloc 的指针传给 kernel 会触发非法地址访问。
    • 统一的错误检查宏(本讲所有代码都使用它):
#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)
代码示例与性能分析
示例 1:Lab 0 —— Device Query(查询并推导本机 GPU 的硬件档案)
// 文件: deviceQuery.cu
// 编译: nvcc -O3 -arch=sm_80 deviceQuery.cu -o deviceQuery
// 运行: ./deviceQuery
//
// Lab 0 参考实现:枚举本机 CUDA 设备,打印设备名、计算能力、SM 数、每 SM 最大线程数、
// 共享内存大小、常量内存大小、warp 大小、全局内存大小,并由 memoryClockRate 与
// memoryBusWidth 推导理论峰值显存带宽;再用 cudaOccupancyMaxActiveBlocksPerMultiprocessor
// 打印不同 blockDim 下的理论占用率,最后给出推荐的启动配置。

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

// 统一的 CUDA 错误检查宏:任何 CUDA API 失败都立刻报出文件、行号与错误串
#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)

// 供占用率查询使用的最小 kernel:不做有意义的工作,只为拿到真实的寄存器/共享内存用量
__global__ void occupancyProbeKernel(const float *in, float *out, int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        out[i] = in[i] * 2.0f + 1.0f;
    }
}

// 由 compute capability 推断架构代号(用于打印可读信息)
static const char *archName(int major, int minor)
{
    if (major == 7 && minor == 0) return "Volta (sm_70)";
    if (major == 7 && minor == 5) return "Turing (sm_75)";
    if (major == 8 && minor == 0) return "Ampere GA100 (sm_80)";
    if (major == 8 && minor == 6) return "Ampere GA10x (sm_86)";
    if (major == 8 && minor == 9) return "Ada Lovelace (sm_89)";
    if (major == 9 && minor == 0) return "Hopper (sm_90)";
    if (major >= 10)              return "Blackwell or newer";
    return "unknown";
}

// 每个 SM 的 FP32 FMA lane 数。NVIDIA 没有提供查询接口,只能按架构查表。
// 这是"估算峰值 FLOPS"的唯一来源,因此下面的 TFLOPS 数字是估算值。
static int fp32LanesPerSm(int major, int minor)
{
    if (major == 7) return 64;                  // Volta / Turing
    if (major == 8 && minor == 0) return 64;    // Ampere GA100 (A100)
    if (major == 8) return 128;                 // Ampere GA10x / Ada
    if (major == 9) return 128;                 // Hopper
    if (major >= 10) return 128;                // Blackwell
    return 64;                                  // 保守默认
}

// 打印"占用率相关参数"表:同时给出运行时 API 结果与手工 min() 估算结果
static void printOccupancyTable(const cudaDeviceProp &p)
{
    cudaFuncAttributes attr{};
    CUDA_CHECK(cudaFuncGetAttributes(&attr, occupancyProbeKernel));

    printf("\n  [占用率相关参数] probe kernel: %d regs/thread, %zu B static smem\n",
           attr.numRegs, attr.sharedSizeBytes);
    printf("    %-10s %-12s %-12s %-11s %-11s %-11s %-11s\n",
           "blockDim", "blocks(API)", "threads/SM", "occupancy",
           "byThreads", "byRegs", "byBlocks");

    const int blockSizes[] = {64, 128, 256, 512, 1024};
    for (int k = 0; k < 5; ++k) {
        int bs = blockSizes[k];
        if (bs > p.maxThreadsPerBlock) {
            continue;
        }

        // 运行时权威答案:把 regs / smem / warp / block 四类上限一起取 min
        int blocksApi = 0;
        CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
            &blocksApi, occupancyProbeKernel, bs, 0));

        // 手工估算:寄存器按 warp 为单位、256 个寄存器为粒度向上取整
        int warpsPerBlock = (bs + p.warpSize - 1) / p.warpSize;
        int regsPerWarp   = attr.numRegs * p.warpSize;
        int regsRounded   = ((regsPerWarp + 255) / 256) * 256;
        int byRegs        = p.regsPerMultiprocessor / (regsRounded * warpsPerBlock);
        int byThreads     = p.maxThreadsPerMultiProcessor / bs;
        int byBlocks      = p.maxBlocksPerMultiProcessor;

        double occ = 100.0 * (double)(blocksApi * bs) /
                     (double)p.maxThreadsPerMultiProcessor;

        printf("    %-10d %-12d %-12d %-10.1f%% %-11d %-11d %-11d\n",
               bs, blocksApi, blocksApi * bs, occ, byThreads, byRegs, byBlocks);
    }
    printf("    说明: blocks(API) 为运行时按四类上限取 min 的结果(权威值);\n");
    printf("          byRegs 为按 256 寄存器分配粒度手工估算的 block 数上限,\n");
    printf("          二者相差 1 属于正常现象,以 API 值为准。\n");
}

int main(int argc, char **argv)
{
    (void)argc;
    (void)argv;

    int deviceCount = 0;
    CUDA_CHECK(cudaGetDeviceCount(&deviceCount));
    printf("Detected %d CUDA capable device(s)\n", deviceCount);
    if (deviceCount == 0) {
        printf("There is no device supporting CUDA.\n");
        return EXIT_SUCCESS;
    }

    for (int dev = 0; dev < deviceCount; ++dev) {
        CUDA_CHECK(cudaSetDevice(dev));

        cudaDeviceProp p{};
        CUDA_CHECK(cudaGetDeviceProperties(&p, dev));

        printf("\n================ Device %d: %s ================\n", dev, p.name);
        printf("  Compute capability (计算能力) : %d.%d   [%s]\n",
               p.major, p.minor, archName(p.major, p.minor));
        printf("  multiProcessorCount (SM 数量) : %d\n", p.multiProcessorCount);
        printf("  warpSize (一个 warp 的线程数) : %d\n", p.warpSize);
        printf("  maxThreadsPerBlock            : %d\n", p.maxThreadsPerBlock);
        printf("  maxThreadsPerMultiProcessor   : %d\n", p.maxThreadsPerMultiProcessor);
        printf("  maxBlocksPerMultiProcessor    : %d\n", p.maxBlocksPerMultiProcessor);
        printf("  regsPerMultiprocessor         : %d  (32-bit)\n", p.regsPerMultiprocessor);
        printf("  regsPerBlock                  : %d\n", p.regsPerBlock);

        printf("\n  --- 存储层次容量 ---\n");
        printf("  sharedMemPerBlock             : %8zu B  (%7.1f KB)\n",
               p.sharedMemPerBlock, p.sharedMemPerBlock / 1024.0);
        printf("  sharedMemPerBlockOptin        : %8zu B  (%7.1f KB)\n",
               p.sharedMemPerBlockOptin, p.sharedMemPerBlockOptin / 1024.0);
        printf("  sharedMemPerMultiprocessor    : %8zu B  (%7.1f KB)\n",
               p.sharedMemPerMultiprocessor, p.sharedMemPerMultiprocessor / 1024.0);
        printf("  totalConstMem (常量内存)      : %8zu B  (%7.1f KB)\n",
               p.totalConstMem, p.totalConstMem / 1024.0);
        printf("  l2CacheSize                   : %8d B  (%7.1f MB)\n",
               p.l2CacheSize, p.l2CacheSize / 1048576.0);
        printf("  totalGlobalMem (全局内存)     : %8.2f GiB\n",
               p.totalGlobalMem / 1073741824.0);

        size_t freeB = 0, totalB = 0;
        CUDA_CHECK(cudaMemGetInfo(&freeB, &totalB));
        printf("  cudaMemGetInfo                : free %.2f GiB / total %.2f GiB\n",
               freeB / 1073741824.0, totalB / 1073741824.0);

        printf("\n  --- 显存带宽 ---\n");
        printf("  memoryClockRate               : %d kHz (%.3f GHz)\n",
               p.memoryClockRate, p.memoryClockRate / 1.0e6);
        printf("  memoryBusWidth                : %d bit\n", p.memoryBusWidth);

        double memClockHz = (double)p.memoryClockRate * 1.0e3;   // kHz -> Hz
        double busBytes   = (double)p.memoryBusWidth / 8.0;      // bit -> Byte
        double bwSimplified = memClockHz * busBytes / 1.0e9;     // 讲义简化公式
        double bwPeak       = 2.0 * memClockHz * busBytes / 1.0e9;  // DDR 修正:每时钟 2 次传输

        printf("  理论带宽 (clock * busWidth / 8)      : %.1f GB/s   <-- 少算一半\n",
               bwSimplified);
        printf("  理论峰值带宽 (clock * busWidth / 8 * 2): %.1f GB/s  <-- GDDR/HBM 为 DDR\n",
               bwPeak);

        printf("\n  --- 峰值算力与机器平衡点(估算)---\n");
        // 注意: clockRate / memoryClockRate 在 CUDA 13 中已弃用,
        //       新代码可改用 cudaDeviceGetAttribute(&clockKHz, cudaDevAttrClockRate, dev)
        int clockKHz = 0;
        CUDA_CHECK(cudaDeviceGetAttribute(&clockKHz, cudaDevAttrClockRate, dev));
        double clockHz = (double)clockKHz * 1.0e3;
        long long lanes = (long long)p.multiProcessorCount *
                          fp32LanesPerSm(p.major, p.minor);
        double peakFlops = 2.0 * (double)lanes * clockHz;   // FMA 记 2 个 FLOP

        printf("  clockRate (查询值,非 boost)  : %.3f GHz\n", clockKHz / 1.0e6);
        printf("  FP32 FMA lane 总数 (查表估算) : %lld\n", lanes);
        printf("  峰值 FP32 (2 * lanes * clock) : %.2f TFLOPS (估算)\n",
               peakFlops / 1.0e12);
        printf("  机器平衡点 (FLOPS / Byte)     : %.2f FLOP/Byte\n",
               peakFlops / (bwPeak * 1.0e9));
        printf("     -> 算术强度低于该值的 kernel 一定受显存带宽限制\n");

        printf("\n  --- 推导量 ---\n");
        int maxWarpsPerSm = p.maxThreadsPerMultiProcessor / p.warpSize;
        printf("  每 SM 最大 warp 数            : %d\n", maxWarpsPerSm);
        printf("  满占用时每线程寄存器预算      : %d regs\n",
               p.regsPerMultiprocessor / p.maxThreadsPerMultiProcessor);
        printf("  满占用时每 block 共享内存预算 : %.1f KB (假设 2 个 block/SM)\n",
               p.sharedMemPerMultiprocessor / 1024.0 / 2.0);

        printOccupancyTable(p);

        printf("\n  --- 推荐的启动配置 ---\n");
        int bs = 256;
        int blocksPerSm = 0;
        CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
            &blocksPerSm, occupancyProbeKernel, bs, 0));
        int gridFull = p.multiProcessorCount * blocksPerSm;
        printf("  blockDim = %d 时每 SM 可驻留 %d 个 block\n", bs, blocksPerSm);
        printf("  满占用 grid = %d 个 block = %d 个线程 = %.2f M 线程\n",
               gridFull, gridFull * bs, gridFull * (double)bs / 1.0e6);
        printf("  grid-stride loop 建议用这个 grid 大小(每个线程循环多次)\n");

        printf("\n  --- 编译目标提示 ---\n");
        printf("  本程序以 -arch=sm_80 编译,设备为 sm_%d%d。\n", p.major, p.minor);
        if (p.major == 8 && p.minor == 0) {
            printf("  完全匹配:可直接运行 sm_80 的 SASS。\n");
        } else if (p.major == 8 && p.minor > 0) {
            printf("  sm_80 的 PTX 会被驱动 JIT 编译成本机 SASS(首次启动略慢)。\n");
        } else {
            printf("  不匹配:请改用 -arch=sm_%d%d 或 -arch=native 重新编译。\n",
                   p.major, p.minor);
        }
    }

    CUDA_CHECK(cudaDeviceReset());
    return EXIT_SUCCESS;
}
  • 【代码做什么?】
    1. cudaGetDeviceCount(&deviceCount) 询问系统里有几张 CUDA 设备;若为 0 则打印提示并正常退出(这一步在 Delta 上如果忘记 srun/salloc 申请 GPU 资源时会返回 0,是最常见的环境问题)。
    2. 对每张设备先 cudaSetDevice(dev)cudaMemGetInfo 等 API 作用于”当前设备”),再 cudaGetDeviceProperties(&p, dev) 一次性取回 cudaDeviceProp 结构体——它包含本讲要求打印的全部字段。
    3. 按四组打印:身份信息namemajor.minor、架构代号)、并行资源multiProcessorCountwarpSizemaxThreadsPerBlockmaxThreadsPerMultiProcessormaxBlocksPerMultiProcessorregsPerMultiprocessor)、存储层次sharedMemPerBlocksharedMemPerBlockOptinsharedMemPerMultiprocessortotalConstMeml2CacheSizetotalGlobalMem)、带宽与算力
    4. 带宽推导memoryClockRate 的单位是 kHz,先换算成 Hz;memoryBusWidth 单位是 bit,除以 8 得到”每时钟传输的字节数”;再乘以 2(GDDR/HBM 是 DDR,每个时钟周期有上升沿和下降沿两次数据传输)。A100 代入:2 × 1.215e9 Hz × (5120/8) B = 2 × 1.215e9 × 640 = 1.5552e12 B/s = 1555.2 GB/s,与硬件档案一致。只写 clock × bus / 8 会得到 777.6 GB/s,正好少一半。
    5. 峰值算力推导峰值 FP32 = 2(FMA)× lane 总数 × clock。A100 代入:2 × 108 × 64 × 1.41e9 = 19.49e12 = 19.5 TFLOPSfp32LanesPerSm() 是按架构查表的(NVIDIA 未提供 lane 数查询接口),而且查到的 clockRate 通常不是 boost 频率,所以这是估算值——但在教学和优化方向判断上完全够用。
    6. 机器平衡点peakFlops / peakBandwidth = 19.5e12 / 1.5552e12 = 12.54 FLOP/Byte。这个数字是 Roofline 模型的折点,任何 AI 低于它的 kernel 都是带宽受限。
    7. 占用率表:对 blockDim = 64/128/256/512/1024 分别调用 cudaOccupancyMaxActiveBlocksPerMultiprocessor(),得到运行时权威的每 SM 常驻 block 数;同时用 byThreads / byRegs / byBlocks 手算一遍,让学生看到”取 min”这条规则如何工作,以及寄存器分配粒度带来的 ±1 差异。
    8. 推荐配置:用 blockDim = 256 算出的 blocksPerSm(A100 上 probe kernel 约 8)乘以 SM 数(108),得到”满占用 grid = 864 个 block = 221 184 个线程”。这个数字就是后续所有 grid-stride kernel 的 grid 大小。
    9. 最后用 major/minor 判断所编译的 sm_80 二进制能否在本机直接运行,还是需要驱动 JIT 重编译。
  • 【并行机制与硬件映射解说】
    • 本程序没有并行 kernel 负载occupancyProbeKernel 只被用来查询资源占用cudaFuncGetAttributes 拿到 numRegs),从未真正启动。真正干活的是主机端的 cudaGetDevicePropertiescudaOccupancyMaxActiveBlocksPerMultiprocessor,它们都是同步的运行时调用
    • block/warp 划分在 Lab 0 中是被”查询”而非被”执行”的对象warpSize = 32 是硬件常量;threadIdx.x / warpSize 就是线程所属的 warp 编号(前提是 block 只用 x 维)。若 blockDim.x = 100,那么一个 block 有 ceil(100/32) = 4 个 warp,最后一个 warp 只有 4 个活跃 lane,其余 28 个 lane 仍占发射槽——block 大小取 32 的整数倍(128、256、512)能避免这种浪费
    • A100 上的常驻量(代入)maxThreadsPerMultiProcessor = 2048warpSize = 322048/32 = 64 个 warp/SM。probe kernel 约 10 个寄存器,byRegs = 65536/((10×32 向上取整到 256)×8 warps/block) = 65536/(512×8) = 16byThreads = 2048/256 = 8byBlocks = 32maxBlocksPerMultiProcessor);三者取 min → 8 个 block/SM = 2048 线程 = 100% 占用率
    • 共享内存无关:probe kernel 的 sharedSizeBytes = 0,所以共享内存不构成限制。若 kernel 使用 48 KB 静态共享内存,则 bySmem = 164/48 = 3,占用率掉到 3×256/2048 = 37.5%
    • 不存在 warp 发散:probe kernel 里的 if (i < n) 在所有线程上取值一致(在 grid 覆盖范围内),因此无发散。仅当 n 不是 blockDim 的整数倍时,最后一个 block 的最后一个 warp 才会发生一次发散。
    • 不存在全局内存合并不良:probe kernel 的访问模式是 in[i],一个 warp 的 32 个线程读连续的 32 个 4 字节字 = 128 字节 = 恰好 1 条 cache line。如果改成 in[i * 2](步长 2),一个 warp 就会横跨 256 字节 = 2 条 cache line,有效带宽打对折。
    • 不存在 bank conflict:完全没有共享内存访问。
    • 不同架构的常驻量对比(同一份 probe kernel、blockDim = 256):
    GPUmaxBlocksPerMultiProcessorbyThreadsbyRegs(以 10 regs 计)实际 blocks/SM占用率
    A100 (sm_80)322048/256 = 865536/(512×8) = 168100%
    H100 (sm_90)322048/256 = 8168100%
    RTX 4090 (sm_89)241536/256 = 665536/(512×6)=216100%
    RTX 2080 Ti (sm_75)161024/256 = 4324100%
  • 【性能优化分析】
    • 占用率(occupancy):probe kernel 在 blockDim = 256 时达到 100%(8 blocks × 256 threads = 2048 = maxThreadsPerMultiProcessor)。这意味着该 kernel 在寄存器维度上没有任何压力(用了 10 个寄存器,预算 32 个)。
      • 计算式:每线程寄存器预算 = regsPerMultiprocessor / maxThreadsPerMultiProcessor = 65536 / 2048 = 32
      • 实测占用率:blocksApi × blockDim / maxThreadsPerMultiProcessor = 8 × 256 / 2048 = 100%
      • 最坏情况:若某 kernel 用了 64 个寄存器,则 byRegs = 65536/((64×32 向上取整到 2048)×8) = 65536/16384 = 4,占用率降到 4×256/2048 = 50%,延迟隐藏能力减半。
    • 算术强度(arithmetic intensity)与该程序的关系:Device Query 本身几乎不搬数据(totalGlobalMemcudaMemGetInfo 都是元数据查询),AI ≈ 0,是纯粹的延迟受限/API 调用受限程序,耗时以微秒计的固定开销为主(每次 CUDA 运行时调用约 1–5 µs)。它的价值不在性能,而在为后续 kernel 提供定标常数
    • Roofline 定标:把本程序打印出的两个数写进笔记本——
      • 峰值带宽 BW = 1555 GB/s(A100)
      • 机器平衡点 AI* = 19.5 TFLOPS / 1555 GB/s = 12.5 FLOP/Byte 有了这两个数,后续任何一个 kernel 都能在 30 秒内判断瓶颈: 预期耗时 ≈ max( 访存字节数 / BW, FLOP 数 / peakFLOPS )
    • 可直接执行的优化方向
      1. 记住本机的”满占用 grid / 满占用线程数”(A100 + 256 线程 = 864 blocks = 221 184 线程),后续 grid-stride kernel 一律用它。
      2. maxThreadsPerBlocksharedMemPerBlockOptintotalConstMem 写进笔记本——它们分别是”每 block 线程数的硬上限”、”每 block 共享内存的硬上限(申请超过 48 KB 必须显式 opt-in)”、”常量内存容量(65 536 B)”,是后面调参时的护栏。
      3. 把理论带宽与实际测得的带宽做对比。Device Query 打印的是理论值 1555 GB/s;真实 kernel 通常只能跑到 85%–93%(1320–1445 GB/s)。用实测值而不是理论值做 Roofline 分析,预测会准确得多。
示例 2:Hello, GPU —— 最小 kernel 与五步模板骨架
// 文件: hello_gpu.cu
// 编译: nvcc -O3 -arch=sm_80 hello_gpu.cu -o hello_gpu
// 运行: ./hello_gpu           (默认 4 个 block x 8 个线程)
//       ./hello_gpu 8 4       (8 个 block x 4 个线程 -> 共 2 个 warp,便于观察 warp 划分)
//       ./hello_gpu 4 32      (4 个 block x 32 个线程 -> 每个 block 恰好 1 个 warp)
//
// 目的:用最少的代码完整展示 CUDA 程序的五步结构,并让每个线程打印自己的身份,
//       从而看清 grid / block / thread 三层索引与 warp 的对应关系。

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

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

// __global__ = kernel:由 host 启动、在 device 上执行、必须返回 void
__global__ void helloKernel(const float *in, float *out, int n)
{
    unsigned int gtid = blockIdx.x * blockDim.x + threadIdx.x;  // 全局线性索引

    // device 端 printf(需要 compute capability >= 2.0)。注意:
    //   - 输出顺序不确定(warp 调度顺序不保证)
    //   - 输出缓冲在 kernel 结束(或缓冲区满)时才刷出
    printf("  [device] blk %u/%u  thr %2u/%-2u  warp %u lane %2u  gtid %2u  %s\n",
           blockIdx.x, gridDim.x,
           threadIdx.x, blockDim.x,
           threadIdx.x / warpSize, threadIdx.x % warpSize,
           gtid,
           (gtid < (unsigned)n) ? "work" : "idle");

    if (gtid < (unsigned)n) {
        out[gtid] = in[gtid] + 1.0f;   // 让每个有效线程真的访问一次全局内存
    }
}

int main(int argc, char **argv)
{
    int nBlocks  = (argc > 1) ? atoi(argv[1]) : 4;
    int nThreads = (argc > 2) ? atoi(argv[2]) : 8;

    if (nBlocks <= 0 || nThreads <= 0 || nThreads > 1024) {
        printf("用法: %s [nBlocks(>0)] [nThreads(1..1024)]\n", argv[0]);
        return EXIT_FAILURE;
    }

    int n = nBlocks * nThreads;
    size_t bytes = (size_t)n * sizeof(float);

    printf("Hello, GPU: %d blocks x %d threads = %d threads = %d warp(s)\n",
           nBlocks, nThreads, n, (n + 31) / 32);
    printf("索引公式: i = blockIdx.x * blockDim.x + threadIdx.x\n\n");

    // ---------- 主机端准备 ----------
    float *h_in  = (float *)malloc(bytes);
    float *h_out = (float *)malloc(bytes);
    for (int i = 0; i < n; ++i) {
        h_in[i]  = (float)i;
        h_out[i] = -1.0f;
    }

    // ---------- 第 1 步:在 device 上分配全局内存 ----------
    float *d_in  = nullptr;
    float *d_out = nullptr;
    CUDA_CHECK(cudaMalloc((void **)&d_in,  bytes));
    CUDA_CHECK(cudaMalloc((void **)&d_out, bytes));

    // ---------- 第 2 步:把输入从 host 拷到 device (H2D) ----------
    CUDA_CHECK(cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice));

    // ---------- 第 3 步:启动 kernel ----------
    dim3 grid(nBlocks, 1, 1);
    dim3 block(nThreads, 1, 1);
    helloKernel<<<grid, block>>>(d_in, d_out, n);

    // kernel 启动是异步的:这一行返回时 kernel 可能还没开始跑。
    // cudaGetLastError() 用来捕获"启动配置非法"这类即时错误。
    CUDA_CHECK(cudaGetLastError());
    // 显式同步:等 device 上所有先前的工作(含 printf 缓冲)完成。
    CUDA_CHECK(cudaDeviceSynchronize());

    // ---------- 第 4 步:把结果从 device 拷回 host (D2H) ----------
    CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));

    // ---------- 第 5 步:释放 device 内存 ----------
    CUDA_CHECK(cudaFree(d_in));
    CUDA_CHECK(cudaFree(d_out));

    // ---------- 主机端自检 ----------
    int bad = 0;
    for (int i = 0; i < n; ++i) {
        if (h_out[i] != h_in[i] + 1.0f) {
            if (bad == 0) {
                printf("首个不匹配: i=%d got %.1f expected %.1f\n",
                       i, h_out[i], h_in[i] + 1.0f);
            }
            ++bad;
        }
    }
    printf("\nhost check: %s (%d/%d 元素正确)\n",
           (bad == 0) ? "PASS" : "FAIL", n - bad, n);

    free(h_in);
    free(h_out);
    return (bad == 0) ? EXIT_SUCCESS : EXIT_FAILURE;
}
  • 【代码做什么?】
    1. 从命令行读 nBlocksnThreads,据此算出总线程数 n 与字节数 bytes。默认 4 × 8 = 32 个线程,正好等于 1 个 warp。
    2. 主机端 malloc 两个数组:h_in[i] = ih_out[i] = -1.0f(哨兵值,用于确认每个元素都被真的写过)。
    3. 第 1 步:两次 cudaMalloc 在显存里开两块 bytes 大小的空间,返回设备指针 d_ind_out。注意 (void **) 强制转换是必需的——cudaMalloc 需要修改指针本身的值,所以传的是指针的地址。
    4. 第 2 步cudaMemcpy(d_in, h_in, bytes, cudaMemcpyHostToDevice) 把输入上行。方向的枚举值写反(写成 DeviceToHost)是经典错误,会报 invalid argument
    5. 第 3 步:构造 dim3 grid(nBlocks,1,1)dim3 block(nThreads,1,1),用 helloKernel<<<grid, block>>>(d_in, d_out, n) 启动。kernel 里:gtid = blockIdx.x * blockDim.x + threadIdx.x 得到全局线性索引(0 到 31),printf 打印 block 编号、block 内线程编号、warp 编号(threadIdx.x / warpSize)、lane 编号(threadIdx.x % warpSize);只有 gtid < n 的线程才写 out[gtid] = in[gtid] + 1.0f
    6. 启动后立刻 cudaGetLastError(),再 cudaDeviceSynchronize()(既等待完成,也确保 device 端的 printf 缓冲被刷出到 stdout)。
    7. 第 4 步cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost) 把结果下行。
    8. 第 5 步:两次 cudaFree 释放显存,最后 free 释放主机内存。
    9. 自检:逐元素对比 h_out[i]h_in[i] + 1.0f,打印 PASS/FAIL 与正确元素个数。若 kernel 的 if (gtid < n) 判断写错,h_out 会保留 -1.0 哨兵值,立刻暴露。
  • 【并行机制与硬件映射解说】
    • warp 划分(用默认 4 block × 8 线程观察):每个 block 只有 8 个线程,8 / 32 = 0.25,因此每个 block 只有 1 个 warp,且该 warp 只有 8 个活跃 lane,其余 24 个 lane 是空的
      • 硬件仍然按 32 个 lane 分配资源和发射槽位:占用率 = 8/32 = 25%(相对 warp 槽位)
      • 这是”小 block 浪费发射槽”的直接演示。用 ./hello_gpu 4 32 运行时,每个 block 恰好 1 个满 warp,效率最高(但 block 数太少,仍无法填满 A100 的 108 个 SM)。
    • 在 A100 上执行 ./hello_gpu 4 8 的真实映射
      • grid 有 4 个 block,A100 有 108 个 SM → 只有 4 个 SM 各拿到 1 个 block,另外 104 个 SM 完全空闲。GPU 的算力利用率 ≈ 4/108 = 3.7%。
      • 每个有活的 SM 上:1 个 warp 被分配到 4 个 warp scheduler 中的一个,该 scheduler 每个周期检查它是否就绪;由于这条 warp 的绝大多数时间在等 printf 的 I/O 和 out[gtid] 的写回,没有任何其他 warp 可以切换——这是延迟隐藏完全失效的极端例子。
    • warp 发散分析:kernel 中的 if (gtid < n) 在默认参数下(n = 32,全部线程都满足)不发散。若运行 ./hello_gpu 4 8 并故意把 n 设成不是 32 的倍数,只有最后一个 warp 会发生一次两路发散,代价约为该 warp 分支路径的额外一遍。
    • 全局内存访问是否合并in[gtid]out[gtid] 都是连续访问。一个满 warp(32 线程)读 in[0..31] = 128 字节 = 恰好 1 条 cache line = 1 次内存事务;写 out[0..31] 同理。完全合并(perfectly coalesced)。相比之下,若写成 in[gtid * 4],一个 warp 会跨越 32 × 16 = 512 字节 = 4 条 cache line,需要 4 次事务换取 128 字节有用数据,有效带宽降到 25%
    • 寄存器使用量与占用率:这个 kernel 大约用 8–12 个寄存器。A100 的预算是每线程 32 个(在 2048 线程满占用下),所以寄存器完全不构成限制。真正的限制是根本没有足够的线程——只有 32 个线程,而 A100 一次可以常驻 108 × 2048 = 221 184 个线程。这就是”小问题在 GPU 上跑不满“的原因:不是 GPU 慢,而是它的并行度没被喂饱。
    • block 调度:4 个 block 会被分配到 4 个不同的 SM(CUDA 运行时倾向于把 block 尽量分散到空闲 SM)。block 之间没有任何同步机制,执行顺序也不保证——printf 输出的顺序在不同次运行中可能不同,这是非确定性(non-determinism)的第一课(讲义 “Determinism and Reproducibility”)。
  • 【性能优化分析】
    • 占用率(occupancy)
      • warp 粒度占用率:每个 block 1 个 warp / 每 SM 最多 64 个 warp。4 个 block 分散到 4 个 SM,每 SM 1 个 warp → 占用率 = 1/64 = 1.6%
      • 线程粒度占用率:8 threads / 2048 threads = 0.39%
      • 结论:这是一个”延迟完全暴露”的程序,其耗时几乎全部是”等 DRAM 与 printf”。
    • 算术强度(arithmetic intensity)
      • 每个元素:读 4 B(in)+ 写 4 B(out)= 8 B;运算 1 次加法 = 1 FLOP。
      • AI = 1 FLOP / 8 B = 0.125 FLOP/Byte
      • A100 平衡点 AI* = 12.5 FLOP/ByteAI / AI* = 0.125/12.5 = 1%带宽受限,且受限于”并行度不足”而非”带宽不足”
    • Roofline 定量分析(代入 n = 32)
      • 访存字节数 = 32 × 8 B = 256 B;访存下界 = 256 B / 1555 GB/s = 0.165 ns(即 0.000000165 ms)。
      • 计算量 = 32 FLOP;计算下界 = 32 / 19.5e12 = 0.0016 ns
      • 两个下界都远小于任何可测量的时间尺度。实际耗时由三部分主导:三次 CUDA API 调用(cudaMalloc ×2 各 5–20 µs)、一次 kernel 启动(3–10 µs)、两次主机-设备往返(微秒级),总计 20–60 µs,是理论下界的 10⁸ 倍以上
      • 这个对比本身就是最重要的教学点:当问题规模太小时,性能由固定开销(launch overhead、API 调用)决定,与 Roofline 无关。后续所有 Lab 的问题规模都被刻意设计得足够大(N = 10⁶–10⁷ 量级),就是为了让 Roofline 分析成立。
    • 可执行的优化方向
      1. blockDim.x 设为 32 的整数倍(推荐 128 或 256),避免空 lane 浪费。
      2. 让 grid 至少覆盖 SM 数 × 每 SM 常驻 block 数(A100 上 864 个 block),否则部分 SM 闲置。
      3. 一次 cudaMalloc 分配一大块内存再手动切分,减少 API 调用次数(每次 cudaMalloc 都是同步操作,代价 5–20 µs)。
      4. 把多次小 cudaMemcpy 合并成一次大的(PCIe 每次传输都有固定启动开销)。
      5. 在真实项目中不要用 device 端 printf 调试大循环——它会串行化输出并严重拖慢 kernel。
示例 3:CUDA 向量加法 —— 全书基准模板(朴素版 vs grid-stride 版)
// 文件: vecadd.cu
// 编译: nvcc -O3 -arch=sm_80 vecadd.cu -o vecadd
// 运行: ./vecadd            (默认 N = 2^24 = 16777216 个 float)
//       ./vecadd 1000003    (故意取非 256 倍数,验证边界处理)
//
// 这是 ECE408 全书使用的"基准模板":完整演示 CUDA 五步结构、错误检查、
// cudaEvent 计时、H2D/kernel/D2H 三段耗时拆分、CPU 参考实现比对、
// 以及"一线程一元素"与"grid-stride loop"两种线程映射的对比。

#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 BLOCK_SIZE 256

// ---------------- 版本 1:一线程一元素(最朴素) ----------------
__global__ void vecAddKernel(const float *A, const float *B, float *C, int n)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {                       // 边界检查:grid 向上取整后必然有越界线程
        C[i] = A[i] + B[i];
    }
}

// ---------------- 版本 2:grid-stride loop(线程数固定,线程处理多个元素) ----------------
__global__ void vecAddKernelGridStride(const float *A, const float *B, float *C, int n)
{
    int stride = gridDim.x * blockDim.x;            // 整个 grid 的线程总数
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += stride) {
        C[i] = A[i] + B[i];
    }
}

static void initData(float *a, float *b, int n)
{
    for (int i = 0; i < n; ++i) {
        a[i] = 1.0f * (float)i;
        b[i] = 2.0f * (float)i;
    }
}

static bool verify(const float *a, const float *b, const float *c, int n)
{
    for (int i = 0; i < n; ++i) {
        float ref = a[i] + b[i];
        float tol = 1e-3f * (1.0f + fabsf(ref));
        if (fabsf(c[i] - ref) > tol) {
            printf("  首个不匹配: i=%d  got %.6f  expected %.6f\n", i, c[i], ref);
            return false;
        }
    }
    return true;
}

static float elapsedMs(cudaEvent_t start, cudaEvent_t stop)
{
    float ms = 0.0f;
    CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
    return ms;
}

int main(int argc, char **argv)
{
    int n = (argc > 1) ? atoi(argv[1]) : (1 << 24);
    if (n <= 0) {
        printf("用法: %s [N>0]\n", argv[0]);
        return EXIT_FAILURE;
    }
    size_t bytes = (size_t)n * sizeof(float);

    printf("N = %d floats\n", n);
    printf("  每数组        = %.2f MiB\n", bytes / 1048576.0);
    printf("  理论访存流量  = %.2f MiB (读 A + 读 B + 写 C, 共 3N*4 字节)\n",
           3.0 * bytes / 1048576.0);

    // ---------- host 准备 ----------
    float *h_A  = (float *)malloc(bytes);
    float *h_B  = (float *)malloc(bytes);
    float *h_C  = (float *)malloc(bytes);   // 版本 1 的结果
    float *h_C2 = (float *)malloc(bytes);   // 版本 2 的结果
    initData(h_A, h_B, n);

    // ---------- 第 1 步:分配 device 内存 ----------
    float *d_A = nullptr;
    float *d_B = nullptr;
    float *d_C = nullptr;
    CUDA_CHECK(cudaMalloc((void **)&d_A, bytes));
    CUDA_CHECK(cudaMalloc((void **)&d_B, bytes));
    CUDA_CHECK(cudaMalloc((void **)&d_C, bytes));

    cudaEvent_t t0, t1, t2, t3;
    CUDA_CHECK(cudaEventCreate(&t0));
    CUDA_CHECK(cudaEventCreate(&t1));
    CUDA_CHECK(cudaEventCreate(&t2));
    CUDA_CHECK(cudaEventCreate(&t3));

    // ---------- 第 2 步:H2D 拷贝 ----------
    CUDA_CHECK(cudaEventRecord(t0));
    CUDA_CHECK(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaEventRecord(t1));
    CUDA_CHECK(cudaEventSynchronize(t1));
    float h2d_ms = elapsedMs(t0, t1);

    // ---------- 第 3 步:启动 kernel(版本 1:一线程一元素) ----------
    int  blockSize = BLOCK_SIZE;
    int  gridSize  = (n + blockSize - 1) / blockSize;   // 向上取整的整数写法
    CUDA_CHECK(cudaEventRecord(t1));
    vecAddKernel<<<gridSize, blockSize>>>(d_A, d_B, d_C, n);
    CUDA_CHECK(cudaEventRecord(t2));
    CUDA_CHECK(cudaEventSynchronize(t2));
    CUDA_CHECK(cudaGetLastError());
    float k1_ms = elapsedMs(t1, t2);

    // ---------- 第 4 步:D2H 拷贝 ----------
    CUDA_CHECK(cudaEventRecord(t2));
    CUDA_CHECK(cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost));
    CUDA_CHECK(cudaEventRecord(t3));
    CUDA_CHECK(cudaEventSynchronize(t3));
    float d2h_ms = elapsedMs(t2, t3);

    bool ok1 = verify(h_A, h_B, h_C, n);

    // ---------- 版本 2:grid-stride,grid 取"刚好占满 GPU" ----------
    int dev = 0, smCount = 0, blocksPerSm = 0;
    CUDA_CHECK(cudaGetDevice(&dev));
    CUDA_CHECK(cudaDeviceGetAttribute(&smCount, cudaDevAttrMultiProcessorCount, dev));
    CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &blocksPerSm, vecAddKernelGridStride, blockSize, 0));
    int residentBlocks = smCount * blocksPerSm;
    int gsGrid = (residentBlocks < gridSize) ? residentBlocks : gridSize;

    CUDA_CHECK(cudaEventRecord(t1));
    vecAddKernelGridStride<<<gsGrid, blockSize>>>(d_A, d_B, d_C, n);
    CUDA_CHECK(cudaEventRecord(t2));
    CUDA_CHECK(cudaEventSynchronize(t2));
    CUDA_CHECK(cudaGetLastError());
    float k2_ms = elapsedMs(t1, t2);

    CUDA_CHECK(cudaMemcpy(h_C2, d_C, bytes, cudaMemcpyDeviceToHost));
    bool ok2 = verify(h_A, h_B, h_C2, n);

    // ---------- 性能报告 ----------
    double traffic = 3.0 * (double)bytes;                       // 总访存字节
    double bw1 = traffic / (k1_ms * 1.0e-3) / 1.0e9;            // GB/s
    double bw2 = traffic / (k2_ms * 1.0e-3) / 1.0e9;
    double gflops1 = (double)n / (k1_ms * 1.0e-3) / 1.0e9;      // 每元素 1 次加法
    double ai = (double)n / traffic;                            // FLOP / Byte
    double h2dBw = 2.0 * (double)bytes / (h2d_ms * 1.0e-3) / 1.0e9;
    double d2hBw = (double)bytes / (d2h_ms * 1.0e-3) / 1.0e9;

    printf("\n================ 性能报告 ================\n");
    printf("[PCIe] H2D   : %8.3f ms   %7.1f GB/s\n", h2d_ms, h2dBw);
    printf("[GPU ] kernel v1 (一线程一元素): %8.3f ms   %7.1f GB/s   %7.2f GFLOP/s\n",
           k1_ms, bw1, gflops1);
    printf("[GPU ] kernel v2 (grid-stride) : %8.3f ms   %7.1f GB/s\n", k2_ms, bw2);
    printf("[PCIe] D2H   : %8.3f ms   %7.1f GB/s\n", d2h_ms, d2hBw);

    double xfer = h2d_ms + d2h_ms;
    printf("\n[总账] 传输合计 %.3f ms, kernel %.3f ms, 传输/kernel = %.1f 倍\n",
           xfer, k1_ms, xfer / k1_ms);

    printf("\n[调度] grid v1 = %d blocks, 每个 block %d 线程\n", gridSize, blockSize);
    printf("       该 GPU 有 %d 个 SM, 每 SM 可常驻 %d 个 block -> 同时可跑 %d 个 block\n",
           smCount, blocksPerSm, residentBlocks);
    printf("       需要 %.2f 个 wave, 尾部浪费约 %.2f%%\n",
           (double)gridSize / residentBlocks,
           100.0 * (1.0 - (double)gridSize /
                    (ceil((double)gridSize / residentBlocks) * residentBlocks)));
    printf("       grid v2 = %d blocks (正好 1 个 wave, 每个线程循环约 %.1f 次)\n",
           gsGrid, (double)gridSize / gsGrid);

    printf("\n[算术强度] AI = %.4f FLOP/Byte  (每元素读 8 B 写 4 B -> 12 B, 做 1 次加法)\n", ai);
    printf("[判定] AI 远低于 A100 的机器平衡点 12.5 FLOP/Byte -> 纯内存带宽受限\n");
    printf("       优化方向是【减少字节搬运】, 而不是【减少指令数】\n");

    printf("\n[自检] version1: %s\n", ok1 ? "PASS" : "FAIL");
    printf("[自检] version2: %s\n", ok2 ? "PASS" : "FAIL");

    // ---------- 第 5 步:释放资源 ----------
    CUDA_CHECK(cudaEventDestroy(t0));
    CUDA_CHECK(cudaEventDestroy(t1));
    CUDA_CHECK(cudaEventDestroy(t2));
    CUDA_CHECK(cudaEventDestroy(t3));
    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));
    free(h_A);
    free(h_B);
    free(h_C);
    free(h_C2);

    return (ok1 && ok2) ? EXIT_SUCCESS : EXIT_FAILURE;
}
  • 【代码做什么?】
    1. main 解析 N(默认 2²⁴ = 16 777 216),计算 bytes = N × 4,打印每数组大小(64 MiB)与理论访存流量(192 MiB,含读 A、读 B、写 C 共 3N 个 float = 201.33 MB)。
    2. host 准备:四个 malloc——h_Ah_B(输入)与 h_Ch_C2(两个版本的输出,便于分别验证);initDataA[i] = iB[i] = 2i
    3. 第 1 步:三次 cudaMalloc 分配 d_Ad_Bd_C,各 N × 4 字节。
    4. 创建 4 个 cudaEvent 作为计时锚点。
    5. 第 2 步:用 cudaEventRecord(t0) / (t1) 夹住两次 cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice),测出 H2D 耗时;cudaEventSynchronize(t1) 保证时间戳已经落实。
    6. 第 3 步gridSize = (N + 255) / 256(整数向上取整,不用浮点 ceil(N/256.0),避免大 N 时的精度问题)。用事件夹住 vecAddKernel<<<gridSize, 256>>>(d_A, d_B, d_C, N),测出 kernel 耗时。kernel 里 i = blockIdx.x * blockDim.x + threadIdx.xif (i < n) 保护后写 C[i] = A[i] + B[i]
    7. 第 4 步:测 D2H 耗时并把 d_C 拷回 h_Cverify() 与 CPU 参考 A[i]+B[i] 逐元素比较(相对容差 1e-3),输出首个不匹配位置。
    8. 版本 2:查询当前设备 smCount(A100 = 108)与 vecAddKernelGridStride 的每 SM 常驻 block 数(约 8),得 residentBlocks = 864;取 gsGrid = min(864, gridSize)。kernel 用 stride = gridDim.x * blockDim.x 的 grid-stride 循环,每个线程处理约 65536 / 864 ≈ 75.85 个元素。
    9. 性能报告:分别算出 H2D/D2H 的实测带宽、两个 kernel 的 GB/s 与 GFLOP/s、传输与计算的耗时比、wave 数量与尾部浪费、算术强度与瓶颈判定。
    10. 第 5 步:销毁事件、cudaFree 三个设备数组、free 四个主机数组。以 ok1 && ok2 作为进程退出码(脚本化测试时非常有用)。
    11. 边界测试./vecadd 1000003(不是 256 的倍数)时,gridSize = (1000003 + 255)/256 = 39073907 × 256 = 1 000 192 > 1 000 003,因此最后 189 个线程是越界的,必须靠 if (i < n) 拦住;verify 会确认结果仍然完全正确。
  • 【并行机制与硬件映射解说】
    • block 内线程如何被划分为 warpblockDim.x = 256 → 每个 block 8 个满 warp(256/32 = 8,无空 lane)。warp 0 = threadIdx.x ∈ [0,31],warp 1 = [32,63],以此类推。i = blockIdx.x*256 + threadIdx.x 使得同一个 warp 的 32 个线程拿到连续的 32 个 i——这正是合并访问的前提。
    • warp 如何被调度到 SM 的 warp scheduler:A100 每个 SM 有 4 个 warp scheduler,每个 scheduler 管理一批 warp,每周期发射 1 条指令。一个 block 的 8 个 warp 会被尽量均匀地分给 4 个 scheduler(每 scheduler 2 个 warp)。由于 A100 每 SM 可驻留 8 个这样的 block(2048 线程 / 64 warp),每个 scheduler 平均管理 64/4 = 16 个 warp,切换余量充足。
    • 一个 SM 上可同时驻留多少 block/warp(A100,代入本 kernel)
      • 线程维度:2048 / 256 = 8 个 block
      • warp 维度:64 / 8 = 8 个 block
      • 寄存器维度:本 kernel 约 10 个寄存器/线程 → regsPerWarp = 10 × 32 = 320,按 256 粒度向上取整到 512byRegs = 65536 / (512 × 8) = 16 个 block
      • 共享内存维度:不使用共享内存,无限制
      • 硬件 block 上限:maxBlocksPerMultiProcessor = 32
      • 取 min → 8 个 block/SM = 2048 线程 = 64 warp = 100% 占用率
      • 全卡同时常驻:108 × 8 = 864 个 block = 864 × 256 = 221 184 个线程
    • 共享内存访问与 bank conflict:本 kernel 完全不使用共享内存,所有数据走寄存器与全局内存,因此 bank conflict 不适用。作为对照:若把 A、B 的一段数据放进 __shared__ float tile[256],那么 tile[threadIdx.x] 的访问模式是”32 个线程访问 32 个连续 float”——A100 的共享内存有 32 个 bank,每个 bank 宽 4 字节(32-bit),地址 addr 落在 bank = (addr / 4) % 32。线程 t 访问 tile[t]bank = t32 个线程命中 32 个不同 bank,零冲突(conflict-free)。反之若写成 tile[threadIdx.x * 32],则 bank = (t*32) % 32 = 0——32 个线程全部命中 bank 0,造成 32 路冲突(32-way bank conflict),共享内存访问被串行化 32 倍,访问时间从 20–30 cycle 恶化到 640–960 cycle。
    • 全局内存访问是否合并(按 32 线程 × 4 字节 = 128 字节事务粒度分析)
  版本 1 的访问模式: i = blockIdx.x*256 + threadIdx.x
  ---------------------------------------------------------------------------
  warp 内 32 个线程的 i 值:  base+0, base+1, base+2, ..., base+31
  A 数组地址:  A + (base+k)*4 字节,k = 0..31
  覆盖的字节区间: [A + base*4, A + base*4 + 128)  = 恰好 128 字节
  ---------------------------------------------------------------------------
  -> 1 次 128 字节内存事务 = 1 条 cache line            【完全合并 coalesced】
  -> 每个 warp 的 C=A+B 需要 3 次事务:读 A(128B)、读 B(128B)、写 C(128B)
  -> 有效带宽利用率 = 100%

  反例(步长访问) i = (blockIdx.x*256 + threadIdx.x) * 4
  ---------------------------------------------------------------------------
  warp 内 32 个线程访问的地址: A + base*16, A + (base+4)*16, ... 步长 16 字节
  覆盖的字节区间: [A + base*16, A + base*16 + 512)  = 512 字节 = 4 条 cache line
  ---------------------------------------------------------------------------
  -> 需要 4 次 128 字节事务,但只有 128 字节是有用数据
  -> 有效带宽利用率 = 128/512 = 25%   【浪费 4 倍带宽】

  反例(广播式越界尾块) N 不是 256 的倍数时的最后一个 block
  ---------------------------------------------------------------------------
  最后一个 warp 中只有部分 lane 满足 i < n,其余 lane 被掩码屏蔽。
  被屏蔽的 lane 不产生访存请求,因此【不浪费带宽】,只浪费发射槽。
  -> 这正是 if (i < n) 边界检查的代价:几乎没有。
  • 寄存器使用量与 warp 发散情况:两个 kernel 都只用约 8–14 个寄存器(grid-stride 版因为多一个 stride 与循环变量,略多 2–4 个),远低于 32 的满占用预算。分支方面:if (i < n) 仅在”最后一个 block 的最后一个 warp”上可能发散一次,占比 1/(N/32/…) 可忽略;grid-stride 版的 for 循环也让同一 warp 内所有 lane 的迭代次数一致(因为 nstride 对所有 lane 相同,只在最后一次迭代可能不同),几乎是完全无发散的代码

  • 【性能优化分析】

    • 占用率(occupancy)定量
      • blockDim = 256occupancy = blocksPerSm × blockDim / maxThreadsPerMultiProcessor = 8 × 256 / 2048 = 100%
      • 线程总数 N = 16 777 216 远大于单次可常驻的 221 184,并行度充足,不存在 Hello 示例那样的”喂不饱”问题。
      • block 数 65536 远大于 864,因此 block 会被反复”退役-补充”,形成 wave(波次)65536 / 864 = 75.85 个 wave。取整到 76 个 wave,尾部波次只填了 85.3%,浪费 1 - 65536/(76×864) = 0.195%——可忽略。这个数字说明:对于这种规整的流式 kernel,朴素的”一线程一元素”完全没有调度缺陷
    • 算术强度(arithmetic intensity)定量
      • 每元素:读 A[i](4 B)+ 读 B[i](4 B)+ 写 C[i](4 B)= 12 B;运算 1 次浮点加法 = 1 FLOP
      • AI = 1 / 12 = 0.0833 FLOP/Byte
      • A100 机器平衡点 AI* = 19.5e12 / 1555e9 = 12.54 FLOP/Byte
      • AI / AI* = 0.0833 / 12.54 = 0.66%离平衡点差 150 倍,是极端的内存带宽受限程序
    • Roofline 模型定量(代入 N = 2²⁴)
      • 访存流量 = 3 × N × 4 B = 3 × 16 777 216 × 4 = 201 326 592 B = 201.33 MB
      • 访存下界 T_mem = 201.33e6 B / 1555e9 B/s = 1.295e-4 s = 129.5 µs
      • 计算量 = N × 1 FLOP = 16.78 MFLOP;计算下界 T_flop = 16.78e6 / 19.5e12 = 0.86 µs
      • T = max(T_mem, T_flop) = 129.5 µs,且 T_mem / T_flop = 150.5
      • 在 Roofline 图上:可达性能 = AI × BW = 0.0833 × 1555 = 129.6 GFLOP/s,而峰值是 19 500 GFLOP/s,只有 0.66%
      • 实测若得到 140 µs,则 实测带宽 = 201.33e6 / 140e-6 = 1438 GB/s = 理论峰值的 92.5%——这已经是这类 kernel 能达到的实践上限(受限于 DRAM 刷新开销、ECC、页激活等)。
    • 瓶颈判定:内存带宽受限(memory-bandwidth-bound)。三条判据同时成立:(1) AI (0.0833) << AI* (12.54);(2) 实测带宽已达理论峰值的 92%,没有提升余地;(3) 减少指令数(如用 float4 向量化)只改变每字节的指令数,不改变总字节数,因此对耗时几乎没有影响。
    • PCIe 传输成本(这一讲最反直觉的量化结论)
      • H2D 传输 2 个数组 = 134.22 MB。用可分页内存(默认 malloc),实测带宽约 6–8 GB/s → 134.22e6 / 7e9 = 19.2 ms
      • 锁页内存cudaHostAlloc),实测约 25 GB/s → 134.22e6 / 25e9 = 5.37 ms(讲义指出锁页内存让 cudaMemcpy 快约 2 倍,与此吻合)。
      • D2H 传输 1 个数组 = 67.11 MB → 锁页下 2.68 ms
      • 对比 kernel 的 0.14 ms传输/计算 = (5.37 + 2.68) / 0.14 ≈ 57 倍(可分页内存下约 130 倍)。
      • 结论:向量加法这类低算术强度 kernel,应用总耗时几乎完全由 PCIe 决定,优化 kernel 本身是徒劳的。
    • 可执行的优化方向(按收益排序)
      1. 提高算术强度:把多个逐元素操作合并成一个 kernel(kernel fusion),例如把 C = A + B 与后续的 D = C * s 合并为 D = (A + B) * s,访存从 5N×4 B 降到 3N×4 B流量减少 40%,耗时同比例下降
      2. 使用向量化访存float4 让每个线程一次搬 16 字节,warp 一次事务仍是 128 字节但只需 8 个线程参与,指令数减少 4 倍;在带宽受限的 kernel 中收益有限,但在指令发射受限时收益明显。
      3. 让数据留在 device 上:把”上传-计算-下载”改成一串只在 device 内存中流转的 kernel,把 PCIe 往返从”每步一次”降到”整个流程一次”。
      4. 重叠传输与计算:用多个 cudaStream_t,让第 k+1 块的 H2D 与第 k 块的 kernel 执行重叠;配合锁页内存可让有效时间趋近 max(传输, 计算) 而非 传输 + 计算
      5. 用锁页内存cudaHostAlloc / cudaFreeHost):直接让 cudaMemcpy 快约 2 倍。
      6. grid-stride 的选择依据:本 kernel 两种版本性能基本相同(都打满带宽)。grid-stride 的真正价值在于:cudaOccupancyMaxActiveBlocksPerMultiprocessor 配合可以让 grid 恰好等于一个 wave,从而消除尾部波次浪费;并且当 N 极大(超过 gridDim.x 的 2³¹-1 上限)时仍然安全;还便于配合持久化 kernel(persistent kernel)与动态负载均衡。
性能优化技巧总结
  1. 先算 Roofline 再动手:用 T ≈ max(字节数/带宽, FLOP数/算力) 估出下界,若实测已接近下界,说明代码没有优化空间——继续调参只是浪费时间。为什么有效:它把”猜测式调优”变成”有上界的收敛判断”。
  2. 算术强度低于机器平衡点(A100 为 12.5 FLOP/Byte)时,一切优化都应指向”减少字节搬运”:融合 kernel、复用数据、降低精度(fp32→fp16)、跳过零元素。为什么有效:带宽是瓶颈时,减少指令数、增加并行度都无法缩短时间。
  3. 让每个 warp 的一次访存恰好覆盖一条 128 字节 cache line:保证 32 个线程访问连续的 32 个 4 字节字。为什么有效:一次事务换回 128 字节有用数据,有效带宽 100%;步长访问会让有效带宽降到 25% 甚至 3%。
  4. blockDim 取 32 的整数倍(128/256/512),并用 Device Query 验证占用率:避免空 lane 浪费,同时确保每 SM 常驻 8 个 block 达到 100% 占用率。为什么有效:占用率直接决定有多少 warp 可用于隐藏 400–800 周期的显存延迟。
  5. cudaOccupancyMaxActiveBlocksPerMultiprocessor() 而不是手算来决定 grid 大小:手算因寄存器分配粒度会偏差 ±1 个 block。为什么有效:运行时 API 反映硬件真实约束。
  6. grid 至少覆盖”SM 数 × 每 SM 常驻 block 数”(A100 + 256 线程 = 864 个 block)。为什么有效:低于这个数量就有 SM 完全闲置;高于它则进入多个 wave,尾部浪费只有约 0.2%。
  7. 优先消灭串行部分,而不是堆并行度:按 Amdahl 定律,s = 10% 时代码最多快 10 倍,P = 64 时效率已跌到 13.7%。为什么有效:串行段的耗时与 P 无关,它会成为不可压缩的下界。
  8. 用锁页内存(cudaHostAlloc)做主机端缓冲区:让 cudaMemcpy 快约 2 倍。为什么有效:DMA 需要物理地址固定的缓冲区,可分页内存必须先经一次 CPU 拷贝到内部暂存区。
  9. 尽量不要在计算流程中间把数据搬回主机:把多次”上传-计算-下载”合并为”一次上传、多次计算、一次下载”。为什么有效:PCIe Gen4 x16 单向约 25–32 GB/s,仅为 A100 显存带宽(1555 GB/s)的 1/50。
  10. 用多 stream 重叠传输与计算:为什么有效:PCIe 引擎与 SM 是独立硬件,单 stream 时它们被迫串行,多 stream 可让总时间趋近 max(传输, 计算)
  11. 共享内存访问要避开 bank conflict:保证同一个 warp 内各线程访问的 bank 不同(例如 tile[threadIdx.x])。为什么有效:32 路冲突会把 20–30 cycle 的共享内存访问拖到 ~700 cycle,与直接访存 DRAM 无异。
  12. 避免 warp 发散:把 threadIdx.x % 2 这类按 lane 奇偶分支改写成”前半 warp 做 A、后半 warp 做 B”,或用分支无关的算术(select/三元运算)。为什么有效:发散让两条路径串行执行,利用率最多掉一半。
  13. 不要在 kernel 里做大量 printf:device 端 printf 会串行化输出并占用缓冲。为什么有效:它把并行的核心变成串行的 I/O。
  14. cudaEvent 分别测量 H2D、kernel、D2H 三段耗时:只有拆开测才知道瓶颈到底在哪一段。为什么有效:向量加法的传输耗时是计算的 57 倍,不测就永远找不到真正的瓶颈。
  15. 为 kernel 设置 __launch_bounds__ 控制寄存器用量:当 nvcc 为了减少指令数而用掉过多寄存器、导致占用率掉到 50% 以下时,限制寄存器数往往更快。为什么有效:它把编译器的优化目标从”单线程最快”改成”整机吞吐最高”,这正是 GPU 的优化目标。
关键要点
  • 2004 年的功耗墙终结了主频红利,硬件转向并行;2007 年的 CUDA 让 GPU 的并行算力变得可编程。 Dennard 缩放的理想模型(f → 1/α、每管功耗 → α²)在 VDD 无法继续下降后失效,芯片设计由此分化为延迟导向的 CPU 与吞吐量导向的 GPU。
  • 并行不改变算法复杂度,只提供受 Amdahl 定律约束的固定倍数收益:Speedup ≤ 1/s 串行占比 25% 时上界是 4×;10% 时用 64 个核只能拿到 8.77×、效率 13.7%。先优化串行瓶颈,再谈并行技巧。
  • GPU 靠”延迟隐藏”而非”降低延迟”取胜:一个 SM 上驻留 64 个 warp(A100),用海量线程的并发掩盖 400–800 周期的显存延迟。 Little 定律给出定量需求:在途字节 = 延迟 × 带宽,A100 需要约 772 KB(每个 SM 约 56 个 warp 的在途 load)才能打满带宽——这就是”GPU 必须用海量线程”的数学解释。
  • 算术强度(FLOP/Byte)与机器平衡点(A100 为 12.5 FLOP/Byte)决定了 kernel 的瓶颈类型。 向量加法的 AI 仅 0.083,比平衡点低 150 倍,理论下界 129.5 µs 完全由 201.33 MB 的访存流量决定;减少字节数(kernel fusion、tiling、向量化)是唯一的优化方向。
  • CUDA 程序的五步模板(分配 → 拷贝入 → 启动 → 拷贝回 → 释放)是所有后续讲座的基准,其中”拷贝入/拷贝回”往往比 kernel 本身贵几十倍。 向量加法在 A100 上 kernel 约 0.14 ms,而 PCIe 往返约 8.05 ms(锁页内存)——Bandwidth 是现代计算机系统的重力
  • 统一着色器架构(2006, G80/Tesla)把分离的顶点/像素着色器合并为通用标量处理器阵列,是 CUDA 得以存在的硬件前提;而”占用率 = 各类资源上限取 min”这条规则从 CC 1.3 到 sm_90 二十年未变。 residentBlocks = min(byThreads, byRegs, bySharedMem, byBlocksPerSM) 是每一次 kernel 调参的出发点。
常见陷阱与注意事项
  • 忘记 cudaDeviceSynchronize() / 在 kernel 结果未就绪时就使用它:现象是读回的数组全是初始值或垃圾数据,或者程序打印的顺序与预期不符 → kernel 启动从 CUDA 1.0 起就是异步的;要在 host 使用结果,必须 cudaDeviceSynchronize()、或用 cudaMemcpy(隐式同步)、cudaStreamSynchronize()cudaEventSynchronize()。注意:cudaMemcpy 会同步,所以”忘了同步”常常被它掩盖,直到你插入计时或并行 stream 才暴露。
  • 忘记 __syncthreads():现象是 block 内线程读到的共享内存值来自”写入之前”的旧值,结果随运行次数变化(非确定性)→ 在”写共享内存 → 读共享内存”之间必须调用 __syncthreads()。同时要清楚它的边界:__syncthreads() 保证同一个 block 内所有线程在 barrier 之前的全局/共享内存访问对该 block 内所有线程可见;但不保证其它 block 的执行进度、不保证数据已写回 DRAM、不保证 printf 的输出顺序。
  • __syncthreads() 放进分支里:现象是程序挂死(hang)或行为不可预测 → __syncthreads() 必须位于所有线程都会执行到的位置。若部分线程因 if 跳过 barrier,硬件会一直等待那些永远不会到达的线程。正确做法是把 barrier 提到条件外面,或者让人为的 barrier 前置于分支。
  • 共享内存 bank conflict:现象是使用了共享内存却比不用还慢 → A100 的共享内存有 32 个 bank、每 bank 4 字节宽。tile[threadIdx.x](bank = t)零冲突;tile[threadIdx.x * 32](bank 恒为 0)是 32 路冲突,访问被串行化 32 倍。解决方式:改访问模式、加 padding(__shared__ float tile[32][33])、或做对角线重排。
  • warp 发散(warp divergence):现象是某个 kernel 的耗时远高于按指令数估算的值 → 同一个 warp 内不同 lane 走不同分支时,硬件串行执行各条路径。典型诱因是 if (threadIdx.x % 2)if (i < n) 之外的按 lane 分支、以及负载方差大的并行源(图的节点度、三角矩阵的行长)。正确做法:让人为分支与 warp 边界对齐(前 16 个 lane 走 A、后 16 个走 B),或改用分支无关的算术。
  • 未检查 CUDA API 返回值:现象是程序”静默失败”——kernel 根本没启动,或 cudaMalloc 失败返回了 nullptr,后续访问产生难以定位的错误 → 所有 CUDA 调用都必须包在 CUDA_CHECK 宏里;kernel 启动后立刻 cudaGetLastError()(因为 <<<>>> 没有返回值)。特别注意:不检查时,某次成功的 API 调用会清空之前的错误码,让你永远看不到真正的失败点。
  • 越界访问:现象是结果局部错误、cuda-memcheck/compute-sanitizer 报 “Invalid global write”,或者在某些运行中”碰巧正确” → gridSize = (N + blockSize - 1) / blockSize 会向上取整,最后一个 block 里有 gridSize*blockSize - N 个线程的索引是越界的,必须在 kernel 里用 if (i < n) 拦住。整数除法写成 N / blockSize(向下取整)则是另一个方向的错误:网格覆盖不足,尾部元素永远不被计算
  • cudaMemcpy 方向写错:现象是 invalid argument 错误,或数据被拷贝到错误方向后结果全错 → 四个参数是 (目的指针, 源指针, 字节数, 方向),方向枚举 cudaMemcpyHostToDevice / DeviceToHost / DeviceToDevice 必须与两个指针的实际位置一致。尤其注意上下行的指针顺序是反的:H2D 时目的在前(d_A, h_A),D2H 时也是目的在前(h_C, d_C),很容易把其中一个写反。
  • 共享内存容量超限:现象是 kernel 启动失败报 invalid argument,或 cudaFuncSetAttribute 报错 → A100 上静态共享内存每 block 上限是 48 KB,要超过它必须改用动态共享内存并调用 cudaFuncSetAttribute(kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, bytes) 显式 opt-in(A100 每 block 最多 163 KB、每 SM 共 164 KB)。同时要意识到:共享内存开得越大,每 SM 能驻留的 block 越少,占用率会同步下降(例如 48 KB/block 时 A100 每 SM 只能放 3 个 block,占用率 37.5%)。
  • host 指针与 device 指针混用:现象是段错误(segmentation fault)、invalid device pointer,或者 kernel 里访问到主机内存导致非法地址错误 → cudaMalloc 返回的指针只能用于 device 代码与 CUDA API,绝不能传给 free()memcpy();反之 malloc 返回的指针传入 kernel 会触发 Illegal address access。命名上养成习惯:device 指针加 _d 后缀,host 指针加 _h 后缀
  • 主机端读回结果前忘记检查 kernel 是否真的跑过:现象是 verify 报 FAIL 但看不出原因 → 先 CUDA_CHECK(cudaGetLastError())CUDA_CHECK(cudaDeviceSynchronize()),两步都用宏包住;然后把 cudaMemcpy 回来的数组与 CPU 参考实现逐元素对比并打印首个不匹配的下标与期望值,这比打印”FAIL”有用一百倍。
  • 对”默认流”的同步语义想当然:现象是在多 stream 程序中,计时变快但结果变错 → kernel 启动是异步的,多个 stream 之间的顺序完全没有保证;默认流(stream 0)有特殊的隐式同步行为,混用默认流与自定义流时尤其容易出错。用 cudaEventRecord + cudaStreamWaitEvent 显式建立依赖关系。
思考题(带答案)

Q1. 某应用有 8% 的代码无法并行化。在 128 个处理器的机器上,理论最大加速比是多少?如果工程师花两个月把这 8% 串行代码优化到只占原始串行时间的 1/4(即优化后串行部分只占并行版本的 2%),加速比会变成多少? 按 Amdahl 定律 Speedup(P) = 1/(s + (1-s)/P)。第一种情况 s = 0.08、P = 128:1/(0.08 + 0.92/128) = 1/(0.08 + 0.0071875) = 1/0.0871875 = 11.47×,效率 11.47/128 = 8.96%。注意这远低于上界 1/0.08 = 12.5×,说明此时还受并行部分限制。第二种情况:串行部分的绝对时间变成原来的 1/4,设原始总时间为 1,则新串行时间 0.02、新并行时间仍为 0.92;新的串行占比 s' = 0.02/(0.02+0.92) = 0.02128Speedup(128) = 1/(0.02128 + 0.97872/128) = 1/(0.02128 + 0.007646) = 34.6×——串行时间减少 4 倍,总加速比提升 3 倍,且上界从 12.5× 抬升到 47×。这正是”先消灭串行瓶颈”的量化依据。

Q2. 一个 kernel 在 A100 上每个线程需要 40 个寄存器、每个 block 256 个线程、使用 0 字节共享内存。(a) 每 SM 最多能驻留多少个这样的 block?(b) 占用率是多少?(c) 如果改用 __launch_bounds__(256, 8) 强制编译器把寄存器压到 32 个以内,占用率变成多少,是否一定更快? (a) 逐项取 min:线程维度 2048/256 = 8 个 block;warp 维度 64/(256/32) = 8 个 block;寄存器维度按 warp 为单位、256 个寄存器为粒度向上取整:regsPerWarp = 40 × 32 = 1280,已对齐到 256 的倍数,byRegs = 65536/(1280 × 8) = 6.4 → 6 个 block;共享内存维度无限制;硬件 block 上限 32。取 min → 6 个 block/SM。(b) 占用率 = 6 × 256 / 2048 = 768/2048 = 37.5%(warp 粒度也是 6×8/64 = 37.5%)。(c) 压到 32 个寄存器后 byRegs = 65536/((32×32)×8) = 65536/8192 = 8 个 block,占用率变成 8×256/2048 = 100%提升 2.67 倍。但不一定更快:编译器为了把寄存器压到 32 个,可能把局部变量溢出到 local memory(本质是显存),每次溢出访问都要 400–800 周期,反而拖慢。判定方法:编译时加 -Xptxas -v 查看 “spill stores/loads” 计数,若为 0 则收益是确定的;若出现大量 spill,应当接受较低占用率,或改写 kernel 降低寄存器需求。

Q3. 在 A100(1555 GB/s、19.5 TFLOPS)上,一个 kernel 处理 N = 2²⁴ 个 float,每个元素读 1 次、写 1 次、做 4 次浮点运算。(a) 它的算术强度是多少?(b) 是计算受限还是带宽受限?(c) 理论最短耗时是多少?(d) 如果实测耗时是理论值的 1.4 倍,最可能的两个原因是什么? (a) 每元素访存 4 B 读 + 4 B 写 = 8 B,运算 4 FLOP,AI = 4/8 = 0.5 FLOP/Byte。(b) 平衡点 AI* = 19500/1555 = 12.5 FLOP/ByteAI = 0.5 << 12.5内存带宽受限。(c) 访存下界 = 2 × 2²⁴ × 4 B / 1555e9 = 134.22e6/1555e9 = 86.3 µs;计算下界 = 2²⁴ × 4 / 19.5e12 = 67.1e6/19.5e12 = 3.44 µs。取 max → 86.3 µs。(d) 最可能的两个原因是:(1) 访问未完全合并或未对齐——若每个线程用步长访问而非连续访问,一个 warp 会占用多条 cache line,有效带宽降到 25%–50%;(2) 占用率不足导致延迟无法被隐藏——按 Little 定律 A100 需要约 56 个 warp 的在途 load 才能打满带宽,若寄存器用量把占用率压到 50%(32 warp)以下,实际带宽就会显著低于峰值。其他可能原因还包括 L2 到 DRAM 的写合并效率、ECC 开销,以及测得的”1.4 倍”里混入了 kernel 启动开销(约 3–10 µs,占 86 µs 的 4%–12%)。