ETHZ-异构系统编程笔记-全-

ETHZ 异构系统编程笔记（全）

001：课程介绍与异构计算概述

在本节课中，我们将学习异构计算系统的基本概念、课程目标以及现代计算系统中各种加速器（如GPU、TPU等）的概览。我们将从经典的弗林分类法开始，探讨单指令多数据（SIMD）架构，并了解其在当前CPU扩展和专用加速器中的应用。

课程概述

欢迎来到“异构系统编程：基于GPU和加速器”课程的第一讲。本课程旨在教授如何为包含GPU和其他加速器的异构计算系统进行编程。我们将介绍异构系统的需求、各种设备类型，并重点学习GPU并行编程。从下一次讲座开始，我们将深入探讨GPU架构、编程模型和并行模式。

异构计算的需求

现代重要工作负载（如机器学习、人工智能、生物信息学、医学成像等）对性能和能效的扩展提出了挑战。这催生了对异构设备的需求。除了传统的中央处理器（CPU），图形处理器（GPU）、现场可编程门阵列（FPGA）以及张量处理单元（TPU）等专用加速器在现代计算系统中变得越来越普遍。

计算机架构分类：弗林分类法

理解异构系统的一个起点是回顾迈克尔·弗林在1966年提出的计算机分类法。他根据设备如何处理数据，将计算机分为四类：

单指令单数据（SISD）：一条指令操作一个数据元素。这是传统的顺序处理器。
单指令多数据（SIMD）：一条指令流同时操作多个数据元素。这是本课程的重点，也是GPU的基础架构之一。
多指令单数据（MISD）：多条指令操作一个数据元素。现实中例子较少，脉动阵列是近似例子。
多指令多数据（MIMD）：多条指令流操作多个数据元素。例如多核处理器。

SIMD架构实例：CPU扩展

即使在CPU内部，也存在异构性。一个典型的例子是CPU中的SIMD指令集扩展。其核心思想是：一条指令同时处理多个数据元素。

例如，我们可以使用一个32位寄存器，将其划分为4个8位部分。通过修改算术逻辑单元（ALU），阻止这4个8位值之间的进位传播，就能用一条指令同时完成4个独立的加法运算。

一个开创性的例子是英特尔在90年代推出的MMX扩展。它主要用于多媒体和图形操作。以下是一个使用MMX指令进行图像合成的伪代码示例，其核心操作是并行比较和筛选像素：

; 假设 MM0 寄存器包含多个像素数据（例如，8个8位像素值）
; MM1 寄存器预先加载了“蓝色”的参考值（例如，8个相同的代表蓝色的8位值）

; 并行比较：生成一个掩码，标记哪些像素是蓝色
PCMPEQB MM1, MM0  ; 比较 MM0 和 MM1，结果掩码存入 MM1

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/ethz-hetsys-prog/img/8653ce9ed13dbbdeb34f33283c0474f2_2.png)

; 使用生成的掩码，并行地从背景图（Y）和前景图（X）中筛选出需要的像素
; （此处省略具体的位操作指令）
; ...

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/ethz-hetsys-prog/img/8653ce9ed13dbbdeb34f33283c0474f2_4.png)

; 合并筛选后的像素，生成最终图像
POR MM2, MM3  ; 合并结果

现代异构计算系统概览

现代片上系统（SoC）或计算系统通常集成多种设备：

GPU（图形处理器）：不仅用于图形，也广泛用于通用并行计算。其内部也有异构性，例如包含用于通用计算的CUDA核心和用于机器学习的专用张量核心。
专用加速器：
- TPU（张量处理单元）：谷歌设计的用于神经网络训练的脉动阵列，专精于矩阵乘法。
- Wafer-Scale Engine（晶圆级引擎）：Cerebras公司推出的巨型机器学习加速器。
- Groq张量流处理器：采用确定性执行模型，专为规则并行负载设计。
- 特斯拉Dojo：特斯拉自研的神经网络训练系统。
近内存/存内处理（PIM）：为了缓解“内存墙”问题，将处理单元放置在内存芯片内部。例如：
- UPMEM的PIM引擎：在DRAM芯片内集成小型处理器（DPU）。
- 三星HBM-PIM：在高带宽内存（HBM）中集成AI处理单元。
- SK海力士GDDR6-AiM：在GDDR6内存中集成加速单元。

设备间通信与一致性

在包含多种设备的系统中，高效的通信和数据一致性至关重要。诸如CCIX、OpenCAPI和最新的Compute Express Link（CXL） 等一致性互连标准，旨在为不同设备提供统一的内存视图，实现更紧密的集成和更细粒度的协作，从而提升性能、能效和可编程性。

课程目标与要求

本课程的主要目标是学会如何通过编程、分析工作负载、提出卸载与调度策略来利用现有的异构设备。

课程关键收获：

加深对计算机体系结构和异构系统的理解。
获得异构架构（如GPU/FPGA）编程的技术技能。
培养批判性思维和分析能力。
熟悉异构计算的前沿研究方向。
提升技术演讲能力。

先修知识： 需要具备数字设计和计算机体系结构的基础，熟悉C/C++编程。对FPGA或GPU编程有基本了解更佳，但最重要的是对计算机架构和解决系统效率问题有浓厚兴趣。

课程安排与资源

形式：讲座将以预录制视频为主，在每周五上午提供。我们会有定期的在线会议讨论项目和答疑。
项目：课程核心是个人实践项目。下周将详细介绍项目选题，之后学生将与指定的导师定期会面，开展项目工作。
资源：所有课程材料（幻灯片、视频、推荐阅读）均可在课程网站找到。也推荐参考过往学期的资料。

下节预告：SIMD处理器与GPU

下一讲我们将深入探讨数据并行的利用。我们将详细讲解SIMD处理器，包括阵列处理器与向量处理器的区别，并阐述现代GPU如何融合这两者的特性。我们也将初步了解如何将计算任务在CPU（处理串行部分）和GPU（处理大规模并行部分）之间进行分配，引导大家建立并行编程思维。

例如，一个简单的数组加法 C[i] = A[i] + B[i] 在GPU上的并行化，其核心思想是将每个循环迭代映射到一个独立的线程上：

// 每个GPU线程执行以下操作
int i = threadIdx.x + blockIdx.x * blockDim.x;
if (i < N) {
    C[i] = A[i] + B[i];
}

总结

本节课我们一起学习了异构计算的基本动机，回顾了弗林分类法，并重点了解了SIMD架构在CPU扩展中的应用。我们还概览了现代计算系统中丰富的异构设备生态，包括GPU、TPU、PIM等，以及连接它们的一致性互连技术。最后，我们明确了本课程的学习目标、要求和安排。希望这为大家后续深入学习异构系统编程打下坚实的基础。

002：SIMD处理与GPU架构 🚀

在本节课中，我们将学习SIMD处理器的基本架构和编程模型，并了解现代图形处理单元如何作为SIMD处理器的一种实现。我们将从SIMD的基本概念开始，逐步深入到向量处理器和GPU的工作原理。

概述

SIMD是一种计算范式，其中单个指令同时对多个数据元素进行操作。这种范式是现代高性能计算，尤其是GPU加速计算的核心。理解SIMD是理解GPU编程和异构系统优化的基础。

SIMD处理：时间与空间

上一节我们介绍了SIMD的基本概念，本节中我们来看看SIMD的两种主要实现方式：在时间上或在空间上执行。

SIMD处理的核心是单个指令操作多个数据元素。这可以通过两种方式实现：

空间并行：使用多个处理单元在同一时间执行同一指令。
时间并行：使用同一个处理单元在连续的时间步上对多个数据元素执行同一指令。

以下是两种经典SIMD处理器的例子：

阵列处理器：属于空间并行。它拥有多个相同的处理单元，每个单元在同一时钟周期内对不同的数据元素执行相同的操作。
向量处理器：属于时间并行。它拥有多个专用的功能单元，一条向量指令中的不同操作（如加载、加法、存储）由不同的单元在流水线上依次处理。

向量处理器详解

现在，让我们更详细地探讨向量处理器，它是理解现代SIMD架构的重要基础。

向量是一维的数字数组。向量处理器的指令直接对向量进行操作，而非标量值。要实现这一点，需要满足几个基本要求：

以下是向量处理器的核心组件：

向量寄存器：用于存储向量。其大小为 n，意味着它可以容纳 n 个 m 位的数据值。总容量为 n * m 位。
向量控制寄存器：
- 向量长度寄存器：定义当前操作的向量长度。
- 向量步长寄存器：定义内存中连续向量元素之间的地址间隔。
- 向量掩码寄存器：用于条件执行，指示对哪些向量元素进行操作。

存储系统：体交织内存

为了能够每个周期加载/存储一个向量元素（即使单次内存访问延迟很高），向量处理器使用体交织内存。

内存被划分为多个可以独立访问的体。数据字以交错方式分布在这些体上。例如，元素0在体0，元素1在体1，依此类推。这样，处理器可以每个周期向不同的体发起访问请求。虽然每个请求仍有延迟，但在初始延迟后，每个周期都可以从不同的体获得一个数据，从而实现了高带宽的向量数据流。

代码向量化

程序要在向量处理器上高效运行，其代码必须是可向量化的。

一个循环可向量化的关键条件是：每次迭代之间相互独立。例如，对于循环 for(i=0; i<50; i++) C[i] = (A[i] + B[i]) >> 1，每次 i 的迭代计算都不依赖于其他 i 的结果，因此它可以被转换为向量指令。

向量化后的伪代码大致如下：

VL = 50        // 设置向量长度为50
VS = 1         // 设置步长为1
V0 = load[A]   // 向量加载A
V1 = load[B]   // 向量加载B
V2 = V0 + V1   // 向量加法
V2 = V2 >> 1   // 向量右移（等价除以2）
store[C] = V2  // 向量存储到C

性能优化技术

在理解了基础执行过程后，我们来看看如何优化向量处理器的性能。

以下是两种关键优化技术：

链接：允许将一个功能单元的结果直接转发给下一个功能单元，无需等待整个向量操作完成。这减少了操作间的空闲时间。
多端口内存体：为每个内存体增加负载/存储端口，允许同时进行多个向量加载/存储操作，使它们能够重叠执行。

通过结合使用链接和多端口内存，可以显著减少总执行周期数。

条件执行与向量掩码

当循环中包含条件语句时，我们需要使用向量掩码来实现条件执行。

向量掩码寄存器为每个向量元素包含一个标志位。在执行一条向量指令时，只有掩码位为“真”的对应元素才会被实际计算和更新。这允许向量处理器高效地处理像 if (A[i] != 0) C[i] = A[i] * B[i] 这样的条件循环。

体冲突与解决方案

当数据访问模式与内存体组织不匹配时，会发生体冲突，降低内存带宽利用率。

体冲突的一个典型例子是：当访问步长与体数量有公因数时（例如，16个体，步长为2），多个请求会指向同一个体，造成排队等待。解决方案包括：

增加体的数量。
为每个体增加更多访问端口。
改变数据布局以匹配访问模式（例如，矩阵转置）。
使用随机或伪随机交织策略将地址映射到体。

现代SIMD：从CPU到GPU

SIMD思想已广泛应用于现代计算设备。

在CPU中，通过SIMD指令集扩展实现，如x86的SSE、AVX指令集，它们提供了越来越宽的向量寄存器。在机器学习加速器中，如谷歌的TPU，也大量使用SIMD处理。

然而，当今SIMD最主要的代表是图形处理单元。GPU本质上是一个大规模并行的SIMD引擎。其核心由许多流多处理器组成，每个SM内部包含多个标量核心、寄存器文件和共享内存（一种体交织内存）。

GPU编程模型：SPMD

GPU的独特之处在于其编程模型与执行模型的分离。

编程模型：程序员使用单程序多数据模型编写代码。程序员创建大量线程，每个线程执行相同的核函数代码，但处理不同的数据。
执行模型：硬件将多个线程（例如32个）动态分组为一个线程束。线程束中的所有线程以锁步方式执行相同的指令，就像一个SIMD单元。这就是SIMT执行模型。

这种设计结合了编程的灵活性和硬件执行的效率。线程可以独立发散，而硬件则通过掩码等技术管理控制流分支。

GPU执行：细粒度多线程

为了进一步隐藏长延迟操作（如访问全局内存），GPU采用了细粒度多线程技术。

GPU的SM会同时管理多个活跃的线程束。当一个线程束因等待内存数据而停滞时，SM会立即切换到另一个就绪的线程束执行。这种快速的上下文切换使得计算单元始终保持忙碌，极大地提高了硬件利用率。

总结

本节课中我们一起学习了SIMD处理的基本原理。我们从阵列处理器和向量处理器的概念出发，探讨了向量化、体交织内存、条件执行等关键技术。随后，我们将这些概念与现代GPU架构联系起来，理解了GPU如何通过SPMD编程模型和SIMT执行模型，将标量线程程序映射到底层的大规模SIMD硬件上高效执行。掌握这些基础知识是后续进行高效GPU编程和优化的关键。

003： GPU软件层次结构

在本节课中，我们将开始学习GPU编程，并重点介绍GPU的软件层次结构。我们将探讨计算如何分配给不同的执行线程，以及这些线程如何被分组以实现更高效的执行并映射到GPU硬件上。

概述

上一讲我们介绍了GPU和SIMD处理器的硬件架构与编程模型。GPU本质上是一种SIMD引擎。即使像NVIDIA H100 GPU这样的架构看起来是同质的，深入其核心会发现不同的执行单元，这让人联想到我们之前讲解的向量处理器和SIMD处理器。GPU核心包含用于32位整数、32位浮点或64位浮点运算的单元，也有用于内存加载/存储的专用单元，以及用于神经网络和机器学习的专用张量核心。GPU核心可以访问不同的内存空间：寄存器、L1数据缓存、共享内存（一种暂存器内存），而在GPU核心之外，还有L2内存和片外DRAM（全局内存）。我们将在下一讲详细讨论GPU内存，本节课我们聚焦于软件层次结构。

同时要记住，GPU内部以线程束为单位执行指令，线程束是基本的SIMD单元。在GPU核心的流水线上，我们采用细粒度多线程执行，这意味着每个周期我们从不同的线程束获取指令。通过这种方式，线程束的执行在GPU流水线上交错进行，这种细粒度多线程执行能够容忍长延迟操作，例如访问片外内存的指令，其延迟可以通过在流水线中执行其他指令来隐藏。

现在，让我们开始讨论GPU编程。在上一讲的幻灯片中，我们引用了Fisher在1983年ISCA上的论文，其中阐述了为SIMD处理器编程的难度。幸运的是，GPU编程相对容易一些。GPU通用编程的兴起始于十多年前，在某种意义上代表了高性能计算的民主化，因为人们可以以相对较低的成本在个人电脑上获得具有极高浮点运算吞吐量的强大GPU。许多工作负载，如矩阵计算、图像处理、神经网络等，都能从这类设备中受益。但我们需要学习如何为这类系统编程，理解新的编程模型，并了解编程时的关键瓶颈以及缓解这些瓶颈的方法。

从CPU到GPU：异构计算的基本流程

异构计算通常从CPU和CPU内存（系统主内存）开始，主要包含三个步骤：

将输入数据从CPU内存传输到GPU内存（通常称为全局内存）。
启动一个在内核（将在GPU上运行的函数）上执行。
内核执行完毕后，CPU获取结果，将其从GPU内存移回CPU内存。

这三个步骤体现在典型的CPU-GPU异构系统编程结构中：我们通常在CPU上运行顺序或适度并行的代码段，而将大规模并行代码段卸载到GPU上执行。通常流程如下：首先，一个或多个线程在CPU上运行代码；在某个时刻，它们会将一些计算任务卸载到GPU设备端，这个计算任务称为内核；内核使用大量线程执行，这些执行线程被分组为线程块（在NVIDIA CUDA术语中）。内核完成后，控制权返回给CPU，可能还会有一些顺序或适度并行的执行，之后CPU可能启动另一个内核。

GPU编程模型：SPMD/SIMT

GPU编程采用单程序多数据编程模型。我们通常为单个标量线程编写代码，但这些线程随后会被分组为SIMD单元（在NVIDIA术语中称为线程束），并在类似SIMD硬件的GPU核心上执行。

在SIMD编程模型中，每个程序或线程处理不同的数据，可以执行不同的控制流路径，但也可以同步。因为底层硬件是SIMD的，当我们把许多线程或过程分组到同一个线程束时，我们希望它们执行基本相同的计算，这样执行速度最快。否则，我们将不得不使用谓词执行、掩码等技术，先为线程束中的部分线程运行计算，再为其他线程运行，这虽然可行且大大方便了编程，但为了编写更优化的程序，我们必须意识到这一点。

CUDA和OpenCL是使用SPMD（也称为SIMT，单指令多线程）的编程模型。它们通常被称为同步编程框架，因为同步执行内核的不同线程或所有线程的方式是终止内核，这被称为全局或栅栏同步。

主机（通常是CPU）负责分配内存、复制数据和启动内核。设备（通常是GPU）则执行内核。

GPU执行层次：网格、线程块与线程

那么GPU如何执行内核呢？GPU运行线程，但单个内核的所有线程被分组为网格（在OpenCL中称为NDRange）。在网格内部，有多个线程块（在NVIDIA术语中简称为块，在OpenCL中称为工作组）。在同一线程块内运行的所有线程可以共享内存和同步，原因是它们将在同一个GPU核心内执行。然后我们还有线程（在OpenCL中称为工作项）。

以下是CUDA程序的典型结构：

包含在CPU上运行的函数原型（如幻灯片中所示的串行函数）。
包含在GPU上运行的内核函数。注意，GPU内核总是以 __global__ 限定符开头。
然后是主机代码，即执行开始的地方。在某个时刻，主机（CPU）会将计算任务卸载到GPU端。

在主机的主函数中，我们首先在设备（GPU）上分配内存，使用 cudaMalloc API 来分配一定字节数的空间给指针 d_in。然后需要将数据从主机传输到设备，即使已经在GPU内存中分配了空间，数据最初仍驻留在CPU内存中，因此我们需要从 h_in 数组复制到 d_in 数组。接着，我们需要定义执行配置，即网格中的块数和每个块的线程数。然后我们可以调用内核，这是将计算任务卸载到GPU的时刻，内核带有参数以及之前定义的执行配置。内核终止后，我们可以将结果从GPU内存传输回CPU内存，再次使用 cudaMemcpy 操作，但这次目标是在CPU内存中分配的 h_out，源是在GPU内存中分配的 d_out。这个复制、启动内核、再复制回CPU的序列可以根据需要重复，例如在多次调用同一内核的迭代算法中，或者在其他具有多个内核卸载到GPU的应用程序中。

在内核内部，我们可以使用自动变量，这些变量透明地分配给GPU核心的寄存器。我们还可以使用共享内存，这是所有在同一线程块、同一GPU核心内运行的线程都可以访问的共享内存空间，我们使用 __shared__ 限定符来标识那些分配在共享内存中的数组。在同一线程块内，线程也可以同步，这被称为块内同步，我们使用 __syncthreads() 指令来实现。

代码示例：向量加法

为了开始更深入地理解如何为GPU系统编程，我们将从一个简单的代码示例——向量加法开始。幻灯片上显示的代码是主机代码，是一个名为 vecAdd 的函数。假设这个函数是从主机CPU的 main 函数调用的。如你所见，在 vecAdd 函数中，我们首先做的是分配GPU内存，为输入数组 A 和 B 以及输出数组 C 分配内存。然后我们将数据从主机复制到设备，如参数所示，我们将所有字节从驻留在CPU内存中的数组 A 复制到驻留在GPU内存中的数组 A_d。我们对数组 B 执行相同的操作。接着，我们将在GPU上启动计算，这是内核调用，我们将在后面的幻灯片中详细查看。最后，当内核终止时，我们将结果从设备内存（GPU内存）复制回主机内存（CPU内存），最后可以释放GPU内存中的这些空间。

向量加法作为我们的第一个编程示例，就是执行两个数组的逐元素加法。在最直接的实现中，我们可以为输入和输出向量中的每个元素分配一个GPU线程。因此，每个线程将只负责一次加法操作：它需要访问内存，读取 A[0] 和 B[0]，将它们相加，然后将结果 C[0] 写入内存。这是线程0执行的计算，线程1执行类似的计算，依此类推。所有这些线程的集合，如前所述，称为网格。但我们需要一种方法将这些线程分配给GPU核心。GPU核心的数量因GPU而异，因此我们将线程分组为线程块。例如，这四个线程可能属于块0、块1、块2和块3。在这个简单的示例中，每个块只有四个线程，通常我们至少会有32个（NVIDIA GPU的线程束大小），更典型的是这个线程束大小的倍数。

启动网格与内核编写

我们如何启动一个网格呢？首先，同一网格中的线程执行相同的函数，即内核。网格可以通过调用内核并使用适当的网格和块大小配置来启动。记住，我们需要先定义执行配置，即每个块的线程数（本例中为512）和总块数，然后我们可以使用这行代码调用内核。

这就是我们在向量加法主机代码示例中缺失的部分：定义执行配置，然后使用这里的执行配置和括号内的参数启动内核。现在这是主机代码，但GPU代码呢？在GPU代码中，我们将有类似这样的内容。在这个例子中，向量加法内核非常简单。注意它是一个内核，因为它以 __global__ 限定符开头。这里我们可以看到输入数组 A、B，输出数组 C 和元素总数 n。如果元素总数是块大小的倍数，我们可以编写像这样简单的代码：对于每个单独的线程，每个线程有自己的线程索引，每个块有自己的块索引和块大小（即 blockDim.x）。对于每个单独的线程，我们获得一个索引 i。使用这个 i，每个线程可以访问内存，从 A 和 B 中获取相应的元素，将它们相加，并将结果存储在输出数组 C 中。如果元素数量不是每个块线程数的倍数，那么我们需要检查边界条件。因此，在内核代码中，我们将不得不添加一个 if 语句，检查 i 是否小于 n（即向量中的元素总数）。这是一个简单的步骤。

多维网格与线程块

这就是我们所说的一维网格，其中线程块被认为是单维的。因为它们只有一个维度，所以当我们定义网格维度（网格中的块数）和块维度（每个块的线程数）时，我们使用 .x。同样，我们使用 blockIdx.x 和 threadIdx.x，因为我们使用的是单一维度。我们可以将这种一维线程块和网格也应用于不需要是一维的数据结构。例如，我们可以将不同的线程和线程块分配给图像（一个二维数据结构）进行计算。如图所示，图像的这部分分配给块0，这部分分配给块1，这部分分配给块2，依此类推。然后，对于每个块，我们将有一定数量的线程（在这个简单示例中是四个）。为了访问图像中的任何元素，我们需要使用这个二维数据结构，需要使用块索引、块维度和线程索引。例如，这里的这个元素是图像的第25个元素或像素，我们通过使用块索引（6）、块维度（4，记住在这个简单示例中我们只有四个线程）和线程索引（1，因为这是分配给块6中线程1的像素）来计算并获得了确切的地址。

但线程块也可以是二维或三维的，甚至更多维，因为当我们处理二维或三维数据结构时，使用多维块可能更方便。因此，在二维块的情况下，我们将有 gridDim.x 和 gridDim.y，如你所见，在这个简单示例中，两者的值都是4。这个块也将有一个二维索引，我们有 blockIdx.x 和 blockIdx.y，对于这个特定的块，x是2，y是1。然后，同样，每个线程块中的单个线程也有相同的x和y索引。现在，我们计算图像中特定像素位置的方式是通过使用块和线程在y维度和x维度上的相应索引和维度来获取其行和列。例如，我们在上一张幻灯片中检查的同一个像素，现在我们看到这个像素在第3行第1列。

总结与下节预告

本节课的内容就到这里。你可以通过阅读推荐教材《大规模并行处理器编程》中与本讲对应的章节（第1、2、3章）来了解更多。你也可以通过访问本幻灯片中的链接，重新观看本讲的更长版本。

在下一讲中，我们将继续讨论GPU编程，这次将重点放在GPU的内存层次结构上，探讨在使用CUDA或OpenCL编程时如何处理GPU内部不同的内存空间，以及不同线程块中的线程如何访问全局内存（片外内存），或者同一线程块内的线程如何使用共享内存共享数据。

希望你觉得本讲内容有趣。非常感谢你的关注，期待在下一讲中与你再见。

004：GPU内存层次结构

在本节课中，我们将要学习GPU的内存层次结构，即线程在GPU编程中可以使用的不同内存空间。我们将快速深入细节，但首先，让我们简要回顾一下上一讲关于GPU编程的内容。

课程回顾与背景

上一讲我们介绍了GPU编程的基础。本课程的参考书是《Programming Massively Parallel Processors》。当使用CPU和GPU编程异构系统时，我们通常使用在CPU上运行的CPU线程和在GPU上运行的GPU内核。CPU线程执行在CPU上运行的顺序或适度并行的代码段，而大规模并行的代码段则在GPU上运行。

我们通常从主机（CPU）开始执行，运行一些串行或适度并行的代码，然后通过启动一个GPU内核将计算卸载到GPU。启动GPU内核时，我们需要首先定义执行配置，通常是每个网格的块数（即总块数）和每个块的线程数。GPU内核执行完毕后，控制权返回给主机，主机可能会执行更多串行代码，最终可能启动另一个内核。

我们的第一个例子是向量加法。以下是向量加法的主机端代码：

// 分配GPU内存
cudaMalloc(&d_A, size);
cudaMalloc(&d_B, size);
cudaMalloc(&d_C, size);

// 将数据从CPU内存复制到GPU内存
cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

// 定义执行配置并启动内核
int threadsPerBlock = 512;
int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);

// 将结果从GPU内存复制回CPU内存
cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

// 释放GPU内存
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);

GPU软件层次结构由线程、线程块和网格组成。线程是最低层级，它们被分组到线程块中，线程块被映射到不同的GPU核心上，所有线程块的集合构成了网格。

在NVIDIA最新的H100架构中，在线程块和网格之间引入了一个新的层级：线程块簇。现在，网格由线程块簇组成，而线程块簇又由线程块组成。线程块簇的优点是，簇内的线程块可以通过H100架构引入的新网络直接通信。

GPU内存架构概述

在回顾了GPU软件层次结构之后，让我们来讨论GPU内存。任何NVIDIA GPU的布局都包含GPU核心（也称为流多处理器或SM）。在GPU核心内部，有不同的内存空间，如寄存器、L1缓存或称为共享内存的暂存器内存。

GPU核心可以访问L2缓存，也可以访问片外内存（在H100上是HBM3内存）。每个SM都有一些控制硬件和CUDA核心，用于执行来自不同线程束的指令。从这些CUDA核心，我们可以访问寄存器。

每个单独的线程都可以访问自己的寄存器，访问速度极快（大约一个周期）。线程也可以访问L1数据缓存或附近的共享内存。共享内存是一个软件管理的缓存，可用于在同一线程块内的线程之间共享数据。访问共享内存和L1缓存的延迟大致相同（大约5个周期）。我们还可以访问常量缓存，这是用于常量数据的另一个缓存。

如果需要访问片外内存（称为全局内存），延迟会高得多，大约数百个周期（例如500个周期）。在全局内存和SM之间是L2缓存。

H100 GPU的新特性包括：可以直接从全局内存复制到共享内存，而无需经过GPU核心中的寄存器和通道；在同一线程块簇内的SM之间可以直接通信。

H100 GPU拥有相当大的L2缓存（约50 MB）和全局内存（HBM3内存，80 GB），其全局内存带宽也非常高，达到3 TB/s，大约是上一代A100 GPU带宽的两倍。

GPU内存空间详解

现在我们知道存在不同的内存空间：寄存器、常量缓存、共享内存、L1缓存、L2缓存、全局内存。在我们的程序中，有不同的方式来使用它们。

我们将使用限定符来定义特定值或数组的分配位置：

自动变量会直接进入寄存器，由编译器分配。
在GPU内核中定义的数组，例如 float temp[N];，将是每个线程的私有数组。如果可能，编译器会将其分配在寄存器中；如果寄存器不足（例如N太大），则必须分配在片外内存中。
当我们想在共享内存中分配时，使用 __shared__ 限定符。
当我们想将数据分配在常量内存中以便从常量缓存访问时，使用 __constant__ 限定符。
当我们想在全局内存中分配时，可以使用 __device__ 限定符。

从运行在GPU核心上的GPU线程的角度来看：

它可以访问自己的寄存器供私有使用。
属于同一线程块的线程可以访问共享内存中的相同数据。
不同线程块中的线程可以从外部片外内存（称为全局内存）以及常量内存中访问相同或不同的数据。常量内存也分配在片外内存中，但可以在每个SM内缓存。

共享内存与数据复用

接下来，我们将重点讨论共享内存，看看不同的线程如何共享分配在共享内存中的数据，以及如何使用共享内存进行通信。

当同一线程块内存在数据复用时，共享内存将特别有用。例如，如果我们需要对图像应用一个滤波器，并将图像的每个像素分配给GPU中的一个线程。在应用滤波器时，我们通常需要读取围绕特定像素的邻近像素，并使用这些邻近像素进行计算。

你可以观察到，如果一个像素被分配给一个线程，而另一个邻近像素被分配给同一线程块中的另一个线程，它们将共享一些都需要读取的值。这就是我们所说的数据复用。

在最直接的代码实现中，我们可能会像下面这样编写高斯滤波的代码：

for (int i = -1; i <= 1; i++) {
    for (int j = -1; j <= 1; j++) {
        sum += image[row + i][col + j] * filter[i + 1][j + 1];
    }
}

但这里存在大量的数据复用，因为多个线程会重复读取相同的像素数据。我们可以利用这种数据复用，采用一种称为分块的技术。

在GPU的特定情况下，我们将把这个数据块放入共享内存中，因为共享内存可以被属于同一线程块的所有线程访问。具体做法是：

在共享内存中定义一些空间。
将整个数据块从全局内存加载到共享内存中。
然后线程从共享内存中的数组访问它们计算所需的像素。

在将数据块加载到共享内存和实际计算之间，我们需要使用 __syncthreads() 指令或函数来同步一个块内的所有线程。这是为了避免当同一线程块内的线程访问共享内存或全局内存时，发生读后写或写后读等危险。

这一点很重要，因为线程块通常有超过32个线程（一个线程束的大小是32）。如果一个线程块有64个线程，就意味着有两个线程束。请记住，GPU核心管道通过从不同线程束发出指令来执行。这意味着，即使指令相同，来自不同线程束的指令也在不同的时间执行。因此，如果有两个线程束正在将数据从全局内存加载到共享内存，我们需要确保所有需要的数据（整个数据块）在开始使用该数据块中的数据之前都已加载完毕。这就是为什么在加载之后，我们需要使用 __syncthreads() 进行同步。

案例研究：矩阵乘法

为了更好地理解如何利用数据复用，我们将分析一个非常重要的内核：矩阵乘法。

在矩阵乘法中，我们对两个矩阵A和B进行行乘以列的操作，得到一个输出矩阵C。最简单的并行化方法是将输出矩阵C的每个元素分配给一个线程。

在GPU中，这些线程被分组到线程块中。因此，我们已经将输出矩阵C划分为多个块，分配给不同的线程块。一个线程块负责计算C的一个子块。

一个特定的线程将读取A的整行和B的整列，执行这两个向量的点积运算，并将结果累加到C的相应输出元素中。

以下是基本的矩阵乘法内核代码：

__global__ void matrixMul(float* A, float* B, float* C, int M, int N, int K) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < M && col < N) {
        float sum = 0.0f;
        for (int i = 0; i < K; i++) {
            sum += A[row * K + i] * B[i * N + col];
        }
        C[row * N + col] = sum;
    }
}

但这里的观察是存在大量的数据复用。为什么？因为计算C的某个子块的四个线程都需要访问矩阵A的同一行的所有元素，以及矩阵B的同一列的所有元素。

我们可以通过执行分块矩阵乘法来利用这种数据复用。在分块矩阵乘法中，我们将输入矩阵A和B划分为块，逐步将块加载到共享内存中，然后使用共享内存中的数据执行计算，因为访问共享内存的速度比访问全局内存快得多（大约快两个数量级）。

负责C的某个块的线程块将执行以下操作：

从内存中读取A的一个块和B的一个块。
将它们放入共享内存。
执行部分矩阵乘法以获得中间结果。
对下一个块重复此过程，直到计算出C的该块的所有元素的最终结果。

以下是分块矩阵乘法的代码框架：

__global__ void matrixMulTiled(float* A, float* B, float* C, int M, int N, int K) {
    __shared__ float sA[TILE_SIZE][TILE_SIZE];
    __shared__ float sB[TILE_SIZE][TILE_SIZE];

    int bx = blockIdx.x, by = blockIdx.y;
    int tx = threadIdx.x, ty = threadIdx.y;

    int row = by * TILE_SIZE + ty;
    int col = bx * TILE_SIZE + tx;

    float sum = 0.0f;
    for (int i = 0; i < (K + TILE_SIZE - 1) / TILE_SIZE; i++) {
        // 协作将A和B的块加载到共享内存
        if (row < M && (i * TILE_SIZE + tx) < K) {
            sA[ty][tx] = A[row * K + i * TILE_SIZE + tx];
        } else {
            sA[ty][tx] = 0.0f;
        }
        if (col < N && (i * TILE_SIZE + ty) < K) {
            sB[ty][tx] = B[(i * TILE_SIZE + ty) * N + col];
        } else {
            sB[ty][tx] = 0.0f;
        }
        __syncthreads(); // 等待块加载完成

        for (int j = 0; j < TILE_SIZE; j++) {
            sum += sA[ty][j] * sB[j][tx];
        }
        __syncthreads(); // 确保块内所有线程用完当前共享内存数据后再加载下一块
    }

    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

请注意，在内层for循环之后也需要 __syncthreads()，因为线程需要等待，以确保在将共享内存中的数据替换为下一个块的内容之前，当前块中的数据已被充分使用。

总结

本节课我们一起学习了GPU的内存层次结构。我们回顾了GPU编程的基本流程和软件层次，并详细介绍了GPU中不同的内存空间，包括寄存器、共享内存、常量缓存、L1/L2缓存和全局内存，以及它们的特点和访问方式。

我们重点探讨了共享内存的作用，它作为线程块内线程间通信和协作的高速暂存器，对于存在数据复用的算法（如图像滤波、矩阵乘法）性能提升至关重要。我们通过高斯滤波和矩阵乘法的例子，深入分析了如何利用分块技术将数据加载到共享内存中，并通过 __syncthreads() 进行同步以避免内存访问冲突。

最后，我们简要提及了NVIDIA H100架构在内存层次方面的新特性，如线程块簇和直接内存复制。要了解更多细节，建议阅读《Programming Massively Parallel Processors》第5章，或观看本讲的长版本视频。

005：GPU性能考量 🚀

在本节课中，我们将学习GPU编程中的关键性能考量因素。我们将探讨影响GPU程序性能的主要瓶颈，并介绍优化内存访问、提高线程利用率和重叠计算与通信等核心技巧。

课程回顾

上一节我们介绍了GPU的内存层次结构，包括寄存器、共享内存和全局内存。我们以矩阵乘法为例，学习了如何使用分块技术，将数据从全局内存加载到共享内存中，以减少全局内存访问延迟并提高数据复用率。

核心代码片段展示了分块加载和同步的过程：

// 从全局内存加载A和B的块到共享内存
__shared__ float As[TILE_SIZE][TILE_SIZE];
__shared__ float Bs[TILE_SIZE][TILE_SIZE];
As[ty][tx] = A[row * N + (bx * TILE_SIZE + tx)];
Bs[ty][tx] = B[(by * TILE_SIZE + ty) * N + col];
__syncthreads();

// 使用共享内存中的块进行计算
for (int k = 0; k < TILE_SIZE; ++k) {
    sum += As[ty][k] * Bs[k][tx];
}
__syncthreads();

GPU性能瓶颈

GPU程序主要有两个性能瓶颈，均与数据移动有关：

CPU与GPU之间的数据传输：在拥有独立GPU的系统中，数据需要通过PCIe等总线传输，其带宽有限。
对全局内存的访问：从全局内存读取数据需要数百个时钟周期，延迟很高。

我们的优化目标就是通过充分的计算来分摊数据移动的成本，并优化内存访问模式。

内存访问优化

GPU是细粒度多线程架构，能够通过交错执行不同线程束的指令来隐藏长延迟操作（如全局内存访问）。为了实现有效的延迟隐藏，我们需要关注以下几个关键点。

占用率

占用率 是指活跃线程束数量与GPU流多处理器所能支持的最大线程束数量之比。高占用率意味着有足够多的线程可以交错执行，从而更好地隐藏内存访问等操作的延迟。

占用率受限于GPU核心的资源，例如：

每个线程块的最大线程数
每个线程使用的寄存器数量
每个线程块使用的共享内存大小

虽然高占用率通常是好的，但并非所有内核都需要100%的占用率才能达到最佳性能，因为过高的占用率可能导致对共享内存等资源的争用。

内存布局与合并访问

在C/C++中，二维数组（如矩阵）通常按行主序存储。这意味着矩阵的第一行所有元素在内存中是连续的，接着是第二行的所有元素，依此类推。

为了高效访问全局内存，我们应该尽量让同一个线程束内的线程访问连续的、物理上相邻的内存地址。这种访问模式称为合并内存访问。

以下是不合并与合并访问的对比：

不合并访问：线程束中的线程访问跨距很大的内存地址（例如，访问矩阵中不同行的相同列）。这会导致多个内存事务，降低有效带宽。
合并访问：线程束中的线程访问连续的内存地址（例如，访问矩阵中同一行的连续列）。这通常只需一次内存事务，能最大化内存带宽。

利用共享内存可以改善合并访问。例如，在分块矩阵乘法中，我们首先以合并的方式将数据块从全局内存加载到共享内存，然后在共享内存中进行可能不合并的访问，从而避免了低效的全局内存访问。

DRAM访问与存储体

DRAM芯片内部由存储单元阵列组成。访问数据时，需要先将一整行数据激活到行缓冲区中，然后才能读写该行中的特定列。

行命中：如果后续访问的数据在同一行（即同一DRAM页）中，则可以直接从行缓冲区读取，速度很快。
行冲突：如果访问的数据在不同行，则需要先关闭当前行，再激活新的一行，这个过程延迟很高。

为了进一步提高带宽，现代内存系统采用了多存储体和多通道技术：

多存储体：将内存地址交错分布到多个可以并行访问的存储体上。当一个存储体在激活新行时，其他存储体可能正在服务数据访问，从而重叠了延迟。
多通道：使用多个独立的内存总线和DRAM阵列，进一步增加总带宽。

然而，如果同一个线程束内的多个线程试图访问同一个存储体或同一个通道内的数据，就会发生存储体冲突或通道冲突，导致访问被序列化，降低性能。可以通过地址哈希或填充等技术来缓解冲突。

填充示例：如果线程访问步长为32（存储体数量），所有线程都会访问存储体0。通过在每个数据元素后添加一个空位（填充），可以改变地址映射，将冲突分散到不同存储体。

// 原始数组，可能引起存储体冲突
__shared__ int array[32 * N];
// 填充后的数组，有助于减少冲突
__shared__ int array_padded[33 * N]; // 每行多一个元素

线程束利用与分支发散

GPU以线程束为单位执行指令，同一个线程束内的32个线程是锁步执行的。这意味着所有线程在同一周期执行相同的指令。

分支发散 发生在同一个线程束内的线程根据数据不同走上不同的执行路径时（例如，有的线程执行if分支，有的执行else分支）。GPU会先执行所有走if分支的线程，再执行所有走else分支的线程，导致线程束的有效利用率降低，执行时间翻倍甚至更多。

优化方法：尽可能重新组织数据或计算，让同一个线程束内的线程执行相同的控制流路径。例如，在向量归约的朴素实现中，大量线程会提前闲置，导致严重发散。通过调整线程与数据元素的映射关系，可以让活跃的线程保持连续，减少发散。

原子操作

原子操作 确保对同一内存地址的“读-修改-写”操作是不可分割的，常用于避免数据竞争，例如在并行归约或直方图计算中更新累加器。

原子操作虽然保证了正确性，但可能导致串行化。当多个线程同时原子访问同一地址时，这些操作必须一个一个顺序执行，成为性能瓶颈。

因此，应谨慎使用原子操作，并考虑是否能用不同的算法（如先局部归约再全局归约）来避免或减少其使用。

异步数据传输与流

为了缓解CPU与GPU间数据传输的瓶颈，可以尝试重叠计算与通信。这通过CUDA流来实现。

CUDA流 是一个操作序列（如内存拷贝、内核执行），其中的操作按顺序发布，但GPU可以并行执行不同流中的操作。

基本思想：将输入数据分成若干独立的块。我们可以在流0中传输数据块0，然后在流0中启动处理数据块0的内核。与此同时，可以在流1中开始传输数据块1。这样，数据块的传输和计算就在时间上重叠了起来。

性能提升估算：

若内核计算是主要耗时部分，重叠通信可隐藏部分通信时间。
若数据传输是主要耗时部分，重叠计算可隐藏部分计算时间。

这种方法在视频处理等应用中非常有效，可以独立处理不同的视频帧块。

总结

本节课我们一起学习了GPU编程的核心性能考量：

数据移动是关键瓶颈：需优化CPU-GPU数据传输和全局内存访问。
利用延迟隐藏：通过保持足够的占用率，让GPU的细粒度多线程架构能够隐藏内存访问延迟。
优化内存访问模式：确保合并内存访问，并利用共享内存进行数据复用和访问模式转换。
最大化硬件利用率：避免线程束分支发散，并谨慎使用可能导致串行化的原子操作。
重叠计算与通信：使用异步传输和CUDA流，将数据分块处理，以隐藏数据传输延迟。

掌握这些性能考量原则，是编写高效GPU程序的基础。在后续课程中，我们将深入更多并行模式及其优化实践。

推荐阅读：《Programming Massively Parallel Processors》第6章及第20章第5节。也可观看本讲加长版视频以获取更多细节。

006：归约

在本节课中，我们将要学习并行编程中的一个核心模式：归约。我们将从归约的基本概念开始，逐步深入到其在GPU上的高效实现，包括如何避免性能瓶颈、利用共享内存和Warp Shuffle指令，以及应用线程协作等优化技术。

概述

归约是一种将一组值合并为单个值的操作，例如求和、求积、求最大值或最小值。它是许多并行程序（如MapReduce模型）的关键原语，也广泛应用于神经网络和机器学习工作负载中。本节课将介绍归约的串行实现、经典的树形归约并行算法，并重点讲解如何在GPU上高效地实现归约，同时考虑内存访问、线程利用率和同步等关键性能因素。

回顾：GPU性能考量

在深入归约之前，让我们简要回顾上一讲关于GPU性能考量的内容。优化GPU代码时，内存访问模式是关键。为了充分利用带宽，我们需要实现合并内存访问，即连续的线程访问连续或相邻的内存地址。当所有线程访问的数据位于同一缓存行或几个连续的缓存行时，可以实现最佳带宽利用。

共享内存是改善合并访问的有效工具。例如，在处理矩阵时，我们可以将数据块加载到共享内存中，以合并的方式访问全局内存，然后从共享内存中读取数据，从而避免对全局内存的非合并访问。

另一个重要考量是线程利用率。当同一Warp内的线程遵循不同的控制流路径时，会发生Warp内部分歧，导致指令执行效率低下。通过重新分配计算任务，可以避免这种分歧，提高Warp利用率。这一点对于归约操作尤为重要。

什么是归约？🔍

归约是一种操作，它将一组值减少为单个值。常见的归约操作包括求和、求积、求最大值和求最小值。这些操作具有结合律、交换律和单位元等数学性质。例如，加法的单位元是0。

归约是并行计算中的关键原语，广泛应用于各种并行程序。MapReduce编程模型就是一个很好的例子。此外，归约也常见于神经网络和其他现代机器学习工作负载中。

串行归约

串行归约的实现非常简单。我们遍历整个输入数组，使用一个累加器逐步合并每个元素。

以下是串行求和的伪代码：

sum = 0; // 单位元
for (i = 0; i < n; i++) {
    sum = sum + array[i];
}

对于一个包含 n 个元素的数组，串行归约需要 n 次迭代。然而，归约操作中的许多步骤是独立的，这为并行化提供了机会。

并行树形归约 🌳

经典的并行归约方法是树形归约。其核心思想是将数组元素配对，由不同的线程或工作单元计算部分和，然后递归地归约这些部分和，直到得到最终结果。

对于一个包含 n 个元素的数组，树形归约只需要 log₂(n) 次迭代，远少于串行归约的 n 次迭代。

在GPU上映射树形归约时，我们通常：

将输入数组划分为多个数据块，分配给不同的线程块。
在每个线程块内部，进一步将数据块分配给不同的Warp。
线程在Warp内执行部分归约，将中间结果存储在共享内存或寄存器中。
由于同一线程块内的Warp并非严格同步执行，需要使用 __syncthreads() 进行块内同步，确保所有线程能看到最新的中间结果。
每个线程块产生一个部分结果。
最后，需要对这些跨线程块的部分结果进行最终归约。

跨线程块同步

在获得每个线程块的部分结果后，需要进行跨线程块同步以完成最终归约。主要有三种方法：

以下是几种实现跨块同步的方法：

CPU端归约：将所有部分结果复制到CPU内存，在CPU上完成最终归约。
启动新内核：在GPU上启动一个新的归约内核，以上一步的部分结果为输入，继续执行归约直到得到单个值。
原子操作：每个线程块使用一个线程，通过原子操作将其部分结果累加到全局内存中的一个原子变量上。

提高Warp利用率

在树形归约的初始实现中，我们可能让偶数编号的线程负责相加连续的元素对。这会导致严重的Warp利用率低下问题，因为在第一次迭代后，就有一半的线程（奇数编号线程）闲置，造成Warp内部分歧。

为了解决这个问题，我们可以采用避免分歧的映射方式。让连续的线程负责处理连续的数据对。这样，在归约的早期阶段，活跃的线程会集中在少数Warp中，而不是分散在所有Warp里，从而最大化Warp的利用率。

利用Warp Shuffle指令 🚀

在归约的最后阶段，当活跃线程仅剩一个Warp时，我们可以使用 Warp Shuffle 指令来进一步提升性能。这些指令允许同一Warp内的线程直接交换寄存器中的数据，无需通过共享内存。

Warp Shuffle指令比先写入共享内存再同步读取的方式更快，延迟更低。它们自Kepler架构引入，有多种变体，例如：

__shfl_sync: 从指定通道（线程）直接复制数据。
__shfl_up_sync: 从ID较低的通道复制数据。
__shfl_down_sync: 从ID较高的通道复制数据。
__shfl_xor_sync: 基于位异或操作交换数据。

在归约中，我们可以在最后一个Warp内使用 __shfl_down_sync 指令来高效地累加部分和。

归约代码示例（含Warp Shuffle）

以下是一个在GPU上实现并行归约的核函数示例，它结合了共享内存树形归约和最后的Warp Shuffle优化：

__global__ void reduction_kernel(int* input, int* partial_sums, int n) {
    __shared__ int s_data[BLOCK_SIZE];
    int tid = threadIdx.x;
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    // 每个线程加载两个元素到共享内存（假设grid size是数组大小的一半）
    if (i * 2 < n) s_data[tid] = input[i*2] + input[i*2+1];
    else s_data[tid] = 0;
    __syncthreads();

    // 在共享内存中进行树形归约
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            s_data[tid] += s_data[tid + s];
        }
        __syncthreads();
    }

    // 最后一个Warp使用Warp Shuffle完成最终归约
    if (tid < 32) {
        int val = s_data[tid];
        for (int offset = 16; offset > 0; offset >>= 1) {
            val += __shfl_down_sync(0xffffffff, val, offset);
        }
        if (tid == 0) {
            partial_sums[blockIdx.x] = val; // 存储该线程块的部分和
        }
    }
}

注意：这是一个简化示例，实际应用中需要考虑数组大小与线程网格的匹配、边界条件等。

现代GPU的Warp级归约指令

自Compute Capability 8.0（如NVIDIA Ampere架构）起，引入了更高效的 Warp级归约内在函数，例如 __reduce_add_sync。这些指令可以直接在硬件上执行Warp内的归约操作，进一步简化代码并提升性能。我们可以用一行指令替代之前需要循环和Shuffle指令的Warp内归约步骤。

线程协作优化 🧵

另一个重要的优化技术是线程协作。其核心思想是让每个线程处理输入数组中的多个元素，而不是只处理一两个。这样，我们可以用更少的线程块覆盖整个数组，减少需要启动的线程块总数和相应的同步开销。

实现方式是：每个线程使用一个循环，以网格总线程数为步长，遍历并累加分配给它的那些元素。每个线程先计算一个本地和，然后再与其他线程在共享内存中汇合，进行后续的树形归约。这种方法在数据量很大时尤其有效，因为它减少了线程块的数量和 __syncthreads() 的调用次数。

总结

本节课我们一起学习了并行编程中的核心模式——归约。我们从归约的基本定义和串行实现开始，深入探讨了在GPU上实现高效归约的多种技术：

理解了树形归约这一并行算法如何将复杂度从O(n)降低到O(log n)。
认识到内存访问模式和Warp利用率对GPU性能至关重要，并学习了如何通过重新映射计算来避免Warp内部分歧。
掌握了使用共享内存作为中间存储来进行块内归约，并使用 __syncthreads() 进行同步。
学习了利用 Warp Shuffle 指令在Warp内高效交换数据，避免额外的共享内存访问。
了解了跨线程块同步的几种策略：CPU归约、启动新内核或使用原子操作。
探讨了线程协作优化，通过让每个线程处理更多工作来减少线程块数量和同步开销。
提及了现代GPU提供的专用Warp归约指令（如 __reduce_add_sync）可进一步简化代码并提升性能。

归约的实现拥有巨大的优化空间，不同的库和编译器会自动探索各种代码变体以适应特定硬件架构。掌握这些基础技术，是编写高性能GPU程序的关键一步。

007：直方图

在本节课中，我们将学习一个新的并行模式：直方图计算。我们将从回顾上一节的内容开始，然后介绍直方图的基本概念、其串行与并行实现，并重点讨论如何使用原子操作来确保并行计算的正确性。最后，我们将探讨一些优化技术，如私有化和合并原子操作，以提高直方图计算的性能。

回顾：归约操作

上一节我们介绍了归约操作。归约是一种将一组值减少为单个值的操作，它具有结合律、交换律和单位元等特性。

我们解释了如何在GPU上实现基于树的归约。通常的做法是将输入数组划分为不同的块，分配给不同的线程块。在每个线程块内部，再将数据分配给不同的线程束。每个线程束会得到一个部分归约结果，这些结果会进一步归约为每个线程块的部分结果。经过一些线程块间的同步后，最终得到完整的归约结果。

我们还讨论了避免线程束分化的映射方法，以确保活跃线程始终属于同一个或少数几个线程束，从而提高计算核心利用率和最终性能。

此外，我们提到了原子操作可以作为实现归约的一种方式。原子操作能保证读-修改-写操作的原子性，可用于共享内存或全局内存。它们对于防止多个线程更新同一内存位置时的数据竞争非常有用。

直方图简介

直方图是一种数据结构，用于从大型数据集中降维并提取显著的特征和模式。它广泛应用于图像处理和欺诈检测等领域。

在基础直方图中，我们将输入值的可能范围划分为若干个“桶”。每个桶关联一个计数器。对于输入中的每个元素，我们检查其值，确定它属于哪个桶，然后递增该桶对应的计数器。

串行直方图计算

串行直方图计算非常简单。我们只需逐个遍历输入元素，并更新对应的直方图桶计数器。

例如，给定一个输入数组和一个直方图，单个线程会迭代地读取每个输入元素，并递增相应的计数器。

以下是串行直方图计算的伪代码：

for (int i = 0; i < input_size; i++) {
    value = input[i];
    histogram[value]++;
}

并行直方图计算

为了提高效率，我们尝试并行实现直方图计算。多个线程可以同时处理输入数组的不同部分。

一个直观的并行内核可能如下所示：

__global__ void histogram_kernel(int *input, int *histogram, int size) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    int idx = tid;
    while (idx < size) {
        int value = input[idx];
        histogram[value]++; // 这里存在数据竞争风险
        idx += stride;
    }
}

在这个内核中，每个线程处理输入数组中相隔stride距离的元素。线程读取输入值后，尝试直接递增直方图中对应的桶。

然而，这个实现存在一个问题。当多个线程同时尝试更新同一个直方图桶时，会发生数据竞争，导致结果错误。

使用原子操作

为了解决数据竞争问题，我们需要使用原子操作来确保对直方图桶的更新是原子的。

以下是使用原子操作的正确并行直方图内核：

__global__ void histogram_kernel_atomic(int *input, int *histogram, int size) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    int idx = tid;
    while (idx < size) {
        int value = input[idx];
        atomicAdd(&histogram[value], 1); // 使用原子加操作
        idx += stride;
    }
}

通过使用atomicAdd，我们保证了即使多个线程同时更新同一个桶，计数器也能被正确递增。但代价是，当发生原子冲突时，这些操作会被序列化，可能增加延迟。

优化：私有化

在处理自然图像等数据时，像素值分布可能具有空间相似性，导致大量线程更新少数几个桶，从而引发频繁的原子冲突。为了优化性能，我们可以采用私有化技术。

私有化是一种优化技术，它维护输出的多个私有副本，待计算完成后再更新全局副本。此技术要求操作满足结合律和交换律，而直方图使用的加法操作正好满足。

私有化的优势在于减少了全局副本的争用。如果输出足够小，我们可以将其放在共享内存等快速存储器中，从而降低原子操作的延迟。

以下是结合了私有化和合并的优化思路：

将图像划分为多个块，每个块分配给一个线程块。
每个线程块在共享内存中维护自己的子直方图。
线程块处理完分配的数据块后，得到完整的子直方图。
最后，通过一次并行归约，将所有子直方图合并到全局内存中的最终直方图。

结合合并技术，可以让每个线程块处理多个图像块，从而减少需要合并的子直方图数量，并减少直方图初始化的开销。

以下是使用私有化和合并的并行直方图内核伪代码结构：

__global__ void histogram_privatized(int *input, int *global_hist, int size) {
    __shared__ int s_hist[NUM_BINS]; // 每个块的私有直方图
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    // 初始化共享内存中的直方图
    if (threadIdx.x < NUM_BINS) s_hist[threadIdx.x] = 0;
    __syncthreads();

    // 计算每个块的子直方图
    int idx = tid;
    while (idx < size) {
        int value = input[idx];
        atomicAdd(&s_hist[value], 1); // 在共享内存中使用原子操作
        idx += blockDim.x * gridDim.x;
    }
    __syncthreads();

    // 将子直方图合并到全局内存
    if (threadIdx.x < NUM_BINS) {
        atomicAdd(&global_hist[threadIdx.x], s_hist[threadIdx.x]);
    }
}

在这个内核中，大部分时间花费在计算共享内存中的子直方图上。虽然最后一步在全局内存中使用原子操作成本较高，但由于只执行一次，对整体性能影响很小。

进一步优化：合并原子操作

我们还可以利用线程束同步原语进一步优化原子操作。例如，使用__shfl、__ballot等函数。

合并原子操作的基本思想是：当同一个线程束中的多个线程需要更新同一个内存位置时，先在线程束内部进行一次局部归约，得到一个总和，然后只让一个线程执行一次原子操作来更新目标内存位置。这样可以显著减少原子操作的次数和冲突。

以下是合并原子操作的简化思想（具体实现更复杂）：

识别出线程束中所有需要更新同一桶的线程。
在这些线程内部，使用线程束洗牌函数进行局部求和归约。
指定一个线程（如lane 0）代表整个线程束，执行一次原子加操作，将归约后的总和写入目标桶。

这种方法可以进一步优化原子操作的使用效率。

总结

本节课我们一起学习了直方图这一并行模式。我们首先回顾了归约操作和原子操作的基础知识。然后，我们介绍了直方图的概念及其串行实现。

在并行实现部分，我们指出了直接更新全局直方图会导致数据竞争，因此必须使用原子操作来保证正确性。接着，我们探讨了通过私有化技术来减少原子冲突的优化方法，即在共享内存中维护每个线程块的子直方图，最后再合并到全局内存。

我们还简要介绍了利用线程束同步原语进行合并原子操作的进一步优化思路。

直方图计算是许多应用中的基础操作，理解其并行实现和优化技巧对于编写高效的GPU程序至关重要。

008：卷积

在本节课中，我们将学习一个新的并行模式：卷积。卷积在信号处理、图像处理和机器学习等领域有广泛应用。我们将从一维卷积的基本概念开始，逐步深入到如何在GPU上高效实现卷积，包括利用常量内存、分块技术以及将卷积转换为矩阵乘法等优化策略。

卷积简介

上一节我们讨论了直方图计算，本节中我们来看看卷积操作。卷积是一种广泛应用于信号处理、图像处理、视频处理和计算机视觉的运算。其基本思想是，将一个滤波器（或称为掩码）应用于输入的每个元素（如信号、图像或视频帧），通过计算该元素及其邻域元素的加权和，得到一个新的输出值。这在图像处理中可用于平滑、锐化、模糊、边缘检测或去噪。如今，卷积在机器学习和人工智能中也至关重要，特别是作为卷积神经网络中的卷积层。

一维卷积示例

让我们从一个常用于音频处理的一维卷积示例开始。假设我们有一个由连续元素组成的音频输入，以及一个通常包含奇数个元素的掩码（本例中为5个元素）。

为了计算输出中的某个元素（例如 P2），我们需要将掩码中心对准输入中对应的元素。然后执行逐元素乘法（部分乘积），最后将这些部分乘积求和。得到的值（本例中为57）存储在输出数组的相应位置。

核心计算过程可以用以下伪代码描述：

对于输出中的每个元素 i：
    sum = 0
    对于掩码中的每个偏移量 j：
        input_index = i + j - (mask_size / 2)
        sum += input[input_index] * mask[j]
    output[i] = sum

边界条件处理

计算卷积时，必须注意边界条件。当计算靠近输入数组边界的输出元素时，掩码可能会超出输入数组的范围。我们需要处理这些“幽灵”元素。

以下是几种常见的边界处理策略：

零填充：假设边界外的元素值为0。
边界值复制：使用边界上的值来填充外部元素。
其他方法：如环绕或镜像。

在零填充的例子中，超出边界的输入元素被视为0，因此对应的部分乘积为0，最终输出值为38。

基础GPU卷积内核

在GPU上实现基础卷积时，可以为输出数组的每个元素分配一个线程。每个线程需要计算其对应输出元素的值。

以下是实现的关键步骤：

计算线程对应的输出索引。
确定需要访问的输入元素的起始地址。
循环遍历掩码的所有元素。
在循环中，检查当前输入索引是否越界，并根据边界处理策略（如零填充）获取输入值。
执行输入值与对应掩码值的乘法，并累加结果。
将最终累加值写入输出数组。

利用常量内存存储掩码

GPU架构中的常量内存是一种缓存型内存，每个GPU核心都有专用的常量缓存。由于卷积操作中的掩码是常数、数据量小且被所有线程访问，因此非常适合存放在常量内存中。当同一个线程束中的所有线程访问相同的常量值时，效率尤其高。

要使用常量内存，需要执行以下步骤：

使用 __constant__ 限定符在全局作用域声明掩码变量。
在主机端，使用 cudaMemcpyToSymbol 函数将掩码数据从主机内存复制到设备的常量内存。

示例代码：

__constant__ float MASK[MASK_SIZE];

// 主机端代码
cudaMemcpyToSymbol(MASK, host_mask, MASK_SIZE * sizeof(float));

使用分块技术优化卷积

之前课程中讨论过的分块（或平铺）技术，同样可以应用于卷积以利用数据局部性。对于一维卷积，我们可以将输入和输出数组划分为多个块，由不同的线程块处理。

每个线程块将其负责的输入数据块（包括计算边界输出所需的额外“光晕”区域）加载到共享内存中。这样，线程在计算卷积时可以从快速的共享内存中读取数据，减少对全局内存的访问。

加载数据到共享内存通常分为三步：

加载左光晕：加载当前数据块左侧边界外的元素。
加载内部元素：加载数据块本身的元素。
加载右光晕：加载当前数据块右侧边界外的元素。

加载过程中需要进行地址计算和边界检查，以确保不会访问无效的输入区域。所有数据加载完毕后，需要调用 __syncthreads() 来确保共享内存中的数据对块内所有线程可见。

二维卷积与机器学习中的应用

卷积在机器学习中，特别是作为卷积神经网络中的层，非常有用。与全连接层相比，卷积层具有只使用局部权重、参数数量相对较少的优势。这使得权重可以轻松放入GPU的片上内存（如共享内存或常量内存），从而让分块技术大显身手。

二维卷积的操作与一维类似，只是掩码在二维平面上移动。在CNN中，输入可能是特征图，通过对局部窗口内的输入值与滤波器权重进行乘加运算，得到输出特征图上的一个点。

一个简化的二维卷积层前向传播内核可能包含以下步骤：

计算当前线程对应的输出特征图位置坐标。
循环遍历所有输入通道。
在每个通道内，循环遍历滤波器（如KxK大小）的每个元素。
执行乘积累加操作。
将最终结果写入输出特征图。

将卷积转换为矩阵乘法

卷积运算可以通过“im2col”（图像到列）等方法转换为矩阵乘法。这种转换非常有益，因为矩阵乘法是计算密集型操作，而GPU（特别是其张量核心）对此类运算进行了深度优化。

基本思想是将输入特征图的局部感受野展开成矩阵的一列，并将多个卷积滤波器展开成矩阵的行。这样，一次大型矩阵乘法就能完成整个卷积层的计算。这种方法不仅提升了性能，还有助于在CNN的不同层之间保持稳定的并行度。

在更现代的GPU架构（如Volta及以后的架构）中，可以使用专门的张量核心来加速这种矩阵乘法，从而获得极高的吞吐量。

总结

本节课我们一起学习了卷积这一重要的并行模式。我们从一维卷积的基本原理和边界条件处理开始，探讨了如何在GPU上实现基础卷积内核。接着，我们介绍了利用常量内存存储掩码以提升访问效率，以及使用共享内存分块技术来优化数据重用。我们还了解了卷积在机器学习中的应用，并学习了通过将卷积运算转换为矩阵乘法来充分利用GPU计算能力的高级优化技术。掌握这些概念和技术，对于在GPU上高效实现图像处理、信号处理和深度学习算法至关重要。

009：高级矩阵乘法分块技术 🧩

在本节课中，我们将学习如何为矩阵乘法应用高级分块技术，特别是结合使用寄存器和共享内存，以实现更高的数据复用和计算性能。

上一节我们介绍了卷积操作及其与矩阵乘法的关系。本节中，我们将深入探讨如何将卷积层高效地转化为矩阵乘法，并应用更精细的分块策略。

概述：从卷积到矩阵乘法

卷积是信号处理、图像处理和卷积神经网络中的核心操作。为了在GPU上高效执行卷积，一种常见策略是将其转化为矩阵乘法，从而利用GPU擅长矩阵运算的特性。

在卷积神经网络的前向传播中，输入特征图与卷积滤波器进行卷积运算，生成输出特征图。通过将输入特征图和滤波器“展开”成矩阵形式，我们可以将卷积层计算转化为一个大型矩阵乘法问题。

矩阵乘法的分块策略回顾

为了高效地在GPU上执行矩阵乘法，我们通常使用分块技术。基本思想是将输出矩阵划分为多个小块，每个线程块负责计算一个输出块。

以下是实现分块矩阵乘法的关键步骤：

从全局内存加载输入矩阵A和B的相应数据块到共享内存。
线程协作，利用共享内存中的数据执行部分矩阵乘法计算。
将部分结果累加到输出矩阵C的对应块中。

这种方法的核心优势在于数据复用：加载到共享内存中的数据可以被线程块内的多个线程多次使用，减少了对全局内存的访问。

联合寄存器与共享内存分块 🚀

为了进一步提升性能，我们可以采用更高级的联合寄存器与共享内存分块技术。这种方法的核心思想是利用不同层级存储器的特性：

寄存器：访问速度极快，但为每个线程私有，需要线程协作才能实现数据共享。
共享内存：访问速度比寄存器慢，但比全局内存快得多，并且对同一线程块内的所有线程可见。

在这种策略中，我们为两个输入矩阵使用不同大小的分块：

将矩阵M的一个垂直分块存储在寄存器中。
将矩阵N的一个水平分块存储在共享内存中。
输出矩阵P的对应分块也存储在寄存器中。

关键计算模式

假设我们为矩阵M设置的分块宽度为 Tile_M，为矩阵N设置的分块宽度为 Tile_N。

每个线程块拥有 Tile_M 个线程。
每个线程负责计算输出分块P中的一行，即 Tile_N 个元素。这实现了线程协作，协作因子为 Tile_N。
从矩阵M加载到寄存器的每个数据元素，可以被重用 Tile_N 次。
从矩阵N加载到共享内存的每个数据元素，可以被重用 Tile_M 次。

优化：平衡线程负载

在基础方案中，可能只有少数线程参与将N的分块加载到共享内存，导致线程束分化，降低效率。

一个更平衡的优化方案是让线程块内的所有线程协作加载一个矩形分块的数据到共享内存。这个矩形分块的尺寸是 Tile_N x K，其中 K = Tile_M / Tile_N。这样，所有线程都参与加载工作，消除了线程束分化。

同时，每个线程可以加载K个矩阵M的元素到其私有寄存器中。此时，每个线程在每个迭代中计算K个步骤，最终完成 Tile_N x Tile_M 个输出元素的计算。

实际应用中的参数选择

在实际内核实现中，参数选择需考虑硬件限制：

Tile_N 的典型值为16。其上限受限于可用寄存器数量，因为输出元素也需存储在寄存器中。
Tile_M 至少为64，即每个线程块至少64个线程，对应输出分块的行数。
变量K受限于每个线程可用的寄存器数量，因为它决定了每个线程能加载的M矩阵元素数量。

总结

本节课我们一起学习了用于矩阵乘法的高级分块技术。我们首先回顾了将卷积层转化为矩阵乘法的动机，然后深入探讨了如何联合使用寄存器和共享内存进行分块计算。这种技术通过线程协作和数据复用，显著提升了计算效率。关键点在于为不同矩阵设计不同大小的分块，并利用所有线程协作加载数据，以实现负载均衡和高内存带宽利用率。

进一步学习资源：

欲深入了解此技术，可参考《大规模并行处理器编程》一书的第7、8、16章。
可观看本课程以往学期的加长版讲座视频获取更多细节。

010：前缀和（扫描）🔢

在本节课中，我们将学习另一个重要的并行模式：前缀和，也称为扫描操作。我们将介绍其基本概念、两种主要类型（包含式与排除式），并深入探讨如何在GPU上高效地实现分层扫描算法。

概述

前缀和（扫描）是一种并行原语，它接收一个输入数组和一个满足结合律的运算符（例如加法、乘法、最大值、最小值），并返回一个输出数组。输出数组的每个元素是输入数组从起始位置到当前（或前一个）位置所有元素递归应用该运算符的结果。它在许多并行算法中作为基础构建块，例如流压缩、分区、排序等。

前缀和的定义与类型

前缀和操作有两种主要形式：排除式扫描和包含式扫描。

排除式扫描：对于输出数组的每个索引 i，其值 Y[i] 是输入数组从 0 到 i-1 所有元素运算的结果。公式表示为：
Y[i] = input[0] op input[1] op ... op input[i-1]
其中 op 是满足结合律的运算符。输出数组的第一个元素通常被初始化为该运算符的单位元（例如，加法为0，乘法为1）。
包含式扫描：对于输出数组的每个索引 i，其值 Y[i] 是输入数组从 0 到 i 所有元素运算的结果。公式表示为：
Y[i] = input[0] op input[1] op ... op input[i]
输出数组的第一个元素就是输入数组的第一个元素。

让我们看一个具体的例子。假设输入数组为 [3, 1, 7, 0, 4, 1, 6, 3]，运算符为加法。

排除式扫描输出：[0, 3, 4, 11, 11, 15, 16, 22]
（第一个元素是单位元0，后续每个元素是前一个输出元素加上当前输入元素）
包含式扫描输出：[3, 4, 11, 11, 15, 16, 22, 25]
（第一个元素是输入的第一个元素3，后续每个元素是前一个输出元素加上当前输入元素）

在GPU上实现分层扫描

在具有多层次并行性的机器（如GPU）上实现扫描，我们采用分层方法。GPU的计算由线程块组织，因此我们可以将输入数组划分给不同的线程块处理。

上一节我们介绍了前缀和的基本概念，本节中我们来看看如何在GPU架构上具体实现它。整体流程分为三个主要步骤：

每个线程块内扫描：每个线程块负责处理输入数组的一个连续块，并计算该块内的包含式扫描，得到中间结果。
扫描部分和：收集每个线程块中间结果的最后一个元素（即该块的部分和），对这些部分和数组执行一次扫描操作。
添加偏移量：将第二步计算得到的对应偏移量，加到每个线程块的中间结果上，得到最终的全局扫描结果。

这种方法需要全局同步。在第一步完成后，所有线程块的部分和必须就绪，才能开始第二步。在GPU上，这通常通过内核终止（结束一个内核，再启动新内核）或使用原子操作来实现。

线程块内的高效扫描：Kogge-Stone算法

现在，我们聚焦于如何在一个线程块内高效地实现包含式扫描。我们将介绍Kogge-Stone算法。

该算法为每个输入元素分配一个线程。它通过迭代方式工作，每次迭代中，每个线程将其当前值与一定“跨度”之前的另一个值进行运算。跨度在每次迭代中翻倍（1， 2， 4， ...），直到覆盖整个数组。

以下是该算法的核心代码逻辑（以加法为例）：

// 假设每个线程处理一个元素，`input` 和 `output` 在全局或共享内存中
int idx = threadIdx.x; // 线程索引
int val = input[idx];  // 当前线程负责的输入值

// 初始化输出为输入值（对于包含式扫描）
output[idx] = val;
__syncthreads(); // 块内同步

// 迭代 log2(blockDim.x) 次
for (int stride = 1; stride < blockDim.x; stride *= 2) {
    int element_to_add;
    if (idx >= stride) {
        element_to_add = output[idx - stride];
    }
    __syncthreads(); // 确保所有线程已完成读取
    if (idx >= stride) {
        output[idx] += element_to_add; // 执行加法操作
    }
    __syncthreads(); // 确保所有线程已完成写入，再进行下一轮
}

为了提升性能，我们可以进行两项优化：

使用共享内存：将数据缓存到速度更快的共享内存中，利用数据重用。
双缓冲技术：使用两个共享内存缓冲区，在每次迭代后交换输入和输出角色，从而消除一次 __syncthreads() 调用。

进一步的优化：Warp级扫描与分层融合

我们可以在线程块内部继续应用分层思想。一个线程块包含多个Warp（通常是32个线程）。我们可以：

每个Warp内扫描：首先，每个Warp使用其内部的线程独立计算一个Warp级的包含式扫描。这里可以使用洗牌指令，这是一种允许同一Warp内线程快速交换寄存器值的硬件原语，速度极快。

// 简化的Warp内扫描函数（加法）
__device__ int warpScan(int val) {
    for (int offset = 1; offset < 32; offset *= 2) {
        int y = __shfl_up_sync(0xffffffff, val, offset); // 洗牌指令
        if (lane_id >= offset) val += y; // lane_id 是线程在Warp内的索引
    }
    return val;
}

扫描Warp的部分和：收集每个Warp扫描结果的最后一个值（即该Warp的部分和），由一个Warp（例如第一个Warp）对这些部分和执行一次扫描。
添加Warp偏移量：将第二步得到的偏移量加到对应Warp的每个元素上，完成线程块内的最终扫描。

这种方法将计算进一步细化，并利用Warp内的高效同步和通信机制。

算法变体与性能分析

根据不同的实现策略，扫描算法在全局内存访问量上有显著差异：

Scan-Scan-Add：先进行块扫描（写回全局内存），再扫描部分和，最后添加偏移。总访问量约为 4N（N为元素总数）。
Reduce-Scan-Scan：先进行块归约（只计算部分和，不写完整中间结果），再扫描部分和，最后在一个内核中完成块扫描并添加偏移。总访问量约为 3N，对于大数组通常性能更好。

使用原子操作进行块间同步

为了进一步减少全局内存访问和内核启动开销，我们可以使用相邻块同步技术。其核心思想是使用全局内存中的原子操作和标志数组，让线程块动态协调，隐式地完成部分和的扫描与传播。

每个线程块计算完自己的部分和后，将其写入一个全局数组。
该线程块的一个领导线程（如0号线程）会轮询检查前一个线程块设置的“完成标志”。
一旦前一个块标志就绪，当前块就读取前一个块的部分和结果，将其与自己的部分和相加，更新自己的全局部分和值。
当前块设置自己的“完成标志”，通知下一个块可以开始工作。

由于GPU调度线程块的顺序不确定，我们需要使用一个全局原子计数器为线程块分配“虚拟ID”，以确保同步顺序的正确性。这种方法可以将总全局内存访问量降低到约 2N，从而获得最佳性能。

前缀和的应用

正如开头所提，前缀和是一个极其有用的并行原语，是许多并行算法的构建基石。以下是部分应用：

流压缩：从一个数组中筛选出满足特定条件的元素。
分区：根据条件将数组元素重新排列。
选择/筛选：选出排名第k的元素。
去重：在排序数组中移除重复项。
基数排序：一种高效的整数排序算法。

在我们的研究中，通过使用基于相邻块同步的高效扫描实现，我们开发了一套原地数据滑动算法（如流压缩、分区等），其输出直接覆盖输入内存空间，节省了大量内存。实验表明，该实现相比标准库（如Thrust）获得了约3倍的加速，几乎达到了内存带宽的理论上限。

总结

本节课中，我们一起学习了并行编程中的核心模式——前缀和（扫描）。我们从其定义和两种类型（包含式与排除式）开始，详细探讨了在GPU上实现高效扫描的分层方法，包括线程块内的Kogge-Stone算法和Warp级扫描。我们还分析了不同算法变体（Scan-Scan-Add， Reduce-Scan-Scan）的性能，并介绍了使用原子操作进行相邻块同步以最大化性能的高级技术。最后，我们了解了扫描在流压缩、分区等多种并行算法中的关键作用。掌握前缀和的实现，是理解许多复杂并行算法的基础。

011：稀疏矩阵

在本节课中，我们将要学习稀疏矩阵及其压缩格式，并以稀疏矩阵向量乘法（SpMV）为例算法进行讲解。

首先，让我们快速回顾一下本课程中已经介绍过的其他并行模式。

我们以归约操作开始，这是一种将一组值缩减为单个值的操作。归约需要具备三个主要属性：结合律、交换律和单位元。它是并行计算中的关键原语。我们讨论了如何在GPU上实现无分支映射以最大化利用率。

接下来是直方图计算，这是一种常用于降低数据维度和提取显著特征的操作。在直方图计算中，我们通常逐个读取输入元素，并递增直方图中的计数器或“桶”。问题在于，当多个线程更新直方图时，需要使用原子操作来确保结果正确，而原子操作会串行化执行。因此，我们需要找到减轻这种串行化的方法，其中一种方法称为“私有化”，其基本思想是在共享内存中放置多个子直方图，每个子直方图由不同的线程块创建，当所有子直方图计算完成后，再合并成最终的直方图。

在后续课程中，我们介绍了卷积。在卷积中，我们对输入的每个元素应用一个滤波器或掩码，并通过加权和来计算输出元素。卷积在信号处理、图像处理、视频处理、计算机视觉以及当今的机器学习和人工智能（例如卷积神经网络）中应用非常广泛。我们讨论了如何执行一维卷积，例如应用一个大小为5的掩码到输入的每个元素，计算部分乘积后进行最终归约。我们还讨论了如何进行二维卷积。在另一讲中，我们探讨了如何将卷积层实现为矩阵乘法，这基本上要求我们将输入特征图和卷积滤波器的矩阵展开并排列成矩阵形式，从而像执行矩阵乘法一样进行计算，这对于GPU来说是非常适合的操作。我们确实讨论了更高级的优化技术，如寄存器平铺和共享内存平铺，以及如何结合它们来高效地实现矩阵乘法。

在上一讲中，我们讨论了前缀和（扫描）操作。这是一种操作，我们取一个输入和一个结合性运算符，通过递归地对输入数组的元素应用该运算符来计算输出元素。请记住，有独占扫描和包含扫描之分。我们讨论了如何实现一种称为包含扫描的层次结构，该结构可以适应具有两级或更多级并行性的系统，例如在GPU中，我们有线程块、线程束和线程。我们首先要做的是在线程块之间划分输入，每个线程块将执行一个块扫描。但在线程块内部，我们也会使用线程束和单个线程。在这种层次化扫描中，有不同的步骤，幻灯片上展示的是“扫描-扫描-相加”实现，我们首先执行每个块的扫描，然后扫描部分和，最后执行相加操作，将偏移量添加到每个对应块的扫描元素中。这是一种实现方式。我们讨论的另一种是“归约-扫描-扫描”实现，它节省了一些内存访问，特别适用于非常大的数组。

今天，我们将讨论稀疏矩阵及其计算。首先，让我们定义什么是稀疏矩阵。我们都知道什么是稠密矩阵，即大多数元素不为零的矩阵。而在稀疏矩阵中，大多数元素为零。这是一个许多元素为零的矩阵，这在许多现实世界的系统中都会出现。但稀疏矩阵也提供了一些机会。首先，不需要为零元素分配空间，这可以节省内存容量。其次，不需要将零从内存加载到处理器，这可以节省内存带宽。最后，不需要对零进行计算，这将节省计算时间。例如，在稀疏矩阵向量乘法中，乘以零是没有意义的，因为我们知道结果将是零。这些都是机会，但要利用这些机会也带来了一些挑战，因为以压缩格式存储需要更复杂的数据访问方式，而如果我们以稠密矩阵存储一切，则访问方式会简单得多。因此，我们今天将对此进行一些讨论，如果你想了解更多，可以观看本讲座的扩展版本。

让我们从一些动机开始。稀疏矩阵在今天非常普遍，你可以在推荐系统、图分析（如广度优先搜索等算法）以及神经网络（包括稀疏深度神经网络或图神经网络）中找到它们。现实情况是，现实世界的矩阵具有很高的稀疏性。如果我们观察连接图，例如Facebook或YouTube的图，我们会发现它们非常稀疏。因此，良好的稀疏矩阵压缩对于实现高效的存储和计算至关重要。确实，稀疏矩阵有许多存储格式，最著名的是坐标格式（COO）和压缩稀疏行格式（CSR），我们今天将简要讨论这两种格式。但还有其他更复杂的格式，可能对GPU上的某些内核或算法带来更好的性能，例如LPAac格式、ELL格式或JDS格式。此外，还有更多格式，例如基于位图的格式，对于某些矩阵来说至少在存储上非常高效。在选择存储格式时，需要考虑一些设计因素。首先是空间效率，即它们消耗多少内存。其次是灵活性，即添加或重新排序元素的难易程度。然后是访问性，即查找所需数据的难易程度。我们将看到这些因素的一些例子。还有内存访问模式，即在像GPU这样的并行系统中，它们是否能够实现合并内存访问。最后是负载均衡，即它们能在多大程度上最小化控制流分支，因为矩阵的稀疏性可能导致不同行具有不同数量的元素，如果我们让不同的线程处理不同的行，这可能会成为分支的来源，我们将在后面看到一些例子。

我们将以SpMV为例。最佳存储格式的选择取决于计算和矩阵特性，主要是稀疏性。如前所述，我们将使用稀疏矩阵向量乘法作为例子来研究不同的格式。在矩阵向量乘法中，我们有一个稀疏矩阵（可以在图的左侧看到）和一个稠密向量，结果将是一个输出稠密向量。

让我们从坐标格式（COO）开始。请注意，这里我们有原始矩阵，这里我们以稀疏压缩格式存储矩阵。注意，对于稀疏矩阵的每个元素，我们存储非零元素及其行索引和列索引。例如，这个值为1的元素具有索引0和0，这就是我们在稀疏矩阵格式中存储的内容。

当使用COO执行矩阵向量乘法（稀疏矩阵向量乘法）时，一种典型的并行化方法是给每个非零元素分配一个线程，就像我们在图中做的那样。然后每个线程必须访问输入向量的对应元素，执行乘法，然后累加。问题在于，我们有多个线程写入同一个输出向量的同一个元素，这将要求我们使用原子操作。例如，你可以在这张幻灯片上看到示例代码。在这个SpMV COO内核中，首先要定义每个线程将要处理的元素的索引。使用这个索引，我们首先访问矩阵（稀疏压缩矩阵）中的行和列索引，然后读取该特定元素的值。使用这个值以及由列索引给出的输入向量中的对应元素，我们执行部分乘法，然后需要累加。如前所述，因为我们有多个线程，或者可能有多个线程在同一行上工作并更新输出向量的同一个元素，我们需要使用原子加法来避免数据竞争。COO有一定的权衡，它有优点也有缺点。在优点方面，它很灵活，因为很容易向矩阵添加新元素。非零元素可以按任何顺序存储，因为每个非零元素都伴随着其行索引和列索引，所以添加新元素非常容易，只需将其放在末尾即可。在访问性方面，给定一个非零元素，很容易找到其行和列。在使用SpMV时，它具有合并内存访问，因为连续的线程访问矩阵的连续元素，这保证了合并内存访问。同时，没有控制流分支，你可以在前面的代码中检查这一点。然而，它也有缺点。首先，如果给定一行或一列，要求找到该行或列的所有非零元素，那将不容易，我们必须扫描整个稀疏矩阵，除非它们已经排序，或者我们在搜索前进行预排序。其次，使用COO的SpMV版本需要使用原子操作，正如我们所见。

我们将简要讨论另一种重要的稀疏矩阵格式，它可能应用更广泛，即压缩稀疏行格式（CSR）。在这种情况下，CSR有一个关键优势，即通常提供高压缩比，并且在许多CPU和GPU的库和框架中广泛使用。在压缩稀疏行格式中，关键是拥有一个行指针数组。我们将同一行的非零元素连续存储，并且有一个索引或指针指向每行的第一个元素。因此，现在我们不再为每个值存储行和列，我们仍然存储列索引，但由于我们基于同一行的所有值将相邻存储的假设来实现CSR格式，我们可以使用这个行指针数组来指示列索引数组和值数组中每行对应元素的起始位置。

现在，当为SpMV实现并行化方法时，通常的做法是给不同的线程分配不同的行。这样，每个线程将负责计算一行元素与输入向量的部分乘积和累加。例如，如果线程从这里开始，在每次迭代中，它们将访问矩阵，读取矩阵的一个元素，然后访问向量，读取向量的对应元素，执行部分乘法并累加，在下次迭代中，它们将访问下一个元素，依此类推。请注意，我们已经可以看到的一件事是，很可能存在线程分支，即线程束间的分支。为什么会这样？因为不同的线程在不同的行上工作，不同的行可能具有不同数量的非零元素，这将成为CSR格式的缺点之一。在这里你可以看到代码。首先，我们根据块索引、块维度和线程索引确定要处理的行。基本上，每个线程处理一行。然后我们遍历该行的所有元素，读取列索引，然后读取值，并根据列索引和稀疏矩阵的值与向量的对应元素执行乘法。之后，我们累加这个sum变量，最后将其写入输出向量。

与COO相比，CSR也有一些权衡，有优点也有缺点。例如，在优点方面，它空间效率高，因为行指针比行索引小。注意，现在我们不需要为每个线程保留行索引，我们只为每行保留一个指向矩阵中该行起始位置的指针。在访问性方面，给定一行，很容易找到所有非零元素，如果我们给定一行，所有非零元素将在同一行中相邻存储。在优点方面，SpMV的实现，正如我们刚刚看到的，不需要原子操作，因为每个线程拥有自己的输出，并分配一个输出元素。在缺点方面，灵活性不太好，因为很难向矩阵添加新元素，我们需要重新排列或重组不同行中的数据。在访问性方面，给定一个非零元素，很难找到其所在行；给定一列，很难找到该列的所有非零元素，我们必须遍历所有行段，查找哪些元素在我们正在寻找的特定列中。正如我们所见，SpMV实现中的内存访问不是合并的，原因在于我们在压缩格式中按行存储，并且将不同的行分配给不同的线程，所以一个线程访问这里，另一个线程访问那里，它们之间有一定的距离，因此不是合并访问。最后，如前所述，它存在控制流分支，因为每行非零元素的数量不同，所以每个线程必须执行的迭代次数也不同。

SpMV通常是一个具有挑战性的算法，难以在并行机器上高效实现。但它还有一个额外的缺点，即它是一个内存受限的计算。内存受限意味着它需要许多内存访问，内存访问非常密集。请注意，正如我们刚刚在SpMV实现中看到的，我们不仅需要访问非零值，而且在CSR格式中，我们还必须使用行指针数组，并且还必须访问列索引数组，因此有许多内存访问，而计算量并不大，因为每个元素只有一次乘法和一次加法。在这个意义上，它是一个内存受限的操作。如果我们使用屋顶线模型进行分析，我们会发现SpMV落在屋顶线的内存受限区域。正因为如此，SpMV被认为是内存处理系统中的一个重要且合适的原语。在内存处理系统中，我们通常有一种称为PIM的内存，其中DRAM存储体或内存阵列与某种计算单元或处理元素（或PIM核心）集成在一起，允许我们在内存附近进行计算。这样，这些核心享有低内存访问延迟和大内存带宽。随着UPMEM的出现，以及三星、海力士等公司的原型，这些PIM系统正在成为现实。在我们的研究中，我们对SpMV进行了广泛的分析，并为真实的处理内存系统创建了一个SpMV库。这个库包含25个SpMV内核，适用于不同的压缩格式、数据类型、数据分区技术、负载均衡技术和同步方法。我们还提供了对首个真实世界处理内存系统（来自UPMEM）上SpMV实现的全面分析，使用了多达26个稀疏矩阵。正如你在这项工作中将看到的，我们探索了不同的分区技术，例如一维分区或二维分区。在一维分区中，我们可以在各个PIM核心上执行完整的计算。在二维分区中，我们需要在CPU上合并结果，但同时我们节省了主内存和PIM内存之间的数据移动，因此这里存在有趣的权衡。在这里，你可以找到包含该库的存储库链接，该库包含所有不同的一维和二维版本，使用不同的稀疏矩阵、不同的负载均衡、同步方法和数据类型。同样，在本课程中，我们提供了一个关于这个SpMV库的简短介绍讲座，如果你有兴趣了解更多关于SpMV及其关键限制，以及如何通过处理内存系统克服这些限制，可以在这里找到链接。

如果你总体上对学习更多关于稀疏矩阵计算感兴趣，我推荐你阅读《大规模并行处理器编程》的第14章，当然还有本讲座的更长版本，在那里我们讨论了更多的压缩格式，包括一些使用位图的新格式，这些格式非常高效，并且允许在SpMV或SpMM和其他稀疏矩阵运算上进行更快的计算。这就是今天的全部内容，非常感谢你的关注。如果你想讨论关于本讲座或课程的任何内容，请告诉我，希望在下一次讲座中见到你。

012：图搜索 🕸️

在本节课中，我们将学习另一种重要的并行模式：图搜索。这是继规约、直方图、卷积、前缀和以及稀疏矩阵之后，我们介绍的又一个核心并行模式。图搜索在社交网络分析、路径规划、电子设计自动化等许多领域都有广泛应用。

概述

在前几节课中，我们学习了多种并行模式。我们首先从规约开始，这是一种将一组值合并为单个值的操作，并讨论了通过发散前映射来最大化资源利用率的实现方法。接着，我们介绍了直方图计算，这是一种用于降低数据维度并提取大型数据集特征的计算。在直方图计算中，我们通常需要使用原子操作，并讨论了通过私有化来缓解串行化的方法。今天，我们将在图搜索的背景下再次讨论私有化技术。

之后，我们学习了卷积操作，它通过对输入的每个元素应用滤波器或掩码来计算加权和，这在机器学习和人工智能的卷积层中非常常见。我们还学习了前缀和，这是一种基于先前输出递归计算输出的操作，并讨论了其分层实现方法。在上节课中，我们讨论了稀疏矩阵及其压缩格式，例如压缩稀疏行格式。由于图通常非常稀疏，这些压缩格式在图搜索中也将至关重要。

图搜索在处理动态数据结构时面临重要挑战，因为数据布局不易实现局部性、合并访问和避免争用。此外，每个计算阶段的并行工作量可能增长或收缩。在本节课中，我们将以广度优先搜索为例，探讨如何应对这些挑战。

图与稀疏矩阵

图与稀疏矩阵表示密切相关，因为图通常表示为邻接矩阵。压缩这种邻接矩阵需要使用稀疏矩阵格式。

例如，我们可以使用压缩稀疏行格式来表示一个无权图的邻接矩阵。该格式包含一个行指针数组，指向列索引数组，列索引数组存储每条边的目标顶点。由于是无权图，非零元素数组中的所有值均为1。

CSR格式示例：

row_ptr: 行指针数组，指向每个源顶点的边在col_idx中的起始位置。
col_idx: 列索引数组，存储每条边指向的目标顶点。
values: 非零值数组（在无权图中通常全为1）。

广度优先搜索

广度优先搜索的目标是确定从源节点到目标节点（或所有目标节点）所需的最小跳数。

例如，在下图中，以节点0为源节点，不同顶点到源节点的距离（跳数）不同：两个顶点距离为1跳，五个顶点距离为2跳，一个顶点距离为3跳。

并行化图处理的方法

并行化图处理主要有两种方法：顶点中心法和边中心法。

顶点中心法

在顶点中心法中，我们为每个顶点分配一个线程来处理相关操作，并通常使用CSR或CSC等压缩格式存储图，因为给定一个顶点，很容易找到其所有邻居。

顶点中心法又分为两种策略：自顶向下和自底向上。

自顶向下策略：从单个源顶点开始，创建一个待访问顶点集合（称为“前沿”）。在每次算法迭代中，处理当前前沿中的顶点，并发现它们的邻居，这些邻居构成下一次迭代的前沿。这是我们本节课将详细讲解的方法。

自底向上策略：同样为每个顶点分配一个线程。在每次迭代中，每个线程检查其对应的顶点，看是否有邻居顶点在前一次迭代中被访问过。如果是，则将该顶点标记为属于当前迭代。这种方法在初始阶段所有线程都处于活动状态，计算量更恒定。

边中心法

在边中心法中，我们为每条边分配一个线程，并对该边执行一些计算，例如查找源顶点和目标顶点。这种方法通常使用坐标格式存储图，因为给定一条边，很容易找到其源和目标顶点。还有一些混合方法，例如为每条边分配线程，但需要查找源顶点和目标顶点的邻居，这时结合使用COO和CSR/CSC格式可能更有效。

自顶向下的BFS算法

让我们聚焦于使用顶点中心法、自顶向下策略的BFS。

我们以源节点0开始，识别并标记所有可以在1、2、3跳内到达的节点。这些节点将在BFS算法的每次连续迭代中被发现。

第一次迭代：我们发现顶点1和2，将它们放入输出队列，该队列将成为下一次迭代的输入队列。
第二次迭代（前沿）：我们处理顶点1和2的邻居，即节点3到7，并将它们加入输出队列。
第三次迭代：最后的前沿将是节点8。

至此，我们得到了以节点0为源的BFS最终结果。如果选择其他节点作为源，前沿会不同，但经过几次迭代后，最终会访问到图中的所有节点。

处理前沿与并行化挑战

当处理一个前沿时，我们需要访问源数组、目标数组以及一个标签数组。在标签数组中，我们将标记每个特定节点是在算法的哪个阶段被发现的。

我们可以尝试简单地并行化这个过程。假设当前输入前沿中有三个节点（5、6、7）。我们可以为队列中的不同节点分配不同的线程并行处理。然而，当线程需要将其发现的邻居加入全局输出队列时，为了避免数据竞争，必须使用原子操作来更新一个全局计数器以获取写入位置。虽然这能保证结果正确，但原子操作会导致串行化，严重损害性能。

解决方案：层次化队列与私有化

我们已经知道如何处理原子操作带来的问题，一种方法就是私有化。我们在直方图计算的课程中解释过，这里将以不同的方式应用它。

我们不使用所有GPU线程共享的单个全局输出队列，而是使用位于共享内存中的本地队列。每当一个线程块中的线程发现新邻居并需要将其加入前沿时，它们会先写入这个私有队列。当计算完成或私有队列满时，我们再通过更新一个全局计数器（或进行一次扫描操作，如果我们知道每个本地队列的确切元素数量）来确定在全局队列中的写入位置。

这就是我们所说的两层层次结构：块队列（通常驻留在共享内存）和全局队列（驻留在全局内存）。我们只在块计算完成时才将块队列中的元素插入全局队列。

在块队列中仍然可能存在冲突，但由于使用了私有化且共享内存速度更快，其危害远小于全局队列中的冲突。我们还可以根据算法复杂度、图的大小以及可用内存（如寄存器和片上内存）向此层次结构添加更多级别，例如线程级队列或线程束级队列。

全局同步与内核执行

我们将启动多次迭代，每次迭代启动一个新的内核。这就需要创建全局屏障，而创建全局屏障的方法就是终止一个内核，然后为算法的下一次迭代调用下一个内核。

这种内核终止和重新启动可能会带来很大开销。因此，人们提出了一些解决方案。

分层内核安排

一种解决方案是分层内核安排，即根据前沿的大小使用不同类型的内核。

对于非常小的前沿：我们可能只需要一个GPU核心和一个线程块。这个线程块可以使用__syncthreads()进行线程块内同步，处理速度快，主要涉及共享内存和全局内存访问。
对于较大的前沿：我们需要启动更大的内核，利用更多GPU核心。当前沿处理完成后，仍需终止内核并重新启动新内核以进行下一次迭代。

在许多情况下，如果内核内部有足够多的计算量，这种内核启动开销是可以接受的。但对于像BFS这样主要涉及内存访问和简单更新的计算，内核启动开销可能仍然显著。

动态并行

使用动态并行可以在不终止父内核的情况下启动新的子内核来处理新发现的邻居，这可能会减少一些开销。我们将在后续课程中详细讨论动态并行。

持久化线程与块间同步

最后一种解决方案是使用持久化线程和块间同步。这种方法结合了上述两种内核的优点，通过避免内核重新启动来减少开销。

其核心思想是：我们不启动足够多的线程块来覆盖整个输入前沿，而是只启动能在现有GPU核心上并发运行的线程块数量。例如，如果两个流式多处理器最多只能同时容纳四个线程块，我们就只启动四个持久化的线程块。

这些线程块会处理当前前沿的一部分，完成后，它们会通过原子操作和全局内存进行同步，然后去处理前沿中尚未被处理的部分。这样，多个迭代可以在同一个内核中完成，无需反复启动和终止内核。

总结

本节课我们一起学习了并行模式中的图搜索，重点以广度优先搜索为例。

我们回顾了图与稀疏矩阵表示的关系，介绍了并行化图处理的顶点中心法和边中心法，并深入探讨了自顶向下的BFS策略。我们分析了简单并行化面临的原子操作串行化挑战，并提出了通过层次化队列和内存私有化来提升性能的解决方案。最后，我们讨论了管理BFS多迭代执行的几种全局同步策略，包括分层内核安排、动态并行和持久化线程。

图搜索是处理不规则、动态数据结构的典型范例，理解其并行化策略对于高效利用GPU等加速器至关重要。

013：动态并行性

概述

在本节课中，我们将学习GPU编程框架中的一个有趣特性：动态并行性。动态并行性提供了一种接口，能以更自然的方式表达动态细化算法。它允许GPU线程在动态发现新工作时启动新的GPU内核。

动态并行性简介

在之前的课程中，我们介绍了多种并行模式。本节中，我们将讨论动态并行性。动态并行性允许GPU线程从设备端启动内核，这被称为设备端内核启动。该特性随Kepler架构出现，在此之前，只能从主机处理器启动内核。

动态并行性的典型用例包括：具有动态负载均衡的工作负载、具有数据依赖执行的工作负载，以及需要递归实现的工作负载。它对于库调用也很有用。动态并行性的最大优势在于提高了程序的可编程性和可维护性。

回想一下，在没有动态并行性的情况下启动内核，总是需要主机处理器的参与来执行内核启动，这在编程上相对繁琐。使用动态并行性，一旦在GPU内部发现新工作，就可以直接启动新的内核。这样，我们只需从主机启动一个初始内核，之后，设备中运行的线程将能够启动其子内核来处理动态发现的工作并并行执行。

同步与内存一致性

从父内核启动子内核时，一个重要的考虑因素是同步和内存一致性的需求。这里需要注意两点：

从父内核到子内核的内存一致性是有保证的。这意味着当子内核开始运行时，父内核需要写入全局内存的所有数据都已就绪。
从子内核到父内核的方向，我们需要使用 cudaDeviceSynchronize 来确保在子内核完成后，父内核才会继续执行。

一个简单示例

为了了解如何使用动态并行性，让我们从一个简单的示例开始。首先看一段没有使用动态并行性的内核代码：

__global__ void kernelWithoutDP(int start, int end) {
    for (int i = start + threadIdx.x; i < end; i += blockDim.x) {
        // 执行某些工作
    }
}

只要所有线程的 start 和 end 值相同，每个线程的迭代次数就相同，负载是均衡的。然而，如果每个线程的 start 和 end 值不同，就会导致负载不均衡，影响性能。

使用动态并行性，我们可以为每个父线程启动一个子内核。子内核中的线程数等于父线程中for循环的迭代次数。这样，即使每个父线程的工作量不同，每个子线程的工作量却是相同的，从而实现了更好的负载均衡。

__global__ void parentKernel() {
    int my_start = ...; // 每个线程自己的起始值
    int my_end = ...;   // 每个线程自己的结束值
    int iterations = my_end - my_start;

    // 启动子内核，线程数等于迭代次数
    childKernel<<<1, iterations>>>(my_start);
}

__global__ void childKernel(int start) {
    int i = start + threadIdx.x;
    // 执行原本在for循环中的工作
}

递归应用：四叉树构建

动态并行性对于实现递归程序也很有用。本节我们来看一个四叉树的例子。四叉树是一种用于在非均匀平面中划分数据的技术，目的是将点聚类，以便后续更好地进行合并访问和计算。

在四叉树中，我们通过递归地将一个2D平面划分为四个象限来分割它，直到每个象限中的点数少于一个阈值（假设阈值为2）。这是一个递归过程。

以下是如何在GPU上实现它的步骤：

主机启动一个仅包含一个线程块的内核，并将整个2D空间分配给该线程块。
该线程块中的线程首先检查当前节点（象限）中的点数是否小于等于阈值，或者是否已达到最大递归深度。如果是，则线程块退出。
否则，线程计算包含所有点的包围盒的中心点，该中心点定义了四个象限的分界。
接着，确定每个象限中的点数。
执行扫描操作，以确定每个象限的点在输出数组中的起始偏移量。
根据计算出的偏移量，将输入缓冲区中的点重新排序到输出缓冲区中，使得属于同一象限的点连续存放。
最后，由线程块中的一个代表线程（例如线程0）启动一个包含四个子线程块的新内核，每个子线程块处理一个象限。

这个过程会递归进行，直到所有象限都满足退出条件。通过动态并行性，每个递归层级的线程块都可以自主启动处理其子象限的新内核。

性能考量与优化

尽管动态并行性可以带来更好的负载均衡，并在可编程性和可维护性方面具有优势，但在实现时也面临一些性能挑战。

主要挑战包括：

如果启动的网格包含的线程数量很少，可能导致GPU资源利用不足。一般建议是，子网格应包含尽可能多的线程块，以充分利用GPU资源。
对于递归实现的并行性，最大嵌套深度受硬件限制，因此只能高效实现相对浅的树结构。此外，建议使用“粗粒度”的树节点，即每个节点部署大量线程，或者分支因子较大（每个父节点有许多子节点），这与本页的一般建议是一致的。

为了优化动态并行性，减少内核启动开销并提高GPU资源利用率，可以采用内核启动聚合技术。

其核心思想是将多个小的内核启动合并为单个更大的内核启动。例如：

Warp级聚合：每个Warp只启动一个内核，由Warp内的线程通信，将各自需要的工作量汇总，然后由一个领导线程启动一个包含足够多线程块的大内核。
Block级聚合：类似地，在线程块级别进行聚合。
内核级聚合：甚至可以从主机端直接启动一个大的内核。

这种方法可以减少内核启动的总次数，从而降低开销，并有可能生成更高效、能更好利用GPU资源的工作负载。

总结

本节课我们一起学习了GPU编程中的动态并行性。我们了解了它的基本概念、如何用于改进负载均衡和实现递归算法（如四叉树构建）。同时，我们也讨论了使用动态并行性时需要注意的性能挑战，并简要介绍了一种通过聚合内核启动来优化性能的技术。动态并行性是一个强大的工具，能够简化某些复杂并行模式的编程，但在使用时需仔细考虑其对性能的影响。

014：协同计算

概述

在本节课中，我们将学习异构系统中的协同计算，也称为细粒度协同或异构执行。我们将重点关注如何在包含不同设备（特别是CPU和GPU）的系统中实现协同计算。实现协同计算的一个必要特性是统一内存，这是一种可以被多种设备同时访问的虚拟内存系统。

统一内存简介

在传统的GPU计算中，我们无法拥有统一的内存视图。通常需要在主机处理器和设备GPU上分别进行内存分配，并使用显式的数据传输操作（如 cudaMemcpy）在设备间移动数据。这种方法自然适用于独立GPU的系统。

以下是一段传统方法的示例代码：

// 主机内存分配
float *h_data = (float*)malloc(N * sizeof(float));
// 设备内存分配
float *d_data;
cudaMalloc(&d_data, N * sizeof(float));
// 数据传输：主机到设备
cudaMemcpy(d_data, h_data, N * sizeof(float), cudaMemcpyHostToDevice);
// 启动GPU内核
kernel<<<...>>>(d_data);
// 同步并传回数据
cudaDeviceSynchronize();
cudaMemcpy(h_data, d_data, N * sizeof(float), cudaMemcpyDeviceToHost);

然而，通过统一内存，我们可以拥有一个统一的虚拟地址空间，该空间在主机和设备之间共享。这一特性始于CUDA 6和Kepler架构，但最初它只是隐藏了双重内存分配和显式数据传输的层。自CUDA 8和Pascal架构起，GPU开始支持页面错误，这意味着GPU内存可以被超额订阅，GPU可以访问CPU内存来获取数据并继续执行，唯一的限制是系统内存大小。

统一内存不仅支持内存超额订阅，还简化了编程。现在，我们不再需要进行双重内存分配，只需使用 cudaMallocManaged 这个新的API在统一内存空间中分配数组。这样，CPU和GPU都可以访问同一块内存，无需进行显式的数据传输。这种方法自然适用于CPU和GPU集成在同一芯片上并共享物理内存的系统，但也可以为独立GPU实现，并且所有新版本的CUDA都支持此功能。

协同计算模式与应用案例

上一节我们介绍了统一内存的概念，本节中我们来看看如何利用这些新特性来实现协同应用程序。我们将讨论几个使用CPU和GPU的协同计算用例。

协同计算应用在CPU和GPU系统中一直是可能的，主要是因为从主机端启动内核是异步的。这意味着CPU在GPU上启动内核后可以继续执行其他任务，然后在某个时刻使用 cudaDeviceSynchronize 来检查GPU内核是否执行完毕，再复制数据。这在统一内存出现之前就已经可以实现。

但有了统一内存，协同计算变得更加容易，并且支持更细粒度的协作。特别是在Pascal架构之后，统一内存支持CPU和GPU内存一致性以及系统范围的原子操作，使得我们可以编写更复杂的协同代码。

以下是使用统一内存和系统原子操作的代码示例：

// 在统一内存中分配输入和输出
cudaMallocManaged(&input, size);
cudaMallocManaged(&output, size);
// 启动GPU内核
kernel<<<...>>>(input, output);
// CPU可以立即继续计算，甚至可以访问GPU正在使用的数据
// 使用系统范围的原子操作防止数据竞争
atomicAdd_system(&shared_var, value);

在GPU端，系统范围的原子操作语法如下，它与普通的原子操作（如 atomicAdd）语义相同，但使用 _system 后缀来标识其作用范围是整个系统：

atomicAdd_system(&data, increment);

现在我们对统一内存及其在细粒度协作中的应用有了更多了解。接下来，我们将简要介绍一些协同模式及其示例。

数据分区模式

数据分区模式是指将数据并行任务划分到不同的设备上执行。例如，设备1可以是CPU，设备2可以是GPU。某些任务在一个设备上执行可能比在另一个设备上更快，但另一个设备可能具有更高的并行度。在执行过程中，我们可能需要在继续执行第二个内核之前进行粗粒度同步。

以下是实现数据分区的一种方式：

将任务划分为多个数据并行子任务。
根据设备特性（如CPU擅长小任务，GPU擅长大规模并行任务）静态或动态地将子任务分配给CPU和GPU。
设备并行执行各自分配的任务。
在必要时进行同步，然后继续执行程序的下一个阶段。

任务分区模式

任务分区模式是指将不同的任务（或内核）分配给不同的设备执行。这可以实现并发执行，例如通过流水线方式，当设备2正在计算当前任务块时，设备1可以计算前一个任务块。

统一内存和系统范围的原子操作还使我们能够实现细粒度任务分区。在这种模式下，我们知道某些子任务在一个设备上运行得更快，因此可以使用细粒度的通信来通知另一个设备上的线程。例如，设备1完成特定子任务（黄色子任务）后，可以通知设备2上的线程，使其可以继续执行任务的第二部分（深绿色子任务）。我们还可以为同一程序中的不同内核应用不同的分区策略。

具体应用案例

案例一：曲面生成（数据分区）

我们以一个在三维空间中生成曲面的算法为例。该算法使用控制点网格，并基于B样条多项式进行参数化有理公式计算。

在这个例子中，我们将曲面划分为多个图块，并将不同的图块静态分配给CPU或GPU。我们可以预先分析CPU和GPU计算这些图块的速度，从而进行静态分区。即使没有统一内存，也可以通过启动CPU线程和GPU内核来实现，并在最后将GPU计算的部分复制回CPU内存。

使用统一内存后，实现变得更加简单：

使用 cudaMallocManaged 分配单个曲面内存。
启动CPU线程和GPU内核。
由于使用统一内存，无需将GPU计算的部分显式复制回CPU内存，这会隐式完成。

在支持系统范围原子操作的更新版本统一内存中，我们甚至可以实施动态图块分配，无需进行离线性能分析。这对于数据依赖性强或不规则、负载不平衡的应用程序效率更高。动态分配需要一个全局变量（图块计数器）来告诉每个GPU线程块或CPU线程下一个要计算的图块，并使用 atomicAdd_system 来更新这个原子变量。

在一项测试中，使用这种协作方式，在Ada GPU上结合CPU核心，相比纯GPU版本获得了高达40%的性能提升。

案例二：广度优先搜索（协同分区）

广度优先搜索算法包含两个不同的内核：一个用于小规模边界，另一个用于大规模边界。对于小规模边界，即使我们优化了内核，GPU资源也可能利用不足。

在NVIDIA Jetson设备上的测试表明，小规模边界在CPU上运行得更快（对应执行初期的红色条），而大规模边界则适合GPU处理。通过协同实现，我们可以根据边界大小动态选择最合适的设备。

没有统一内存时，我们需要显式地进行内存拷贝，并在主机代码中使用循环来检查边界大小并决定启动CPU线程还是GPU内核。

有了统一内存，编程更简单，无需显式内存拷贝。而有了支持系统范围原子操作的统一内存，我们甚至可以在开始时同时启动CPU线程和GPU内核，让它们持续运行，根据边界大小通过系统原子操作进行通信和同步，共同处理同一个边界。

对于某些图，这种方法的加速比可达34-39%，优于纯GPU版本。当然，对于非常大的图，其边界通常也很大，大部分时间还是会使用GPU。

案例三：RANSAC算法（细粒度任务分区）

RANSAC算法用于寻找符合特定观测值的最佳模型，可应用于物体检测等领域。该算法包含多个迭代，直到找到最佳模型。每个迭代包含一个拟合阶段（顺序计算）和一个评估阶段（数据并行计算）。

由于迭代是基于随机采样的，因此各迭代之间相互独立。我们可以并行化所有迭代，在CPU线程和GPU线程块中计算它们，并使用原子操作进行细粒度通信。这样，CPU和GPU可以协作处理不同的迭代，或者协作处理同一迭代的不同部分，从而加速整体计算。

总结与资源

本节课我们一起学习了异构系统中的协同计算。我们首先介绍了实现协同计算的关键——统一内存，它简化了内存管理并支持更细粒度的设备协作。接着，我们探讨了数据分区、任务分区和细粒度任务分区等协同模式，并通过曲面生成、广度优先搜索和RANSAC算法三个具体案例，说明了如何在实际应用中实现CPU与GPU的协同工作。

更多示例可以在CHAI基准测试套件中找到，该套件包含了离散和统一内存版本、CUDA/OpenCL实现，甚至GPU模拟器版本。所有课程材料均公开提供。

posted @ 2026-03-26 08:20 布客飞龙II 阅读(0) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

人最大的痛苦就是说一些自己都不相信的话。