贝鲁特大学-CMPS297S-GPU-编程笔记-全-

贝鲁特大学 CMPS297S GPU 编程笔记（全）

GPU计算：P1：引言 🚀

在本节课中，我们将要学习GPU计算的基本概念，了解其诞生的历史背景和设计哲学。我们将从计算机技术的发展趋势讲起，探讨为何并行计算成为主流，并深入对比CPU与GPU这两种不同设计取向的处理器。

摩尔定律的兴衰 💡

上一节我们介绍了课程主题，本节中我们来看看推动计算技术发展的核心动力——摩尔定律。

摩尔定律并非物理定律，而是英特尔联合创始人戈登·摩尔提出的一个预测。他预测单位面积上的晶体管数量每18到24个月会翻一番。这个预测在过去的半个世纪里基本成立，成为了半导体行业发展的目标。

然而，这一趋势正在走向终结。晶体管尺寸已缩小到物理极限，难以继续微缩。更重要的是，在2005年左右，处理器频率的提升也停滞了。这是因为更快的开关速度意味着更高的功耗，而功耗会产生热量。当热量超过冷却技术的处理能力时，处理器就无法在更高频率下稳定运行。

在2005年之前，如果你购买一台新电脑，即使运行相同的程序，速度也会因为频率提升而自动变快。2005年之后，这种“免费的性能午餐”结束了。单线程性能虽然仍在提升（得益于编译器优化和更复杂的硬件技术，如分支预测、乱序执行），但提升幅度已大不如前。

并行计算成为主流 ⚙️

由于单线程性能提升放缓，但晶体管数量仍在增长，硬件设计者开始将这些额外的晶体管用于增加处理器核心数量。这就是并行计算成为主流的转折点。

硬件的变化驱动了软件的变革。为了获得更好的性能，软件开发人员必须重写程序，使其能够并行运行。免费的午餐结束了，性能提升必须通过并行编程来实现。

两种设计哲学：延迟导向 vs. 吞吐量导向 🏎️🚌

在深入处理器细节前，我们先了解两种通用的系统设计思路：延迟导向设计和吞吐量导向设计。

延迟导向设计旨在最小化完成单个任务所需的时间。一个典型的例子是小汽车，它追求以最快速度将一个人从A点送到B点。

吞吐量导向设计则旨在最大化在给定时间内完成的任务数量。一个典型的例子是公交车，它虽然运送每个乘客的时间更长，但一次可以运送很多人。

以下是这两种设计思路的对比：

延迟导向：优化单个任务的速度。
吞吐量导向：优化单位时间内完成的任务总数。

这两种设计哲学同样适用于处理器设计。

CPU与GPU的架构对比 🖥️🎮

上一节我们介绍了两种设计哲学，本节中我们来看看它们如何具体体现在CPU和GPU的架构差异上。

CPU是延迟导向设计的典范，而GPU则是吞吐量导向设计的代表。这种根本差异影响了它们各个组件的设计。

以下是CPU与GPU在几个关键组件上的设计差异：

1. 算术逻辑单元

CPU：拥有少量但非常强大的ALU。这些ALU占用大量芯片面积，经过高度优化，旨在将一次算术运算的延迟降到最低。
GPU：拥有大量小型、简单的ALU。每个ALU执行一次运算的延迟较长，但通过大规模并行和流水线技术，可以获得极高的总体算术吞吐量。

2. 缓存

CPU：配备大容量缓存。目的是将访问慢速DRAM内存的高延迟操作，转化为访问快速缓存的低延迟操作，从而降低平均内存访问时间。
GPU：缓存容量较小。这意味着缓存命中率较低，更多需要访问高延迟的内存。但节省下来的芯片面积可以用于放置更多的ALU来计算。

3. 控制逻辑

CPU：采用复杂技术（如分支预测、数据前递、乱序执行）来减少控制冒险和数据冒险，优化单线程的执行延迟。这些技术需要额外的硬件支持。
GPU：控制逻辑相对简单。为了容纳更多计算单元，GPU需要容忍更长的分支延迟和数据冒险。

4. 多线程与延迟隐藏

CPU：由于操作延迟相对较低，只需要适度的多线程（如每个核心2个线程，即超线程技术）来填充流水线中的空闲周期，隐藏短延迟。
GPU：由于ALU、内存访问和控制决策的延迟都很高，为了隐藏这些长延迟、避免流水线停滞，GPU采用了大规模的多线程技术。单个GPU核心上可以同时驻留并调度数十个线程，确保在任何周期都有可执行的指令。

5. 时钟频率

CPU：通常具有较高的时钟频率。
GPU：通常运行在较低的时钟频率下。这是因为GPU已经集成了大量计算单元，功耗很高，降低频率有助于控制总功耗和散热。

GPU的起源与发展 📈

GPU最初是为图形处理而设计的，这也是其名称（图形处理器）的由来。图形工作负载（如渲染数百万像素）天生具有大规模并行性，且像素间通常相互独立，这正好契合吞吐量导向的设计理念。

2007年是一个关键转折点。英伟达发布了CUDA平台，它提供了一套编程接口，让开发者能够直接编写通用目的程序在GPU上运行，而无需再将计算任务“伪装”成图形渲染任务。这标志着GPU计算时代的开始，GPU开始被广泛用于科学计算、深度学习等非图形领域。

数据显示，自CUDA发布以来，GPU在通用计算领域的应用呈爆炸式增长，包括开发者数量、学术论文、高校课程以及高性能计算系统的采用率。如今，世界排名前十的超级计算机中，有多台都采用了GPU作为加速器。

为何是GPU取得成功？🎯

并行处理器有多种实现方式，为何GPU能脱颖而出？核心原因在于市场规模和成本摊销。

芯片设计制造成本极高，需要巨大的销量来分摊成本。当通用并行计算需求兴起时，GPU已经拥有了一个庞大而稳定的市场——游戏产业。这使得GPU厂商能够基于现有规模，相对容易地改造架构以适应通用计算，而无需像一家新公司那样，从零开始为某个细分市场（如科学计算）设计专用芯片，并面临市场狭小、成本难以回收的风险。

值得注意的是，游戏仍是GPU最大的营收来源。因此，GPU架构的任何改进都必须权衡对所有主要市场（尤其是游戏）的影响。一项能大幅提升科学计算性能但可能损害游戏性能的改进，很可能不会被采纳。

在本课程中，我们将聚焦于如何使用GPU进行通用目的编程，重点学习如何为科学计算等数据中心负载编写高效的GPU程序。

总结 📚

本节课我们一起学习了：

摩尔定律的终结是推动计算架构向多核并行转变的根本原因。
延迟导向设计（如CPU）追求最快完成单个任务，而吞吐量导向设计（如GPU）追求单位时间内完成最多任务。
CPU与GPU在ALU、缓存、控制逻辑、多线程程度和时钟频率上存在显著差异，这些差异源于其不同的设计目标。
GPU起源于图形处理，因其与生俱来的并行性而采用吞吐量导向设计。CUDA的发布使其能够方便地用于通用计算，从而开启了GPU计算时代。
GPU的成功得益于其早已建立的巨大市场规模（主要是游戏），这使得其在向通用计算领域扩展时具有成本优势。

对应教材章节：请参考课程指定教材《Programming Massively Parallel Processors: A Hands-on Approach》的相关引言和背景章节。

GPU计算：P2：数据并行编程

概述

在本节课中，我们将要学习数据并行编程的基本概念，并了解如何在GPU上实现它。我们将以向量加法为例，介绍CUDA编程模型的基础知识，包括内存管理、内核启动和线程组织。

课程回顾

上一节课我们简要介绍了GPU计算，探讨了其兴起的原因和历史背景。我们提到，尽管摩尔定律预测晶体管数量会持续增长，但自2005年起，处理器频率的提升遇到了瓶颈，这主要是由于功耗墙的限制。因此，为了从硬件中获得更多性能，行业转向了增加核心数量。

在2005年，我们见证了处理器核心数量的爆炸式增长，以及并行处理器，特别是作为大规模并行处理器的GPU的出现。我们讨论了面向延迟的设计和面向吞吐量的设计之间的区别。CPU是面向延迟处理器的例子，而GPU是面向吞吐量处理器的例子。这种差异体现在架构设计上：CPU拥有少量但功能强大的ALU，可以快速完成单个操作；而GPU拥有许多小型ALU，每个操作时间较长，但可以同时执行更多操作。

我们还介绍了CUDA，它使人们能够以通用方式对GPU进行编程，从而推动了使用CUDA的应用程序数量以及GPU在全球超级计算机中提供的计算能力的爆炸式增长。GPU在从具有大量并行性的应用程序中提取性能方面已被证明非常有用。

数据并行与任务并行

有两种主要的并行类型：任务并行和数据并行。它们是两种从应用程序中提取并行性的不同方法，并且可以同时使用。

任务并行是指对相同或不同数据执行的不同操作可以并行进行。例如，在使用Microsoft Word时，文本编辑和后台运行的拼写检查器可以并行运行。这被称为任务并行。

数据并行是指对大量数据的不同部分执行相同的操作可以并行进行。一个明显的例子是屏幕渲染，需要计算每个像素的值。计算不同像素值的操作非常相似，但每个操作都在不同的输入数据上进行，以计算不同的输出数据。这被称为数据并行。

通常，任务并行只能释放适度的并行性，因为一个应用程序中同时执行的不同类型的任务数量通常不会很多。而数据并行则有可能释放出大量的并行性，因为它处理的是海量数据。为了增加程序中的数据并行性，你不需要编写更多代码，只需在更大的数据集上运行相同的程序即可。这就是为什么数据并行非常适合GPU的原因。

向量加法示例

为了展示数据并行性以及如何在GPU上实现它，我们将从一个简单的计算示例开始：向量加法。你可以将向量加法视为数据并行编程的“Hello World”。

向量加法就是将两个向量相加。给定输入向量X和Y，将它们相加得到输出向量Z。具体操作是：将X的第一个元素与Y的第一个元素相加，得到Z的第一个元素；将X的第二个元素与Y的第二个元素相加，得到Z的第二个元素，依此类推。

如果要在CPU上顺序编程实现，我们只需编写一个循环来遍历这些向量并进行加法运算。

以下是CPU上顺序实现的代码示例：

void vecAddCPU(float* x, float* y, float* z, int n) {
    for (int i = 0; i < n; ++i) {
        z[i] = x[i] + y[i];
    }
}

GPU系统架构

在讨论如何在GPU上并行执行向量加法之前，我们需要了解包含GPU的系统的组织方式。

一个典型的系统有一个CPU，它可以访问一些主内存。在CUDA和GPU编程的上下文中，我们将CPU称为主机，将主内存称为主机内存。

当我们使用GPU计算时，系统中还会有一个GPU。GPU被称为设备。GPU可以访问自己的内存，我们称之为全局内存。CPU和GPU拥有独立的内存，它们不能直接访问彼此的内存，除非使用像统一虚拟内存这样的高级功能。

如果要在GPU上对数组进行加法运算，我们需要将数据从CPU的主内存传输到GPU的全局内存。这通常通过PCIe或NVLink等互连技术完成。

将计算从CPU卸载到GPU的典型操作序列如下：

在GPU上分配内存。
将数据从CPU复制到GPU。
在GPU上执行计算。
将结果从GPU内存复制回主内存。
释放GPU上的内存。

CUDA内存管理

现在，让我们看看如何为向量加法示例执行这五个步骤。首先从分配内存开始。

在GPU上分配内存，可以使用CUDA函数cudaMalloc。它接受两个参数：一个设备指针（用于存放分配内存的地址）和要分配的数据大小。

以下是如何使用cudaMalloc在GPU上分配内存的代码示例：

float *x_d, *y_d, *z_d;
cudaMalloc((void**)&x_d, n * sizeof(float));
cudaMalloc((void**)&y_d, n * sizeof(float));
cudaMalloc((void**)&z_d, n * sizeof(float));

分配内存后，我们也需要释放内存。对应的函数是cudaFree，它接受一个指向要释放数组的指针。

cudaFree(x_d);
cudaFree(y_d);
cudaFree(z_d);

cudaMalloc和cudaFree的返回类型都是cudaError_t，这是一个错误代码，用于错误检查。

数据传输

下一步是将数据复制到GPU以及从GPU复制回来。这可以通过函数cudaMemcpy完成。它接受目标地址、源地址、要复制的数据大小以及复制方向。常见的复制方向是cudaMemcpyHostToDevice（从CPU到GPU）和cudaMemcpyDeviceToHost（从GPU到CPU）。

以下是如何复制数据的代码示例：

// 将数据从主机复制到设备
cudaMemcpy(x_d, x, n * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(y_d, y, n * sizeof(float), cudaMemcpyHostToDevice);

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/8d35d538943219e48d439ec958f88b53_44.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/8d35d538943219e48d439ec958f88b53_46.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/8d35d538943219e48d439ec958f88b53_47.png)

// ... 执行GPU计算 ...

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/8d35d538943219e48d439ec958f88b53_49.png)

// 将结果从设备复制回主机
cudaMemcpy(z, z_d, n * sizeof(float), cudaMemcpyDeviceToHost);

启动GPU内核进行计算

数据就位后，我们需要在GPU上运行加法代码。我们希望并行执行向量加法，而不是顺序执行。GPU拥有大量可以并行运行的线程。

我们可以为数组中的每个元素启动一个GPU线程，每个线程负责添加一对元素并存储结果。这样，我们就将一个GPU线程映射到一个向量元素。

在GPU上，线程被组织成一个称为网格的数组。网格中的线程进一步分组为线程块。同一块中的线程可以以不同块中的线程无法实现的方式进行交互和协作。

要启动一个网格，我们需要指定网格中的块数以及每个块中的线程数。所有在同一个网格中启动的线程都执行相同的函数，这个函数被称为内核。

启动网格是通过调用一个特殊函数（内核）并告诉它网格大小（块数）和块大小（每个块的线程数）来完成的。

以下是如何启动内核的代码示例。我们假设有一个名为vecAddKernel的内核函数来执行并行向量加法。

// 定义每个块的线程数
const int threadsPerBlock = 512;
// 计算需要的块数（向上取整）
int blocks = (n + threadsPerBlock - 1) / threadsPerBlock;

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/8d35d538943219e48d439ec958f88b53_73.png)

// 启动内核
vecAddKernel<<<blocks, threadsPerBlock>>>(x_d, y_d, z_d, n);

实现GPU内核

内核类似于C或C++函数，但在其声明前有关键字__global__，表示这是一个内核。内核内部使用特殊的关键字来区分不同的线程。

gridDim.x：网格中x维度的块数。
blockIdx.x：线程所在块在网格x维度上的索引。
blockDim.x：每个块在x维度上的线程数。
threadIdx.x：线程在其块内x维度上的索引。

为了确定每个线程应该处理哪个数组元素，线程需要计算其在网格中的全局索引。全局索引可以通过以下公式计算：
int i = blockIdx.x * blockDim.x + threadIdx.x;

然后，每个线程使用其全局索引i来访问数组元素并执行加法：z[i] = x[i] + y[i];。这种并行编程方法称为单程序多数据，即多个线程执行相同的程序，但每个线程操作不同的数据集。

由于我们启动的线程总数（blocks * threadsPerBlock）可能大于数组大小n，我们需要添加边界检查，确保只有索引i小于n的线程才执行操作。

以下是完整的vecAddKernel实现：

__global__ void vecAddKernel(float* x, float* y, float* z, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        z[i] = x[i] + y[i];
    }
}

代码编译

要编译包含CUDA代码的文件，我们使用NVIDIA CUDA编译器nvcc。nvcc会将代码分离为主机（CPU）代码和设备（GPU）代码。主机代码被编译为在CPU上运行，设备代码（内核）被编译为一种称为PTX的虚拟指令集，然后在运行时即时编译为特定的GPU汇编代码。

在代码中，我们使用关键字来区分函数应该在何处运行：

__host__：默认，表示函数在CPU上运行。
__global__：表示内核，从CPU调用，在GPU上执行。
__device__：表示函数在GPU上执行，只能从GPU代码调用。
__host__ __device__：表示函数同时为CPU和GPU生成代码，可从两者调用。

性能考量与同步

默认情况下，内核调用是异步的。CPU调用内核后，会立即继续执行后续代码，而不会等待GPU内核完成。这允许CPU和GPU并行工作。如果需要等待内核完成，可以调用cudaDeviceSynchronize()。

在测量GPU内核执行时间时，通常需要在启动内核前开始计时，在调用cudaDeviceSynchronize()后停止计时，以确保内核已完成。

错误检查

CUDA函数通常返回cudaError_t类型的错误代码。良好的实践是检查这些返回值以确保操作成功。可以使用cudaGetLastError()来获取最后一个错误的信息。

cudaError_t error = cudaGetLastError();
if (error != cudaSuccess) {
    // 处理错误
    printf("CUDA error: %s\n", cudaGetErrorString(error));
}

总结

本节课中我们一起学习了数据并行编程的基础知识。我们以向量加法为例，介绍了如何在GPU上使用CUDA进行编程。关键步骤包括：在GPU上分配内存、在主机与设备间传输数据、启动指定网格和块大小的内核、在内核中计算线程的全局索引并执行并行操作，最后处理结果并释放内存。我们还了解了代码编译过程、主机/设备函数关键字、内核异步执行与同步，以及基本的错误检查方法。通过本课，你已经掌握了编写简单CUDA程序的基本技能。

GPU计算：P3：多维网格与数据 🚀

在本节课中，我们将要学习如何在GPU上处理多维数据。我们将探讨如何创建多维线程网格，以及如何高效地存储和访问多维数据。通过三个具体的例子——RGB图像转灰度图、图像模糊和矩阵乘法——我们将掌握这些核心概念。

上一节我们介绍了数据并行编程的基本概念，并使用向量加法作为例子。本节中我们来看看如何处理更复杂的多维数据。

多维网格与数据布局

许多应用中的数据是多维的，例如图像（二维）或矩阵。为了高效地处理这类数据，CUDA支持创建多维的线程网格。这简化了线程索引的计算，使代码更直观。

创建多维网格

在CUDA中，我们使用 dim3 类型来定义网格和线程块的维度。dim3 是一个包含x、y、z三个整数的结构体。

以下是配置一个二维网格的示例代码：

// 定义每个线程块的维度：32x32 个线程
dim3 numThreadsPerBlock(32, 32); // z维度默认为1

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/efca0f9f3448fd2be3ffe506fa0fb9d9_3.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/efca0f9f3448fd2be3ffe506fa0fb9d9_5.png)

// 计算网格中线程块的数量
int width = 1024; // 图像宽度
int height = 768; // 图像高度
dim3 numBlocks((width + numThreadsPerBlock.x - 1) / numThreadsPerBlock.x,
               (height + numThreadsPerBlock.y - 1) / numThreadsPerBlock.y);

// 启动内核
myKernel<<<numBlocks, numThreadsPerBlock>>>(...);

核心概念：dim3 类型用于定义多维网格和线程块。numBlocks 和 numThreadsPerBlock 的x、y分量分别对应网格的宽度和高度。

线程索引计算

在内核函数中，每个线程可以通过内置变量确定自己在多维网格中的全局位置。

以下是计算线程对应行和列的代码：

__global__ void myKernel(...) {
    // 计算线程对应的输出行和列
    unsigned int row = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned int col = blockIdx.x * blockDim.x + threadIdx.x;
    // ... 后续处理
}

核心概念：blockIdx 和 threadIdx 是 dim3 类型，分别表示线程块在网格中的索引和线程在线程块中的索引。blockDim 表示线程块的维度。

多维数据的内存布局

在C/C++中，动态分配的多维数组在内存中通常以行主序方式连续存储。这意味着数据的第一行先被存储，接着是第二行，依此类推。

给定行号 row 和列号 col，在一维数组中的索引计算公式为：

公式：index = row * width + col

其中 width 是二维数据的宽度（列数）。

示例一：RGB图像转灰度图 🖼️➡️⚫⚪

在这个例子中，我们将彩色图像（每个像素包含红、绿、蓝三个通道）转换为灰度图像（每个像素一个强度值）。

并行化策略

一个直观的并行化方法是为图像中的每个像素分配一个线程。每个线程负责读取对应像素的RGB值，计算加权平均值，并将结果写入输出数组。

以下是RGB转灰度内核的实现：

__global__ void rgbToGrayKernel(unsigned char* red, unsigned char* green,
                                unsigned char* blue, unsigned char* gray,
                                int width, int height) {
    // 1. 计算线程对应的像素位置
    unsigned int row = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned int col = blockIdx.x * blockDim.x + threadIdx.x;

    // 2. 边界检查：确保线程在图像范围内
    if (row < height && col < width) {
        // 3. 计算一维数组索引
        unsigned int idx = row * width + col;

        // 4. 执行RGB到灰度的转换（加权平均）
        // 常见权重：红色30%，绿色60%，蓝色10%
        gray[idx] = (unsigned char)(0.3f * red[idx] +
                                    0.6f * green[idx] +
                                    0.1f * blue[idx]);
    }
}

核心概念：边界检查至关重要，因为启动的线程总数可能略多于实际像素数（由于网格/块尺寸取整）。只有位于有效图像范围内的线程才应执行计算。

示例二：图像模糊 🔍➡️🌫️

图像模糊操作中，每个输出像素是其周围一个区域（例如3x3、5x5）内输入像素的平均值。

并行化策略

我们同样为每个输出像素分配一个线程。每个线程需要读取其对应输入像素及其周围多个像素的值，计算平均值，然后写入输出。

以下是图像模糊内核的实现：

#define BLUR_SIZE 1 // 模糊半径，例如1表示3x3区域

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/efca0f9f3448fd2be3ffe506fa0fb9d9_91.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/efca0f9f3448fd2be3ffe506fa0fb9d9_93.png)

__global__ void blurKernel(unsigned char* image, unsigned char* blurred,
                           int width, int height) {
    // 计算线程对应的输出像素位置
    int outRow = blockIdx.y * blockDim.y + threadIdx.y;
    int outCol = blockIdx.x * blockDim.x + threadIdx.x;

    // 边界检查：输出位置必须在图像内
    if (outRow < height && outCol < width) {
        int sum = 0;
        int pixelCount = 0;

        // 遍历以(outRow, outCol)为中心的模糊区域
        for (int inRow = outRow - BLUR_SIZE; inRow <= outRow + BLUR_SIZE; ++inRow) {
            for (int inCol = outCol - BLUR_SIZE; inCol <= outCol + BLUR_SIZE; ++inCol) {
                // 关键：对每次输入访问进行边界检查！
                if (inRow >= 0 && inRow < height && inCol >= 0 && inCol < width) {
                    int idx = inRow * width + inCol;
                    sum += image[idx];
                    pixelCount++;
                }
            }
        }
        // 计算平均值并写入输出
        int outIdx = outRow * width + outCol;
        blurred[outIdx] = (unsigned char)(sum / pixelCount);
    }
}

重要规则：在并行编程中，每次内存访问都应有对应的边界检查。对于模糊操作，即使输出像素在边界内，其读取的周围输入像素也可能越界，必须进行检查。

示例三：矩阵乘法 ➕✖️➡️📈

矩阵乘法 C = A × B 中，输出矩阵C的每个元素是矩阵A的一行与矩阵B的一列的点积。

并行化策略

最直接的并行化方法是为输出矩阵C的每个元素分配一个线程。每个线程负责计算一个点积。

以下是矩阵乘法内核的实现（假设A、B、C都是N x N矩阵）：

__global__ void matrixMulKernel(float* A, float* B, float* C, int N) {
    // 计算线程对应的输出元素位置
    unsigned int row = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned int col = blockIdx.x * blockDim.x + threadIdx.x;

    // 边界检查
    if (row < N && col < N) {
        float sum = 0.0f;
        // 执行点积计算
        for (int i = 0; i < N; ++i) {
            sum += A[row * N + i] * B[i * N + col];
        }
        // 将结果写入输出矩阵
        C[row * N + col] = sum;
    }
}

性能提示：这个基础实现（通常称为“朴素”矩阵乘法）在GPU上已经能带来显著加速（例如数十到上百倍），因为它暴露了巨大的数据并行性（N²个独立的点积计算）。后续课程将介绍如何通过共享内存等技术进一步优化。

总结与关键要点 🎯

本节课中我们一起学习了GPU上处理多维数据的核心方法。

以下是关键知识点总结：

多维网格：使用 dim3 类型定义网格(numBlocks)和线程块(numThreadsPerBlock)的维度，简化了多维数据的处理。
索引计算：线程通过 blockIdx, threadIdx, blockDim 等内置变量计算其在多维数据中的全局位置（行、列）。
数据布局：理解C语言中多维数据以行主序连续存储至关重要，索引转换公式为 index = row * width + col。
边界检查：这是并行编程中最容易出错的部分之一。必须确保每个线程的内存读写操作都在数组有效范围内。规则是：每次内存访问都应有对应的索引边界检查。
并行策略：对于像图像处理（像素级操作）和矩阵乘法（元素级操作）这类问题，常见的模式是“一个输出元素一个线程”。
性能：即使使用这些基础的多维网格技术，也能在GPU上获得相对于CPU的显著性能提升，因为GPU能够同时启动成千上万个线程来利用数据并行性。

通过RGB转灰度、图像模糊和矩阵乘法这三个例子，我们实践了从问题分析、并行策略制定到CUDA内核实现的全过程。掌握这些基础是进行更复杂GPU编程的基石。

GPU计算：04：GPU架构

概述

在本节课中，我们将暂时告别代码编写，深入探讨GPU的硬件架构。我们将了解GPU如何组织其计算单元（流式多处理器SM），线程块和线程如何被调度到这些硬件上执行，以及SIMD（单指令多数据）执行模型如何影响程序性能。理解这些底层原理对于编写高效的CUDA程序至关重要。

回顾：多维数组与数据并行

上一节我们学习了如何创建多维网格，存储和访问多维数据，并编写了三个处理二维数据的并行程序示例。

以下是上一节的核心内容总结：

RGB转灰度图像：我们为输出图像的每个像素分配一个线程。这展示了如何创建二维网格（包含x和y维度的块），以及块本身也可以是二维的线程数组。我们使用了blockIdx、threadIdx和blockDim的x、y分量。
数据布局：在C语言中动态分配多维数组时，需要将其线性化为行优先的一维数组进行存储。
边界检查：处理多维网格时，需要在多个维度上进行边界检查。
图像模糊：此例的边界检查更为复杂，因为输入和输出的边界条件不同。
矩阵乘法：我们为输出矩阵的每个元素分配一个线程，该线程顺序遍历矩阵A的对应行和矩阵B的对应列进行计算。

GPU架构基础

现在，让我们开始了解GPU的硬件架构。

GPU由多个流式多处理器组成，我们简称其为SM。每个SM包含多个核心，这些核心是执行算术运算的单元。同一个SM内的核心共享控制逻辑（用于取指、译码等）和一些内存资源。此外，所有的SM都能访问同一块全局内存，也就是我们在第一讲中从CPU复制数据到GPU时所用的内存。

以本课程使用的Volta V100 GPU为例，它拥有80个SM，每个SM有64个核心，总计5120个核心。

线程块到SM的映射

我们编写的CUDA内核配置了包含线程块和线程的网格。那么，这些线程块和线程是如何在GPU上运行的呢？

其工作方式是，线程以块为粒度被分配到SM上。这意味着：

同一个线程块内的所有线程将被分配到同一个SM上执行。
一个SM可以同时容纳多个线程块。
但一个线程块不能被拆分到多个SM上。

线程在执行时需要资源，例如寄存器和一些线程特定的控制数据。由于SM的硬件资源是有限的，因此它能同时支持的线程和线程块数量也是有限的。这在一定程度上解释了为什么线程块中的线程数量存在上限（例如1024）。如果一个网格启动的线程块数量超过了GPU能同时执行的数量，多余的线程块将排队等待，直到有SM上的资源被释放。

线程块内的协作与同步

由于同一个线程块内的线程被分配到同一个SM上，它们能够以不同线程块间无法实现的方式进行协作。这也是我们将网格划分为线程块的重要原因之一。

以下是线程块内协作的主要方式：

屏障同步：线程可以调用__syncthreads()函数。这会使得块内的所有线程在代码中的某个点相互等待，直到所有线程都到达该屏障后，才能继续执行。这有助于协调线程间的工作。
共享内存：每个SM上都有一块共享内存，同一线程块内的线程可以快速访问和共享数据，而不同块间的线程则不能。这是一种高效的协作方式。

将同一线程块的线程分配到同一个SM上，使得上述协作（如同步）能够更高效地实现，因为它们无需跨SM进行通信。

一个重要的原则是：线程块必须整体地分配到SM上。也就是说，SM必须有足够的资源容纳块内的所有线程，才能开始执行该块。这是为了避免死锁。例如，如果只分配了部分线程，而这些线程在屏障处等待尚未分配的其他线程，就会导致程序无法继续执行。

线程块间的独立性与可扩展性

与线程块内的紧密协作相反，不同线程块之间不应进行同步或直接协作。CUDA编程模型刻意不允许块间同步（没有内置的跨块同步原语）。

这种设计带来了一个关键优势：线程块之间是相互独立的。这意味着：

线程块可以以任意顺序执行。
它们可以并行执行，也可以顺序执行。

这种独立性实现了透明可扩展性。同一段CUDA代码可以在拥有不同数量SM（即不同硬件并行度）的设备上运行。在SM较少的设备上，线程块可能顺序执行；在SM较多的设备上，更多的线程块则可以并行执行，从而自动利用更多的硬件资源，而无需修改代码。

程序员不应编写试图让不同线程块同步的代码（例如通过全局内存进行自旋锁等），因为这破坏了块间的独立性假设。如果尝试同步的线程块并未被同时调度（例如一些在排队），就极有可能导致死锁。

Warp：SM上的调度单元

线程块被分配到SM后，还会被进一步划分为更小的单元，称为Warp。Warp是SM上进行调度的基本单位。

关于Warp需要了解的是：

Warp的大小是设备特定的，但至今一直是32个线程。
一个包含1024个线程的块会被分成32个Warp（32线程/个）。
一个包含64个线程的块会被分成2个Warp。

SIMD执行与Warp

Warp的特殊之处在于，其内部的线程按照SIMD模型一起被调度执行。

SIMD代表单指令，多数据。这意味着：

一个Warp中的所有线程在同一个时钟周期内执行相同的指令。
但每个线程操作的是不同的数据。

这种模式的优势在于摊销控制开销。只需要一套取指、译码和分发逻辑，就可以驱动32个核心同时工作，极大地节省了硬件资源，使得我们可以将更多的晶体管用于增加计算核心，而非控制单元。

下图展示了Volta V100 GPU中一个SM的详细架构。可以看到，64个FP32核心被组织成4个处理块，每个处理块有16个核心，并共享一个指令分发单元和一个寄存器文件。一个Warp（32线程）会在这样一个处理块上执行。

控制流分歧

SIMD模型也带来了一个显著的劣势：控制流分歧。

当Warp中的线程遇到条件分支（如if语句）时，如果部分线程走then路径，另一部分走else路径，就会发生分歧。由于Warp必须执行相同的指令，硬件会如何处理呢？

硬件会串行化所有不同的执行路径：

首先，所有线程一起执行then路径，但只有条件为真的线程是活跃的，条件为假的线程被禁用（核心闲置）。
然后，所有线程再一起执行else路径，此时之前活跃的线程被禁用，之前禁用的线程变为活跃。
最后，所有线程在分支后重新汇合，继续同步执行。

这导致了SIMD效率的降低，即部分计算核心在某个时刻处于闲置状态。分歧越严重，效率损失越大。

另一个例子是循环次数不同的线程。所有线程必须执行完循环次数最多的那个线程所需的迭代后，才能一起退出循环，这也会造成大量闲置。

程序员可以通过重构算法或数据布局（例如，预先对数据进行排序，让具有相似执行路径的线程聚集在同一个Warp内）来减轻控制流分歧的影响。

延迟隐藏与多线程

GPU核心可能会遇到长延迟操作，例如未命中缓存的内存访问或复杂的多周期算术运算。为了避免核心流水线停滞，GPU SM采用了基于Warp的多线程技术。

其工作原理如下：

当正在执行的Warp遇到长延迟操作时，它会被换出核心。
调度器立即从已分配到该SM的、就绪的Warp中挑选一个换入核心执行。
当第一个Warp的延迟操作完成（如数据从内存返回），它重新变为就绪状态，等待被调度执行。

通过这种方式，SM的流水线始终保持忙碌，从而隐藏了操作延迟，提高了整体吞吐量。

占用率

为了有效隐藏延迟，我们希望在每个SM上同时驻留尽可能多的线程（即多个Warp）。占用率衡量了SM上实际活跃的Warp（或线程）数量与硬件支持的最大数量之比。

通常，更高的占用率有助于更好地隐藏延迟。然而，占用率可能受到以下因素限制：

线程块大小：如果线程块大小设置不当，可能无法充分利用SM的线程容量。
每个SM的线程块数量上限：SM能同时容纳的块数量有限。
寄存器使用量：每个线程使用的寄存器数量。如果线程使用过多寄存器，SM的寄存器文件可能无法支持最大数量的线程。
共享内存使用量：每个线程块使用的共享内存量。

例如，在V100上（最大2048线程/SM，最大32块/SM）：

选择256线程/块：需要8个块即可达到2048线程，且未超过32块上限，可实现高占用率。
选择32线程/块：需要64个块才能达到2048线程，但受限于32块上限，最多只能有1024线程，占用率仅为50%。
选择768线程/块：2048无法被768整除，只能容纳2个块（1536线程），剩余线程槽位浪费，占用率也非最优。

因此，选择线程块大小时需要考虑目标GPU的硬件限制，以优化占用率。CUDA提供了cudaGetDeviceProperties API，用于查询设备的这些属性。

总结

本节课我们一起深入学习了GPU的架构。我们了解了GPU由多个SM组成，线程块如何映射到SM上执行，以及Warp作为基本调度单元遵循SIMD模型所带来的优势（控制开销摊销）和挑战（控制流分歧）。我们还探讨了如何通过基于Warp的多线程来隐藏操作延迟，以及占用率的概念及其对性能的影响。理解这些硬件特性是后续进行CUDA程序性能分析和优化的基础。

GPU计算：第5讲：内存与分块

在本节课中，我们将学习GPU的内存架构、CUDA编程模型中的内存访问方式，以及一个重要的内存优化技术——分块（Tiling）。我们将以矩阵乘法为例，详细讲解如何利用共享内存来减少全局内存访问，从而提升程序性能。

概述：性能指标与计算瓶颈

处理器设计者通常通过一些指标来告知用户处理器的性能。最值得注意的两个指标是：

峰值浮点运算速率（FLOPS）：处理器每秒能执行的浮点运算次数。这反映了处理器核心的计算能力。
峰值内存带宽：处理器的内存每秒能向核心提供的数据字节数。这反映了内存系统的数据传输能力。

以V100 GPU为例，其峰值浮点运算速率约为14 TeraFLOPS（14万亿次/秒），峰值内存带宽约为900 GB/s（9000亿字节/秒）。这些是理论峰值，实际程序性能通常低于此值，但它们为评估代码性能提供了重要参考。

根据程序的特点，其性能可能受限于计算能力或内存带宽：

计算受限（Compute Bound）：程序性能受限于处理器的浮点运算速率。此时，计算核心始终满负荷工作，而内存系统相对空闲。
内存受限（Memory Bound）：程序性能受限于内存带宽。此时，计算核心经常因等待数据而空闲，内存系统成为瓶颈。

为了充分利用GPU强大的计算能力，我们希望程序是计算受限的。这引出了一个关键指标：期望的计算与全局内存访问比。其公式为：
期望计算/内存访问比 = 峰值FLOPS / 峰值内存带宽
对于V100，这个比值约为 14 TeraFLOPS / 900 GB/s ≈ 15.6 次操作/字节。这意味着，平均每从内存加载1字节数据，需要执行约15.6次浮点运算，才能让计算核心完全“吃饱”，达到峰值性能。

接下来，我们通过两个例子来分析这个比值。

向量加法的计算/内存比

向量加法的内核代码如下：

z[i] = x[i] + y[i];

每个线程执行1次浮点加法，但需要从全局内存加载两个float类型的数据（共8字节）。因此，其计算/内存比为 1次操作 / 8字节 = 0.125 次操作/字节。这个值远低于V100所需的15.6，因此向量加法是典型的内存受限型应用。这也是为什么之前我们看到，即使GPU有大量ALU，向量加法的加速比也并不惊人的原因。

矩阵乘法的潜力与问题

我们之前实现的朴素矩阵乘法内核代码如下：

for (int k = 0; k < N; ++k) {
    sum += A[row * N + k] * B[k * N + col];
}

每个线程在循环的每次迭代中执行1次乘法和1次加法（共2次浮点运算），并加载2个float（共8字节）。因此，其计算/内存比为 2次操作 / 8字节 = 0.25 次操作/字节，虽然比向量加法高，但仍远低于目标。

然而，矩阵乘法本身具有很高的数据复用潜力。计算一个N×N的矩阵乘法，总共需要执行约2N³次浮点运算，但理论上只需加载2 * N² * 4字节的输入数据（假设每个数据只加载一次）。因此，其潜在的计算/内存比可达 (2N³次操作) / (8N²字节) = 0.25N 次操作/字节。当N较大时，这个值可以很高。

问题在于我们朴素的实现方式没有利用这种复用。例如，计算输出矩阵中同一行的不同元素时，需要重复加载矩阵A的同一行数据；计算同一列的不同元素时，需要重复加载矩阵B的同一列数据。这造成了大量冗余的全局内存访问。

上一节我们介绍了性能瓶颈的概念，并看到矩阵乘法有巨大的优化空间。本节中，我们将深入GPU内存架构，并学习如何通过“分块”技术来挖掘这种潜力。

GPU内存架构与CUDA内存模型

在深入优化之前，我们需要了解GPU的内存是如何组织的，以及CUDA编程模型为我们提供了哪些内存管理工具。

GPU内存层次结构

GPU的内存是一个层次化结构，访问速度和容量各不相同：

全局内存（Global Memory）：所有流多处理器（SM）共享的大容量、高延迟（约数百周期）内存。主机与设备之间的数据拷贝通常发生在这里。
L2缓存（L2 Cache）：位于芯片上，所有SM共享，用于缓存全局内存数据。
流多处理器（SM）内部：
- 寄存器（Registers）：线程私有，访问速度极快（通常单周期）。
- L1缓存/共享内存（L1 Cache / Shared Memory）：在物理上通常是统一的资源，可配置分配。访问延迟约数个周期。
  - L1缓存：由硬件自动管理，缓存频繁访问的全局内存数据。
  - 共享内存：由程序员显式管理，是块内线程可共享的快速内存。
- 常量缓存（Constant Cache）：用于缓存只读的常量数据。

CUDA编程模型中的内存

在CUDA编程模型中，我们可以通过特定的限定符来指定变量的存储位置和生命周期：

__device__：变量位于全局内存。作用域为整个网格（所有线程），生命周期为整个应用程序。
__constant__：变量位于常量内存。作用域为整个网格，生命周期为整个应用程序。
__shared__：变量位于共享内存。每个线程块拥有该变量的一个独立副本，块内所有线程共享此副本。其作用域为线程块，生命周期为块执行期间。
局部变量：在kernel函数内声明，无特殊限定符。通常存储在寄存器中，每个线程拥有私有副本，作用域和生命周期限于线程内部。

了解内存模型后，我们就可以利用共享内存来优化程序了。接下来，我们将聚焦于共享内存，并学习如何用它来实现分块优化。

分块优化与共享内存实践

分块（Tiling）是一种经典的内存优化技术，其核心思想是：将输入数据分割成小块（Tile），先将这些小块从慢速的全局内存加载到快速的共享内存中，然后让线程块内的所有线程协作，从共享内存中重复访问这些数据，从而减少对全局内存的访问次数。

矩阵乘法的分块策略

回顾我们之前对矩阵乘法的并行化：将输出矩阵划分为小块，每个线程块负责计算一个小块，块内的每个线程负责计算小块中的一个元素。

在朴素实现中，计算同一输出块的所有线程，会重复加载矩阵A的相同行和矩阵B的相同列。分块优化旨在改变这一点：

将输入矩阵也分块：我们将矩阵A和B也划分为与输出块尺寸相匹配的小块。
协作加载：线程块内的所有线程协作，将当前计算所需的一小块A和一小块B从全局内存加载到共享内存中。
共享内存计算：所有线程同步，确保数据加载完毕。然后，每个线程从共享内存中读取数据，进行部分点积计算。
循环推进：处理完当前小块后，线程再次同步，然后协作加载下一对小块，重复上述过程，直到计算完所有数据。

这样，每个线程对全局内存的访问次数从原来的O(N)次减少到O(N / TileWidth)次，而增加的是对共享内存的访问。由于共享内存速度快得多，因此能显著提升性能。

分块矩阵乘法代码实现

以下是分块矩阵乘法内核代码的关键步骤解析。我们假设分块大小TILE_DIM为32，线程块大小也是32x32。

步骤1：声明共享内存数组

__shared__ float As[TILE_DIM][TILE_DIM];
__shared__ float Bs[TILE_DIM][TILE_DIM];

每个线程块将拥有As和Bs这两个共享内存数组的独立副本，用于存储当前正在处理的输入数据块。

步骤2：循环遍历所有数据块

for (unsigned int tile = 0; tile < (N / TILE_DIM); ++tile) {
    // 加载数据块到共享内存
    // 从共享内存计算部分点积
}

外层循环遍历所有需要的数据块对。

步骤3：协作加载数据到共享内存

// 每个线程加载一个元素到共享内存中
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_39.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_41.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_42.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_44.png)

As[threadIdx.y][threadIdx.x] = A[row * N + (tile * TILE_DIM + threadIdx.x)];
Bs[threadIdx.y][threadIdx.x] = B[(tile * TILE_DIM + threadIdx.y) * N + col];

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_46.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_48.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_50.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_52.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_54.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_56.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_58.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/3b83afc778a0465a3d1b03e67ec73347_60.png)

__syncthreads(); // 确保所有线程都完成加载

线程根据其在线程块内的位置（threadIdx），负责加载输入块中对应位置的一个元素。
__syncthreads()是一个屏障同步函数，调用它的所有线程会在此等待，直到同一线程块内的所有线程都执行到此位置，才能继续向下执行。这确保了在计算开始前，共享内存中的数据已准备就绪。

步骤4：从共享内存计算部分点积

float sum = 0.0f;
for (unsigned int i = 0; i < TILE_DIM; ++i) {
    sum += As[threadIdx.y][i] * Bs[i][threadIdx.x];
}
__syncthreads(); // 确保所有线程都完成对此数据块的计算

每个线程利用已加载到共享内存中的小块数据，计算其负责的输出元素的部分和。
计算完成后再次同步，这是为了防止某些线程过快地进入下一次循环加载新数据，而其他线程还在使用共享内存中的旧数据。

步骤5：循环结束后存储结果

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

在完成所有数据块的计算后，将累加和写入全局内存中的输出矩阵。

边界条件处理

在分块实现中，边界条件处理需要格外小心。因为线程块的大小是固定的（如32x32），但当矩阵维度不是分块大小的整数倍时，位于边缘的线程块可能只有部分线程参与有效计算。
然而，在加载输入块时，即使某些线程不参与最终输出计算，它们也可能需要参与数据加载工作，以确保共享内存中被有效线程需要的数据被正确加载。因此，不能简单地用if (row < N && col < N)来禁用所有无效线程，而需要更精细地控制不同阶段（加载、计算、存储）哪些线程是活跃的。

性能提升与注意事项

应用分块优化后，矩阵乘法的性能通常能得到显著提升（例如2倍或更多），因为它大幅减少了昂贵的全局内存访问。

使用共享内存时还需注意：

资源限制：每个SM的共享内存容量有限。使用过多共享内存可能会降低占用率（Occupancy），即同时活跃的线程块数量，从而影响利用线程隐藏延迟的能力。
动态共享内存：除了示例中的静态声明，共享内存大小也可以在内核启动时动态指定，这为处理不同大小的分块提供了灵活性。

CPU上的分块优化

有趣的是，分块优化并非GPU专属，在CPU上同样有效。CPU虽然缺乏程序员可管理的“共享内存”，但其大容量缓存（L1、L2、L3）同样受益于分块技术。通过将循环重构成对数据块的遍历，可以提升数据的时间局部性，使得正在处理的数据更有可能驻留在高速缓存中，减少对主内存的访问。

CPU上的分块矩阵乘法代码结构会包含多层循环：外层循环遍历输出行块、输出列块和输入块，内层循环遍历块内的行、列以及点积计算。其逻辑与GPU内核中的循环结构有清晰的对应关系。

总结

本节课我们一起深入学习了GPU内存相关的关键知识：

性能指标与瓶颈：理解了峰值FLOPS、内存带宽、计算受限与内存受限的概念，以及期望的计算/内存访问比。
内存架构与模型：了解了GPU的层次化内存结构（全局内存、共享内存、寄存器等）及其在CUDA编程模型中的对应（__device__, __shared__等限定符）。
分块优化：掌握了利用共享内存进行分块优化的核心思想。通过让线程块内的线程协作，将数据从全局内存加载到共享内存并复用，从而减少全局内存访问，提升性能。我们以矩阵乘法为例，详细分析了代码实现步骤，包括共享内存声明、协作加载、同步以及边界条件处理。
优化通用性：认识到分块是一种通用的数据局部性优化技术，在CPU上通过改善缓存利用率也能带来性能提升。

通过本课的学习，你不仅知道了如何编写更高效的CUDA内核，也加深了对计算机体系结构中“内存墙”问题以及通过数据复用提升性能这一普适思想的理解。

GPU计算：P6：性能考量 🚀

在本节课中，我们将继续探讨GPU计算的性能考量。我们将回顾之前提到的内存优化，并引入两个新的关键优化概念：内存合并与线程粗化。理解这些概念对于编写高效的GPU程序至关重要。

概述 📋

上一节我们介绍了GPU架构中的内存层次结构以及CUDA编程模型，并学习了如何使用共享内存分块技术来优化内存访问，并将其应用于矩阵乘法示例。

本节中，我们将首先深入了解DRAM架构的工作原理，这将引出内存合并优化的重要性。接着，我们将探讨线程粒度与线程粗化，这是一种通过让单个线程处理更多工作来减少冗余、提升性能的技术。

DRAM架构回顾 🧠

为了理解内存合并，我们需要先了解DRAM（动态随机存取存储器）的基本工作原理。

一个DRAM单元可以逻辑上视为一个存储电荷的电容器和一个三态器件。当启用三态器件时，电容器放电以读取存储的值。

一个DRAM阵列由多个连接到列线的DRAM单元组成。在任意时刻，只能读取一行中的单元。我们通过行地址来激活特定行，该行所有单元的值会被读出并存储到感应放大器和列锁存器中。

整个被读取的行称为一个“突发传输”。内存地址分为行地址和列地址。行地址用于选择激活哪一行，列地址则通过多路复用器从已读出的整行数据（突发）中选择我们请求的特定部分。

从DRAM阵列读取数据、进行感应放大并写回（因为读取会破坏数据）的过程是缓慢的。而使用多路复用器从已读出的突发中选择数据则相对较快。

关键点：访问同一突发内的不同数据速度很快，因为只需改变多路复用器的选择。而访问不同突发中的数据则需要重新执行耗时的行激活和读取过程。

内存合并优化 ⚡

基于DRAM的工作原理，我们引出了内存合并的概念。

当同一个Warp中的线程访问连续的、位于同一DRAM突发内的内存位置时，这些访问可以被合并，并由一次DRAM事务完成服务。这被称为内存合并。

反之，如果同一个Warp中的线程访问的内存位置分散在不同的DRAM突发中，则这些访问无法合并，需要发起多次内存事务。这被称为内存发散，会显著降低性能，因为Warp必须等待所有独立的内存加载完成才能继续执行。

在我们之前编写的代码中，内存访问模式天然地支持了合并：

向量加法：线程0到31访问数组x和y中连续的元素0到31，这些元素在内存中相邻，很可能位于同一突发内。
矩阵乘法：在分块和非分块版本中，通过合理的线程块组织（Warp包含连续的threadIdx.x线程），对矩阵A和B的全局内存访问也是合并的。

确保内存访问模式是合并的，是GPU编程的一项基本优化。后续课程中，我们会遇到需要手动调整以实现内存合并的例子。

多DRAM存储体与延迟隐藏 🏦

现代内存系统通常由多个DRAM存储体（Bank）组成。这带来了并行访问的可能性，有助于隐藏内存访问延迟。

如果只有一个DRAM存储体，即使在充分利用突发传输的情况下，在读取完一个突发后，仍需等待下一个突发的读取过程，存在空闲时间。

拥有多个DRAM存储体时，当一个存储体正在服务数据时，其他存储体可以并行地进行下一轮耗时的行激活和读取操作。这样，当处理器处理完一个存储体的数据时，另一个存储体的数据可能已经准备就绪，从而有效地“隐藏”了内存访问延迟。

为了有效利用多存储体并实现延迟隐藏，我们需要同时发起大量的内存访问请求。这再次强调了高占用率的重要性。高占用率不仅能隐藏核心流水线的延迟，还能通过提供足够多的并发内存请求来隐藏访问DRAM的延迟。

线程粒度与线程粗化 🧵

到目前为止，我们的并行化策略是让线程尽可能“细粒度”，即每个线程负责最小的并行工作单元（如一个向量元素、一个输出像素、一个输出矩阵元素）。

细粒度线程的优势在于为硬件提供了尽可能多的线程，以充分利用所有计算资源，并实现透明可扩展性：硬件可以根据自身资源多少，自动并行化或串行化这些线程。

然而，细粒度线程也可能带来冗余。例如，在分块矩阵乘法中，水平相邻的两个线程块会加载相同的A矩阵块。如果这两个线程块能在不同SM上并行执行，这种冗余是换取并行性可接受的代价。但如果由于硬件资源有限（例如SM数量少），这两个线程块最终被硬件串行执行，那么这种冗余加载就变成了不必要的开销。

线程粗化优化正是为了解决这个问题。其核心思想是让一个线程（或线程块）顺序处理多个输出单元，从而重用已加载的数据或已进行的计算，减少并行化带来的冗余开销。

以下是线程粗化的实现思路，以矩阵乘法为例：

定义粗化因子（例如 COARSE_FACTOR = 4），表示每个线程负责的输出元素数量。
调整线程块网格维度，将X方向的线程块数量除以粗化因子。
在线程内部，使用数组（大小为粗化因子）来存储多个部分和。
在加载输入块的主循环中，外层循环遍历A的块，内层循环遍历B的多个相邻块（数量等于粗化因子）。这样，加载一次的A块可以被重用于计算多个输出列。
最终，线程将其负责的多个结果写回全局内存。

线程粗化减少了并行化的开销（如冗余内存加载、计算或同步），但同时也增加了每个线程的资源使用（如寄存器、本地内存）。如果粗化过度，可能导致占用率下降，反而损害性能。此外，最优的粗化因子通常依赖于具体的硬件设备，这在一定程度上牺牲了透明可扩展性。

优化清单与权衡 ✅

至此，我们学习了一系列常见的GPU优化技术：

调整资源使用以最大化占用率
最小化控制发散以提高SM效率
确保内存合并访问模式
使用共享内存分块以捕获数据重用
应用线程粗化以减轻并行化开销

需要注意的是，这些优化之间有时存在权衡。例如：

最大化占用率可能加剧缓存抖动。
使用过多共享内存会限制占用率。
线程粗化减少冗余工作但增加单线程资源消耗，可能限制占用率。

因此，在应用优化时，诊断性能瓶颈至关重要。我们需要明确是内存带宽、计算能力、控制发散还是其他因素限制了程序性能，从而有针对性地选择和应用优化策略，找到最佳的平衡点。

总结 🎯

本节课我们一起深入探讨了GPU性能优化的两个重要方面。
首先，我们从DRAM架构出发，理解了内存合并的原理与重要性，它确保了Warp内线程的高效内存访问。
其次，我们引入了线程粗化的概念，通过让线程处理更多工作来减少并行化带来的冗余开销，但需要注意其与占用率之间的权衡。
最后，我们回顾了目前已学的优化技术清单，并强调了根据实际性能瓶颈进行针对性优化的重要性。掌握这些核心概念，将为后续学习更复杂的并行模式打下坚实基础。

GPU计算：P7：性能剖析 🚀

在本节课中，我们将学习如何剖析CUDA应用程序，以了解其性能表现并识别瓶颈所在。

概述

上一节我们介绍了DRAM存储体、内存合并与分散、线程束合并等概念，并探讨了不同优化策略（如资源调优、内存访问模式、共享内存分块、线程合并）之间的权衡关系。为了在这些相互竞争的优化中找到最佳平衡点，关键在于识别应用程序的性能瓶颈。本节我们将学习如何使用NVIDIA提供的剖析工具来诊断瓶颈。

性能瓶颈与优化权衡

一个瓶颈是指限制应用程序性能的约束条件。不同的应用程序消耗不同的硬件资源（如计算单元、内存带宽、线程槽位、共享内存等），当其中一种资源先于其他资源耗尽时，便形成了瓶颈。瓶颈取决于具体的应用程序和硬件设备。

我们之前学到的优化策略，通常是以一种资源换取另一种资源。例如：

共享内存分块：使用更多共享内存，以减少全局内存带宽的使用。
线程合并：减少冗余工作（如数据加载、计算或同步开销），但可能增加每个线程所需的资源，从而限制占用率。

因此，在应用优化之前，必须准确诊断当前的瓶颈。如果应用程序的瓶颈是共享内存限制导致的低占用率，那么增加共享内存使用量的优化（如更大的分块）反而会进一步降低性能。正确的做法可能是减少共享内存的使用。

使用 `nvprof` 进行基础剖析

nvprof 是CUDA SDK中自带的命令行剖析工具。最简单的使用方式是运行 nvprof [你的可执行文件]。

以下是一个向量加法示例的剖析步骤：

编译程序：确保代码中包含了 cudaProfilerStop() 调用，以便在程序退出前保存所有剖析信息。
```
nvcc -o vector_add vector_add.cu
```
运行 nvprof：
```
nvprof ./vector_add
```
该命令会输出一个表格，显示应用程序中各个活动（如内核执行、内存拷贝、API调用）所花费的时间、调用次数、占总时间的百分比等信息。这比在代码中手动插入计时器更方便。

生成剖析文件以供深入分析

nvprof 可以将详细的剖析数据导出到文件，以便用可视化工具进行深入分析。

以下是生成不同类型剖析文件的命令：

生成时间线文件：记录各活动的执行时间线。
```
nvprof -o timeline.prof ./vector_add
```
生成性能指标文件：收集诸如占用率、缓存命中率、内存带宽等硬件计数器数据。使用 --metrics all 收集所有可用指标（注意：这会运行内核多次，以收集不同组的指标，因此总执行时间会变长）。
```
nvprof --metrics all -o metrics.prof ./vector_add
```

使用 NVIDIA Visual Profiler (nvvp) 进行可视化分析

nvvp 是一个图形化剖析工具，可以更直观地分析 nvprof 生成的 .prof 文件。

启动 nvvp：
```
nvvp &
```
导入剖析文件：在nvvp中创建新会话，选择“导入”由 nvprof 生成的数据。分别指定时间线文件（如 timeline.prof）和指标文件（如 metrics.prof）。

分析时间线

时间线视图直观展示了GPU上各种活动（如内存拷贝、内核执行）的起止时间和重叠情况。这有助于理解应用程序的整体执行流程和时间分布。

分析内核性能

在nvvp中，可以针对特定的内核进行深入分析。

检查GPU使用情况：工具会给出高层次建议，例如指出内存拷贝效率低、内核并发性不足等。
执行内核分析：
- 利用率分析：显示计算单元和内存系统的利用率。例如，在向量加法中，内存系统利用率很高（>80%），而计算单元利用率很低（~10%），这清楚地表明该内核是内存带宽受限型。
- 指令构成：显示内核执行中算术指令、控制流指令和内存指令的比例，这与代码逻辑相符（向量加法包含较多加载/存储指令和少量算术、分支指令）。

针对性瓶颈分析

根据初步分析，nvvp会推荐进行更具体的分析：

内存带宽分析：展示各级内存（共享内存、L2缓存、全局内存等）的利用率。对于内存带宽受限的内核，你会看到设备内存（全局内存）利用率接近饱和，而缓存利用率很低。
计算分析：展示各功能单元（如加载存储单元、单精度浮点单元）的利用率以及各类指令的执行计数。对于计算受限的内核，这里会显示计算单元的高利用率。
延迟分析：重点分析占用率。工具会显示当前占用率、理论最大占用率以及限制占用率的资源（如线程块大小、共享内存用量）。它还提供了交互式图表，可以预测更改线程块大小或共享内存使用量对占用率的影响。
PC采样（如果设备支持）：可以揭示线程束停顿的主要原因（例如，等待内存依赖、执行依赖、同步等）。这对于理解底层延迟来源非常有帮助。

案例分析：改变线程块大小的影响

让我们通过一个实验来验证占用率对性能的影响。将向量加法的线程块大小从512改为32。

修改并编译代码：将 blockSize 设置为32。

重新剖析：

nvprof -o timeline_smallblock.prof ./vector_add_smallblock
nvprof --metrics all -o metrics_smallblock.prof ./vector_add_smallblock

在nvvp中分析：导入新的剖析文件。你会发现：
- 内核执行时间显著增加（性能下降）。
- 占用率分析显示，理论最大占用率降至50%（受限于每个SM的最大线程块数量），实际占用率约为35%。
- 延迟分析会明确指出“GPU利用率可能受限于线程块大小”，并高亮显示 Threads Per Block 是限制因素。

这个案例表明，不恰当的线程块配置会导致GPU硬件资源未被充分利用，从而成为性能瓶颈。

总结

本节课我们一起学习了GPU性能剖析的核心方法。我们了解到识别性能瓶颈对于实施有效优化至关重要。通过使用 nvprof 命令行工具和 nvvp 可视化剖析器，我们可以：

量化应用程序中各个部分的时间消耗。
判断内核是受限于内存带宽、计算能力还是其他资源。
分析占用率，并确定限制占用率的因素（如寄存器用量、共享内存、线程块大小）。
通过PC采样理解线程束停顿的原因。

掌握这些剖析技能，能够帮助你在优化CUDA程序时做出数据驱动的决策，避免盲目应用可能适得其反的优化策略，从而更高效地提升程序性能。

GPU计算：第8讲：卷积

在本节课中，我们将开始学习并行模式，第一个要讨论的模式是卷积。在此之前，我们先快速回顾一下之前课程的内容。

课程回顾

到目前为止，我们已经学习了GPU计算的基础知识。我们探讨了单线程性能因频率停滞和功耗墙而停滞不前，这推动了并行计算的发展。我们比较了CPU（面向延迟的设计）和GPU（面向吞吐量的设计）的架构差异。

CPU（面向延迟）：专注于让单个任务尽可能快，使用强大的ALU、大容量缓存、复杂的控制流（如分支预测、数据转发）和少量多线程来隐藏剩余延迟。
GPU（面向吞吐量）：拥有大量较简单的ALU、较小的缓存和更简单的控制逻辑，通过海量线程并行执行来隐藏操作延迟，从而获得高吞吐量。

我们了解了典型的GPU系统架构：CPU拥有主内存，GPU拥有设备内存。通常的流程是在GPU上分配内存，将数据从CPU复制到GPU，在GPU上执行内核（访问设备内存），然后将结果复制回CPU，最后释放GPU内存。我们也提到了统一内存等更便捷的方式。

我们从简单的向量加法内核开始，学习了数据并行性，以及如何组织网格（Grid）、线程块（Block）和线程（Thread）。我们学习了如何计算线程索引，以便每个线程对不同的数据执行相同的操作。

随后，我们学习了网格可以组织成多维数组（如2D、3D网格），以及多维数据通常按行主序（Row-Major Order）存储。我们通过多个例子（如RGB转灰度、图像模糊、矩阵乘法）练习了如何将多维索引转换为一维索引以访问动态分配的数组。

在介绍了CUDA编程模型的基础后，我们转向GPU架构。GPU由多个流式多处理器（SM）组成，每个SM包含多个共享内存和控制的核，所有SM共享全局内存。网格在SM上以线程块为单位进行调度，同一块内的所有线程在同一SM上同时执行。线程块在SM上进一步划分为线程束（Warp，大小为32），线程束是调度单位，遵循SIMD模型（单指令多数据），即同一线程束内的所有线程执行相同的指令。

这引出了控制发散（Control Divergence）的问题：如果同一线程束内的线程需要执行不同的控制路径（如if-else），硬件会串行化执行所有路径，导致部分线程闲置，降低硬件利用率。我们需要尽量减少这种情况。

我们还讨论了延迟隐藏（Latency Hiding）。通过在SM上调度远多于核心数量的线程（即高占用率），当一个线程束遇到长延迟操作（如内存访问）时，可以切换到另一个就绪的线程束执行，从而隐藏延迟。占用率受限于每个SM的最大线程数、每个块的线程数、寄存器数量和共享内存大小等因素。

接下来，我们研究了GPU内存架构。同一SM上的线程可以访问共享内存（Shared Memory）。我们学习了如何使用共享内存来缓存将被重用的数据，例如在矩阵乘法中，通过线程块协作将输入矩阵的块（Tile）加载到共享内存，减少对全局内存的访问次数，提高计算与全局内存访问的比率。

我们探讨了DRAM的组织方式。访问DRAM阵列较慢，但一旦开始突发传输（Burst），访问突发内的数据就很快。因此，最好让线程访问同一DRAM突发内的数据，以减少访问DRAM阵列的次数。多个DRAM库（Bank）可以并行工作，进一步隐藏延迟。高占用率有助于产生大量内存访问请求，从而更好地隐藏内存访问延迟。

基于以上，我们总结了一个常见优化清单：

调整资源以最大化占用率。
最小化控制发散。
确保合并的内存访问模式。
使用共享内存分块（Tiling）以捕获数据重用。
使用线程粗化（Coarsening）以减轻并行化开销。
私有化（Privatization，后续会讲）。

选择哪种优化取决于应用程序的性能瓶颈。优化通常是用一种资源换取另一种资源，因此需要针对瓶颈进行优化。CUDA提供了性能分析工具来帮助识别瓶颈。

以上回顾结束了GPU计算基础部分。接下来，我们将进入课程的新部分——并行模式。我们将花大量时间讨论这些模式，并希望在每个模式中介绍新的架构特性或优化技术。今天，我们从卷积开始。

什么是卷积？ 🌀

卷积是一种运算。在2D卷积中，输出的每个元素都是其对应输入元素及其邻域元素的加权和。

我们之前见过的图像模糊（Image Blur）就是卷积的一个特例，其中所有权重都相同（计算平均值）。在更一般的卷积形式中，每个邻域输入元素可以有不同的权重。

这些权重由一个卷积掩码（Convolution Mask）（有时被称为卷积核，但为避免与CUDA内核函数混淆，这里称为掩码）决定。掩码包含一组权重，将这些权重应用于输入邻域，即可计算出加权和作为输出。

卷积在信号处理、图像处理、视频处理等领域有广泛应用。它通常用于将信号（一维、二维如图像像素等）转换为更理想的值。例如：

高斯模糊：一种更复杂的模糊，距离中心越远的像素权重越低。
图像锐化
边缘检测

卷积所实现的变换取决于掩码中的权重。本节课我们以2D卷积为例进行讲解，但卷积也有一维和三维形式。

如何并行化卷积？ ⚙️

基于我们目前的知识，一个简单的并行化方法是：为每个输出元素分配一个线程。

这个线程负责遍历对应的输入邻域元素和掩码，计算加权和。这是一种直接的方法，但不是唯一的方法。对于典型的、输出尺寸较大的工作负载（如图像处理），这种方法通常就足够了。

我们注意到掩码的一些特点：

掩码通常很小（例如5x5）。
掩码是常量（在整个内核执行期间不变）。
网格中的所有线程都会访问相同的掩码。

因此，一个优化思路是：将掩码存储在常量内存（Constant Memory）中，以实现更快的访问。

常量内存介绍 🧠

在CUDA编程模型中，线程可以访问寄存器（私有）、共享内存（同一线程块内共享）和全局内存（所有线程可访问）。此外，所有线程还可以访问常量内存。

常量内存具有以下特点：

速度快：访问速度比全局内存快。
容量小：通常只有64KB。
只读：GPU内核执行期间不能写入，但可以从CPU复制数据到其中。

将常量数据（如卷积掩码）放入常量内存的优势在于硬件可以为其建立高效的缓存（常量缓存）。由于数据是只读的，这种缓存无需支持写回、脏位跟踪和缓存一致性协议，实现更简单高效。其容量小也有助于降低缓存未命中率。

每个SM都有自己的常量缓存。虽然程序员也可以手动将常量数据加载到共享内存，但使用常量内存由硬件管理，更加方便。

在CUDA中使用常量内存 💻

以下是如何在代码中声明和使用常量内存。

首先，在全局作用域声明一个常量内存数组：

__constant__ float mask_c[MASK_DIM][MASK_DIM];

这里使用 __constant__ 限定符，并遵循命名约定（例如加 _c 后缀表示常量内存）。

然后，在主机代码中，使用 cudaMemcpyToSymbol 函数将数据从CPU复制到GPU的常量内存：

cudaMemcpyToSymbol(mask_c, mask, MASK_DIM * MASK_DIM * sizeof(float));

其中 mask 是主机上的掩码数组。

实现卷积内核 🛠️

现在，让我们编写卷积内核。思路与图像模糊内核类似：每个线程计算一个输出元素。

首先，计算当前线程对应的输出行列索引：

int out_row = blockIdx.y * blockDim.y + threadIdx.y;
int out_col = blockIdx.x * blockDim.x + threadIdx.x;

接着，进行边界检查，确保线程在输出图像范围内：

if (out_row < height && out_col < width) {
    // 计算逻辑
}

在计算逻辑中，初始化累加器 sum，然后循环遍历掩码的每个元素：

float sum = 0.0f;
for (int mask_row = 0; mask_row < MASK_DIM; ++mask_row) {
    for (int mask_col = 0; mask_col < MASK_DIM; ++mask_col) {
        // 计算对应的输入元素索引
        int in_row = out_row - MASK_RADIUS + mask_row;
        int in_col = out_col - MASK_RADIUS + mask_col;
        // 边界检查：确保输入索引有效
        if (in_row >= 0 && in_row < height && in_col >= 0 && in_col < width) {
            // 累加：输入值 * 掩码权重
            sum += input[in_row * width + in_col] * mask_c[mask_row][mask_col];
        }
    }
}
// 将结果写入输出
output[out_row * width + out_col] = sum;

这里，MASK_RADIUS 是掩码半径（例如，对于5x5掩码，半径为2）。输入索引通过输出索引减去半径再加上掩码内的偏移量得到。边界检查确保不会访问输入图像之外的内存。

关于边界检查引起的控制发散：只有图像边缘的少数线程会经历发散，对整体性能影响不大，通常是可接受的。

编译运行此代码，与CPU版本相比，可以获得显著的加速比。性能良好的原因包括：高度并行、数据重用（邻近线程访问重叠的输入数据）、掩码访问快（常量内存）、以及全局内存访问模式良好（合并访问）。

卷积的数据重用与分块优化 🧩

观察卷积的数据访问模式可以发现显著的数据重用。处理相邻输出元素的线程会访问大量重叠的输入数据。

为了捕获这种重用，我们可以应用之前学过的优化：共享内存分块（Shared Memory Tiling）。

基本思想是：一个线程块负责计算一个输出块（Tile）。该输出块对应的输入块比输出块更大（因为卷积需要邻域信息）。具体来说：

输入块尺寸 = 输出块尺寸 + 2 * 掩码半径

或

输入块尺寸 = 输出块尺寸 + 掩码尺寸 - 1

实现策略如下：

启动与输入块尺寸对应的线程块（即“过度配置”线程）。
所有线程协作，将整个输入块加载到共享内存中（每个线程加载一个元素）。
然后，只有对应于输出块的那些线程（即输入块内部的一个子集）参与计算，从共享内存中读取数据并与掩码进行卷积运算。
对于加载时越界的“幽灵（Ghost）元素”，一个常见的处理技巧是直接在共享内存中将其置为0，这样在计算时就不需要再进行边界检查，简化了逻辑。

这种分块优化可以显著减少对全局内存的访问次数。

计算与全局内存访问比率分析 📊

我们来分析分块优化对“计算操作与全局内存访问字节数比率”的改善。

非分块版本：

每个线程加载 M^2 个输入元素（M为掩码尺寸），每个元素4字节（float），故全局内存访问量为 4 * M^2 字节。
每个线程进行 M^2 次乘加运算，即 2 * M^2 次浮点操作。
比率为 (2 * M^2) / (4 * M^2) = 0.5 操作/字节。

分块版本（以线程块为单位分析）：

设输入块尺寸为 T，输出块尺寸为 T - M + 1。
一个线程块的计算操作：输出块有 (T - M + 1)^2 个元素，每个元素需要 M^2 次乘加，即 2 * M^2 * (T - M + 1)^2 次操作。
一个线程块的全局内存加载：加载整个输入块 T^2 个元素，即 4 * T^2 字节。
比率为 [2 * M^2 * (T - M + 1)^2] / (4 * T^2) = 0.5 * M^2 * (1 - (M-1)/T)^2。

当 M=5, T=32 时，比率约为 9.57 操作/字节，相比非分块版本的 0.5 操作/字节，有近19倍的提升。这清晰地展示了分块优化通过数据重用来改善计算内存比的强大效果。

需要注意的是，如果处理的是图像（像素通常为8位字节），使用 char 类型而非 float，计算内存比会有所不同（访问量更小，但操作也可能是整数运算）。优化需根据具体数据类型和应用需求进行。

总结 📝

本节课我们一起学习了：

卷积的定义及其作为并行模式的应用。
并行化卷积的一种简单方法：每个输出元素一个线程。
常量内存的特性、优势（只读、硬件管理缓存、快速访问）及其在CUDA中的使用方法（__constant__ 声明，cudaMemcpyToSymbol 复制）。
实现了一个基本的卷积内核，并讨论了边界处理。
分析了卷积中的数据重用模式，并引入了共享内存分块作为优化方向，以显著提升计算与内存访问的比率。
通过公式对比了分块与非分块实现的性能潜力。

在接下来的作业中，你将有机会实践分块卷积的实现。更多细节可以参考教材第7章。

下节预告：在下一讲中，我们将继续探索其他并行模式及相关的GPU优化技术。

GPU计算：P9：模板计算 🧮

在本节课中，我们将要学习一种新的并行模式——模板计算。我们将探讨其基本概念、与卷积运算的异同，并学习如何利用CUDA的3D线程网格、共享内存和寄存器平铺等技术来高效地实现和优化3D模板计算。

概述

上一节我们介绍了卷积运算，它引入了常量内存的使用，并且每个输出元素都是对应输入元素及其周围元素的加权和。本节中我们来看看模板计算，它是卷积运算的一个特例，通常用于结构化网格上的计算，其中每个网格点的值基于其邻居子集的值来计算。

什么是模板计算？

模板计算模式指的是一类在结构化网格上进行的计算，其中网格点的值基于该点邻居的一个子集来计算。

例如，在2D网格中，一个五点模板意味着某个输出网格点的值是其自身输入值以及上下左右四个邻居输入值的组合。在3D中，一个七点模板则包含在X、Y、Z三个维度上的邻居。

为了便于说明，我们将在2D和3D示例之间切换。通常，这些网格点被方便地存储为多维数组。

与卷积的异同

模板计算与卷积非常相似，实际上是卷积的一个特例。主要区别在于：

卷积通常访问输入元素周围的整个块（例如3x3的卷积核）。
模板通常只访问特定方向的直接邻居（例如上下左右），这为某些优化提供了可能，而这些优化在通用卷积中可能无法实现。

并行化策略

最细粒度和直观的并行化方式是为每个输出网格点分配一个线程，这与卷积类似。

对于2D模板，我们使用2D线程网格。对于3D模板，我们自然需要使用3D线程网格。

边界条件简化

为了简化边界条件的处理，我们做一个贯穿课程的假设：只计算内部的输出值，不计算边界上的值。这确保了我们在访问邻居输入值时，这些值始终在有效范围内。这个假设通常是合理的，因为在许多实际应用中，边界元素可能由其他GPU或进程负责计算。

基础3D模板实现

现在，让我们使用3D线程网格来实现一个基础的3D模板计算。

内核配置

我们配置一个3D的线程块（例如 8x8x8，共512个线程），并计算一个3D的网格块数量来覆盖整个输出空间。输出和输入都是 N x N x N 的立方体。

dim3 blockDim(BLOCK_SIZE, BLOCK_SIZE, BLOCK_SIZE);
dim3 gridDim((N + BLOCK_SIZE - 1) / BLOCK_SIZE,
             (N + BLOCK_SIZE - 1) / BLOCK_SIZE,
             (N + BLOCK_SIZE - 1) / BLOCK_SIZE);
stencil_kernel<<<gridDim, blockDim>>>(d_in, d_out, N);

内核函数步骤

以下是内核函数的主要步骤：

计算线程对应的输出索引：每个线程根据其3D的块索引和线程索引，计算出它负责的3D输出数组中的位置 (i, j, k)。
```
int i = blockIdx.z * blockDim.z + threadIdx.z;
int j = blockIdx.y * blockDim.y + threadIdx.y;
int k = blockIdx.x * blockDim.x + threadIdx.x;
```
边界检查：确保只计算内部的输出元素（忽略 i, j, k 为 0 或 N-1 的边界）。
```
if (i >= 1 && i < N-1 && j >= 1 && j < N-1 && k >= 1 && k < N-1) {
    // 进行计算
}
```

执行模板计算：对于内部的每个输出元素，计算其值为对应输入元素及其六个邻居（X、Y、Z方向各两个）的线性组合。这里使用常量 C0 和 C1 作为权重。

int idx = i * N * N + j * N + k; // 3D数组转1D线性索引
out[idx] = C0 * in[idx]
           + C1 * (in[idx - 1] + in[idx + 1]     // X方向邻居
                 + in[idx - N] + in[idx + N]     // Y方向邻居
                 + in[idx - N*N] + in[idx + N*N]); // Z方向邻居

这个基础版本直接访问全局内存，性能有提升空间。接下来，我们探讨如何优化。

优化：共享内存平铺

与卷积类似，处理相邻输出元素的线程块会加载许多相同的输入元素，存在数据复用。因此，第一个优化思路是使用共享内存平铺。

平铺的挑战

输入瓦片和输出瓦片的尺寸不同。对于3D七点模板，输入瓦片在每个维度上都比输出瓦片大2（多一层“边界”）。因此，我们需要启动足够多的线程来加载整个输入瓦片，但在计算时，只使用其中一部分线程来计算内部的输出瓦片。

实现思路

定义瓦片尺寸：BLOCK_SIZE 对应输入瓦片尺寸。输出瓦片尺寸则为 BLOCK_SIZE - 2。
调整网格配置：计算网格大小时，使用输出瓦片尺寸，以确保有足够的块覆盖所有输出元素。
加载输入瓦片到共享内存：
- 每个线程负责将全局内存中的一个输入元素加载到共享内存的3D数组中。
- 加载时需要检查该输入索引是否在全局数组边界内。
同步线程：确保所有线程完成数据加载。
计算输出：
- 只有位于线程块内部（非边界）的线程才参与计算。
- 计算时，从共享内存中读取数据，而不是全局内存。

然而，初步测试可能发现性能提升并不明显，甚至可能更慢。这是因为对于较小的瓦片尺寸（如8x8x8），边界元素的比例较高，数据复用的收益被共享内存同步和索引计算的开销所抵消。

计算与内存访问比分析：

基础版本：每个线程进行8次操作（2次乘法，6次加法），加载7个浮点数（28字节）。计算内存访问比约为 8 ops / 28 bytes ≈ 0.29 ops/byte。
平铺版本（瓦片尺寸T）：每个线程块进行 8 * (T-2)^3 次操作，加载 T^3 * 4 字节数据。比值约为 (T-2)^3 / T^3。当 T=8 时，比值约为 0.84 ops/byte，虽有提升但不显著。

结论：增大瓦片尺寸 T 可以提高内部元素的比例，从而提升计算内存访问比。但直接增大 T 会受到两个硬件限制：

线程块最大线程数（通常为1024）。
共享内存容量有限，大瓦片会降低占用率（Occupancy）。

进一步优化：线程协作与寄存器平铺

为了处理更大的输出瓦片而不增加线程数或共享内存使用，我们引入两种技术：线程协作和寄存器平铺。

线程协作

线程协作允许一个线程处理输出瓦片中的多个元素。这样，我们可以用固定数量的线程（例如32x32）处理一个更大的3D输出瓦片（例如32x32x32），只需让每个线程在Z维度上循环处理多个平面。这减少了对边界元素的冗余加载。

寄存器平铺

观察3D七点模板的计算特点：当处理一个输出平面时，每个线程需要：

当前平面的数据：被该平面内多个相邻线程共享。
前一个平面和后一个平面的数据：每个线程只需要自己对应的那个元素，不被其他线程共享。

因此，我们可以优化内存层次的使用：

将当前平面的数据块放入共享内存，供线程间共享访问。
将前一个和后一个平面中每个线程需要的那一个数据元素，保存在该线程的寄存器中。
随着计算在Z维度上推进，寄存器中的数据在迭代间移动和更新（例如，后一个平面变成当前平面，当前平面变成前一个平面），仅在需要被共享时才写入共享内存。

这种方法被称为寄存器平铺，它显著减少了对共享内存的需求，从而允许使用更大的瓦片尺寸或提高占用率。

寄存器平铺的启示：我们之前在矩阵乘法中已经无形中使用了寄存器平铺——每个线程将累加结果（C矩阵的一个元素）存储在自己的寄存器中。模板计算中的显式应用加深了我们对这一优化技术的理解。

总结

本节课中我们一起学习了模板计算这一重要的并行模式。

核心概念：模板计算是卷积的特例，用于基于邻居更新网格点。我们以实现3D七点模板为例。
基础实现：使用3D线程网格，每个线程处理一个输出点，并直接访问全局内存。
共享内存平铺：利用数据复用，将输入瓦片加载到共享内存。但受限于瓦片尺寸和硬件限制，收益可能有限。
高级优化：
- 线程协作：让一个线程处理多个元素，以增大有效瓦片尺寸，减少边界开销。
- 寄存器平铺：根据数据访问模式，将仅限线程私有的数据保存在寄存器中，将需要线程间共享的数据保存在共享内存中，优化内存层次使用，提升性能。

通过结合共享内存平铺、线程协作和寄存器平铺，我们可以构建出高性能的模板计算内核，充分挖掘GPU的并行计算潜力。

GPU计算：P10：归约

在本节课中，我们将要学习一种新的并行模式——归约。我们将了解什么是归约，如何从简单的串行实现开始，逐步设计出高效的并行GPU算法，并探讨如何通过优化内存访问、减少线程发散和使用线程粗化等技术来提升性能。

回顾：模板计算

上一节我们介绍了模板计算。模板计算是一种在网格上进行的计算，网格中某点的输出值依赖于该点及其邻居点的输入值。我们以一个二维模板为例，通常将网格点存储为二维数组。对于每个输出点，其结果是基于对应输入点及其相邻网格点的输入值计算得出的。

我们还研究并实现了三维模板计算，这是我们首次使用三维线程网格的经验。三维模板与二维类似，只是除了X和Y维度的邻居外，还使用了Z维度的邻居。并行化模板的一种直接方法是为每个输出元素分配一个线程，这可以通过使用三维线程网格来实现。

我们实现了这样一个版本的模板计算，但随后观察到模板计算存在显著的数据重用。如果一个线程块负责一片输出元素，那么这个线程块会共同使用一片输入元素。其中一些输入元素会被多个输出元素使用。由于这种数据重用，我们通常通过共享内存来利用它。我们学习了如何将输入数据块加载到共享内存中，然后用于计算输出元素。

接下来，我们分析了计算与内存访问的比率，发现通过分块，我们能够改善这个比率。并且，分块尺寸越大，计算内存比就越好。因此，我们希望从小分块尺寸过渡到大分块尺寸。直观上看，更大的分块尺寸意味着相对于边界元素，内部元素的数量更多。在模板计算中，边界元素不像卷积那样被重用。因此，更大的分块尺寸允许我们有更少的边界元素和更多的内部元素，从而带来更多的数据重用。

然而，增加分块尺寸（尤其是三维块）的问题在于，我们很快就会超过加载输入块和计算输出块所需的最大每块线程数，以及存储整个输入块可用的共享内存量。即使没有超过共享内存量，使用过多的共享内存也会损害我们的占用率。我们处理这个问题的方法是使用线程粗化，以便在不增加线程数量的情况下处理更大的输出块。我们没有启动三维块来处理整个立方体，而是启动一个二维线程块。这个块将顺序遍历立方体，在每一步中，它只需要输入块中的三个平面。这样做通过粗化减少了线程数量，同时也减少了所需的共享内存，因为我们只需要在任何时间点存储三个平面。

我们观察到的最后一点是，在任何时间点，当我们使用这三个平面（前一个平面、当前平面和下一个平面）时，只有当前平面中的元素实际上被计算不同输出的线程共享。但仔细观察，下一个平面和上一个平面中的元素只由加载它们的线程使用。基于这个观察，我们提出，与其将所有三个平面都存储在共享内存中，不如将下一个平面加载到寄存器中。当它变成当前平面时，我们将寄存器移动到共享内存；当当前平面变成前一个平面时，我们将共享内存移回寄存器。我们称这种将数据块存储在寄存器中而非共享内存的能力为寄存器分块。我们之前实际上已经见过寄存器分块，在矩阵乘法和卷积中，我们一直对输出块应用寄存器分块，只是不那么明显。在模板计算中，寄存器分块的概念变得更加明显，因为同一个数据块有时在共享内存中，有时在寄存器中。

以上就是对上节课内容的快速回顾。

引入归约模式

本节中，我们来看看一种新的并行模式——归约。

归约是一种将一组输入值减少为一个输出值的操作。我们有一组输入值，并将它们归约为一个输出值。归约操作可以是求和、求积、求最小值或求最大值。通常，归约操作是一个满足结合律和交换律的操作，并且具有明确定义的单位元。求和、求积、求最小值、求最大值都满足结合律和交换律，并且都有明确定义的单位元。单位元本质上是累加器的初始值，或者是在没有元素可归约时得到的值。对于求和，单位元是0；对于求积，单位元是1；对于求最小值，单位元是正无穷大；对于求最大值，单位元是负无穷大。

在本讲座中，我们将以求和作为归约的例子。然而，我们将讨论的所有内容同样适用于其他形式的归约，如求积、最小值和最大值。有时人们也称归约为折叠，折叠比归约更通用，但如果你熟悉折叠，归约基本上就是一种折叠。

串行归约实现

首先，让我们看看如何实现串行归约。

如果要实现串行求和归约，可以简单地使用一个循环。将和初始化为0，然后从i = 0循环到n（n是输入元素的数量），执行 sum += input[i]。通常，归约的模式是：有一个累加器，将其初始化为单位元，然后遍历数组元素，对累加器和输入应用某个满足结合律和交换律的操作F，以获得累加器的新值。

如你所见，归约具有相当的顺序性。我们有一个循环，并且存在循环携带依赖。这与向量加法不同，向量加法的所有循环迭代都是独立的。这里，循环迭代之间存在依赖关系。那么，如何并行化这样的操作呢？

并行归约：树形结构

我们可以将其分割，让多个线程计算多个部分，然后将它们相加。例如，对于求和，每个线程可以求和一部分，但最终每个线程会得到一个部分和，我们需要将它们组合在一起。假设每个线程处理一个元素，如何组合它们呢？答案是利用树形结构。

并行归约可以按如下方式进行：如果有一个包含八个元素的数组，可以创建四个线程，每个线程并行地相加两个元素。这样就将需要求和的元素数量减少了一半。现在有四个元素，可以启动两个线程来并行相加这四个元素，得到两个部分和。最后，使用一个线程来并行相加这两个元素，得到最终结果。

我们称之为归约树。它的计算模式类似于树形结构，但我们并没有在内存中实际构建树数据结构。你会注意到，每一步每个线程相加两个元素。如果有n个元素，完全并行地归约需要多少步？答案是log₂(n)。因为从n个元素开始，然后有n/2个部分和，接着是n/4个部分和，依此类推，直到完成。每一步，线程数量减半。如果有n个元素，第一步有n/2个线程，下一步有n/4个线程，再下一步有n/8个线程。每一步都有一半的线程退出。这种模式称为归约树。

需要注意的一点是，当以这种方式归约元素时，执行后续归约的线程必须确保前一步的加法已经完成。这意味着需要能够在线程之间进行同步。之前我们使用同步是为了让从共享内存加载数据的线程相互等待，确保所有人都加载完毕后再开始使用数据。在这里，我们看到同步的不同用途：一个线程依赖于两个值，它需要等待其他线程确保这些值可用。这与共享内存的上下文略有不同，但概念相似。

通常，我们让所有这些线程并行执行归约，然后在它们之间进行同步。接着，让负责这两个元素的线程执行归约，再进行同步。每一步都需要一次同步。

GPU上的同步与分段归约

在GPU上，哪些线程可以相互同步？同一个线程块内的线程可以同步，但不同线程块之间的线程不能同步。因此，如果在一个线程块内进行归约，我们可以做到这一点。但跨多个线程块进行归约树则有些棘手。我们通常在GPU上执行一种称为分段归约的操作。

由于线程必须在步骤之间同步，而我们无法跨块同步线程，因此我们将进行分段归约。每个线程块将其负责的输入段归约，产生一个部分和。然后，我们稍后再将这些部分和归约在一起。换句话说，将输入分成若干段，每个线程块在其段上执行一个需要同步的归约树。当线程块完成后，每个线程块将得到一个部分和，并将这些部分和存储在一个部分和数组中。最后，对这个部分和数组进行高层次归约。

有几种方法可以归约部分和。通常，部分和数组比原始输入数组小得多。一种方法是启动一个新的归约内核，对部分和数组应用同样的归约树，重复此过程直到部分和数组只有一个线程块大小，此时得到的总部分和就是最终总和。另一种方法是使用原子操作将这些部分和原子地累加到一个累加器中。第三种方法是，如果部分和数量足够少，可能不值得并行化部分和的加法，可以直接在CPU上相加。

今天我们将重点关注归约树本身，即如何在单个线程块内执行归约树。其余部分（如部分和的最终归约）将在主机端处理。

线程块内的归约树实现

现在，让我们看看如何在单个线程块内实现归约树。假设这是线程块负责的段，一种方法是为每隔一个输入元素分配一个线程。每个线程将其值与右侧相邻的值相加，这是第一步，现在值的数量减少了一半。

在下一步，需要将这个元素和这个元素相加，这个元素和这个元素相加，依此类推。只需要一个线程来相加这两个元素，而之前负责这个元素的线程不再需要，因此它将退出。使用这个线程来相加这两个元素，而那个线程将不执行任何操作。以此类推，在下一步，只需要这个线程来相加这两个元素，其余三个线程不执行任何操作。继续这个过程，最后一步，还有两个值需要相加，只有一个线程执行加法，其余线程不执行任何操作。这是在GPU上实现归约树的一种方式。

显然，当这些线程执行加法时，负责后续相加的线程必须等待前一步的线程完成计算，以提供所需的值。因此，我们将让线程执行两个元素的加法，然后进行线程同步，接着再执行两个元素的加法，再进行同步，依此类推。

现在，让我们在代码中实现这个逻辑。

基础并行归约代码实现

以下是我们将实现的基础归约内核代码的步骤和逻辑。

首先，每个线程块需要确定它负责的段。注意，一个线程块负责的元素数量实际上是线程数量的两倍，因为在第一次迭代中，每个线程将相加两个元素。我们不希望为线程块负责的每个元素都启动一个线程，因为那样从一开始就会有一半的线程不执行任何操作。

为了找到线程块负责的段，我们需要为每个线程块跳过两倍于线程数量的元素。段起始索引 segment_start 是 blockIdx.x * (blockDim.x * 2)。

接下来，每个线程需要确定它在该段内负责哪个元素。线程0负责元素0，线程1负责元素2，线程2负责元素4，依此类推。通常，线程索引为 threadIdx.x 的线程负责段内索引为 2 * threadIdx.x 的元素。该元素在全局数组中的绝对索引 i 是 segment_start + 2 * threadIdx.x。

现在，每个线程知道了它在输入中负责的元素，我们需要一个循环来迭代这些步骤以执行归约树。一个方便的循环变量是步长。在第一步，每个线程相加右侧距离为1的元素；第二步，距离为2；第三步，距离为4；第四步，距离为8。步长从1开始，每次迭代翻倍，直到达到线程块维度（即线程数量）。因此，循环可以写为 for (int stride = 1; stride <= blockDim.x; stride *= 2)。

在每次迭代中，每个线程将把它拥有的元素与右侧距离为步长的元素相加：input[i] += input[i + stride]。但是，并非每个线程在每次迭代都执行这个加法。当步长为1时，所有线程都执行；步长为2时，每隔一个线程执行（线程0, 2, 4, 6）；步长为4时，每四个线程执行（线程0, 4）；步长为8时，只有线程0执行。一般来说，只有线程索引是步长倍数的线程才执行计算。因此，需要一个保护条件：if (threadIdx.x % stride == 0)。

在执行加法后，所有线程必须同步，然后才能进行下一步。使用 __syncthreads()。

完成所有步骤后，每个线程块将在块内线程0的位置得到一个部分和。然后，只有线程0将其结果存储到部分和数组中：if (threadIdx.x == 0) partial_sums[blockIdx.x] = input[i]。

这就是我们的基础归约代码。

性能分析与初始优化

运行此代码，在CPU上耗时21毫秒，在GPU上内核执行耗时3.8毫秒。但如果包括所有数据拷贝时间，GPU和CPU的总时间几乎相同。归约不是像矩阵乘法或卷积那样高度并行的操作，由于其顺序性，我们不期望获得那么大的性能提升。然而，如果数据已经在GPU上，那么在GPU上进行归约仍然是值得的。

现在，让我们讨论一下这个代码有哪些不优化的地方，以及如何进一步优化。

首先，存在大量的内存访问，我们反复访问全局内存。
其次，线程数量不断减少，线程在退出。虽然理论上可以用其他工作填充这些空闲线程周期，但问题在于这些退出的线程仍然会占用执行资源，因为GPU以线程束为单位执行。如果同一个线程束中的一些线程执行计算而另一些不执行，就会导致控制流发散，造成资源浪费。
此外，内存访问模式不连续。我们曾讨论过内存合并访问，希望同一个线程束中的线程访问全局内存中相邻的数据以实现合并访问。但在这里，线程访问的是间隔一个元素的数据，然后访问右侧距离步长的数据。访问模式不连续，并且随着步长增加而恶化。

我们可以通过重新排列线程访问数据的方式来解决合并访问和控制流发散问题。

优化：改善合并访问与减少线程发散

我们不再为数组中每隔一个的元素分配线程，而是将块内的八个线程分配给数组的前八个初始元素。在第一次迭代中，它们都将这八个元素与右侧距离为块大小的八个元素相加，并将结果存储在这里。现在，当这些线程访问左侧元素时，访问是连续的，因为它们是相邻的；访问右侧元素时也是连续的。在下一步，上半部分线程退出，下半部分线程保持活动状态，并以这种方式执行加法。在后续步骤中，依此类推。

这种方法同时解决了合并访问和控制流发散问题。原因是，现在线程同时访问的数据是相邻的，并且只有末尾的线程退出，而不是中间的线程退出。如果前四个线程构成一个线程束（实际是32个，这里为说明而缩小），那么在第一次迭代后，这四个线程将一起退出，而另外四个线程将进行计算。因此，要么整个线程束都在计算，要么整个线程束都已退出，没有控制流发散。

让我们来实现这种方法。

在代码中，每个线程现在负责的元素是 segment_start + threadIdx.x，而不是 segment_start + 2 * threadIdx.x。循环步长现在从 blockDim.x 开始，然后变为4、2、1。循环条件为 for (int stride = blockDim.x; stride >= 1; stride /= 2)。在每次迭代中，只有线程索引小于步长的线程才执行计算：if (threadIdx.x < stride)。然后执行 input[i] += input[i + stride]，接着同步，最后线程0存储结果。

我们重构了归约树以消除控制流发散并获得更好的合并访问。运行优化后的代码，性能从3.79毫秒提升到2.75毫秒，显示出显著的改进。

进一步优化：使用共享内存

接下来，我们可以使用共享内存进行优化。这里并没有真正的数据重用，因为没有任何一个特定值被多次使用。然而，包含这些值的内存位置实际上被重用了。优化方法是先将输入加载到共享内存中，然后在共享内存上执行归约树。这样可以避免修改原始输入数组。

具体来说，在第一步，线程加载它们的两个初始元素，相加后将结果存储在共享内存数组中。我们只需要与线程数量一样多的共享内存（即块大小）。初始加载来自全局内存，随后的步骤都在共享内存中进行。当完成后，最终的线程将共享内存中的值存入全局内存。这减少了访问全局内存的次数。

寄存器比共享内存更快，但我们不能使用寄存器，因为一个线程写入的值需要被另一个线程读取。寄存器是线程私有的，如果希望同一个块内的多个线程读写同一个值，必须将其放在共享内存中。

让我们实现这个优化。

首先，在共享内存中分配一个数组：__shared__ float input_s[blockDim.x];。在第一步，每个线程加载它负责的元素和步长距离的元素，相加后存入共享内存：input_s[threadIdx.x] = input[i] + input[i + stride];。然后进行同步 __syncthreads()。

现在，在共享内存上继续执行归约树。循环步长从 blockDim.x / 2 开始到1：for (int stride = blockDim.x/2; stride >= 1; stride /= 2)。在循环内，执行 if (threadIdx.x < stride) input_s[threadIdx.x] += input_s[threadIdx.x + stride]; 然后同步。

最后，线程0将共享内存中的最终结果存入全局部分和数组：if (threadIdx.x == 0) partial_sums[blockIdx.x] = input_s[0];。

运行此代码，性能从2.75毫秒提升到2.15毫秒。

高级优化：线程粗化

最后，我们将应用线程粗化。线程粗化的思想是，如果硬件要序列化线程块，我们不如自己在线程块内序列化工作，而不是让硬件来序列化。这样，每个线程可以负责更多的输入元素，从而减少同步和控制流发散的开销。

在归约中，并行化的代价是同步和控制流发散。通过线程粗化，我们可以让每个线程块负责更多的数据，从而减少需要同步和存在控制流发散的步骤数量。

具体来说，我们定义一个粗化因子（例如4）。这意味着每个线程块负责的输入段大小是原来的粗化因子倍。每个线程现在负责 2 * coarsening_factor 个初始元素。线程将循环遍历这些元素，将它们累加到一个寄存器变量中，然后将累加结果存入共享内存。这个初始累加阶段没有控制流发散，也不需要同步，因为每个线程处理的是自己的数据段。只有在将结果存入共享内存后，才执行需要同步和控制流发散的归约树步骤。

通过粗化，我们将初始的多个无发散、无同步的步骤，替换了原来可能被硬件序列化的多个块的归约树步骤，从而减少了总体同步和发散开销。

在代码中，我们定义粗化因子 COARSENING_FACTOR。每个线程块负责的段大小变为 blockDim.x * 2 * COARSENING_FACTOR。每个线程通过一个循环累加它负责的多个元素：for (int tile = 0; tile < COARSENING_FACTOR * 2; ++tile) sum += input[i + tile * blockDim.x];。然后将 sum 存入共享内存 input_s[threadIdx.x] = sum;，接着执行共享内存上的归约树。

由于每个线程块处理更多数据，我们需要启动更少的线程块。在主机代码中，计算块数量时应将每块元素数乘以粗化因子。

运行使用粗化因子的代码，性能得到进一步提升。可以尝试不同的粗化因子，性能会先提高，达到某个点后，由于过度序列化并干扰了硬件的透明可扩展性，性能开始下降。

边界条件处理

最后，在处理边界条件时，最后一个线程块可能包含超出数组范围的元素。在从全局内存加载数据的初始阶段，需要确保加载的索引在边界内，否则可以加载0。在共享内存中执行归约树时，由于共享内存大小固定且由活动线程填充，通常不需要额外处理边界。

总结

本节课中，我们一起学习了归约这一并行模式。我们从串行归约出发，探讨了如何利用树形结构实现并行归约。重点学习了在GPU上实现高效归约的关键技术：通过重新组织线程访问模式来改善内存合并访问并减少控制流发散；利用共享内存减少全局内存访问次数；以及应用线程粗化来减少同步开销和线程发散，从而提升整体性能。我们还简要讨论了边界条件的处理方法。归约是许多并行算法的基础组件，掌握其优化技巧对GPU编程至关重要。

GPU计算：第11讲：扫描（Kogge-Stone算法）🚀

概述

在本节课中，我们将学习一种新的并行模式——扫描。我们将重点介绍一种实现并行扫描的具体方法：Kogge-Stone算法。我们将从扫描的基本概念开始，逐步深入到其并行实现、优化技巧，并分析其性能特点。

回顾：归约操作

上一节我们介绍了归约操作。归约是一种将一组输入值通过某种运算符（如求和、求积、求最小值、求最大值）合并为单个值的操作。该运算符需要满足结合律和交换律，并具有一个单位元。

一个顺序求和归约的伪代码如下：

sum = 0; // 单位元
for (int i = 0; i < n; i++) {
    sum = sum + input[i]; // 运算符
}

我们学习了如何并行化归约，核心思想是使用归约树。通过将线程组织成树状结构，可以在 log n 步内完成 n 个元素的归约。然而，这需要线程块内的全局同步。由于无法在GPU上跨线程块进行全局同步，我们通常采用分段归约策略：每个线程块进行本地归约，产生部分和，然后通过启动新内核、原子操作或在CPU上累加这些部分和来得到最终结果。

我们还探讨了优化方法，例如通过线程粗化来减少同步开销和控制流分叉，从而提升性能。

扫描操作简介

本节中，我们来看看扫描操作。扫描操作接收一个输入数组 X[0...n-1] 和一个满足结合律的运算符（如加法），并返回一个输出数组 Y[0...n-1]。

扫描有两种形式：

包含性扫描：输出 Y[i] 是输入中从 X[0] 到 X[i]（包含 X[i]）所有元素的组合结果。
- 公式：Y[i] = X[0] op X[1] op ... op X[i]
排他性扫描：输出 Y[i] 是输入中从 X[0] 到 X[i-1]（不包含 X[i]）所有元素的组合结果。Y[0] 是单位元。
- 公式：Y[0] = identity; Y[i] = X[0] op X[1] op ... op X[i-1] (对于 i > 0)

示例（加法运算符）：

输入数组: [3, 6, 7, 4, 2, 1, 5, 9]
包含性扫描结果: [3, 9, 16, 20, 22, 23, 28, 37]
排他性扫描结果: [0, 3, 9, 16, 20, 22, 23, 28]

顺序扫描的伪代码如下：

// 包含性扫描
output[0] = input[0];
for (int i = 1; i < n; i++) {
    output[i] = output[i-1] op input[i];
}

// 排他性扫描
output[0] = identity;
for (int i = 1; i < n; i++) {
    output[i] = output[i-1] op input[i-1];
}

分段扫描策略

与归约类似，并行扫描也需要线程间的同步。由于无法跨线程块同步，我们采用分段扫描（或称为层次化扫描）策略。

以下是分段扫描的步骤：

本地扫描：每个线程块负责输入数组的一个分段，并在其内部执行并行扫描。
收集部分和：每个线程块将其分段内所有元素的组合结果（即该分段的最后一个扫描结果）存储到一个部分和数组中。
扫描部分和：对这个部分和数组本身执行一次扫描（可以在GPU上递归调用扫描内核，或在CPU上执行）。
更新结果：每个线程块将上一步得到的、对应其之前所有分段的扫描结果，加到其本地扫描结果的每个元素上，从而得到最终的全局扫描结果。

这种方法允许各个线程块独立、并行地处理数据分段，最后通过组合部分结果来获得完整答案。

Kogge-Stone并行扫描算法

现在，我们聚焦于如何在单个线程块内高效地实现并行扫描。我们将使用Kogge-Stone算法来实现包含性扫描。

该算法的核心思想是通过 log n 步迭代，逐步构建扫描结果。假设有8个元素（x0 到 x7），每个元素由一个线程负责。

算法步骤：

第1步（跨度=1）：每个线程（除了第一个）将其元素与左边相距1个位置的元素相加。结果：x1 包含了 x0+x1，x2 包含了 x1+x2，依此类推。
第2步（跨度=2）：每个线程（索引 >=2）将其当前值（已是部分和）与左边相距2个位置的元素（部分和）相加。结果：x2 包含了 x0+x1+x2，x3 包含了 x0+x1+x2+x3。
第3步（跨度=4）：每个线程（索引 >=4）将其当前值与左边相距4个位置的元素相加。经过 log2(8)=3 步后，所有线程都得到了正确的包含性扫描结果。

关键点：

每一步的跨度翻倍：1, 2, 4, ...
只有线程索引大于等于当前跨度的线程才执行加法操作。
总共需要 log2(n) 步。

基础GPU内核实现

以下是基于Kogge-Stone算法的包含性扫描基础GPU内核代码框架。我们首先在全局内存中操作。

__global__ void scan_kernel(float* input, float* output, float* partial_sums, int n) {
    // 1. 计算全局索引
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= n) return;

    // 2. 将输入数据拷贝到输出数组（原地操作准备）
    output[i] = input[i];
    __syncthreads(); // 确保所有数据已就位

    // 3. Kogge-Stone扫描循环
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
        // 读取所需的数据
        float val = 0;
        if (threadIdx.x >= stride) {
            val = output[i - stride];
        }
        __syncthreads(); // 所有线程完成读取

        // 执行加法并写回
        if (threadIdx.x >= stride) {
            output[i] += val;
        }
        __syncthreads(); // 所有线程完成写入，为下一步做准备
    }

    // 4. 存储本线程块的部分和（由最后一个线程负责）
    if (threadIdx.x == blockDim.x - 1) {
        partial_sums[blockIdx.x] = output[i];
    }
}

注意：上述代码存在竞态条件。在读取 output[i - stride] 和写入 output[i] 之间，其他线程可能已经修改了 output[i - stride] 的值。因此，我们需要将读取和写入分离，并用同步点隔开，如代码中所示。

优化1：使用共享内存

在全局内存中进行多次读写操作效率低下。我们可以利用共享内存（更快）来存储中间数据。

优化步骤：

将数据从全局内存加载到共享内存缓冲区。
在共享内存中执行所有Kogge-Stone扫描步骤。
将最终结果从共享内存写回全局内存输出数组。

这种优化显著减少了访问延迟，提升了内核性能。

优化2：双缓冲

观察发现，每次迭代需要两次 __syncthreads()：一次确保所有读取完成，一次确保所有写入完成。双缓冲技术可以消除其中一次同步。

以下是双缓冲的原理：

我们分配两个共享内存缓冲区：buffer1 和 buffer2。
在每次迭代中，一个作为输入缓冲区，另一个作为输出缓冲区。
每个线程：
- 从输入缓冲区读取数据。
- 如果线程索引 >= 跨度，则进行计算并将结果写入输出缓冲区。
- 否则，直接将输入缓冲区的值复制到输出缓冲区（保持数据完整）。
迭代结束后，交换两个缓冲区的角色（输入变输出，输出变输入）。

因为读写操作发生在不同的物理缓冲区，所以不存在对同一内存位置的读写冲突，从而移除了读写操作之间的同步屏障。每次迭代仅需一次同步来确保所有线程都完成了对输出缓冲区的写入，然后才能交换缓冲区进行下一步。

双缓冲在减少同步开销方面带来了显著的性能提升。

实现排他性扫描

基于已实现的包含性扫描，实现排他性扫描非常简单。核心思想是数据偏移。

方法：

在将数据加载到共享内存时，每个线程加载其前一个输入元素（input[i-1]）。
对于线程块中的第一个线程（threadIdx.x == 0），则加载单位元（例如加法中的0）。
对这样偏移后的共享内存数组执行包含性扫描。
得到的结果就是原始输入数组的排他性扫描结果。

此外，需要注意部分和的计算和后续添加步骤需要相应调整，以确保最后一个元素被正确计入部分和，并且分段组合逻辑正确。

工作效率分析

一个并行算法如果其执行的总工作量（如浮点加法次数）与对应的顺序算法相同，则称其为工作高效的。

顺序扫描：对于 n 个元素，需要 n-1 次加法操作。工作量是 O(n)。
Kogge-Stone并行扫描：分析其操作次数。每一步中，执行加法的线程数递减。总加法次数约为 n * log2(n) - (n - 1)。工作量是 O(n log n)。

因此，Kogge-Stone算法不是工作高效的。它在并行资源充足时可以通过并行性快速完成，但如果硬件资源有限（例如并行度不足），它可能因为做了更多额外工作而比顺序算法更慢。

总结

本节课中我们一起学习了：

扫描操作的定义、包含性与排他性的区别及其顺序实现。
在GPU上实现扫描的分段策略，以解决跨线程块同步的限制。
Kogge-Stone并行扫描算法的原理，它通过 log n 步迭代在线程块内完成扫描。
基础的GPU内核实现，并解决了其中的竞态条件问题。
两个关键优化：使用共享内存减少全局内存访问延迟，以及使用双缓冲技术减少同步开销。
如何通过数据偏移，基于包含性扫描轻松实现排他性扫描。
分析了Kogge-Stone算法的工作效率，指出其 O(n log n) 的工作量并非最优。

下次课我们将探讨另一种工作更高效的并行扫描算法，并对比它们的性能。

GPU计算：第12讲：扫描（Brent-Kung算法） 🚀

在本节课中，我们将学习一种不同于上次讨论的并行扫描算法——Brent-Kung并行扫描算法。我们还将探讨线程粗化在扫描中的应用，以提升算法的工作效率。

概述 📋

上一讲我们介绍了扫描模式，并重点讲解了用于并行扫描的Kogge-Stone方法及其优化技术——双缓冲。扫描模式是指，给定一个数组，输出数组的每个元素将是其前面所有元素的组合（对于包含性扫描，还包括当前元素；对于排他性扫描，则不包括）。我们以加法为例，但该模式同样适用于其他满足结合律的操作符，如乘法、最小值、最大值等。

为了实现并行扫描，我们像处理归约一样对输入进行分段。这是因为并行扫描需要同步，而不同线程块之间无法实现同步。因此，我们将输入分段，每个线程块扫描自己的段，将部分和存储在一个子数组中，然后扫描这个部分和数组。最后，每个线程块将扫描后的部分和（来自其前面的块）加到其内部元素上，从而得到最终的全局扫描数组。

我们用于单个线程块内扫描的算法是Kogge-Stone并行算法。该算法让每个线程加上与其相距特定步长的元素，步长从1开始，每次翻倍，直到覆盖整个数组。由于存在读写冲突，我们需要在读和写之间进行同步。为了减少同步开销，我们引入了双缓冲技术，即使用输入缓冲区和输出缓冲区，并在每次迭代后交换它们。

我们还分析了Kogge-Stone算法的工作效率。虽然顺序扫描执行n次加法，但Kogge-Stone方法在log n步中，每步执行n - 2^步长次操作，总计约O(n log n)次操作。这表明该算法的工作效率不高，因为它比顺序算法执行了更多冗余操作。如果线程块因资源有限而被串行化，这种低工作效率可能会影响性能。

因此，本节课我们将探讨Brent-Kung扫描算法，并学习如何应用线程粗化来优化扫描。

Brent-Kung并行扫描算法 🔄

Brent-Kung并行包含性扫描算法分为两个主要步骤。

第一步：归约阶段

第一步类似于一个归约树。我们首先将相邻的两个元素相加。在这个阶段，我们为每两个元素分配一个线程，因为每一步的加法操作数量在减少。接着，我们将相距两个步长的元素相加，然后是相距四个步长的元素，依此类推，直到完成整个数组的归约树。这个阶段被称为归约阶段。

第二步：后归约阶段

完成归约阶段后，数组中部分扫描值已经就绪，但其他值尚未完成。例如，x4和x5的和缺少x0到x3的部分。在第二阶段，我们将完成这些剩余值的计算。我们通过将左侧相距一个步长的值相加来完成扫描。例如，将x0到x3加到x4和x5上，得到x0到x5的和。类似地，我们将x0和x1加到x2上，得到x0到x2的和，以此类推。这个阶段被称为后归约阶段。

算法对比与分析

以下是Kogge-Stone算法与Brent-Kung算法的对比：

Kogge-Stone算法：步数少（log n步），但总操作数多（O(n log n)）。
Brent-Kung算法：步数多（2 log n - 1步），但总操作数少（O(n)次操作），因此具有更高的工作效率。

具体分析如下：

归约阶段：执行log n步，操作数为n - 1。
后归约阶段：执行log n - 1步，操作数为n - 2 log n - 1。
总计：步数为2 log n - 1，操作数为2n - log n - 2，即O(n)次操作。

因此，Brent-Kung算法虽然步骤更多，但更高效。

算法实现与优化 ⚙️

在实现Brent-Kung算法时，我们应用了共享内存和线程重索引等优化技术。

共享内存的使用

与Kogge-Stone算法类似，使用共享内存可以提高数据重用性。此外，对于Brent-Kung算法，全局内存访问可能无法合并，而使用共享内存可以确保全局内存加载是合并的，从而进一步提升性能。

双缓冲的必要性

在Kogge-Stone算法中，由于同一迭代中有线程读写同一内存位置，我们需要双缓冲来避免冲突。然而，在Brent-Kung算法中，同一迭代中没有线程同时读写同一数据元素，因此不需要双缓冲，可以使用单个缓冲区。

控制流分歧与线程重索引

在Kogge-Stone算法中，控制流分歧不是问题，因为每一步中活跃的线程是连续的。但在Brent-Kung算法中，如果为每个输入值分配一个线程，会导致严重的控制流分歧，因为并非所有输入值在每一步都被处理。

为了避免控制流分歧，我们采用线程重索引技术。我们不再为每个线程固定分配特定的数据元素，而是在每次迭代中根据步长重新分配线程负责的数据元素，确保活跃的线程尽可能连续，从而最大化线程束的执行效率。

线程重索引的关键是计算每个线程在给定步长下负责的数据元素索引。我们推导出索引公式为：
index = (thread_idx + 1) * 2 * stride - 1
其中，thread_idx是线程索引，stride是当前步长。在实现时，需要确保索引在数组边界内。

排他性扫描的实现

对于排他性扫描，有两种实现方式：

将其表述为包含性扫描，即将输入元素左移一位。
修改Brent-Kung算法的后归约阶段，专门实现排他性扫描。

第二种方法通过巧妙地移动和添加元素，在归约阶段后，将最后一个元素（全段和）保存为块的部分和并替换为零，然后在后归约阶段通过一系列操作将零和部分和移动到正确位置，最终得到排他性扫描结果。这种方法在作业中需要实现。

线程粗化优化 🧵

线程粗化是一种优化技术，用于在过度并行化导致工作效率下降时，减少性能损失。在扫描的上下文中，我们可以通过分段扫描应用线程粗化。

具体方法是：在每个线程块内，让每个线程顺序扫描一个子段（例如8个元素），产生该子段的部分和。然后，线程块并行扫描这些线程部分和数组。最后，每个线程将对应的扫描部分和（来自前面的线程）加到其子段的元素上。

这样，我们利用顺序扫描的高工作效率处理子段，仅对部分和进行并行扫描，从而在保持并行性的同时提高了整体算法的工作效率。

以下是应用线程粗化后的性能对比示例：

未使用线程粗化：约4.5毫秒
使用线程粗化（因子为8）：约1.5毫秒
性能提升了约3倍。

总结 📝

本节课我们一起学习了Brent-Kung并行扫描算法。该算法通过归约和后归约两个阶段，以更多的步数换取了更少的操作总数，从而实现了更高的工作效率。我们探讨了其实现中的关键优化技术，包括共享内存的使用和线程重索引以避免控制流分歧。同时，我们分析了排他性扫描的两种实现方式。

最后，我们介绍了线程粗化优化技术，它通过让线程顺序处理数据子段来减少并行算法带来的工作效率损失，并在实际代码中展示了显著的性能提升。

通过对比Kogge-Stone和Brent-Kung算法，我们理解了在并行计算中工作效率与并行度之间的权衡，以及根据具体硬件和问题规模选择合适算法的重要性。

GPU计算：第13讲：直方图 📊

概述

在本节课中，我们将学习一个新的并行模式——直方图。我们将介绍硬件提供的一个新特性：原子操作，并探讨一种新的优化类别：私有化。我们将从直方图的基本概念开始，逐步深入到其并行实现、遇到的问题以及相应的优化策略。

回顾：扫描操作

上一节我们介绍了扫描操作，特别是用于并行扫描的Brent-Kung方法，以及应用于扫描的线程粗化技术及其优势。

我们看到了实现Brent-Kung并行包含性扫描的模式。该模式的第一部分涉及归约步骤，第二部分涉及后归约步骤，该步骤会回溯并更新那些不正确的值，从而得到一个完全扫描的数组。

比较Kogge-Stone和Brent-Kung方法，我们发现Kogge-Stone方法在更少的步骤内完成，而Brent-Kung方法需要更多步骤。然而，Brent-Kung方法在第一步和总体上执行的操作更少，这使其具有更高的工作效率。

总的来说，Kogge-Stone方法在 log n 步内完成，需要 O(n log n) 次操作，而Brent-Kung方法需要 2 log n - 1 步，但只执行 O(n) 次操作，因此比Kogge-Stone方法更高效。这里的工作效率指的是并行算法相对于顺序实现所执行的操作数量。

我们还对扫描应用了各种优化，包括使用共享内存。共享内存除了数据重用外，另一个优势是它实现了合并访问，因为我们在扫描中的访问原本不是合并的。通过以合并的方式加载到共享内存，然后在共享内存中工作，我们避免了内存合并的问题。

我们不需要双缓冲，因为在扫描树的工作方式中，同一迭代中没有元素被不同线程读写。我们还研究了如何通过在每个迭代中重新索引线程，而不是在整个执行过程中为线程分配固定元素，来最小化控制流分歧。此外，我们还探讨了如何使用Brent-Kung方法进行独占扫描。第一种方法是简单地将其表述为独占扫描，就像我们在Kogge-Stone方法中所做的那样。第二种方法是使用不同的后归约步骤，该步骤在归约后保存所有元素的总和，将其替换为零，然后通过另一种结构最终将零传播到开头，计算并更新值，从而得到独占扫描的数组。

我们没有为这个方法编写代码，因为你们应该在相关的作业中实现它。

最后，我们讨论了工作效率，并指出尽管Brent-Kung方法具有更高的理论工作效率，但在实践中，考虑到非活动线程时，其实际资源消耗更接近 O(n log n)。实际上，Brent-Kung方法虽然更高效，但在我们的特定实现中表现相似甚至更差。

我们还讨论了扫描的线程粗化。当我们并行化扫描时，如果工作确实并行运行，我们会牺牲工作效率，但这是值得付出的代价。然而，如果硬件最终将我们的线程块串行化，那么为了并行化而付出的低效率代价就不值得了。因此，我们可以应用线程粗化，让一些线程对其自己的段执行顺序扫描，然后对部分和进行并行扫描，最后将它们重新分配到每个线程拥有的段中。通过这样做，我们提高了算法的工作效率，因为每个线程执行的顺序扫描通常比并行扫描更高效。

这是对上节课内容的快速回顾。

直方图简介

本节我们将讨论一个新的并行模式：直方图。通过直方图，我们将引入硬件提供的一个新特性：原子操作，并探讨一种新的优化类别：私有化。

什么是直方图？

直方图近似表示数据集的分布。其方法是将输入数据集中值可能出现的范围划分为多个区间（有时也称为桶），并计算落在每个区间内的值的数量。

直方图的一个常见应用（当然不是唯一的地方）是颜色直方图。例如，如果有一幅图像，图像由像素组成，这些像素可能有一个取值范围（例如，使用8位整数表示每个像素值，因此像素值可以从0到255）。颜色直方图可以计算落在每个值范围内的像素数量，从而为我们提供图像中颜色分布、亮度或暗度以及明暗变化程度的某种摘要。

顺序实现直方图

顺序实现直方图的方法如下：我们可以简单地使用一个循环，从0遍历到图像的像素数（宽度乘以高度）。对于图像中的每个像素，我们读取像素值。假设图像存储为一维数组 image，其大小为 width * height。对于每个像素，我们读取其值，然后递增 bins 数组中对应索引的值。例如，如果像素值为0，则递增 bins[0]；如果像素值为10，则递增 bins[10]；如果像素值为128，则递增 bins[128]；如果像素值为255（8位整数的最大值），则递增 bins[255]。最终，bins 数组中将包含每个值对应的像素数量。

以下是顺序实现的伪代码：

unsigned char* image; // 输入图像数组
int num_pixels = width * height;
int bins[256] = {0}; // 初始化所有区间计数为0

for (int i = 0; i < num_pixels; i++) {
    unsigned char pixel_value = image[i];
    bins[pixel_value]++; // 递增对应区间的计数
}

并行实现直方图

如何并行执行直方图操作？一种可能的方法是：将图像分配给多个线程，每个线程负责更新自己的区间计数，然后合并这些计数。

本着透明可扩展性的精神，我们可以从尽可能细的粒度开始，为图像中的每个输入像素分配一个线程。然后，每个线程将找到对应的区间并更新它。

以下是实现此方法的初始GPU内核代码框架：

__global__ void histogram_kernel(unsigned char* image, int* bins, int num_pixels) {
    int i = blockIdx.x * blockDim.x + threadIdx.x; // 计算线程负责的像素索引
    if (i < num_pixels) { // 边界检查
        unsigned char b = image[i]; // 读取像素值
        bins[b]++; // 递增对应区间的计数
    }
}

然而，这个实现存在一个问题：多个线程可能同时访问并尝试更新 bins 数组中的同一个内存位置，导致数据竞争。

数据竞争与原子操作

数据竞争问题

当多个线程在没有协调的情况下并发访问同一内存位置，并且至少有一个访问是写操作时，就会发生数据竞争。数据竞争可能导致不可预测的程序输出。

在我们的例子中，bins[b]++ 这个操作不是单一的原子操作。它由多个指令组成：加载 bins[b] 的旧值，将该值加1，然后存储新值。当多个线程同时执行这个操作序列时，指令的交错执行会影响最终结果的正确性。

例如，假设 bins[b] 初始为0，两个线程A和B都试图递增它。如果线程A完全在线程B之前运行，或者线程B完全在线程A之前运行，最终 bins[b] 会变为2，这是正确的结果。然而，如果线程A加载了旧值0，计算新值为1，但在存储之前，线程B也加载了旧值0（此时仍为0），然后线程A存储了1，线程B计算新值为1并存储1，那么最终 bins[b] 只变为1，而不是2。这就是数据竞争导致的问题。

解决方案：原子操作

为了避免数据竞争，对同一内存位置的并发读-修改-写操作需要以互斥的方式进行排序。

在CPU上，一种常见的方法是使用互斥锁（mutex）。然而，在GPU上使用锁可能导致死锁，特别是在同一warp内的线程尝试获取同一个锁时，因为warp遵循SIMD模型，线程必须一起执行指令。

GPU提供了更好的解决方案：原子操作。原子操作是在GPU上执行读-修改-写操作的单个ISA指令。硬件保证在此操作完成之前，没有其他线程可以访问该内存位置。如果对同一内存位置有并发的原子操作，硬件会将这些操作序列化。

CUDA提供了多种原子操作，例如 atomicAdd、atomicSub、atomicMin、atomicMax 等。对于直方图，我们使用 atomicAdd。

atomicAdd 函数的原型如下：

int atomicAdd(int* address, int val);
unsigned int atomicAdd(unsigned int* address, unsigned int val);
unsigned long long int atomicAdd(unsigned long long int* address, unsigned long long int val);
float atomicAdd(float* address, float val);
double atomicAdd(double* address, double val);

它接受一个指向内存位置的指针和一个要加的值，将该值加到内存位置的值上，并返回旧值。

我们可以将之前的 bins[b]++ 替换为原子操作：

__global__ void histogram_kernel(unsigned char* image, int* bins, int num_pixels) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < num_pixels) {
        unsigned char b = image[i];
        atomicAdd(&bins[b], 1); // 使用原子操作安全地递增
    }
}

这样，即使多个线程同时更新同一个区间，也能保证计数的正确性。

原子操作可以在全局内存和共享内存上执行。

性能优化：私有化

私有化的概念

尽管原子操作解决了正确性问题，但性能可能并不理想。全局内存的原子操作延迟很高，特别是当多个线程频繁竞争更新同一内存位置时。

私有化是一种优化技术，用于减少对共享输出的竞争。其核心思想是：不为所有线程维护一个全局的直方图副本，而是为每个线程块创建一个私有的直方图副本。线程块内的线程更新自己的私有副本，从而将竞争限制在块内。当线程块完成计算后，再将私有副本的结果合并到全局副本中。

这样做的好处是：

减少全局竞争：线程块内的竞争发生在私有副本上，减少了全局内存原子操作的次数和序列化开销。
利用更快的内存：如果私有副本足够小，可以将其放入共享内存中。共享内存的访问速度远快于全局内存，原子操作在共享内存上的开销也更低。

私有化的实现步骤

以下是实现私有化的基本步骤：

在共享内存中声明一个私有直方图数组。
初始化共享内存中的直方图为零。
每个线程读取其负责的像素，并使用原子操作更新共享内存中的私有直方图。
使用同步（如 __syncthreads()）确保所有线程完成更新。
最后，由线程块内的线程合作，将私有直方图中的非零值通过原子操作添加到全局直方图中。

伪代码框架如下：

__global__ void histogram_privatized_kernel(unsigned char* image, int* global_bins, int num_pixels) {
    __shared__ int shared_bins[256]; // 在共享内存中声明私有直方图
    // 初始化共享内存直方图为零
    for (int i = threadIdx.x; i < 256; i += blockDim.x) {
        shared_bins[i] = 0;
    }
    __syncthreads();

    // 每个线程处理多个像素（线程粗化）
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    for (int i = tid; i < num_pixels; i += stride) {
        unsigned char b = image[i];
        atomicAdd(&shared_bins[b], 1); // 更新共享内存中的私有副本
    }
    __syncthreads();

    // 将私有副本合并到全局内存
    for (int i = threadIdx.x; i < 256; i += blockDim.x) {
        if (shared_bins[i] != 0) {
            atomicAdd(&global_bins[i], shared_bins[i]);
        }
    }
}

结合线程粗化

我们可以进一步结合线程粗化来优化。线程粗化是指减少线程块的数量，让每个线程处理更多的输入元素。

这样做的好处是：

减少私有副本数量：更少的线程块意味着更少的私有直方图副本，从而减少了最终合并到全局直方图所需的原子操作次数。
提高资源利用率：当硬件资源有限，线程块可能被串行执行时，减少线程块数量可以降低并行化开销，提高工作效率。

在实现时，需要确保线程以合并访问的方式加载输入数据，以保持内存访问效率。

总结

本节课我们一起学习了直方图这一并行模式。我们从直方图的基本概念和顺序实现开始，探讨了其并行化方案。

我们遇到了并行化中的关键挑战：数据竞争。为了解决这个问题，我们引入了GPU的原子操作，它能够以硬件保证的互斥方式安全地更新共享内存位置。

为了提高性能，我们深入探讨了私有化这一优化技术。通过为每个线程块创建私有直方图副本，我们显著减少了对全局内存的竞争。结合共享内存存储私有副本，可以进一步降低访问延迟。此外，线程粗化技术可以与私有化结合，通过减少线程块数量来降低全局合并的开销，从而在特定硬件条件下提升整体效率。

这些概念和优化策略是构建高效GPU程序的重要组成部分。

GPU计算：第14讲：合并模式 🧩

在本节课中，我们将要学习一个新的并行计算模式：合并。我们将探讨如何将两个已排序的数组合并成一个新的有序数组，并学习如何在GPU上高效地实现这一操作。

回顾：直方图模式

上一节我们介绍了直方图模式。在实现直方图时，我们引入了两个新概念：原子操作 和 私有化 优化。

直方图用于近似数据的分布，它统计数据集中落入一系列“箱子”的值的数量。例如，一个颜色直方图统计图像中每个像素值（0-255）出现的次数。

当多个线程并发地更新同一个直方图箱子时，会出现数据竞争问题。数据竞争发生在多个线程同时访问同一内存位置，且至少有一个访问是写操作时。这可能导致更新丢失。

在CPU上，我们通常使用互斥锁来解决这个问题。但在GPU上，锁可能导致死锁，尤其是在同一个warp内的线程之间。

因此，GPU提供了原子操作，例如 atomicAdd，它通过一条指令完成“读取-修改-写入”操作，保证了操作的原子性。

然而，让所有线程直接更新全局直方图会导致大量原子操作，性能低下。我们采用了私有化优化：每个线程块拥有直方图的私有副本，线程在私有副本上更新，最后再将非零的箱子提交到全局直方图。如果直方图足够小，私有副本可以放在共享内存中，以降低原子操作的延迟。此外，我们还应用了线程协同，让每个线程块处理更大的输入段，从而减少线程块数量和最终的全局原子操作总数。

介绍合并模式

本节中，我们来看看合并模式。合并是一种操作，它接收两个已排序的列表，并将它们组合成一个单一的、有序的列表。

例如，给定两个有序数组 A = [1, 2, 4, 6, 8, 9] 和 B = [2, 3, 5, 7, 9]，合并后的结果 C = [1, 2, 2, 3, 4, 5, 6, 7, 8, 9, 9]。

顺序合并实现

在讨论并行化之前，我们先实现一个顺序合并算法。以下是其核心逻辑：

void mergeSequential(int* A, int m, int* B, int n, int* C) {
    unsigned int i = 0; // 数组A的索引
    unsigned int j = 0; // 数组B的索引
    unsigned int k = 0; // 数组C的索引

    // 合并两个数组，直到其中一个耗尽
    while (i < m && j < n) {
        if (A[i] <= B[j]) {
            C[k++] = A[i++];
        } else {
            C[k++] = B[j++];
        }
    }
    // 将剩余元素复制到C中
    while (i < m) C[k++] = A[i++];
    while (j < n) C[k++] = B[j++];
}

该算法使用三个索引 i, j, k 分别遍历数组A、B和输出数组C。它比较 A[i] 和 B[j]，将较小的值放入 C[k]，并递增相应的索引。当一个数组耗尽后，将另一个数组的剩余元素直接复制到C中。

并行化合并的思路

合并操作看起来是顺序的，如何并行化呢？关键在于对输出数组进行分区，而不是对输入数组分区。

如果我们均匀地划分输入数组A，那么A中的一小段可能对应B中很长的一段，导致线程间工作量不均衡。相反，如果我们均匀地划分输出数组C，并让每个线程负责合并写入C的某一段，那么每个线程的工作量（写入C的元素数量）是相同的，从而获得更好的负载均衡。

因此，并行合并的策略是：

将输出数组C划分为大小相等的段。
每个线程负责将对应C段所需的元素从A和B中合并出来。

核心挑战在于：给定一个输出索引K（C段的起始位置），如何找到在输入数组A和B中对应的起始索引I和J？ 我们称I和J为K的协同秩。

计算协同秩

我们的目标是：给定K，找到对应的I和J。我们有以下关系：

公式：K = I + J

因为C中K之前的元素数量，等于A中I之前的元素数量加上B中J之前的元素数量。因此，如果我们能找到I，就可以通过 J = K - I 计算出J。

确定I的边界

首先，我们需要确定I的可能取值范围：

I必须是A的有效索引：0 <= I <= m（I可以等于m，表示A中所有元素都已用完）。
由于 J = K - I 且 0 <= J <= n，我们可以推导出：K - n <= I <= K。

综合这两个条件，I的边界为：

下界：low = max(0, K - n)
上界：high = min(m, K)

公式：low = max(0, K - n), high = min(m, K)

在边界内搜索正确的I

在边界 [low, high] 内，我们需要找到满足合并条件的I。正确的I满足以下两个条件：

A[I-1] <= B[J] （如果I>0且J<n）
B[J-1] <= A[I] （如果J>0且I<m）

这意味着A中最后一个被取出的元素应小于等于B中第一个将被取出的元素，反之亦然。

由于数组A和B都是有序的，我们可以在边界 [low, high] 内进行二分查找来定位满足条件的I。在二分查找的每一步，我们检查当前猜测值 I_guess：

如果 A[I_guess-1] > B[J_guess]，说明猜测的I太大，应调低上界。
如果 B[J_guess-1] > A[I_guess]，说明猜测的I太小，应调高下界。
否则，猜测正确。

以下是协同秩函数的核心代码框架：

__device__ unsigned int coRank(int* A, int m, int* B, int n, unsigned int K) {
    unsigned int low = (K > n) ? (K - n) : 0;
    unsigned int high = (m < K) ? m : K;

    while (low < high) {
        unsigned int I = (low + high) / 2; // 猜测的I
        unsigned int J = K - I; // 对应的J

        // 检查边界条件，然后进行比较
        bool tooHigh = (I > 0 && J < n && A[I-1] > B[J]);
        bool tooLow = (J > 0 && I < m && B[J-1] > A[I]);

        if (tooHigh) {
            high = I; // 猜测太大，降低上界
        } else if (tooLow) {
            low = I + 1; // 猜测太小，提高下界
        } else {
            return I; // 找到正确的I
        }
    }
    return low;
}

基本的并行合并内核

有了协同秩函数，我们可以实现基本的并行合并内核。每个线程需要：

计算自己负责的输出段起始索引 K。
调用 coRank 函数找到对应的输入起始索引 I 和 J。
计算下一个线程的K值 (K_next)，以确定本线程输入段的大小 (I_next - I 和 J_next - J)。
调用顺序合并函数，合并 A[I : I_next] 和 B[J : J_next] 到 C[K : K_next]。

以下是内核的简化结构：

__global__ void parallelMergeBasic(int* A, int m, int* B, int n, int* C) {
    unsigned int elementsPerThread = 6; // 每个线程处理的元素数
    unsigned int K = (blockIdx.x * blockDim.x + threadIdx.x) * elementsPerThread;

    if (K >= m + n) return; // 边界检查

    // 1. 找到本线程输出段在输入数组中的起始位置
    unsigned int I = coRank(A, m, B, n, K);
    unsigned int J = K - I;

    // 2. 找到下一个线程的K，以确定本段大小
    unsigned int K_next = min(K + elementsPerThread, m + n);
    unsigned int I_next = coRank(A, m, B, n, K_next);
    unsigned int J_next = K_next - I_next;

    // 3. 执行顺序合并
    mergeSequential(&A[I], I_next - I, &B[J], J_next - J, &C[K]);
}

优化：使用共享内存改善合并访问

基本的并行合并存在性能问题：

线程发散：每个线程独立进行二分查找，路径长度不同。
非合并内存访问：二分查找导致对全局内存的随机访问，且线程间的访问模式不连续。

优化思路是在线程块级别进行协作：

块级协同秩：仅让线程块内的一个线程（如thread 0）计算整个线程块输出段对应的输入段边界 (I_block, J_block, I_next_block, J_next_block)。
合并加载：线程块内所有线程协作，将 A[I_block : I_next_block] 和 B[J_block : J_next_block] 连续地、合并地加载到共享内存中。
共享内存内操作：每个线程在共享内存中的数组副本上，进行二分查找和顺序合并。这样，所有的随机访问都发生在低延迟的共享内存中。
合并存储：所有线程协作，将共享内存中合并好的结果段，连续地、合并地写回全局内存。

这种方法将昂贵的非合并全局内存访问，转换为了快速的共享内存访问，并充分利用了内存带宽。

优化的内核结构概述如下：

计算线程块负责的输出段范围 K_block 和 K_next_block。
由单个线程计算块级的协同秩 I_block, J_block, I_next_block, J_next_block，并存入共享内存。
同步线程。
所有线程协作，将A和B的对应段加载到共享内存。
同步线程。
每个线程在共享内存的数据上，计算自己负责的子段的协同秩并进行合并。
所有线程协作，将共享内存中的结果写回全局内存。

总结

本节课中我们一起学习了合并模式。我们从顺序合并算法出发，探讨了如何在GPU上并行化这一操作。关键点在于对输出数组进行分区，并通过计算协同秩来找到每个输出段对应的输入段。我们实现了基本的并行合并内核，并分析了其性能瓶颈。最后，我们介绍了一种利用共享内存和线程块内协作的优化策略，通过将数据加载到共享内存来将非合并的全局内存访问转化为快速的共享内存访问，从而显著提升性能。

合并是许多算法（如排序、数据库操作）的基础组件，掌握其高效的GPU实现对于高性能计算至关重要。

GPU计算：15：排序

在本节课中，我们将学习两种重要的并行排序算法：基数排序和归并排序。我们将重点探讨基数排序的原理、并行化方法及其优化策略，并简要介绍如何将之前学习的并行归并模式应用于归并排序。

概述：从归并到排序

上一节我们介绍了并行归并模式。我们提到，有序归并是将两个有序列表合并成一个有序列表的操作。为了并行化此操作，一种方法是将要生成的输出列表划分为相等的段，并为每个线程分配其中一段。然后，每个线程需要在两个输入数组A和B中找到对应的输入段，并顺序合并这两个输入段。并行归并的关键挑战在于每个线程如何在A和B中找到对应的输入段。

为了解决这个问题，我们引入了协同秩的概念。给定一个指向输出数组C中元素（或段起始点）的索引K，关键挑战是找到指向A和B中对应元素的索引I和J。我们称I为K在A中的协同秩，J为K在B中的协同秩。我们观察到J就是K减去I，因此真正的挑战在于找到I。我们通过为I设定边界，然后在此边界内进行二分搜索来确定I。判断猜测正确、过高或过低的条件基于一个观察：A中I之前的元素应小于B中J之后的元素；同样，B中J之前的元素应小于A中I之后的元素。

每个负责一个输出段的线程都会进行二分搜索来确定I，然后相应地找到J，最后将两个输入段顺序合并到输出段中。此操作的内存访问不是合并的，因为二分搜索涉及对数组的随机访问，并且在顺序合并阶段，每个线程合并的是输入数组中的不同段，访问模式也不连续。

为了优化，我们让块中的一个线程确定整个块对应的输入段（I_block和J_block），然后块内的线程协作将这些段以合并访问的方式加载到共享内存中。接着，所有线程在共享内存中进行各自的二分搜索和顺序合并。最后，我们将合并后的结果从共享内存以合并访问的方式写回全局内存。这样，我们通过合并访问加载到共享内存，在共享内存中进行非合并访问的操作，再以合并访问写回全局内存，从而优化了内存访问。

我们还讨论了线程粗化的代价。并行化归并的代价是每个线程都需要进行自己的二分搜索。线程越多，二分搜索操作就越多。通过为每个线程分配多个输出元素（即一个输出段），我们已经应用了线程粗化，从而分摊了二分搜索的成本。

本节中，我们将开始讨论排序算法，主要关注基数排序，并简要介绍归并排序。

基数排序原理 🔢

基数排序是一种排序算法，其工作原理是基于基数（或进制）将待排序的键值分配到桶中。如果输入键值在某种位置记数系统中以某个基数为表示，那么我们就使用该基数，并迭代键值的每一位数字。每次迭代，我们都根据当前处理的数字将键值分配到对应的桶中。

对键值的每一位数字重复此分配过程。一个重要点是，每次迭代时，每个桶内应保持前一次迭代产生的顺序。这很关键，我们稍后会解释原因。在计算机中实现基数排序时，我们通常喜欢使用2的幂作为基数，因为这简化了对二进制数的处理。每次迭代处理键值中一个固定的比特位切片。

为了简单和便于说明，我们将从基数为2（即1比特）的情况开始。然后在本节后续部分再扩展到更大的基数。

单比特基数排序演示

假设我们有一个数组需要排序。我们需要将数组的键值按比特位划分，并每次迭代处理一个比特位（因为我们使用的是1比特基数）。每次迭代，我们将根据当前处理的比特位将输入分配到桶中。

首先，我们查看每个数字的最低有效位。有些是0，有些是1。我们将根据这个比特位将键值分离到两个桶中：将所有比特位为0的放在一起，将所有比特位为1的放在后面。完成这一步后，数组仅按最低有效位排序。

接着，我们重复这个过程，但每次查看不同的比特位，从低有效位向高有效位进行。在下一次迭代中，我们查看第二个比特位，并根据该比特位是0还是1来分离键值。这里的关键是，在分离时，我们在每个桶内保持了原始顺序。这样做的原因是，如果我们打乱顺序，就会破坏前一次迭代得到的结果。前一次迭代已按第一个比特位排序，第二次迭代在按第二个比特位排序的同时，保持了第一次迭代的顺序。这确保了键值现在按较低的两位比特排序。

然后我们对第三个比特位重复此过程，依此类推，直到处理完所有比特位。最终，我们得到一个完全排序的数组。

基数排序是一种非比较排序算法。它通常要求键值具有固定大小，以便可以划分为数字位。如果键值不满足此条件，则需要使用基于比较的排序算法，如我们稍后将看到的归并排序。

并行化基数排序的关键步骤 ⚙️

现在，我们讨论如何并行化基数排序。核心问题在于：如何根据一个特定的比特位分离键值，同时保持顺序？换句话说，如何获取一个数组，将所有当前比特位为0的元素放在左边，为1的元素放在右边，并保持它们原有的相对顺序？

关键在于为每个元素找到其在输出（排序后）数组中的目标索引。

计算目标索引

让我们分析如何计算目标索引。

对于比特位为0的元素：其目标索引等于它左边0的个数。左边0的个数又等于该元素左边的总元素数减去左边1的个数。而左边的总元素数就是该元素的索引位置。因此，目标索引(0) = 元素索引 - 左边1的个数。
对于比特位为1的元素：其目标索引等于数组中0的总数加上该元素左边1的个数。0的总数等于数组总大小减去1的总数。因此，目标索引(1) = 数组大小 - 1的总数 + 左边1的个数。

由此可见，无论是0还是1，计算目标索引都需要知道每个元素左边1的个数。此外，对于比特位为1的元素，还需要知道整个数组中1的总数。

使用扫描模式

那么，如何找到每个元素左边1的个数呢？这正是独占扫描模式可以解决的问题。

因此，并行化基数排序（单次迭代，处理一个比特位）的逻辑步骤如下：

提取比特位：每个线程负责一个输入元素，提取出当前要处理的比特位（0或1），形成一个比特位数组。
执行独占扫描：线程协作对这个比特位数组执行独占扫描操作。结果数组的每个位置就存储了对应输入元素左边1的个数。同时，扫描过程也能得到整个数组中1的总数。
计算目标索引：每个线程根据其元素的比特位是0还是1，应用上述公式，利用元素索引、左边1的个数、数组大小和1的总数来计算目标索引。
分散写入：每个线程将其负责的输入元素写入到输出数组的对应目标索引处。

优化：改善合并访问 🚀

上述基本并行方法存在效率问题。在最后一步“分散写入”中，每个线程写入全局内存的位置是分散的，导致内存访问不是合并的，这会严重影响性能。

为了优化合并访问，我们可以利用共享内存。思路是：先在线程块内部进行排序和分离，然后再以合并访问的方式将整个桶写入全局内存。

以下是优化步骤：

局部排序：每个线程块将其负责的输入元素加载到共享内存中。块内线程协作，执行上述的扫描和分离步骤，但这次是在共享内存中局部进行。这样，每个线程块在共享内存中将其元素分离到0桶和1桶中。
计算全局桶地址：每个线程块需要知道它的0桶和1桶应该写入全局输出数组的哪个起始位置。为此，每个块将其局部0的个数和1的个数写入一个全局数组。然后，对这个全局数组执行一次独占扫描。扫描结果告诉每个块，在它之前所有块的0的总数（即其0桶的起始地址），以及全局0的总数加上在它之前所有块的1的总数（即其1桶的起始地址）。
合并写入：每个线程块将其共享内存中的0桶和1桶分别整体地、连续地写入到全局内存中计算好的起始地址处。由于块内线程写入的是连续的内存区域，因此访问是合并的。

这种优化将非合并的分散写入，转变为先进行局部（共享内存）非合并操作，再进行全局合并写入，从而显著提高了内存访问效率。

基数选择与线程粗化 ⚖️

之前我们一直使用1比特基数。对于一个32位的键，需要32次迭代。为了减少迭代次数，我们可以使用更大的基数，例如2比特、4比特等。更大的基数意味着更少的迭代，但同时也意味着更多的桶（例如2比特有4个桶，3比特有8个桶）。

更多的桶会恶化我们刚才优化的合并写入效果，因为一个线程块可能需要将数据写入全局内存中更多个不连续的位置。

因此，基数的选择需要在迭代次数和合并访问行为之间取得平衡。我们希望基数足够大以减少迭代，但又不能太大以免严重损害合并访问性能。

另一种优化合并访问的策略是线程粗化。之前我们让每个线程处理一个元素。线程粗化是指让每个线程处理多个元素。这样，每个线程块处理的总元素数增加，意味着每个局部桶的大小也增加。当块将大桶写入全局内存时，合并访问的效果会更好，因为连续写入的数据块更长。

线程粗化通过让每个线程承担更多工作，分摊了扫描等操作的代价，并提高了内存访问的局部性。

归并排序简介 🤝

基数排序并非适用于所有类型的数据（例如可变长度键）。我们还需要一种适用于并行化的基于比较的排序算法，其中一种流行的方法是归并排序。

归并排序是一种分治算法：

分解：将列表递归地分成两半，直到子列表足够小（例如，只有一个元素）。
解决：排序这些子列表（对于很小的列表，排序很简单）。
合并：将已排序的子列表两两合并，形成更大的有序列表，直到最终合并成完整的有序列表。

并行化归并排序

我们可以利用之前学过的并行归并模式来并行化归并排序：

在底层，可以使用简单的排序方法（甚至可以在线程块内部使用基数排序）对小数组进行排序。
在合并阶段，可以将不同的合并任务分配给不同的线程块并行执行。例如，在第一次合并阶段，多个线程块可以同时合并不同的相邻有序子数组对。
每个合并操作本身也可以并行化，正如我们在并行归并模式中学到的那样：将输出段分配给不同的线程块，每个块通过协同秩计算找到对应的输入段并进行合并。

随着合并的进行，需要合并的有序段越来越大，但段的数量越来越少。因此，并行性从“跨多个合并操作”逐渐转向“在单个大型合并操作内部”。在后期阶段，我们使用更多的线程块来并行化单个大型合并操作。

总结 📚

本节课我们一起学习了两种重要的并行排序算法。

我们深入探讨了基数排序，它是一种非比较排序算法。我们解释了其按比特位迭代分配的原理，并详细分析了如何并行化其核心步骤——根据特定比特位分离键值并保持顺序。关键点在于使用独占扫描来计算每个元素的目标索引。我们还讨论了如何通过共享内存局部排序和全局桶地址计算来优化全局内存的合并访问。此外，我们探讨了基数大小选择的权衡以及线程粗化对性能的积极影响。

最后，我们简要介绍了归并排序，说明了如何利用之前学习的并行归并模式，通过分治和并行合并来实现高效的并行排序。

理解这些排序算法的并行化策略及其优化，对于在GPU上高效处理大规模数据至关重要。

GPU计算：第16讲：稀疏矩阵计算（COO与CSR格式） 🧮

在本节课中，我们将要学习稀疏矩阵计算，并重点研究两种存储格式：坐标格式（COO）和压缩稀疏行格式（CSR）。我们将以稀疏矩阵-向量乘法（SPMV）作为案例，探讨不同格式如何影响并行计算的性能、内存访问模式和负载均衡。

什么是稀疏矩阵？ 🤔

在开始讨论存储格式之前，我们首先需要理解什么是稀疏矩阵。

一个稠密矩阵是指矩阵中大多数元素都是非零的矩阵。相反，一个稀疏矩阵则是指矩阵中包含大量零元素的矩阵。在实际应用中，许多系统产生的矩阵都是稀疏的，这意味着矩阵中绝大部分元素都是零。

利用矩阵的稀疏性，我们可以获得多方面的优势：

节省内存：我们可以压缩矩阵，只存储非零值，从而减少内存占用。
节省存储空间：将矩阵存储到文件时，可以占用更少的磁盘空间。
节省内存带宽：计算时无需加载零值，减少了内存访问量。
节省计算时间：无需对零值进行计算，提高了计算效率。

因此，稀疏矩阵计算的核心挑战之一就是设计高效的存储格式，在节省内存的同时，还能优化计算性能。

稀疏矩阵存储格式的设计考量 ⚖️

设计稀疏矩阵存储格式时，需要考虑多个因素：

空间效率：格式能节省多少内存。
灵活性：格式是否便于添加、重新排序或删除矩阵中的值。
可访问性：格式是否便于以特定方式（如按行或按列）访问所需数据。
内存访问模式：格式是否有利于实现内存合并访问（特别是在GPU上）。
负载均衡：格式是否有助于平衡并行工作线程的工作量，在GPU上这关系到控制流发散的程度。

不同的格式在这些考量上各有优劣。在本课程中，我们将以稀疏矩阵-向量乘法（SPMV） 作为具体应用场景，来分析和比较不同格式。

SPMV的运算可以表示为：
y = A * x
其中 A 是一个稀疏矩阵，x 是一个稠密输入向量，y 是结果稠密输出向量。

接下来，让我们看看第一种存储格式。

坐标格式（COO） 📍

COO格式是一种直观的存储方式。对于矩阵中的每一个非零元素，我们存储三个信息：其值、行索引和列索引。

具体来说，我们会创建三个数组：

values[]: 按顺序存储所有非零元素的值。
row_indices[]: 存储每个非零元素对应的行索引。
col_indices[]: 存储每个非零元素对应的列索引。

例如，对于一个稀疏矩阵，其COO表示如下：

值:     [1, 7, 5, 3, 9, 2, 8]
行索引: [0, 0, 1, 1, 1, 2, 2]
列索引: [0, 1, 0, 2, 3, 1, 2]

使用COO格式实现SPMV

在COO格式下，实现SPMV最自然的并行化策略是为每一个非零元素分配一个线程。

每个线程的执行步骤如下：

根据线程ID获取对应的非零元素索引 i。
从 row_indices[i] 和 col_indices[i] 中读取该元素的行号 row 和列号 col。
从 values[i] 中读取该元素的值 val。
从输入向量 x 中读取 x[col]。
计算 val * x[col]，并使用原子操作将结果累加到输出向量 y[row] 中。

需要使用原子操作的原因是，同一行中的多个非零元素可能被不同的线程处理，而这些线程会尝试同时更新同一个输出位置 y[row]。

以下是该内核函数的简化代码示例：

__global__ void spmv_coo_kernel(COOMatrix coo, float* x, float* y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x; // 线程对应的非零元素索引
    if (i < coo.num_nonzeros) {
        int row = coo.row_indices[i];
        int col = coo.col_indices[i];
        float val = coo.values[i];
        atomicAdd(&y[row], val * x[col]); // 使用原子加操作
    }
}

COO格式的优缺点分析

上一节我们介绍了COO格式的基本原理和SPMV实现，本节我们来总结其优缺点。

优点：

灵活性高：非零元素可以任意顺序存储，添加新元素非常容易（只需追加到数组末尾）。
可访问性：给定一个非零元素，可以立即获取其行和列坐标，便于跨非零元素并行。
内存合并访问：线程按顺序访问 values、row_indices、col_indices 数组，访问模式是合并的。
无控制流发散：每个线程处理一个非零元素，工作量相同。

缺点：

空间效率较低：需要为每个非零元素存储行和列两个索引。
可访问性局限：给定一个行号，很难快速找到该行所有的非零元素（需要搜索）。
需要原子操作：多个线程可能更新同一输出位置，必须使用原子操作，这会带来性能开销。

为了克服COO格式需要原子操作的缺点，我们引入了第二种格式。

压缩稀疏行格式（CSR） 📊

CSR格式旨在优化按行访问的模式。它将同一行的所有非零元素连续存储在一起。

CSR格式使用三个数组：

values[]: 按行主序存储所有非零元素的值（第一行的非零值，接着是第二行的，依此类推）。
col_indices[]: 存储每个非零元素对应的列索引。
row_ptr[]: 存储每行非零元素在 values 和 col_indices 数组中的起始位置。row_ptr[i] 指向第 i 行的起始，row_ptr[i+1] 指向第 i 行的结束，因此该数组长度为 行数 + 1。

沿用之前的例子，其CSR表示如下：

值:        [1, 7, 5, 3, 9, 2, 8]
列索引:    [0, 1, 0, 2, 3, 1, 2]
行指针:    [0, 2, 5, 7] // 第0行从索引0开始，有2个元素；第1行从索引2开始，有3个元素...

使用CSR格式实现SPMV

在CSR格式下，实现SPMV最自然的并行化策略是为输出向量 y 的每一个元素（即矩阵的每一行）分配一个线程。

每个线程的执行步骤如下：

根据线程ID获取负责的行号 row。
通过 row_ptr[row] 和 row_ptr[row+1] 确定该行非零元素在数组中的起止位置。
循环遍历该行所有的非零元素：
- 读取列号 col = col_indices[k] 和值 val = values[k]。
- 计算 val * x[col]，并累加到一个本地寄存器中。
循环结束后，将累加和写入 y[row]。由于每个输出位置仅由一个线程写入，因此不需要原子操作。

以下是该内核函数的简化代码示例：

__global__ void spmv_csr_kernel(CSRMatrix csr, float* x, float* y) {
    int row = blockIdx.x * blockDim.x + threadIdx.x; // 线程负责的行号
    if (row < csr.num_rows) {
        float sum = 0.0f;
        int row_start = csr.row_ptr[row];
        int row_end = csr.row_ptr[row + 1];
        for (int j = row_start; j < row_end; j++) {
            int col = csr.col_indices[j];
            float val = csr.values[j];
            sum += val * x[col];
        }
        y[row] = sum; // 无需原子操作
    }
}

CSR格式的优缺点分析

上一节我们介绍了CSR格式和基于它的SPMV实现，本节我们来分析其利弊。

优点：

空间效率高：相比COO，节省了存储行索引的空间（仅用 row_ptr 数组）。
可访问性：给定行号，可以快速定位并遍历该行所有非零元素，便于跨行并行。
无需原子操作：每个输出元素由独立线程计算和写入。

缺点：

灵活性低：添加或删除非零元素可能涉及大量数据移动和 row_ptr 更新，成本高昂。
可访问性局限：给定一个非零元素，很难找到其行号（需在 row_ptr 中搜索）；给定列号，也很难找到该列所有非零元素。
内存访问未合并：不同行的线程访问 values 和 col_indices 数组的偏移量不同，访问模式是分散的，不利于合并。
存在控制流发散：不同行的非零元素数量不同，导致线程的工作量不同，运行时间不一致，造成GPU利用率下降。

此外，还存在一种压缩稀疏列格式（CSC），它是CSR的转置版本，将非零元素按列连续存储。CSC格式便于按列访问，但在SPMV场景下不如CSR常用。

总结 📝

本节课我们一起学习了稀疏矩阵计算的基础和两种经典存储格式。

我们首先了解了稀疏矩阵的概念及其优势。然后，我们深入探讨了坐标格式（COO） 和压缩稀疏行格式（CSR）。COO格式简单灵活，易于实现跨非零元素的并行，但需要原子操作。CSR格式空间效率更高，且避免了原子操作，但导致了未合并的内存访问和控制流发散。

每种格式都在空间、灵活性、可访问性和性能之间做出了不同的权衡。在实践中，没有一种格式在所有情况下都是最优的，需要根据具体的矩阵特性和计算任务进行选择。

在下一讲中，我们将继续探索其他稀疏矩阵存储格式，如ELLPACK（ELL）和锯齿状对角线存储（JDS），它们旨在改善CSR在内存合并和负载均衡方面的不足。

GPU计算：第17讲：稀疏矩阵计算（ELL与JDS格式）🎯

在本节课中，我们将继续学习稀疏矩阵计算，重点介绍两种新的稀疏矩阵存储格式：ELL格式和JDS格式。我们将探讨它们的设计原理、在稀疏矩阵向量乘法中的应用，并分析其性能优劣。

上一讲我们介绍了稀疏矩阵计算，并重点讨论了COO和CSR存储格式。本节中，我们来看看另外两种旨在优化GPU上稀疏矩阵向量乘法的格式：ELL和JDS。

回顾与概述 📋

稀疏矩阵中大部分元素为零，我们可以利用这一特性，不存储零元素以节省内存，也不加载和计算它们以节省内存带宽和计算时间。选择存储格式是优化稀疏矩阵计算的关键考虑因素，因为它决定了我们可以采用的计算方法。

我们主要关注稀疏矩阵向量乘法。评估存储格式时，我们考虑以下几个因素：

空间效率：节省了多少存储空间。
灵活性：添加或重新排序元素的难易程度。
可访问性：使用该格式可以轻松访问哪些信息。
内存访问模式：是否支持合并内存访问。
负载均衡：在GPU上下文中，是否有助于最小化线程束分化。

ELL格式详解 🧱

ELL格式旨在改善CSR格式在SPMV中的内存合并访问问题。

格式原理

ELL格式的核心思想是对行进行填充，使得所有行具有相同数量的元素（包括非零元和填充元素），然后以列主序存储这个填充后的二维数组。

假设我们有一个稀疏矩阵，处理步骤如下：

按行分组：像CSR一样，将每行的非零元分组。
填充行：用无效的填充元素（如0或特殊标记）将每行填充至相同的长度，该长度等于矩阵中最大非零元数量。
列主序存储：将填充后的二维数组按列优先的顺序扁平化存储为一维数组。这意味着内存中首先存储所有行的第一个元素，然后是所有行的第二个元素，依此类推。

公式：对于存储在列主序中的二维数组，元素 (row, col) 在一维数组中的索引 i 可以通过以下公式计算：
i = col * num_rows + row
其中，col 是元素在填充行中的列索引（从0开始），num_rows 是矩阵的总行数。

SPMV并行化策略

在ELL格式上执行SPMV的典型策略是为每一行分配一个线程。每个线程负责计算其对应行的输出向量元素。

代码：线程的计算循环核心逻辑如下（伪代码）：

int row = threadIdx.x + blockIdx.x * blockDim.x; // 线程负责的行
if (row < num_rows) {
    float sum = 0.0f;
    for (int iter = 0; iter < num_nonzeros_per_row[row]; ++iter) { // 仅遍历实际非零元
        int idx = iter * num_rows + row; // 计算在扁平化数组中的索引
        int col = ell_col_indices[idx];
        float val = ell_values[idx];
        sum += input_vector[col] * val;
    }
    output_vector[row] = sum;
}

性能与特性分析

与CSR和COO格式相比，ELL格式具有以下特点：

内存合并访问：由于线程按行组织，且数据按列主序存储，当线程访问其行的第一个元素时，这些访问是连续的，从而实现了良好的内存合并。后续迭代中，线程组同步移动到下一列元素，访问模式同样连续。
空间效率：较差。由于填充，可能会浪费大量存储空间，特别是当某些行的非零元数量远多于其他行时。
灵活性：较好。只要某一行未达到最大非零元数量（即还有填充位），就可以较容易地添加新元素（替换填充位）。删除元素也可通过将其置为填充位实现。
可访问性：优秀。结合了COO和CSR的优点：
- 给定行，可以找到该行的所有非零元（通过固定步长 num_rows 跳跃）。
- 给定非零元在数组中的索引，可以轻松找到其行索引（行索引 = 索引 % num_rows）和列索引。
控制流分化：与CSR类似。虽然通过跟踪每行实际非零元数可以避免处理填充元素，但不同行的线程仍然执行不同次数的循环迭代，导致线程束分化。

过渡：ELL格式解决了CSR的内存合并问题，同时避免了COO所需的原子操作。然而，它引入了空间浪费，并且仍然存在负载不均衡的问题。接下来，我们将看到一种旨在解决负载均衡问题的格式。

混合ELL-COO格式 🤝

为了缓解ELL格式因少数长行导致大量空间浪费的问题，可以采用一种混合策略：ELL + COO。

工作原理

设定一个阈值（例如，大多数行拥有的非零元数量）。
使用ELL格式存储每行前 阈值 个非零元。对于非零元数量少于阈值的行，剩余位置用填充元素补足。
所有行中超出阈值的额外非零元，则使用COO格式存储。

优势

更好的空间效率：减少了不必要的填充。
完全的灵活性：现在可以向任何行添加元素。如果ELL部分有空间，则添加进去；否则，添加到COO部分。
保留ELL优点：对于大多数行（短行），仍然享有ELL带来的良好内存合并访问特性。

JDS格式详解 🧩

Jagged Diagonal Storage格式的主要设计目标是最小化SPMV中的控制流分化（线程束分化）。

格式原理

JDS格式通过两个关键步骤实现负载均衡：

按行长度排序：首先，根据每行非零元的数量（行长度）对矩阵的行进行降序排序。最长的行排在最前面。
锯齿对角线存储：将排序后矩阵的非零元按“锯齿对角线”顺序存储。即，先存储所有行的第一个非零元，然后存储所有行的第二个非零元，以此类推。这类似于没有填充的ELL列主序存储。

此外，还需要两个辅助数组：

置换数组：记录新行顺序对应的原始行索引，用于将计算结果写回正确的输出位置。
迭代指针数组：记录每个“对角线”（即所有行的第k个元素）在值数组和列索引数组中的起始位置。

SPMV并行化策略

在JDS格式上执行SPMV，同样为排序后的每一行分配一个线程。

代码：核心计算逻辑如下（伪代码）：

int sorted_row = threadIdx.x + blockIdx.x * blockDim.x; // 线程负责的排序后行索引
if (sorted_row < num_rows) {
    int original_row = jds_permutation[sorted_row]; // 找到原始行号
    float sum = 0.0f;
    int offset = sorted_row; // 在当前对角线内的偏移
    for (int diag = 0; diag < max_nonzeros_per_row; ++diag) {
        if (offset >= jds_iter_ptr[diag+1]) break; // 如果偏移超出下一个对角线的起始点，说明该行在本对角线已无元素
        int idx = jds_iter_ptr[diag] + offset; // 计算全局索引
        int col = jds_col_indices[idx];
        float val = jds_values[idx];
        sum += input_vector[col] * val;
        // 准备下一个对角线的偏移：仍然是 sorted_row，因为行顺序不变
    }
    output_vector[original_row] = sum; // 结果写入原始行位置
}

在实际高效实现中，循环条件通常直接使用每行非零元数量，而非迭代指针判断。

性能与特性分析

控制流分化：优秀。由于行按长度排序，相邻线程处理的行长度相似。在计算过程中，短行会先完成计算，线程从尾部开始依次退出，使得活跃线程保持连续，极大减少了线程束内分化。
内存合并访问：优秀。每个计算迭代中，所有活跃线程访问的数据在内存中是连续的（因为它们访问同一个“对角线”），实现了完全合并访问，且没有无效的填充数据加载。
空间效率：较好。没有填充开销，但需要额外的置换数组和迭代指针数组。
灵活性：很差。添加或删除元素会改变行的长度和顺序，需要重新排序并更新辅助数组，成本很高。
可访问性：
- 给定排序后的行索引，容易遍历其非零元。
- 给定原始行索引，需要查询置换数组或其逆数组才能找到其在排序后的位置，稍显麻烦。
- 给定一个非零元，难以直接确定其原始行号。
- 给定列号，难以找到该列的所有非零元。

过渡：JDS格式通过牺牲灵活性和增加一些簿记开销，换取了优异的负载均衡和内存访问性能。它通常是针对固定稀疏矩阵进行大量SPMV运算时的优选格式之一。

总结 📝

本节课中我们一起学习了两种重要的稀疏矩阵存储格式：

ELL格式：通过对齐行长度并按列主序存储，有效改善了内存合并访问。它易于修改，但空间效率可能较低，且存在负载不均衡问题。其性能通常优于CSR，但可能不如更专门的格式。
JDS格式：通过按行长度排序并采用锯齿对角线存储，旨在最小化线程束分化，同时保持良好的内存合并访问。它在SPMV上常能获得最佳性能，但矩阵修改成本极高。

选择哪种格式取决于具体应用场景：

如果矩阵需要频繁修改，COO或ELL可能更合适。
如果追求极致的SPMV性能且矩阵稳定，JDS或混合ELL-COO是强有力的候选。
需要在性能、空间和灵活性之间取得平衡时，CSR是一个稳健的默认选择。

理解这些格式的权衡有助于我们为特定的计算任务选择最合适的稀疏矩阵表示方法。

GPU计算：第18讲：图处理

在本节课中，我们将学习如何使用GPU进行图处理。我们将从回顾稀疏矩阵计算开始，然后探讨如何将图表示为稀疏矩阵，并介绍两种主要的图处理并行化方法：顶点中心法和边中心法。最后，我们将以广度优先搜索为例，详细分析这两种方法的具体实现和性能特点。

回顾：稀疏矩阵格式

上一讲我们介绍了稀疏矩阵计算，这是关于稀疏矩阵计算的两讲系列中的第二讲。在第一讲中，我们讨论了COO和CSR格式。上一讲我们讨论了另外两种格式：ELL格式和JDS格式。

ELL格式的做法是，我们像CSR一样对矩阵的非零元素进行分组，但不同之处在于，我们添加了填充，使得所有行（包括填充）都具有相同数量的非零元素。然后我们以列优先顺序存储这个表格。这为我们提供了一个包含非零值和列索引的数组。这种格式的优点是，当每个线程访问其负责的行的第一个元素以执行SPMV时，线程访问的是相邻的数据元素，因此内存访问是合并的。

然而，ELL格式的一个问题是填充带来的开销。如果某些行特别长，我们最终会添加大量填充，占用大量空间。为了减少填充，我们可以使用混合的ELL-COO格式。我们为大多数行保留足够多的列，但对于那些非常长的行，我们使用额外的COO格式来单独处理其多余的非零元素。这样，我们既获得了ELL格式在内存合并方面的好处，又通过单独处理超长行的非零元素来最小化填充。

之后，我们研究了JDS格式，其目的是最小化控制流分歧。与CSR和ELL类似，JDS也对行的非零元素进行分组。为了最小化控制流分歧，我们对行进行排序，并保留一个置换向量来记录每个新行对应的原始行索引。通过排序，当线程处理相邻行时，它们将处理长度相似的行，这有助于减少控制流分歧。为了实现良好的内存合并，我们以列优先顺序存储这些行。这样，当所有线程访问其对应行的第一个元素时，访问是合并的。此外，同一warp中的线程处理的行长度相似，这使我们也能在使用JDS格式时保持较低的控制流分歧。

图的表示

现在，我们来讨论图处理。首先，我们需要了解如何表示图。图由顶点和连接顶点的边组成。逻辑上，我们可以使用邻接矩阵来表示图。邻接矩阵的每一行代表一个源顶点的邻接列表，每一列代表目标顶点。矩阵中的“1”表示从行顶点到列顶点存在一条边。

邻接矩阵通常非常稀疏。我们可以使用之前讨论过的稀疏矩阵存储格式来表示图。我们将重点关注无权图，这意味着邻接矩阵中的所有值都是1，因此我们不需要存储值，只需跟踪是否存在非零元素。我们还将重点关注无向图，这意味着矩阵是对称的。但请注意，我们今天讨论的许多概念同样适用于有权图和有向图。

我们可以使用COO格式来表示图。在COO表示中，我们有一个边列表，包含每条边的源顶点和目标顶点。对于有权图，我们还需要第三个数组来存储每条边的权重。

我们也可以使用CSR或CSC格式来表示图。在CSR表示中，我们有源指针数组，指向目标顶点数组中每个源顶点邻居的起始位置。由于我们处理的是无向图，矩阵是对称的，因此CSR和CSC表示是等价的。如果图是有向的，则需要同时保留CSR和CSC表示。

选择哪种存储格式取决于图的属性和我们要使用的并行化方法，这会影响处理图的并行化策略。

图处理的并行化方法

接下来，我们讨论如何并行化图处理。主要有两种方法：顶点中心法和边中心法。

顶点中心法为图中的每个顶点分配一个线程。每个线程对其负责的顶点执行某些操作，这些线程并行运行，可能需要进行通信和全局同步。如果希望线程能够遍历某个顶点的所有邻居，通常使用CSR或CSC格式，因为它们便于查找给定顶点的所有邻居。在某些情况下，也可以使用ELL或JDS格式来进一步优化访问模式。

边中心法为图中的每条边分配一个线程。每个线程对其负责的边执行某些操作，这通常涉及查看边的源顶点和目标顶点。对于这种方法，通常使用COO格式，因为它便于查找给定边的源顶点和目标顶点。

在某些算法中，可能需要同时使用CSR和COO表示，这称为混合表示。例如，在三角形计数等子图算法中，给定一条边，需要找到其端点的邻居，这就需要同时使用COO（查找边的端点）和CSR（查找端点的邻居）。不过，我们今天要学习的广度优先搜索不需要这种混合表示。

广度优先搜索

现在，我们以广度优先搜索为例，具体应用顶点中心法和边中心法。BFS的目标是找到从某个源顶点到图中每个顶点的距离或层级。

BFS的工作原理是：从源顶点开始，标记其为第0层。然后，在每次迭代中，访问上一层级的所有顶点的邻居，并将未访问过的邻居标记为当前层级。重复此过程，直到没有新的顶点被访问。

顶点中心法：自顶向下方法

在顶点中心的自顶向下方法中，每次迭代为每个顶点启动一个线程。线程检查其负责的顶点是否属于上一层级。如果是，则遍历该顶点的所有邻居，并将未访问的邻居标记为当前层级。

这种方法之所以称为“自顶向下”，是因为如果我们将BFS视为构建一棵BFS树，那么线程就像是树中每个父顶点，负责找到其所有子顶点。

以下是使用CUDA实现自顶向下BFS内核的伪代码：

__global__ void bfs_top_down(CSRGraph graph, int* level, int* visited_flag, int current_level) {
    unsigned int vertex = blockIdx.x * blockDim.x + threadIdx.x;
    if (vertex >= graph.num_vertices) return;

    if (level[vertex] == current_level - 1) {
        for (unsigned int edge = graph.src_ptr[vertex]; edge < graph.src_ptr[vertex + 1]; ++edge) {
            unsigned int neighbor = graph.dst[edge];
            if (level[neighbor] == UINT_MAX) {
                level[neighbor] = current_level;
                *visited_flag = 1;
            }
        }
    }
}

这种方法在具有高度数顶点的图上可能会遇到严重的控制流分歧，因为不同顶点的邻居数量差异很大。

顶点中心法：自底向上方法

在顶点中心的自底向上方法中，每次迭代同样为每个顶点启动一个线程。但线程检查的是其负责的顶点是否未被访问。如果未被访问，则遍历该顶点的所有邻居，检查是否有邻居属于上一层级。如果有，则将该顶点标记为当前层级，并可以提前退出循环。

这种方法之所以称为“自底向上”，是因为线程像是BFS树中每个潜在的子顶点，通过查看其父顶点来判断自己是否属于当前层级。

以下是使用CUDA实现自底向上BFS内核的伪代码：

__global__ void bfs_bottom_up(CSRGraph graph, int* level, int* visited_flag, int current_level) {
    unsigned int vertex = blockIdx.x * blockDim.x + threadIdx.x;
    if (vertex >= graph.num_vertices) return;

    if (level[vertex] == UINT_MAX) {
        for (unsigned int edge = graph.src_ptr[vertex]; edge < graph.src_ptr[vertex + 1]; ++edge) {
            unsigned int neighbor = graph.dst[edge];
            if (level[neighbor] == current_level - 1) {
                level[vertex] = current_level;
                *visited_flag = 1;
                break;
            }
        }
    }
}

自底向上方法在高度数图上通常表现更好，因为线程一旦找到符合条件的邻居就可以提前退出，减少了冗余工作。但在初始迭代中，当已访问顶点很少时，自底向上方法会检查大量未访问顶点的所有邻居，效率较低。

方向优化BFS

方向优化BFS结合了自顶向下和自底向上方法的优点。它在初始迭代中使用自顶向下方法，当已访问顶点数量增多时，切换到自底向上方法。这种策略通常能获得最佳性能。

边中心法

在边中心法中，每次迭代为图中的每条边启动一个线程。线程检查其负责的边的源顶点是否属于上一层级，且目标顶点未被访问。如果是，则将目标顶点标记为当前层级。

以下是使用CUDA实现边中心BFS内核的伪代码：

__global__ void bfs_edge_centric(COOGraph graph, int* level, int* visited_flag, int current_level) {
    unsigned int edge = blockIdx.x * blockDim.x + threadIdx.x;
    if (edge >= graph.num_edges) return;

    unsigned int vertex = graph.src[edge];
    unsigned int neighbor = graph.dst[edge];

    if (level[vertex] == current_level - 1 && level[neighbor] == UINT_MAX) {
        level[neighbor] = current_level;
        *visited_flag = 1;
    }
}

边中心法没有循环，每个线程只处理一条边，因此控制流分歧较少。它在高度数图上通常表现良好。

数据结构的影响

最佳并行化方法的选择取决于图的结构。

高度数图：例如社交网络图，其中存在“名人”顶点，拥有大量连接。自底向上顶点中心法和边中心法在处理这类图的负载不均衡方面表现更好。
低度数图：例如道路网络图，每个顶点的连接数有限。自顶向下顶点中心法在这类图上通常更快，因为每次迭代只需要处理少量已访问顶点的邻居。

因此，没有一种BFS算法在所有数据集上都是最优的，需要根据具体的数据集特性来选择算法。

BFS与SPMV的相似性

BFS与稀疏矩阵-向量乘法在数据访问模式上非常相似。以自底向上BFS为例：

“为每个顶点分配线程”类似于SPMV中“为每一行分配线程”。
“遍历顶点的每条边”类似于SPMV中“遍历行的每个非零元素”。
“查看邻居的层级”类似于SPMV中“查看输入向量在列索引处的值”。
“更新当前顶点的层级”类似于SPMV中“更新输出向量在行索引处的值”。

事实上，许多图问题都可以用稀疏线性代数计算来表述。这样做的优点是可以利用成熟且高度优化的稀疏线性代数库。缺点是，这种表述方式并不总是解决图问题最高效的方法。

冗余工作与前沿优化

我们目前讨论的方法在每次迭代中都会检查每个顶点或每条边，以判断它们是否与当前迭代相关。这种方法的优点是易于实现、并行度高且线程间无需同步。缺点是会做大量不必要的工作，因为许多线程会发现其顶点或边不相关并立即退出。

在下一讲中，我们将探讨另一种方法：前沿优化。我们将尝试只处理与当前迭代相关的顶点或边，通过构建“前沿”来仅启动相关顶点或边的线程。这将减少冗余工作，但需要更多的线程间同步来构建这些前沿。我们将在下次课中探讨这种权衡。

本节课中，我们一起学习了图的稀疏矩阵表示、顶点中心与边中心并行化方法，并深入分析了它们在广度优先搜索中的具体应用与性能表现。请阅读教材第12章的相关章节，我们将在下节课继续讲解该章的后续内容。

GPU计算：19：图处理（第二部分）

在本节课中，我们将继续学习图处理。上一讲我们介绍了图的基本表示方法和广度优先搜索的不同并行实现。本节我们将重点探讨如何通过基于“前沿”的方法来减少冗余计算，并学习相关的性能优化技术。

上一节我们介绍了顶点中心化和边中心化的BFS实现，并观察到这些方法在每次迭代中都会检查所有顶点或边，导致大量线程空转。本节中，我们来看看如何通过只处理“前沿”中的顶点来避免这种冗余。

基于前沿的BFS

我们的目标是：在每一层，只检查属于上一层的顶点，而不是检查图中的所有顶点。具体方法是，在每一层，我们将被访问的顶点放入一个队列中。在下一层，我们只处理这个队列（即“前沿”）中的顶点。

这种方法的优点是减少了冗余计算。但缺点是需要同步：所有线程都试图向一个共享队列中添加元素，这会产生竞争。

让我们通过一个例子来理解这个过程。假设我们从源节点（第0层）开始。初始前沿队列只包含源节点。

从第0层到第1层：我们只启动一个线程来处理源节点。该线程访问其所有邻居，将未访问的邻居标记为第1层，并加入下一前沿队列。
从第1层到第2层：我们为第1层前沿队列中的每个顶点启动一个线程。每个线程访问其邻居，将未访问的邻居标记为第2层，并加入新的前沿队列。
以此类推，直到某一层的前沿队列为空，算法结束。

前沿BFS的实现与性能

以下是基于前沿的顶点中心化BFS核心逻辑的伪代码表示。每个线程处理前沿队列中的一个顶点。

// 伪代码：前沿BFS内核
for each vertex v in previous_frontier (assigned to a thread):
    for each neighbor u of v:
        // 使用原子比较交换操作，确保只有一个线程能成功标记该邻居
        old_level = atomic_compare_and_swap(level[u], UNVISITED, current_level)
        if (old_level == UNVISITED):
            // 成功标记，将邻居加入下一前沿
            index = atomic_add(¤t_frontier_counter, 1)
            current_frontier[index] = u

在这个实现中，我们使用了两个关键的原子操作来解决竞争条件：

atomic_compare_and_swap：用于安全地标记顶点为已访问。只有当顶点之前是未访问状态时，才将其设置为当前层号。这确保了同一个邻居不会被多个线程重复加入前沿队列。
atomic_add：用于安全地获取在全局前沿队列中的写入位置，避免多个线程写入同一位置。

我们对比了原始的顶点中心化自上而下BFS（检查所有顶点）和这个前沿版本在道路网络图上的性能。原始版本耗时约66毫秒，而前沿版本仅需约14.8毫秒，性能提升显著，这得益于冗余计算的减少。

优化：队列私有化

然而，前沿方法引入了一个新的瓶颈：所有线程都在竞争同一个全局计数器（current_frontier_counter）来分配队列中的位置，导致原子操作串行化，成为性能瓶颈。

以下是解决此问题的优化思路：队列私有化。其核心思想是让每个线程块维护一个私有的本地队列（可置于共享内存中以降低延迟）。线程块内的线程向本地队列添加元素，仅需在共享内存上进行原子操作，竞争大大减少。当线程块处理完所有顶点后，再一次性将其本地队列的内容合并到全局队列中。

具体步骤如下：

每个线程块在共享内存中分配一个本地队列 local_queue 和一个本地计数器 local_counter。
线程在处理顶点时，使用原子操作将邻居添加到 local_queue，并更新 local_counter。
如果 local_queue 已满，则溢出的顶点直接通过全局原子操作添加到全局队列。
处理结束后，线程块内同步。然后，由一个线程（如 threadIdx.x == 0）通过一次全局原子操作，为整个线程块的本地队列内容在全局队列中预留连续空间。
该线程块的所有线程协作，将 local_queue 中的数据高效地、合并地写入到全局队列的预留位置。

通过私有化，我们将大量细粒度的全局原子操作，转化为线程块内共享内存上的原子操作，以及最后少量粗粒度的全局原子操作和合并的全局内存写入，从而有效缓解了竞争并提升了内存访问效率。

应用此优化后，性能从14.8毫秒进一步提升到约14.67毫秒。提升幅度取决于具体硬件和图结构，在某些架构或图上可能效果更明显。

优化：减少内核启动开销

我们目前的方法是为图的每一层都启动一次GPU内核。内核启动本身有一定开销，并且我们需要在内核之间拷贝前沿计数器等数据以判断算法是否结束。

我们注意到，在图搜索的开始和结束阶段，前沿中的顶点数量通常很少，可能少到只需一个线程块就能处理。对于这些连续的小层，我们可以将它们合并到同一个内核中执行。

具体做法是：在一个内核中，一个线程块先处理完当前层的前沿，然后通过 __syncthreads() 在块内同步，接着继续处理下一层的前沿，只要下一层的前沿大小也能被一个线程块容纳。这样就避免了为这些小层单独启动内核的开销。

这种优化进一步减少了总的内核启动次数，带来了额外的性能提升。

总结

本节课中我们一起学习了图处理算法的进一步优化。

我们首先引入了基于前沿的BFS，通过只处理活跃顶点来消除冗余计算，并使用原子操作解决同步问题。
接着，我们探讨了队列私有化优化，通过将全局队列竞争转化为线程块内的本地操作，有效缓解了原子操作的瓶颈。
最后，我们学习了减少内核启动开销的方法，通过合并处理连续的小规模前沿层来提升整体效率。

这些技术展示了在图算法这类不规则、数据依赖强的计算任务中，如何通过巧妙的并行策略和同步原语来挖掘GPU的并行潜力。

GPU计算：第20讲：线程束内同步

在本节课中，我们将学习一种高级的GPU编程技术：线程束内同步。我们将探讨如何利用同一线程束内线程的特殊关系，通过数据洗牌和投票指令，实现比传统的共享内存和块内同步更高效的线程间协作。

概述

上一节我们回顾了课程至今的内容，从GPU计算基础到各种并行模式的应用。从本节开始，我们将进入课程的新阶段，探讨一些高级编程特性，这些特性可以帮助你进一步优化项目中的应用程序。今天，我们将重点讨论线程束内同步。

线程束是SM中的基本调度单位，通常包含32个线程。由于线程束内的线程遵循SIMD模型，它们几乎同时执行相同的指令，这使得它们之间的同步比跨整个线程块的同步（如 __syncthreads()）更加高效。CUDA为此提供了两类内置函数：用于在线程间直接共享寄存器数据的洗牌指令，以及用于让线程就某个条件进行集体表决的投票指令。

线程束洗牌指令

线程束洗牌指令允许同一线程束内的线程直接共享寄存器中的数据，而无需通过共享内存。这比使用共享内存和 __syncthreads() 更快，因为它避免了全局同步的开销，并且寄存器访问速度远快于共享内存。

CUDA提供了几种洗牌指令变体：

__shfl_sync：从线程束中指定的通道（即线程）复制数据。
__shfl_up_sync：从ID较低的线程复制数据。
__shfl_down_sync：从ID较高的线程复制数据。
__shfl_xor_sync：通过按位异或操作确定源线程。

应用示例：归约优化

我们曾在归约模式中看到，线程需要读取其他线程在前一次迭代中产生的值。当迭代进行到只剩下一个线程束时，所有读写操作都发生在这个线程束内部。此时，我们可以用洗牌指令替代共享内存和同步操作。

以下是优化后的归约核函数代码片段，展示了如何结合使用共享内存和洗牌指令：

// 第一部分：在共享内存中进行归约树操作，直到步长等于线程束大小
for (int stride = blockDim.x / 2; stride >= WARP_SIZE; stride >>= 1) {
    if (threadIdx.x < stride) {
        sdata[threadIdx.x] += sdata[threadIdx.x + stride];
    }
    __syncthreads();
}

// 第二部分：将数据加载到寄存器，然后使用洗牌指令在单个线程束内完成归约
if (threadIdx.x < WARP_SIZE) {
    // 加载最后两个值到寄存器
    int sum = sdata[threadIdx.x] + sdata[threadIdx.x + WARP_SIZE];
    // 使用洗牌指令继续归约
    for (int stride = WARP_SIZE / 2; stride > 0; stride >>= 1) {
        sum += __shfl_down_sync(0xffffffff, sum, stride);
    }
    // 线程0存储最终结果
    if (threadIdx.x == 0) {
        output[blockIdx.x] = sum;
    }
}

通过这种优化，我们显著提升了归约操作的性能，因为它减少了耗时的全局同步操作，并利用了更快的寄存器间通信。

线程束投票指令

线程束投票指令使得线程束内的线程能够就某个条件进行集体表决。这在需要根据线程的局部状态做出集体决策时非常有用。

主要的投票指令包括：

__all_sync：检查参与的所有线程是否都满足某个条件（谓词非零）。
__any_sync：检查是否有至少一个参与线程满足条件。
__ballot_sync：返回一个掩码，指示哪些参与线程满足条件。
__activemask：返回当前执行点上活跃线程的掩码。

应用示例：队列插入优化

考虑一个内核，其中多个线程根据条件判断是否要将数据插入一个全局队列。原始实现中，每个满足条件的线程都会原子性地递增全局队列计数器，这会造成严重的竞争。

我们可以优化为：让一个线程束内的线程协作，只进行一次原子操作来为该线程束所有需要插入的线程预留空间。

以下是优化步骤和关键代码：

确定活跃线程和领导者：使用 __activemask() 找出当前条件分支中活跃的线程，并选择第一个活跃线程作为领导者。
```
unsigned int active = __activemask();
unsigned int leader = __ffs(active) - 1; // __ffs 返回从1开始的索引
```

计算需要插入的元素数量：使用 __popc 计算活跃掩码中置位比特的数量。
```
unsigned int num_active = __popc(active);
```

领导者线程执行原子操作：只有领导者线程执行一次原子加操作，为整个线程束预留空间。

unsigned int queue_index;
if ((threadIdx.x % WARP_SIZE) == leader) {
    queue_index = atomicAdd(&queue_size, num_active);
}

广播起始索引：领导者线程将其获得的队列起始索引广播给线程束内的所有其他活跃线程。
```
queue_index = __shfl_sync(active, queue_index, leader);
```

计算每个线程的偏移并存储：每个线程计算在已预留空间内的个人偏移量（即它前面有多少个活跃线程），然后进行存储。

unsigned int thread_mask = 1 << (threadIdx.x % WARP_SIZE);
unsigned int prev_mask = thread_mask - 1;
unsigned int prev_active = active & prev_mask;
unsigned int offset = __popc(prev_active);

queue[queue_index + offset] = my_value;

在这个特定例子中，由于操作本身很简单，额外的协作开销可能抵消了减少原子操作带来的收益。然而，在更复杂的、包含更多计算上下文的场景中，这种线程束级别的协作可以显著减少全局原子操作的竞争，从而提升性能。

总结

本节课我们一起学习了GPU编程中的线程束内同步技术。我们首先介绍了线程束洗牌指令，它允许线程直接共享寄存器值，并以归约优化为例展示了其应用。接着，我们探讨了线程束投票指令，它支持线程进行集体表决，并以优化队列插入为例说明了其使用方法和协作流程。

这些高级特性允许程序员更精细地控制线程间的协作，特别是在线程束内部，从而有可能实现更高的性能。理解并恰当应用这些指令，是进行高性能GPU编程的重要技能。

GPU计算：21：固定内存与流

在本节课中，我们将要学习两个重要的CUDA概念：固定内存（Pinned Memory）和流（Streams）。我们将探讨如何通过使用固定内存来加速主机与设备之间的数据传输，以及如何利用流来实现数据传输与内核执行的并行化，从而更高效地利用GPU硬件资源。

回顾：线程束内同步

上一节我们介绍了线程束（Warp）内同步。我们提到，可以利用同一线程束内线程之间的特殊关系进行快速同步，CUDA为此提供了内置函数。线程束同步原语主要有两大类：

以下是两类主要的线程束同步原语：

洗牌指令：允许同一线程束内的线程共享数据，它们可以直接共享彼此的寄存器值，而无需通过共享内存。
投票指令：线程可以对某个条件进行“投票”，这对于检查哪些线程的谓词评估为真，或者检查线程束中哪些线程在某个分支上处于活动状态非常有用。

我们首先研究了洗牌指令。洗牌指令有多种类型，例如直接从一个线程获取数据，或者向上/向下洗牌。还有一种“异或洗牌”在某些情况下也很有用。

我们将洗牌指令应用于归约操作。在之前的实现中，当归约进行到一个线程束时，我们仍然通过共享内存并使用__syncthreads()进行同步。但我们意识到，实际上可以仅使用洗牌指令在线程束内完成归约树，而无需经过共享内存，也无需调用__syncthreads()。

接着，我们研究了线程束投票指令。投票指令有不同变体，例如检查所有线程是否满足某个谓词，或者是否有任一线程满足谓词。ballot指令可以检查哪些线程的谓词评估为真，哪些为假，并返回一个掩码。activemask指令则告诉我们哪些线程在特定分支上是活动的。

线程束投票函数的一个应用场景是，当需要线程束内的线程进行通信，并让一个线程代表其他线程执行操作时。例如，如果我们要向队列中添加元素并递增计数器，与其让每个线程都去原子递增计数器，不如让线程束中的一个“领导者”线程代表所有线程执行原子递增操作。然后，每个线程可以根据领导者线程分配的结果，计算出自己在队列中应占用的位置。

我们看到了实现此过程的步骤：首先分配一个领导者线程，然后计算需要向队列中添加多少元素，接着领导者线程执行原子操作，最后将结果广播给其他线程，每个活动线程再据此确定存储自己结果的位置。

我们使用了activemask来找出哪些线程是活动的，并使用__ffs（查找首个置位位）内部函数在活动掩码中找到第一个置位位，从而确定领导者线程。我们通过对活动掩码进行__popc（人口计数）来计算需要分配多少队列空间。__ffs和__popc本身并非线程束同步原语，它们是每个线程可以独立调用的内部函数，但在使用线程束内部函数时非常有用，因为我们可以将它们应用于活动掩码来完成各种有趣的操作。

在确定了要添加到队列的元素数量后，我们让领导者线程执行原子操作，然后使用__shfl_sync让线程束中的每个线程从领导者线程获取值，以确定其应存储结果的位置。最后，我们使用一些位操作技巧，让每个活动线程计算出自己在队列中的具体位置。

这就是我们上一节讨论的内容。在开始今天的新内容之前，关于线程束内部函数还有什么问题吗？

固定内存

今天我们将转向一个新主题：固定内存。我们还将讨论流。这是两个相关的主题，我们将从固定内存开始。

首先，什么是固定内存？回想一下我们之前的向量加法代码。我们分配了向量内存，复制源向量，调用加法内核，然后复制回结果向量，最后释放内存。

当我们编译并运行这段代码时，可以看到执行时间的分配情况。CPU执行耗时76毫秒，而GPU总耗时169毫秒。但请注意，GPU内核执行时间实际上只有3.2毫秒。GPU总时间主要由什么主导？是复制时间。将两个向量复制到GPU花了58毫秒，将结果复制回来花了94毫秒。

显然，当向量已经在GPU上时，在GPU上执行加法要快得多。但如果向量在CPU上，从GPU获得的加速可能无法抵消复制数据的开销。我们今天要问的问题是：能否优化这种复制？能否减少每次向GPU复制数据以及从GPU复制数据回来的时间？

直接内存访问

为了讨论如何优化复制，我们首先需要了解调用cudaMemcpy时发生了什么。这个过程涉及直接内存访问。

直接内存访问允许硬件单元在不涉及CPU的情况下访问内存。当我们在CPU和GPU之间复制数据时，实际上使用了DMA。这意味着，当数据从主内存复制到GPU时，CPU并不参与实际的复制过程。CPU不会从主内存加载数据然后存储到GPU内存。实际情况是，存在一个DMA引擎。当CPU调用cudaMemcpy时，它只是发起调用，然后可以继续执行其他任务，而DMA引擎负责从CPU的主内存读取数据并将其放入GPU内存。

使用DMA相对于让CPU参与复制有什么优势？DMA是专用硬件，可能复制得更快。此外，CPU可以在内存复制进行的同时去处理其他有用的工作，从而更好地利用资源。

那么使用DMA有什么缺点呢？当CPU访问主内存时，它使用虚拟地址并进行地址转换。而DMA引擎没有进行地址转换的能力，它使用物理地址来读取主内存并写入GPU内存。

DMA引擎使用物理地址的缺点是什么？如果DMA正在读取的内存被操作系统换出，DMA无法检测到这一点。CPU可以检测到，因为如果某个页面被换出，在进行地址转换时会发生页面错误，然后操作系统会将其换入。但对于DMA，当它从某个物理页面读取时，如果该物理页面被操作系统交换出去，DMA引擎无法知晓，它会继续读取同一个物理页面，而此时该页面可能已经包含了新的数据。这显然是我们不希望发生的。

因此，为了避免数据损坏，操作系统必须确保任何被DMA访问的页面都不会被交换出去。我们需要锁定这些页面，防止操作系统交换它们。换句话说，这些页面需要被“固定”。这就是固定内存。任何被DMA访问的页面都需要被页面锁定，或者说，需要被固定。我们必须将其分配为固定内存，这样它就不会被换出。

CUDA内存复制的工作原理

现在我们知道cudaMemcpy使用DMA，而DMA要求页面被固定。那么cudaMemcpy是如何工作的呢？当我们使用malloc分配数据时，并没有做任何特殊处理来锁定它，因此我们假设这些数据可以被换出。

实际上，cudaMemcpy的工作方式如下：当从主机复制到设备时，CUDA运行时会使用一个预先分配的固定内存缓冲区。CPU首先将我们想要复制的数据从普通缓冲区复制到这个固定内存缓冲区。然后，DMA引擎从固定内存缓冲区将数据复制到GPU。当从设备复制回主机时，过程类似：DMA从设备内存复制数据到一个固定内存缓冲区，然后CPU再从这个固定内存缓冲区复制数据到由malloc分配的CPU缓冲区。

由此可见，每次调用cudaMemcpy时，我们并不是只执行一次复制，而是执行了两次复制：一次在普通CPU缓冲区和固定内存缓冲区之间，另一次在固定内存缓冲区和设备内存之间。这显然带来了额外的开销。

优化复制：使用固定内存

那么，如何优化复制，使其更快呢？如果我们不想进行两次复制，而只想进行一次复制，我们可以对我们的缓冲区做什么？我们可以直接固定它们。我们可以在CPU上拥有一个包含原始数据的固定内存缓冲区，然后直接从该缓冲区进行DMA复制，而不是先从一个非固定缓冲区复制到固定内存缓冲区，再由DMA从固定内存缓冲区复制。

CUDA提供了相应的API来直接分配和释放固定内存中的主机数据：

cudaMallocHost：类似于cudaMalloc，但在固定内存中分配主机数据。
cudaFreeHost：释放之前用cudaMallocHost分配的固定内存数据。

既然cudaMallocHost能让复制快这么多，我们是否应该将所有数据都分配在固定内存中？这是一个好主意吗？当我们固定一个页面时，我们是在告诉操作系统不能换出这个页面。如果我们分配了过多的固定内存，就会限制操作系统在交换页面时的灵活性，从而可能降低整体系统性能。因此，固定内存应谨慎使用，只应用于那些需要频繁在CPU和GPU之间传输的数据。

在现代拥有大量内存的系统中，如果你的GPU显存较小（例如4GB），而主机内存很大（例如128GB），那么使用较多的固定内存可能问题不大，因为你需要固定的内存量相对于主机总内存来说并不大。但总的来说，需要小心，不要将所有东西都分配在固定内存中，只应在复制开销确实很重要的地方使用。

性能对比

让我们修改代码，使用cudaMallocHost和cudaFreeHost来分配固定内存，看看性能提升。

修改后运行代码，可以看到GPU时间从169毫秒大幅改善到80毫秒。原因是复制时间减少了：从58毫秒降至43毫秒，从94毫秒降至20毫秒。有趣的是，CPU时间也改善了。这是因为现在CPU页面也被固定了，这意味着在CPU上进行计算时，不会发生页面错误，所有页面都驻留在内存中。

关于固定内存还有什么问题吗？

流

接下来，我们将讨论流。在此之前，我想谈谈系统架构层面的并行性。我们已经讨论了很多关于数据并行性的内容，即让不同线程处理不同数据。今天我想讨论的是任务并行性，更具体地说，是一种流水线并行。

在一个典型的系统中，我们有CPU和GPU。CPU有主内存，GPU有设备内存。通常，我们先将数据复制到GPU，然后在GPU上计算，最后将结果复制回来。但实际上，系统能够同时执行以下操作：在GPU上执行内核、从主机向设备复制数据、从设备向主机复制数据。执行GPU计算的硬件、从主内存复制到GPU内存的硬件、以及从GPU内存复制回主内存的硬件，这是三种不同的硬件资源。

我们可以同时使用这三种资源。但到目前为止，我们的做法是顺序执行：复制到GPU、在GPU上运行、复制回来。每次只使用其中一种资源，没有充分利用所有三种资源。我们的目标是尝试同时利用它们，以尽可能充分地利用硬件。

以向量加法为例，我们之前是：复制A和B到设备，执行网格计算A+B并存到C，然后复制C回主机。在水平时间轴上，当复制A和B到设备时，运行GPU的部分和从设备复制回主机的资源是空闲的。当在GPU上执行时，负责在主机和设备之间复制的硬件是空闲的。当从GPU复制回CPU时，GPU本身和负责从主机复制到设备的硬件是空闲的。

我们如何能更好地利用这些硬件资源呢？我们可以采用流水线技术。与其复制所有的A和B，然后计算，再复制所有的C回来，不如将A、B和C分成若干段。我们可以复制第一段A和B到设备，当第一段在设备上就绪后，执行网格计算A1+B1，同时，我们可以开始复制第二段A2和B2到设备。然后，当计算完成且第二段复制也完成时，我们可以将A1+B1的结果从设备复制回C1，同时，可以执行网格计算A2+B2，并开始复制第三段A3和B3到设备，如此继续。

这样，我们将数组分成段，每一步都在并行地进行：复制一个段、计算另一个段、复制回第三个段。我们正在重叠不同段的复制和计算。

流的引入

为了实现这种流水线，我们需要能够进行异步内存复制。我们需要能够启动一个复制操作，然后无需等待它完成，就可以继续执行其他操作，比如启动内核或另一个复制操作。这就引入了流的概念。

在CUDA中，流是操作序列。放入同一个流中的任务按顺序执行。而放入不同流中的任务可以并行执行。默认情况下，如果不指定流，所有操作都会进入默认流并被序列化。

如果我们希望不同的内存复制和内核执行能够并行，就需要将它们放入不同的流中。此外，我们还需要使用异步内存复制，以便主机在发起复制后可以继续执行，而不必等待复制完成。

使用流的API

以下是相关的API调用：

cudaStreamCreate：创建一个流，并提供一个指向cudaStream_t类型对象的指针来存储流句柄。
cudaMemcpyAsync：异步内存复制API，与cudaMemcpy类似，但多了一个参数用于指定流。
内核启动配置：在内核启动配置中，可以在网格维度和块维度之后，可选地指定动态共享内存大小（默认为0），以及第四个可选参数——流。

代码实现：使用流的向量加法

我们希望将输入分成若干段，然后为每个段在不同的流中执行异步复制到设备、内核启动和复制回主机操作。

首先，我们创建多个流（例如32个）。然后，计算段的大小，并循环遍历每个段。对于每个段，我们确定其起始和结束位置。接着，我们使用cudaMemcpyAsync将数据复制到设备，使用修改后的内核启动（指定流和段参数）执行计算，再使用cudaMemcpyAsync将结果复制回主机。所有这些操作都指定了同一个流，因此对于同一个段，这些操作是顺序的。但不同段之间的操作，由于在不同的流中，可以并行执行。

在启动所有异步操作后，我们调用cudaDeviceSynchronize等待所有流中的操作完成。由于操作是异步和并行的，我们无法再单独测量每次复制或内核的时间，只能测量整个流水线的总时间。

性能分析与可视化

运行修改后的代码，GPU总时间从80毫秒进一步降低到60毫秒。整个复制加内核的时间为49毫秒，这非常接近单独的复制到设备时间（43毫秒），表明我们成功地将内核执行和复制回的时间隐藏在了复制到设备的时间之后。

为了更直观地理解，我们可以使用性能分析器（如nvprof）来生成时间线。在不使用流的版本中，时间线显示所有操作（复制到设备、内核执行、复制回主机）都在默认流中顺序执行。而在使用流的版本中，时间线显示这些操作被分割成许多小块，并且复制到设备、内核执行和复制回主机的操作在不同的流中相互重叠，形成了流水线。这清楚地展示了我们如何通过重叠操作来隐藏延迟。

总结

本节课中，我们一起学习了两个关键的CUDA优化技术。

首先，我们探讨了固定内存。通过使用cudaMallocHost直接在固定内存中分配主机数据，可以避免cudaMemcpy内部的额外复制步骤，从而显著减少主机与设备之间的数据传输时间。但需注意，过度使用固定内存可能会限制操作系统的虚拟内存管理，应谨慎用于频繁传输的数据。

其次，我们学习了流。通过创建多个流并将计算任务分段，我们可以利用cudaMemcpyAsync和指定流的内核启动，实现数据传输与内核执行的流水线并行。这允许我们同时利用系统的复制和计算硬件资源，将内核执行和一部分数据传输时间隐藏在主要的数据传输时间之后，从而进一步提升整体性能。

通过结合使用固定内存和流，我们可以更高效地管理GPU计算任务，尤其是在数据量较大、需要在主机和设备之间频繁移动数据的应用中。

GPU计算：第22讲：动态并行

在本节课中，我们将要学习CUDA编程中的动态并行技术。动态并行允许GPU上的线程启动新的内核网格，这对于处理具有嵌套并行性的应用程序非常有用。

上一节我们介绍了固定内存和流的概念，本节中我们来看看动态并行。

什么是动态并行？

动态并行指的是在GPU上执行的线程能够启动新的网格在GPU上运行的能力。这意味着，如果一个线程在执行过程中发现更多可以并行处理的工作，它可以动态地启动一个新的内核来处理这些工作。

为何使用动态并行？

使用动态并行主要有两个原因：

避免与CPU同步：如果计算由多个内核组成，我们可以直接从GPU启动它们，无需CPU介入。
处理嵌套并行性：许多应用程序具有多层次的并行性。例如，一个并行任务单元内部可能包含更多可并行的工作。当内部并行工作的数量在编译时未知时，动态并行尤其有用。

以下是两种典型的嵌套并行应用场景：

不规则嵌套工作：每个线程内部的并行工作量不同。例如，在图算法中，每个顶点（线程）的邻居数量不同，因此每个线程需要启动不同数量的子线程来并行访问其邻居。
递归深度未知：每个线程可能根据运行时信息决定是否进一步递归（即启动更多线程）。例如，在四叉树/八叉树空间划分或快速排序等分治算法中，是否进一步划分取决于当前分区数据的属性。

动态并行示例：优化广度优先搜索

我们将以广度优先搜索为例，展示如何使用动态并行进行优化。在传统的BFS实现中，每个线程负责处理一个顶点，并顺序遍历其所有邻居。

// 传统BFS内核代码片段（顺序处理邻居）
__global__ void bfs_kernel(...) {
    int tid = ...; // 线程ID对应一个顶点
    int start = csr_graph.offsets[tid];
    int num_neighbors = csr_graph.offsets[tid+1] - start;

    for (int i = 0; i < num_neighbors; i++) { // 顺序循环
        int neighbor = csr_graph.edges[start + i];
        // 访问邻居并可能加入下一层前沿
        if (atomicCAS(&level[neighbor], -1, current_level) == -1) {
            int pos = atomicAdd(&current_frontier_size, 1);
            next_frontier[pos] = neighbor;
        }
    }
}

我们可以使用动态并行来并行化每个线程内部的邻居遍历循环。父线程根据其邻居数量启动一个子网格，每个子线程处理一个邻居。

// 使用动态并行的BFS内核代码片段
__global__ void bfs_parent_kernel(...) {
    int tid = ...;
    int start = csr_graph.offsets[tid];
    int num_neighbors = csr_graph.offsets[tid+1] - start;

    if (num_neighbors > THRESHOLD) { // 仅当邻居数较多时启动子网格
        dim3 threads_per_block(1024);
        dim3 num_blocks((num_neighbors + threads_per_block.x - 1) / threads_per_block.x);
        bfs_child_kernel<<<num_blocks, threads_per_block>>>(..., start, num_neighbors, ...);
    } else {
        // 邻居数较少，仍顺序处理
        for (int i = 0; i < num_neighbors; i++) {
            // ... 处理邻居的代码 ...
        }
    }
}

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-cmps297s-gpu-comp/img/d4537e5532734bbdad21de8bd7c07dd7_17.png)

// 子内核，每个线程处理一个邻居
__global__ void bfs_child_kernel(..., int start, int num_neighbors, ...) {
    int i = threadIdx.x + blockIdx.x * blockDim.x; // 子线程索引
    if (i < num_neighbors) {
        int neighbor = csr_graph.edges[start + i];
        // 访问邻居并可能加入下一层前沿
        if (atomicCAS(&level[neighbor], -1, current_level) == -1) {
            int pos = atomicAdd(&current_frontier_size, 1);
            next_frontier[pos] = neighbor;
        }
    }
}

性能考量与优化

直接让每个线程都启动子网格可能导致性能下降，原因包括：

启动队列溢出：GPU有未决启动数量的限制（默认2048）。超过此限制可能导致错误或程序挂起。可以使用 cudaDeviceSetLimit(cudaLimitDevRuntimePendingLaunchCount, new_limit) 提高限制。
流竞争：默认情况下，同一线程块内的线程共享一个默认流，它们发起的网格会串行执行。可以通过编译器标志 --default-stream per-thread 为每个线程创建独立的流，以增加并行性。
小网格开销：启动一个只有少量线程的网格，其开销可能超过并行带来的收益，并且会降低SM的利用率。

针对以上问题，我们可以进行优化：

设置启动阈值：仅为拥有大量邻居（例如超过1200个）的顶点启动动态并行。对于邻居数少的顶点，仍在父线程中顺序处理。这能显著减少启动数量并避免小网格。
聚合启动（高级优化）：让一个线程（如一个Warp或一个块的主线程）收集多个线程的工作，然后代表它们发起一个更大的、合并后的网格。这能进一步减少启动次数和调度开销。

其他应用：卸载驱动代码

动态并行另一个用途是将控制循环从CPU卸载到GPU。例如，在BFS中，原本由CPU循环调用每一层的内核。现在，可以启动一个单独的GPU线程来执行这个循环，逐层启动内核。

// 在GPU上驱动BFS各层计算的内核
__global__ void bfs_driver_kernel(...) {
    *current_frontier_size = 0; // GPU上直接初始化
    while (*previous_frontier_size > 0) {
        // 动态启动处理当前层的内核
        bfs_parent_kernel<<<...>>>(...);
        cudaDeviceSynchronize(); // 等待当前层内核完成
        // 准备下一层迭代
        *previous_frontier_size = *current_frontier_size;
        *current_frontier_size = 0;
        swap(previous_frontier, next_frontier);
        current_level++;
    }
}

这样做的好处是释放了CPU，使其可以同时执行其他任务。但需要注意的是，单个控制线程在CPU上可能运行得更快，因此除非需要CPU同时处理其他工作，否则此优化可能不会提升整体性能。

重要注意事项

内存可见性：
- 父线程在启动子网格前对全局内存的写入，对子网格中的线程是可见的。
- 子网格对全局内存的写入，只有在父线程调用 cudaDeviceSynchronize() 等待子网格完成后，才对父线程可见。
共享内存和本地内存：不要将指向共享内存或线程本地内存的指针传递给动态启动的子内核。因为子网格可能在不同的SM上执行，无法访问父线程所在块的共享内存或父线程的本地内存。
嵌套深度限制：硬件对动态并行的嵌套深度有限制（当前硬件通常为24层）。即一个网格启动另一个网格，后者再启动一个网格……这样的链不能超过24层。

本节课中我们一起学习了CUDA动态并行的概念、应用场景、实现方法以及关键的优化技术和注意事项。动态并行是处理不规则和递归并行模式的强大工具，但需要谨慎使用以避免性能陷阱。

GPU计算：第23讲：综合专题与课程总结 🎓

在本节课中，我们将快速浏览课程中未深入探讨的剩余主题，包括多GPU编程、互连技术、统一虚拟寻址、零拷贝内存、事件、张量核心、常用库、其他编程接口、其他GPU硬件，以及关于CPU与GPU性能对比的说明。本节课将作为课程的总结。

多GPU编程

上一节我们介绍了单GPU编程模型。本节中我们来看看当系统拥有多个GPU时，如何协调它们的工作。

多GPU编程主要有两种场景：

单节点多GPU：一个CPU连接多个GPU。
多节点多GPU：多个CPU（可能各自连接多个GPU）通过网络互联。

对于单节点多GPU，可以使用CUDA API进行管理。

cudaGetDeviceCount()：获取节点上的GPU数量。
cudaSetDevice(int device)：选择要操作的特定GPU。

通常，我们会结合使用多CPU线程（如OpenMP）来管理多GPU，每个CPU线程负责操作一个GPU。

对于多节点多GPU，通常使用MPI（消息传递接口）库来协调跨网络的不同CPU进程。每个MPI进程（或称为rank）可以控制一个或多个GPU。

跨GPU（尤其是跨网络节点）的数据共享是一个关键挑战。数据需要在GPU内存、CPU内存和网络之间移动，这可能导致高昂的通信开销。一个常见的优化策略是重叠计算与通信。例如，在进行模板计算时，可以在计算内部区域的同时，通过网络交换边界（halo）数据。

互连技术

GPU与CPU之间通过互连技术进行通信。了解不同的互连方式对性能有重要影响。

以下是两种主要的互连技术：

PCIe：一种通用、广泛支持的互连标准，用于连接CPU、GPU、FPGA等多种设备。
NVLink：NVIDIA推出的专用高速互连技术，带宽高于PCIe，但通常只支持特定代的NVIDIA GPU和部分CPU。

选择互连技术时，需要在通用性（PCIe）和性能（NVLink）之间做出权衡。

统一虚拟寻址

之前我们假设CPU内存和GPU内存拥有各自独立的虚拟地址空间。现代GPU系统支持统一虚拟寻址。

在UVA系统中，CPU和GPU共享一个统一的虚拟地址空间，但物理内存仍然是分开的。关键特性是：CPU和GPU的地址范围是互斥的。通过查看指针的值，运行时系统就能判断它指向的是主机内存还是设备内存。

这带来的主要好处是简化了编程。例如，在执行内存拷贝时，不再需要显式指定方向（如cudaMemcpyHostToDevice），可以使用cudaMemcpyDefault，运行时系统会根据指针值自动判断。

// 非UVA时需要指定方向
cudaMemcpy(dev_ptr, host_ptr, size, cudaMemcpyHostToDevice);
// UVA时可以使用Default
cudaMemcpy(dev_ptr, host_ptr, size, cudaMemcpyDefault);

零拷贝内存

传统上，GPU线程无法直接访问CPU内存，需要显式拷贝。零拷贝内存允许GPU线程直接访问固定的主机内存。

当GPU线程访问零拷贝内存时，数据会在需要时按需从CPU内存拷贝到GPU内存，而不是由程序员预先进行批量拷贝。

要使用零拷贝内存，必须使用cudaHostAlloc()分配固定内存。在不支持UVA的系统上，还需要使用cudaHostGetDevicePointer()获取一个GPU可用的指针。在支持UVA的系统上，可以直接将主机指针传递给GPU内核。

零拷贝内存的优势包括：

简化编程：无需管理显式的数据拷贝。
潜在的性能收益：如果只访问大数据集的一小部分，可以避免不必要的批量拷贝；同时，按需拷贝可以与计算重叠，实现某种自动的流水线操作。

然而，其性能并不总是优于精心手动管理的批量拷贝，因为按需拷贝可能无法充分分摊每次数据传输的固定开销。是否使用需根据具体应用场景（如数据访问模式、计算通信比）进行测试。

事件

为了测量性能，我们之前常在每次内核启动或内存拷贝后使用cudaDeviceSynchronize()进行同步。但这会阻碍异步执行和并发。CUDA事件提供了更精细的计时和同步机制。

事件可以插入到CUDA流中，用于记录特定时间点，而无需阻塞整个设备。

以下是使用事件的核心操作：

cudaEventCreate(): 创建事件对象。
cudaEventRecord(cudaEvent_t event, cudaStream_t stream = 0): 在指定流中记录事件。
cudaEventSynchronize(cudaEvent_t event): 主机等待特定事件完成（而不是等待流中所有操作）。
cudaEventElapsedTime(float* ms, cudaEvent_t start, cudaEvent_t end): 计算两个事件间的时间间隔。
cudaEventDestroy(): 销毁事件。

使用事件可以在不干扰异步操作的情况下，精确测量内核或内存拷贝的执行时间。

张量核心

从Volta架构（如V100）开始，NVIDIA GPU引入了张量核心。它们是专门用于加速矩阵乘累加运算的硬件单元。

每个张量核心能在单个指令中完成一个小型矩阵（如4x4）的 D = A * B + D 运算。这在深度学习等以矩阵运算为核心的工作负载中能带来巨大的性能提升。使用张量核心通常需要通过特定的API或库（如cuBLAS）进行调用。

CUDA库

在课程中，我们手动实现了许多并行原语（如归约、扫描、矩阵乘法）以理解其原理。但在实际开发中，应优先使用高度优化的CUDA库。

以下是一些重要的CUDA库：

Thrust / CUB: 提供并行原语，如归约、扫描、排序、筛选等。
cuBLAS: 提供高度优化的密集线性代数运算，如向量加法、矩阵-向量乘法、矩阵-矩阵乘法等。
cuSPARSE: 提供稀疏线性代数运算，支持多种稀疏矩阵存储格式。
cuDNN: 深度神经网络库，被TensorFlow、PyTorch等主流框架用于GPU加速。
其他: 还有用于图形处理的nvGRAPH，用于快速傅里叶变换的cuFFT，用于信号/图像处理的NPP等。

使用这些库可以大幅提高开发效率和应用程序性能。

其他编程接口

CUDA并非编程GPU的唯一方式。根据可移植性和易用性需求，还有其他选择。

以下是几种其他编程接口：

OpenCL: 开放标准，支持跨厂商GPU（如AMD、Intel）、CPU甚至FPGA，可移植性更强，但代码可能比CUDA更复杂。
OpenACC: 指令制导编程模型。通过向C/C++/Fortran代码中添加编译指令（类似OpenMP），编译器会自动将代码并行化并移植到GPU上，简化了编程。
C++ AMP: 一种C++库和语言扩展，允许使用C++语法（如lambda表达式）编写GPU代码。

其他GPU硬件

本课程主要关注NVIDIA的独立GPU。但GPU市场还有其他参与者。

主要的GPU硬件包括：

NVIDIA GPU: 课程重点，广泛用于高性能计算和深度学习。
AMD GPU: 如Radeon系列，同样用于游戏和高性能计算。
Intel GPU: 常见于集成显卡，近年来也推出了独立GPU产品。
Arm Mali GPU: 常见于移动和嵌入式设备。
集成GPU: CPU和GPU集成在同一芯片上，共享物理内存，消除了PCIe数据传输开销。

关于CPU-GPU对比的说明

在课程中，为了突出GPU的并行能力，我们通常将优化的GPU代码与单线程、未向量化的CPU代码进行对比。这夸大了GPU的加速比。

为了进行公平的性能比较，应该将GPU代码与充分并行化且向量化的CPU代码进行对比。现代CPU同样拥有多核、多线程以及SIMD向量单元。本课程未深入探讨CPU并行化，是为了专注于GPU编程原理，但大家在评估性能时应意识到这一点。

总结

本节课中我们一起快速回顾了GPU计算课程的剩余主题。我们简要介绍了多GPU编程的挑战、互连技术、统一虚拟寻址和零拷贝内存如何简化编程、使用事件进行精细计时、专用于矩阵运算的张量核心、一系列可提高开发效率的CUDA库、除CUDA外的其他编程接口、以及不同的GPU硬件生态。

希望本课程为你打下了坚实的GPU并行编程基础。理解这些核心概念和优化原则，将有助于你开发高效的应用或进一步探索更高级的主题。感谢大家参与本课程！

GPU计算：P24：矩阵乘法的高级优化 🚀

在本节课中，我们将学习如何对GPU上的矩阵乘法运算进行高级优化。我们将从分析基础实现的性能瓶颈开始，逐步探讨如何通过共享内存分块、寄存器分块、向量化访存、软件流水线等一系列技术，将矩阵乘法从内存带宽受限的计算转变为计算资源受限的计算，从而最大化GPU的利用率。

算术强度分析

上一节我们介绍了矩阵乘法的基本并行方法。本节中，我们来看看如何分析其性能潜力。

矩阵乘法的核心操作是计算输出矩阵 C 的每个元素，它是输入矩阵 A 的一行和矩阵 B 的一列的点积。一个简单的并行策略是为每个输出元素分配一个线程。

基础实现的算术强度分析：
每个线程需要从全局内存加载两个单精度浮点数（共8字节），执行一次乘法和一次加法（共2次浮点运算）。因此，算术强度为：
算术强度 = 浮点运算次数 / 访存字节数 = 2 / 8 = 0.25 次浮点运算/字节
这个值远低于现代GPU（如H100）达到计算资源饱和所需的约20次浮点运算/字节，因此该实现是内存带宽受限的。

然而，矩阵乘法本质上具有成为计算密集型运算的潜力。对于一个完美的实现（每个输入值仅加载一次），计算一个 N x N 矩阵乘法的算术强度为：
完美算术强度 = (2 * N^3) / (12 * N^2) = N/6 次浮点运算/字节
当 N=1024 时，强度可达约171次浮点运算/字节，远高于计算饱和阈值。因此，我们的优化目标就是尽可能接近这个理想强度。

共享内存分块

为了提升算术强度，我们首先引入共享内存分块技术。

其核心思想是：将一个线程块负责计算的输出矩阵分块，并将计算该输出块所需的输入矩阵也相应分块。线程块先将一对输入块从全局内存加载到共享内存中，然后所有线程再利用共享内存中的数据计算各自负责的输出元素。

共享内存分块的算术强度分析：
假设我们使用 BM x BN 大小的输出块，以及 BM x BK 和 BK x BN 大小的输入块。

访存量：加载两个输入块，共 2 * BM * BK * 4 字节。
计算量：输出块包含 BM * BN 个点积，每个点积涉及 BK 次乘加（2次浮点运算），总计算量为 2 * BM * BN * BK 次浮点运算。
算术强度：(2 * BM * BN * BK) / (8 * BM * BK) = BN / 4 次浮点运算/字节。

分析表明，算术强度与输出块的宽度 BN 成正比，而与输入块的宽度 BK 无关。这意味着：

增大输出块尺寸（BM 和 BN）能直接提升算术强度。
我们可以保持较小的 BK 以节省共享内存，从而将共享内存资源用于增大输出块。

典型的配置如：BM = BN = 128, BK = 8。此时算术强度为 128 / 4 = 32 次浮点运算/字节，足以使计算在高端GPU上达到计算资源饱和。

然而，增大输出块意味着每个线程需要负责更多输出元素，这些累加器需要存储在寄存器中。同时，输入块存储在共享内存中。GPU上这两种资源都是有限的，因此我们需要在资源约束下进行优化。

内核实现框架

在深入优化细节之前，我们先看一下高级优化内核的代码框架。

以下是实现上述分块策略的伪代码框架：

// 每个线程块负责一个 BM x BN 的输出块
// 每个线程负责一个 TM x TN 的输出子块（存储在寄存器中）
// 循环遍历所有输入块对 (BK x BK)
for (int k = 0; k < K; k += BK) {
    // 1. 协作加载：将全局内存中的 A 块和 B 块加载到共享内存中
    load_tile_from_global_to_shared(A_global, A_shared, ...);
    load_tile_from_global_to_shared(B_global, B_shared, ...);
    __syncthreads();

    // 2. 计算贡献：每个线程从共享内存读取数据，计算对其 TM x TN 输出子块的贡献
    //    结果累加到线程的寄存器数组中
    matrix_multiply_accumulate(A_shared, B_shared, reg_C, ...);
    __syncthreads();
}
// 3. 写回结果：每个线程将其寄存器中的 TM x TN 结果写回全局内存
write_tile_from_registers_to_global(reg_C, C_global, ...);

接下来，我们将分解这个框架中的关键部分。

关键组件实现

以下是实现上述框架所需的关键辅助函数及其要点。

声明与初始化寄存器中的输出块

我们希望每个线程的 TM x TN 输出块完全存储在寄存器中，以获得最快的访问速度。在CUDA中，编译器可以将小的、索引恒定的局部数组提升到寄存器。

实现方法：

// 在函数内声明一个固定大小的数组
float reg_C[TM][TN];
// 通过完全展开循环来初始化和访问它，确保索引是编译期常量
#pragma unroll
for (int i = 0; i < TM; ++i) {
    #pragma unroll
    for (int j = 0; j < TN; ++j) {
        reg_C[i][j] = 0.0f;
    }
}

循环展开后，所有数组访问都变成了对固定寄存器的直接访问，避免了局部内存的开销。

从全局内存加载数据块到共享内存

加载函数需要将一个大块（例如128x8）的数据从全局内存高效、合并地转移到共享内存。由于线程块中的线程数（如256）可能少于块中的元素数（如1024），每个线程需要加载多个元素。

初始实现思路：
将输入块划分为多个子块，每个线程依次加载每个子块中的一个元素。这能保证合并访问，但每个线程需要发出多条加载指令。

计算贡献：矩阵乘法核心

当输入数据已在共享内存中后，每个线程需要计算其负责的输出子块。这通过一个三层嵌套循环实现，遍历输出子块的行、列以及共享内存中当前输入块的深度维度（BK）。

将结果从寄存器写回全局内存

计算完成后，每个线程需要将其寄存器中的 TM x TN 结果写回全局内存中的对应位置。类似加载过程，需要处理边界并尽可能实现合并写入。

核心优化技术

我们已经有了一个利用共享内存和寄存器的基础框架。本节中，我们来看看一系列能进一步提升性能的核心优化技术。

向量化加载与存储

现代GPU支持向量化内存事务，允许单个线程单条指令加载或存储多个连续的数据元素（如4个单精度浮点数，共16字节）。

向量化加载的优势：
在加载输入块到共享内存时，使用向量化加载可以让每个线程用更少的指令加载更多数据。这提高了向内存系统发出请求的速率，对于后续要讨论的低占用率情况尤为重要。

向量化存储的优势：
在将寄存器中的结果写回全局内存时，向量化存储有两个好处：

减少指令数：将多个标量存储合并为一条向量存储指令。
改善合并访问：线程写入连续的内存地址，提高了内存事务的效率。例如，每个线程存储其4x4子块中的一行（4个连续元素），同一线程束（Warp）内相邻线程存储的行是连续的，从而实现了良好的合并。

寄存器分块

在基础实现中，线程在计算输出子块时，会反复从共享内存中读取相同的输入值。共享内存虽然比全局内存快，但在计算核心全力工作时，它仍可能成为瓶颈。

寄存器分块思想：
将计算过程进一步分块。每个线程不是一次性计算整个 TM x TN 输出子块，而是分“条”进行：

从共享内存中加载一个 TM x 1 的A条和一个 1 x TN 的B条到寄存器中。
执行这两个向量的外积，更新整个 TM x TN 输出子块（所有数据都在寄存器中）。
重复步骤1和2，遍历完共享内存块的所有“条”。

这种方法将共享内存访问次数降到最低，计算完全在寄存器中进行，形成了密集的乘加指令流，极大地提升了计算单元的利用率。

优化存储的合并访问

即使使用向量化存储，如果线程负责的输出子块在内存中不连续，写入仍可能无法完全合并。

解决方案：重组线程块的工作分配。
我们不直接让一个线程负责一个大的、连续的输出子块（如8x8），而是让一个线程束（Warp） 共同负责一个较大的输出区域（如64x32）。然后，将这个Warp输出区域划分为多个象限，每个线程负责每个象限中的一个小块（如4x4）。

这样，当线程束内的所有线程写入同一个象限时，它们写入的4x4小块在内存中是连续排列的，从而实现了完美的合并存储。

消除共享内存体冲突

共享内存被组织成多个体（Bank）。当线程束中的多个线程同时访问同一个体的不同地址时，会发生体冲突，导致串行访问，降低性能。

问题来源：
在从共享内存加载数据条到寄存器时，如果数据在内存中的布局导致线程束内不同线程访问的地址间隔是体数量的整数倍（如32），就会发生体冲突。

解决方案：共享内存填充。
在声明共享内存数组时，人为地增加每一行的宽度（例如，实际宽度为8，但声明为9）。这种填充改变了数据在体中的映射关系，将原本可能冲突的访问分散到不同的体上。

处理低占用率

由于大量使用寄存器（用于输出子块、输入条以及循环展开产生的临时变量），高级矩阵乘法内核的寄存器使用量非常高，接近每个线程255个寄存器的上限。这严重限制了每个流多处理器（SM）上可同时驻留的线程数（占用率），可能低至12.5%。

低占用率带来两个问题：

隐藏计算延迟能力差：硬件调度器可选择的就绪线程束少。
隐藏内存延迟能力差：同时发出的内存请求少，难以饱和内存带宽。

为了应对低占用率，我们采用了以下策略：

激进的循环展开：增加指令级并行，减少线程束停顿。
向量化加载/存储：让有限的线程能以更高的速率发出内存请求。
软件流水线和线程束专业化：专门用于重叠计算和内存操作。

软件流水线

在基础循环中，加载下一对输入块的操作，必须等待当前输入块的计算完成（因为共享内存缓冲区被复用）。这是一种“假依赖”。

解决方案：双缓冲。
为输入块分配两个共享内存缓冲区。当使用缓冲区0中的数据进行计算时，可以异步地将下一对输入块加载到缓冲区1中。这消除了“计算完成”和“开始下一次加载”之间的依赖，使得加载和计算可以重叠。

软件流水线实现：
重构循环，使得在迭代 i 中，同时进行：

使用缓冲区 i%2 中的数据进行计算。
将下一对输入块加载到缓冲区 (i+1)%2 中。
这样，计算和内存传输指令在同一个迭代中交错，编译器可以更好地调度它们以隐藏延迟。

线程束专业化

软件流水线依赖编译器静态地交错指令。为了更动态地利用硬件调度器，我们可以使用线程束专业化。

思想：
将线程块中的线程束分为两组：

计算线程束：专门负责执行矩阵乘法的计算部分。
加载/存储线程束：专门负责在全局内存和共享内存之间传输数据。
硬件调度器会自动在这些不同任务的线程束之间进行切换，从而实现计算和内存传输的动态重叠，更有效地隐藏延迟。

利用专用硬件单元

现代GPU提供了专用于加速特定操作的硬件单元。

张量内存加速器

这是一种专用的数据搬移引擎，可以异步地将大块数据在全局内存和共享内存之间传输。

优势：

释放计算核心：计算核心不再需要执行加载/存储指令来进行数据搬移，可以专注于计算。
节省寄存器：数据搬移不再经过计算核心的寄存器文件，腾出了更多寄存器用于计算。
更好的重叠：专有引擎与计算核心并行工作，能更彻底地隐藏内存访问延迟。

张量核心

张量核心是执行小型矩阵乘加运算的专用硬件单元（例如，计算 D = A * B + C，其中A、B、C、D是小矩阵）。它们比通用的FP32核心快得多，并且功耗更低。

像cuBLAS、cuDNN和CUTLASS这样的高性能库，在内部已经集成了所有这些优化，并自动利用张量核心。因此，在实际应用中，应优先使用这些库，而不是自己实现。

总结与课程回顾

本节课中我们一起学习了GPU上矩阵乘法的高级优化技术。

我们从分析基础实现的算术强度出发，认识到其内存带宽受限的本质。为了将其转化为计算密集型任务，我们引入了共享内存分块来重用数据，并推导出增大输出块尺寸是提升性能的关键。

为了实现大尺寸输出块，我们采用了寄存器分块技术，将计算核心转移到更快的寄存器中进行。为了高效处理数据，我们使用了向量化加载/存储来提升内存吞吐，并重组了线程工作分配以实现合并访问。我们还通过共享内存填充消除了体冲突。

面对因高寄存器使用导致的低占用率挑战，我们采用了软件流水线和双缓冲来重叠计算与内存操作，并介绍了线程束专业化的概念以动态调度任务。最后，我们了解了如何利用张量内存加速器和张量核心这类专用硬件单元来获得终极性能。

虽然在实际开发中应直接使用优化库（如cuBLAS、CUTLASS），但理解这些底层优化原理至关重要。它们不仅是矩阵乘法的核心优化手段，其思想（如分块提升数据重用、平衡资源使用、重叠计算与访存）也适用于优化其他复杂的GPU计算任务。

posted @ 2026-03-26 12:21 布客飞龙III 阅读(66) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

贝鲁特大学-CMPS297S-GPU-编程笔记-全-

贝鲁特大学 CMPS297S GPU 编程笔记（全）

GPU计算：P1：引言 🚀

摩尔定律的兴衰 💡

并行计算成为主流 ⚙️

两种设计哲学：延迟导向 vs. 吞吐量导向 🏎️🚌

CPU与GPU的架构对比 🖥️🎮

GPU的起源与发展 📈

为何是GPU取得成功？🎯

总结 📚

GPU计算：P2：数据并行编程

概述

课程回顾

数据并行与任务并行

向量加法示例

GPU系统架构

CUDA内存管理

数据传输

启动GPU内核进行计算

实现GPU内核

代码编译

性能考量与同步

错误检查

总结

GPU计算：P3：多维网格与数据 🚀

多维网格与数据布局

创建多维网格

线程索引计算

多维数据的内存布局

示例一：RGB图像转灰度图 🖼️➡️⚫⚪

并行化策略

示例二：图像模糊 🔍➡️🌫️

并行化策略

示例三：矩阵乘法 ➕✖️➡️📈

并行化策略

总结与关键要点 🎯

GPU计算：04：GPU架构

概述

回顾：多维数组与数据并行

GPU架构基础

线程块到SM的映射

线程块内的协作与同步

线程块间的独立性与可扩展性

Warp：SM上的调度单元

SIMD执行与Warp

控制流分歧

延迟隐藏与多线程

占用率

总结

GPU计算：第5讲：内存与分块

概述：性能指标与计算瓶颈

向量加法的计算/内存比

矩阵乘法的潜力与问题

GPU内存架构与CUDA内存模型

GPU内存层次结构

CUDA编程模型中的内存

分块优化与共享内存实践

矩阵乘法的分块策略

分块矩阵乘法代码实现

边界条件处理

性能提升与注意事项

CPU上的分块优化

总结

GPU计算：P6：性能考量 🚀

概述 📋

DRAM架构回顾 🧠

内存合并优化 ⚡

多DRAM存储体与延迟隐藏 🏦

线程粒度与线程粗化 🧵

优化清单与权衡 ✅

总结 🎯

GPU计算：P7：性能剖析 🚀

概述

性能瓶颈与优化权衡

使用 nvprof 进行基础剖析

生成剖析文件以供深入分析

使用 NVIDIA Visual Profiler (nvvp) 进行可视化分析

分析时间线

分析内核性能

使用 `nvprof` 进行基础剖析