UIUC-ECE408-应用并行编程笔记-全-

UIUC ECE408 应用并行编程笔记（全）

001：课程介绍与GPU计算基础

概述

在本节课中，我们将学习应用并行编程课程的基本信息、课程结构、学习目标，并初步了解GPU（图形处理器）计算的基础知识，包括其历史背景、设计哲学以及与CPU（中央处理器）的关键区别。

课程介绍

欢迎来到应用并行编程课程。本课程旨在教授如何为大规模并行处理器（主要是GPU）进行编程，以充分发挥其高性能计算潜力，同时保证代码的功能性和可维护性。

本课程由伊利诺伊大学于2007年首创，由Wen-mei Hwu教授和NVIDIA首席架构师David Kirk共同讲授，是首批教授GPU编程的课程之一。

课程技术主题

我们将涵盖以下技术主题：

  • 并行编程基础
  • 并行算法的原理与模式
  • 编程API、工具和技术
  • 处理器架构及其约束
  • 关键应用案例（尤其是AI领域）

教学团队与资源

  • 讲师：Steve Lumetta教授。办公室时间为周二和周四上午9点至11点（讲座之间的同一时间段）。
  • 助教：Hanwen。办公室时间为周六。
  • 课程页面：所有公告、讲义和代码链接都将发布在课程网页上。
  • Piazza：用于课程讨论和问答。
  • Canvas：仅用于成绩发布。
  • GitHub：用于实验和项目作业的提交。
  • 教材：《Programming Massively Parallel Processors: A Hands-on Approach》（第四版）。

课程评分构成

课程总成绩由以下部分构成：

  • 考试（50%）：两次考试，各占25%。考试将在Zoom上进行，需要开启摄像头。
  • 实验（25%）：共9个实验（含第0个环境设置实验），每个实验权重相等。
  • 项目（25%）：一个大型个人项目，评估功能、性能、演示和代码风格。

迟交政策与学术诚信

  • 每个学生有两次延期提交实验或项目前两个检查点的机会，每次最多延期48小时，但必须在原截止日期前通过Google表单申请。
  • 我们将去掉最低分的实验+测验组合。
  • 我们鼓励讨论课程概念，但严禁复制代码、在考试中作弊、规避实验要求或使用AI生成代码。违反学术诚信政策将导致严重后果（首次违规该次考试记零分，再次违规课程不及格）。

实验环境设置：Delta超级计算机

我们将使用位于NCSA的Delta超级计算机进行实验和项目。

Delta配置简介

Delta节点的主要配置如下：

  • GPU：每个节点有4个NVIDIA A40 GPU，通过NVLink高速互连。每个GPU拥有48GB DDR6显存。
  • CPU：每个节点有1个AMD Milan CPU，包含64个核心，主频约2.5 GHz。
  • 内存：每个节点有256 GB CPU内存。

环境设置步骤

  1. 申请Delta账户：通过课程页面链接申请ACCESS ID，建议使用伊利诺伊大学NetID。
  2. 设置GitHub仓库：通过课程页面链接创建个人课程代码仓库。
  3. 配置身份验证：建议使用个人访问令牌（PAT）而非SSH密钥在Delta上访问GitHub。
  4. 完成实验0：在截止日期前（6月20日周五晚8点，美国中部夏令时），登录Delta，编译并运行示例代码，保存输出结果用于后续测验。

注意：Delta计划于6月19日周三进行维护，届时将无法访问。

GPU计算基础

上一节我们介绍了课程的基本设置，本节中我们来看看GPU为何成为现代计算的核心，以及它与CPU的根本区别。

计算范式的转变

在20世纪，计算主要用于辅助测量、控制和理解物理过程。进入21世纪，计算本身成为了发现和创新的主要驱动力，例如计算成像、药物设计和AI模型。

这一转变部分源于“后登纳德缩放”时代的技术转折。登纳德缩放定律指出，随着晶体管尺寸缩小，芯片功耗可以保持稳定。然而，约在2004年，业界为追求更高性能而大幅提升时钟频率，导致了“功耗墙”问题。此后，行业转向通过增加核心数量（而非提升单核频率）来提升性能。

CPU与GPU的设计哲学

CPU和GPU针对不同的计算模式进行了优化：

  • CPU：面向延迟优化
    • 目标：尽可能快地执行单个任务（线程）。
    • 手段：高时钟频率、大容量缓存、复杂的控制逻辑（如分支预测）、数据前推。
    • 公式：追求最小化 任务完成时间。
  • GPU：面向吞吐量优化
    • 目标：在单位时间内完成尽可能多的任务（线程）。
    • 手段：适中时钟频率、大量计算核心、简单控制逻辑、高能效ALU、深度流水线。
    • 公式：追求最大化 单位时间内完成的任务数量。

为何需要两者协同？

一个成功的应用策略是混合使用CPU和GPU：

  • 将顺序执行、对延迟敏感的部分放在CPU上运行。
  • 将数据并行、计算密集的部分放在GPU上运行。
    程序员需要判断代码的哪部分适合GPU，并管理CPU与GPU之间的数据移动，因为两者拥有独立的内存空间。

并行编程的挑战

编写高效的GPU程序并非简单地将代码移植过来，我们需要考虑：

  1. 负载均衡：所有并行线程应尽可能同时完成，避免某些线程长时间运行而拖慢整体进度。
  2. 内存带宽：计算能力（FLOPS）的增长速度远快于内存带宽（Byte/s）的增长。这形成了“内存墙”，程序员必须精心设计数据访问模式以隐藏内存延迟。
  3. 可扩展性与可维护性：算法应能有效利用不断增加的计算核心，并且代码应易于长期维护。

总结

本节课我们一起学习了应用并行编程课程的框架、评分政策和使用Delta超级计算机的流程。我们探讨了GPU兴起的历史背景，并深入理解了CPU（延迟优化）与GPU（吞吐量优化）两种根本不同的处理器设计哲学。认识到需要根据任务特性（顺序执行 vs. 数据并行）合理分配计算资源，并初步了解了编写高效并行程序所面临的主要挑战（如负载均衡和内存带宽限制）。在接下来的课程中，我们将开始学习CUDA编程模型，亲手将理论应用于实践。

002：CUDA 编程模型入门 🚀

在本节课中，我们将学习数据并行计算的基本概念，以及 CUDA C 编程接口的核心特性。我们将从 GPU 编程面临的挑战开始，逐步深入到 CUDA 的执行模型、线程组织、内存管理以及如何编写和编译一个简单的 CUDA 程序。

回顾：GPU 编程的挑战

上一节我们介绍了在高吞吐量、大规模并行系统（如 GPU）上进行编程时需要面对的几个核心问题。本节中，我们来看看这些挑战的具体表现。

以下是 GPU 编程中需要解决的关键问题：

  • 可扩展性：算法需要能够随着数据量的增加而提供更多的并行性。二次方或更复杂度的算法通常不适合 GPU。
  • 规律性：每个线程或数据实例的处理方式应尽可能相同，以利于并行化。
  • 负载均衡：每个数据块的处理时间应大致相同，以避免部分线程空闲等待。
  • 全局内存带宽：GPU 的计算能力增长速度远超内存带宽，因此程序性能常受限于内存访问速度。
  • 串行化：当大量线程试图同时访问同一资源（如内存地址）时，会产生排队等待，导致并行效率下降。

为什么学习 CUDA 编程？

掌握如何为这些高吞吐量实体编写高效代码至关重要。CUDA 的创始人认为，教授人们编写可扩展、可移植且能跨越近 20 年 GPU 硬件演进的代码，是其对社会最大的贡献之一。可扩展的算法和库是我们能留下的最佳遗产，因为软件的生命周期往往远超预期。

CUDA 执行模型简介

CUDA 采用一种集成的编程模型，程序包含在 CPU（主机）上运行的代码和在 GPU（设备）上运行的代码。编译器通过代码中的注解来区分这两部分。

  • 主机代码：处理串行或适度并行的部分，在 CPU 上执行。
  • 设备代码：处理高度并行的部分，在 GPU 上执行，采用 SPMD 模型。

什么是 SPMD？

SPMD 代表 单程序多数据。在 GPU 上，一个被称为 内核 的函数会被成千上万个线程同时执行。所有线程执行相同的代码（单程序），但处理不同的数据（多数据）。

执行流程如下：

  1. 主机代码运行。
  2. 遇到 内核启动 时，成千上万的线程在 GPU 上开始执行该内核函数。
  3. 所有线程执行完毕后，控制权返回主机，继续执行后续代码。
  4. 此过程可以反复进行。

CUDA 线程组织：网格与线程块

一个 CUDA 内核由一组线程执行，这组线程被称为 网格。网格是一个三维的 线程块 数组，而每个线程块本身又是一个三维的 线程 数组。这种多维结构便于处理图像、矩阵等多维数据。

每个线程都有唯一的标识符，用于决定它处理哪部分数据。

网格维度

网格的维度由 gridDim 这个结构体的 x、y、z 成员表示，它们分别代表网格在三个方向上的线程块数量。

线程块维度与索引

每个线程块的维度由 blockDim 的 x、y、z 成员表示，即每个线程块在各个方向上的线程数量。
线程块在网格中有唯一的 blockIdx。
线程在线程块内有唯一的 threadIdx。

一个线程的唯一全局标识 由 blockIdx 和 threadIdx 共同决定。

线程协作

  • 线程块内协作：线程可以通过共享内存、原子操作和屏障同步进行紧密协作。
  • 线程块间协作：不同线程块间的协作能力较弱，主要通过全局内存进行通信，且需要特别注意同步和一致性问题。

实战：向量加法

让我们通过一个经典的例子——向量加法，来具体了解 CUDA 编程步骤。

目标：计算 C[i] = A[i] + B[i]，其中 i 从 0 到 n-1。

主机端（CPU）C 代码

void vecAddHost(float* A_h, float* B_h, float* C_h, int n) {
    for (int i = 0; i < n; i++) {
        C_h[i] = A_h[i] + B_h[i];
    }
}

CUDA 实现步骤

将一个顺序的 C 函数转化为 CUDA 程序需要三个主要步骤：

  1. 在设备上分配内存并传输数据
  2. 在 GPU 上执行内核（进行计算）
  3. 将结果拷贝回主机并释放设备内存

步骤 1：设备内存分配与数据传输

由于 GPU 拥有独立的内存空间，必须显式地进行内存管理和数据搬运。

  • 分配设备内存：使用 cudaMalloc。
    float *A_d, *B_d, *C_d;
cudaMalloc((void**)&A_d, n * sizeof(float));
cudaMalloc((void**)&B_d, n * sizeof(float));
cudaMalloc((void**)&C_d, n * sizeof(float));
  • 拷贝数据到设备：使用 cudaMemcpy，方向为 cudaMemcpyHostToDevice。
    cudaMemcpy(A_d, A_h, n * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(B_d, B_h, n * sizeof(float), cudaMemcpyHostToDevice);
  • 释放设备内存：计算完成后，使用 cudaFree。
    cudaFree(A_d); cudaFree(B_d); cudaFree(C_d);

步骤 2：编写与启动内核

内核函数 是一个带有 __global__ 限定符的特殊函数，它将在 GPU 上执行。

__global__ void vecAddKernel(float* A_d, float* B_d, float* C_d, int n) {
    // 计算当前线程的全局索引
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    // 确保索引在数组范围内
    if (i < n) {
        C_d[i] = A_d[i] + B_d[i];
    }
}

内核启动 使用特殊的 <<<...>>> 语法配置执行参数。

// 假设每个线程块有 256 个线程
int threadsPerBlock = 256;
// 计算需要的线程块数量，并向上取整以确保覆盖所有数据
int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
// 启动内核
vecAddKernel<<<blocksPerGrid, threadsPerBlock>>>(A_d, B_d, C_d, n);

步骤 3：取回结果

使用 cudaMemcpy 将结果从设备内存拷贝回主机内存，方向为 cudaMemcpyDeviceToHost。

cudaMemcpy(C_h, C_d, n * sizeof(float), cudaMemcpyDeviceToHost);

扩展到多维：图像处理示例

向量加法是一维的。对于像图像处理这样的二维问题，CUDA 的多维线程组织显得更加自然。

灰度转换

将彩色（RGB）图像转换为灰度图像。每个输出像素独立计算，公式通常为：灰度值 = 0.299*R + 0.587*G + 0.114*B。

内核函数示例：

__global__ void colorToGrayscale(unsigned char* rgbImage,
                                 unsigned char* grayImage,
                                 int width, int height) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    if (col < width && row < height) {
        int grayOffset = row * width + col;
        int rgbOffset = grayOffset * 3; // 每个 RGB 像素占 3 个字节
        unsigned char r = rgbImage[rgbOffset];
        unsigned char g = rgbImage[rgbOffset + 1];
        unsigned char b = rgbImage[rgbOffset + 2];
        grayImage[grayOffset] = 0.299f*r + 0.587f*g + 0.114f*b;
    }
}

内核启动（使用二维网格和线程块）：

dim3 blockDim(16, 16); // 16x16 的线程块
dim3 gridDim((width + blockDim.x - 1) / blockDim.x,
             (height + blockDim.y - 1) / blockDim.y);
colorToGrayscale<<<gridDim, blockDim>>>(d_rgbImage, d_grayImage, width, height);

图像模糊（卷积）

图像模糊需要每个输出像素是其周围一个区域（如 3x3、5x5）内像素的平均值。这引入了数据共享和边界处理的问题。

  • 数据共享：相邻的线程需要读取重叠的输入像素区域。
  • 边界处理：对于图像边缘的像素，其卷积核会超出图像边界，需要特殊处理（如忽略界外像素或进行填充）。

内核中需要额外的条件判断来处理边界，确保不访问非法内存地址。

CUDA 函数限定符

CUDA 扩展了 C/C++，引入了特殊的函数限定符来指示函数的执行位置和调用方式：

  • __global__：声明一个内核函数。从主机调用，在设备上执行。必须返回 void。
  • __device__：声明一个设备函数。从内核或其他设备函数调用，在设备上执行。
  • __host__：声明一个主机函数（默认）。在主机上执行。可与 __device__ 联用，使同一个函数同时为主机和设备编译。

编译 CUDA 程序

CUDA 程序使用 nvcc 编译器进行编译。nvcc 的工作流程如下：

  1. 将代码分离为主机代码和设备代码。
  2. 主机代码交给系统的 C/C++ 编译器（如 gcc）编译。
  3. 设备代码被编译为一种称为 PTX 的中间汇编语言。
  4. 在程序运行时，GPU 驱动程序会根据具体的 GPU 架构，将 PTX 代码即时编译为优化的机器码。

重要概念与注意事项

  1. 线程块大小选择：通常选择较大的线程块（如 256、512）以更好地利用 GPU 资源。最大值是 1024。
  2. 线程块执行顺序：线程块以任意顺序执行，且彼此独立。不要在内核代码中假设线程块之间的执行顺序或完成状态。
  3. 内存一致性：不同线程块之间通过全局内存进行通信时，没有严格的读写一致性保证，需要程序员使用原子操作或同步原语来管理。
  4. 资源限制：每个线程块可用的共享内存、寄存器等资源是有限的，这会影响活跃线程块的数量和性能。

总结

本节课中，我们一起学习了 CUDA 并行编程的基础。我们从 GPU 编程的挑战出发，深入探讨了 CUDA 的 SPMD 执行模型，理解了网格、线程块和线程的多维层次组织。我们通过向量加法的实例，掌握了 CUDA 编程的核心流程：设备内存管理、内核编写与启动、数据传输。接着，我们将知识扩展到图像灰度转换和模糊等多维应用，并学习了如何处理边界条件。最后，我们介绍了 CUDA 的函数限定符、编译过程以及编写高效 CUDA 程序需要注意的关键点，如线程块执行顺序和内存一致性。这些概念是构建更复杂、高性能 GPU 应用程序的基石。

003：CUDA线程组织与内存模型 🚀

概述

在本节课中，我们将学习CUDA编程中线程的多维逻辑组织方式，以及CUDA的内存层次结构模型。我们将深入探讨如何通过合理的线程块配置来优化性能，并理解不同内存类型（寄存器、共享内存、全局内存）的特性、访问速度和使用方法。这些知识是编写高效GPU程序的基础。

多维线程组织回顾

上一节我们介绍了CUDA的网格-线程块-线程三级层次结构。本节中，我们来看看如何利用这种结构进行更复杂的计算。

一个CUDA内核启动时，会创建一个由线程块组成的网格。每个线程块本身又是一个由线程组成的多维数组。虽然逻辑上是三维的，但任何维度的尺寸都可以是1，因此我们常用1D或2D结构。

每个线程通过其块索引 (blockIdx) 和线程索引 (threadIdx) 的组合来获得一个唯一的标识符。内核代码可以利用这些标识符来计算数据数组中的唯一索引。

核心概念公式：
对于一个2D网格和2D线程块，计算全局线程索引的常用方法是：

int col = blockIdx.x * blockDim.x + threadIdx.x; // 列索引
int row = blockIdx.y * blockDim.y + threadIdx.y; // 行索引
int global_index = row * total_width + col; // 线性化索引

线程块执行与调度

线程块被分配给流多处理器执行。一旦分配，整个线程块将在同一个SM上运行，不会被拆分。

线程束与SIMD执行

在SM内部，线程块中的线程被分组为大小为32的线程束。线程束是SM的基本调度单元。

  • 线程束形成：将线程块中所有线程的(threadIdx.z, threadIdx.y, threadIdx.x)按x变化最快、z变化最慢的规则线性化，然后每32个线程分为一个线程束。
  • SIMD执行：当SM选择一个线程束执行时，该线程束中的所有32个线程执行相同的指令。这被称为单指令多数据执行模型。

零开销线程束调度

SM采用零开销线程束调度。这意味着SM可以在每个时钟周期无延迟地在所有就绪的线程束之间切换。这种机制对于隐藏内存访问等长延迟操作至关重要。

控制流与分支发散

由于线程束以SIMD方式执行，当线程束内的线程遇到条件分支（如if-else）并走向不同路径时，就会发生分支发散。

如何处理分支发散

GPU使用谓词执行来处理分支发散：

  1. 首先，所有线程执行if的then代码块，但只有条件为真的线程会保留执行结果。
  2. 然后，所有线程执行else代码块，只有条件为假的线程会保留执行结果。

性能影响：分支发散会导致SM串行执行所有分支路径，显著降低性能。最坏情况是线程束中只有一个线程活跃，而其他31个线程空闲。

如何避免分支发散

  • 确保线程束内控制流一致：尽量让同一个线程束内的所有线程执行相同的控制路径。
  • 调整数据访问模式：例如，通过将线程索引除以线程束大小（如 threadIdx.x / warpSize），可以使同一个线程束内的线程获得相同的商值，从而走向相同的分支。

线程块大小的选择策略

选择线程块大小时，需要考虑多种硬件约束，目标是最大化SM上可同时调度的线程束数量，以隐藏延迟。

以下是选择线程块大小需要考虑的几个关键因素：

  1. 每个SM的最大线程数：限制了可以同时驻留的线程总数。
  2. 每个SM的最大线程块数：限制了可以同时调度的线程块数量。
  3. 每个线程块的最大线程数：当前GPU通常为1024。
  4. 共享内存大小：每个线程块可用的共享内存是有限的。
  5. 寄存器文件大小：每个线程可用的寄存器数量受限于SM的总寄存器数和活跃线程数。

选择策略：通常需要在上述约束间进行权衡。一个常见的经验法则是，线程块中的线程数最好是线程束大小（32）的倍数，并且足够大以充分利用SM资源，但又不能太大以至于限制了活跃线程块的数量。

CUDA内存模型

现在，我们转向CUDA的内存层次结构。理解不同内存的带宽、延迟和用途对于优化程序性能至关重要。

内存层次概述

CUDA提供了多种内存空间，供程序员显式或隐式地使用：

  1. 寄存器：速度最快，延迟约1周期。每个线程私有，用于存储自动变量（非数组的局部变量）。生命周期与线程相同。
  2. 共享内存：速度很快，延迟约5-10周期。每个线程块共享，用于线程块内线程间的通信。生命周期与线程块相同。
  3. 全局内存：速度慢，延迟约500-1000周期。所有线程可读写，容量大（GB级别）。数据在主机与设备间传输。生命周期与应用程序相同。
  4. 常量内存：只读，速度较快（有缓存）。所有线程可读，由主机初始化。生命周期与应用程序相同。

在代码中使用不同内存

  • 寄存器：在内核函数中声明非数组的局部变量（如 int temp; float sum = 0.0f;）。
  • 共享内存：使用 __shared__ 限定符声明变量（如 __shared__ float tile[TILE_WIDTH][TILE_WIDTH];）。
  • 全局/常量内存：在全局作用域使用 __device__ 或 __constant__ 限定符声明变量。

案例研究：矩阵乘法

我们将以矩阵乘法为例，展示如何将并行化思想应用于实际问题，并分析其性能瓶颈。

算法描述

给定两个宽度为 Width 的方阵 M 和 N，计算输出矩阵 P = M * N。矩阵乘法的定义是：
P[row][col] = sum over k ( M[row][k] * N[k][col] )
其中 row 和 col 的范围是 0 到 Width-1。

并行化策略

一个直观的并行化策略是每个线程计算输出矩阵P的一个元素。

  • 线程组织：使用2D网格和2D线程块。网格大小覆盖整个输出矩阵P，线程块大小（如16x16）是一个可调参数。
  • 内核伪代码：
    __global__ void matrixMulKernel(float* M, float* N, float* P, int Width) {
    // 计算该线程负责的P矩阵的行和列
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < Width && col < Width) { // 边界检查
        float pValue = 0;
        for (int k = 0; k < Width; ++k) {
            pValue += M[row * Width + k] * N[k * Width + col];
        }
        P[row * Width + col] = pValue;
    }
}

性能瓶颈分析

上述简单内核的性能受限于全局内存带宽。

  • 计算与内存访问比：每个线程计算一个P元素需要进行Width次内积迭代。每次迭代需要从全局内存读取M和N各一个元素（共8字节），并执行一次乘法和一次加法（2次浮点运算）。
  • 算术强度：计算访存比 = 2次浮点运算 / 8字节 = 0.25 FLOP/Byte。
  • 瓶颈：现代GPU的峰值计算能力远高于其内存带宽。例如，一个具有1000 GFLOP/s计算能力和150 GB/s内存带宽的GPU，其理论最大性能受限于内存带宽：150 GB/s * 0.25 FLOP/Byte = 37.5 GFLOP/s，仅为峰值性能的3.75%。

结论：简单的全局内存访问模式无法充分利用GPU的强大计算能力。我们需要通过数据重用和使用更快的内存（如共享内存）来减少对全局内存的访问。

总结

本节课我们一起学习了CUDA线程组织的核心概念，包括线程束、分支发散和线程块大小选择策略。我们还深入探讨了CUDA的内存层次模型，了解了寄存器、共享内存和全局内存的特性与用法。最后，我们以矩阵乘法为例，分析了其并行化实现和由全局内存带宽导致的性能瓶颈。在接下来的课程中，我们将学习如何利用共享内存来优化类似矩阵乘法的计算密集型内核，从而显著提升程序性能。

004：共享内存与分块矩阵乘法

在本节课中，我们将学习如何利用GPU的共享内存来优化矩阵乘法运算。我们将分析全局内存访问的瓶颈，并引入“分块”算法，通过将数据块加载到更快的共享内存中来减少全局内存访问次数，从而显著提升计算性能。

概述：全局内存瓶颈与共享内存的引入

上一节我们分析了简单矩阵乘法内核的性能，发现其受限于全局内存带宽。本节中，我们将探讨如何利用共享内存来缓解这一瓶颈。

简单矩阵乘法内核的每个线程需要从全局内存中读取 2 * width 个浮点数（来自矩阵 M 和 N），并写入1个浮点数（到矩阵 P）。这导致计算访存比很低，约为 2 次浮点运算 / 8 字节数据 = 0.25 FLOPS/字节。在带宽为150 GB/s的GPU上，理论性能上限仅为37.5 GFLOPS。

共享内存的访问延迟远低于全局内存（约5个周期 vs 500个周期），但容量较小（例如48KB）。为了利用它，我们需要将大矩阵 M 和 N 分解成更小的块（称为“分块”或“瓦片”），分阶段加载到共享内存中进行计算。

分块矩阵乘法的核心思想

分块算法的核心是数据复用。在标准矩阵乘法 P = M * N 中：

  • 矩阵 M 的每一行被用于计算 P 的对应整行。
  • 矩阵 N 的每一列被用于计算 P 的对应整列。

简单内核中，同一个 M 或 N 元素会被不同线程重复从全局内存中读取 width 次。分块算法旨在将一个数据块一次性读入共享内存，让线程块内的所有线程共享该数据块，从而将全局内存访问次数减少约 tile_width 倍。

以下是实现此想法的步骤概述：

  1. 分配共享内存：为 M 和 N 的数据块分配共享内存空间。
  2. 协作加载：线程块内的所有线程协作，将全局内存中的一个 M 块和一个 N 块加载到共享内存中。
  3. 屏障同步：确保所有线程完成数据加载后，再开始计算。
  4. 计算部分结果：线程从共享内存（而非全局内存）读取数据，计算部分内积。
  5. 屏障同步：确保所有线程完成当前块的计算后，再覆盖共享内存以加载下一个数据块。
  6. 循环：重复步骤2-5，直到处理完所有需要的 M 和 N 的数据块。
  7. 写入结果：将最终结果写回全局内存中的 P。

内核代码解析：从简单版本到分块版本

首先，回顾一下简单的矩阵乘法内核：

__global__ void MatrixMulKernel(float* d_M, float* d_N, float* d_P, int width) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (row < width && col < width) {
        float p_value = 0;
        for (int k = 0; k < width; ++k) {
            p_value += d_M[row * width + k] * d_N[k * width + col]; // 全局内存访问
        }
        d_P[row * width + col] = p_value; // 全局内存访问
    }
}

现在，我们来看分块版本的内核。以下代码假设矩阵宽度 width 是块大小 TILE_WIDTH 的整数倍（边界处理稍后讨论）。

__global__ void MatrixMulTiledKernel(float* d_M, float* d_N, float* d_P, int width) {
    // 1. 分配共享内存
    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    // 计算线程对应的输出位置
    int bx = blockIdx.x; int by = blockIdx.y;
    int tx = threadIdx.x; int ty = threadIdx.y;
    int col = bx * TILE_WIDTH + tx;
    int row = by * TILE_WIDTH + ty;

    float p_value = 0;

    // 循环遍历所有需要的分块
    for (int m = 0; m < (width / TILE_WIDTH); ++m) {
        // 2. 协作加载分块到共享内存
        // 加载 M 的分块
        Mds[ty][tx] = d_M[row * width + (m * TILE_WIDTH + tx)];
        // 加载 N 的分块
        Nds[ty][tx] = d_N[(m * TILE_WIDTH + ty) * width + col];

        // 3. 屏障同步：等待块内所有线程完成加载
        __syncthreads();

        // 4. 计算部分内积（从共享内存读取）
        for (int k = 0; k < TILE_WIDTH; ++k) {
            p_value += Mds[ty][k] * Nds[k][tx];
        }

        // 5. 屏障同步：确保所有线程用完当前共享内存数据后再加载下一块
        __syncthreads();
    }

    // 7. 将最终结果写回全局内存
    if (row < width && col < width) {
        d_P[row * width + col] = p_value;
    }
}

关键点解析：

  • __shared__：声明变量位于共享内存中，该内存区域在线程块内共享。
  • 索引计算：加载 M 和 N 时，需要根据当前分块索引 m 进行偏移。M 的列偏移和 N 的行偏移都与 m * TILE_WIDTH 相关。
  • __syncthreads()：CUDA内置的屏障同步函数。调用该函数的线程将等待，直到同一线程块内所有活跃线程都调用此函数后，才会继续执行。这对于保证共享内存读写顺序的正确性至关重要。
  • 共享内存使用：计算时，线程从 Mds 中读取一行，从 Nds 中读取一列，在共享内存中进行乘加运算。

性能提升分析

假设我们使用 TILE_WIDTH = 16。

  • 数据复用：每个从全局内存加载到共享内存的 M 或 N 元素，会被线程块内的16个线程各使用一次。全局内存访问次数减少了约16倍。
  • 计算访存比：现在，对于每个加载的 M 和 N 元素（共8字节），我们执行了 TILE_WIDTH (16) 次乘加运算（32次浮点运算）。计算访存比提升为 32 FLOPS / 8 字节 = 4 FLOPS/字节。
  • 理论性能上限：在之前150 GB/s带宽的GPU上，理论性能上限变为 150 * 4 = 600 GFLOPS，更接近GPU的峰值算力（例如1000 GFLOPS）。

使用更大的分块（如32x32）可以进一步提高复用因子，但会受到共享内存容量和SM内最大线程数的限制。

线程同步与屏障 (__syncthreads)

正确使用 __syncthreads() 是编写分块算法的关键。以下是关于CUDA屏障的重要规则：

  1. 作用范围：仅在同一线程块内同步线程。不同块间的线程无法直接通过此函数同步。
  2. 必须全部执行：线程块内所有仍有代码要执行的活跃线程都必须到达同一个 __syncthreads() 调用点。
  3. 静态调用点：所有线程必须从代码的同一个静态位置调用 __syncthreads()。不能根据条件分支让部分线程从一个点调用，另一部分从另一个点调用。
  4. 内存操作保证：在调用 __syncthreads() 之前发起的对全局内存和共享内存的读写操作，保证在该屏障对所有线程完成之后是可见的（即已完成）。

何时需要屏障？

  • 写-读 hazard：当一个线程写入共享内存后，另一个线程要读取该位置时，需要在写之后、读之前同步。
  • 读-写 hazard：当所有线程用完共享内存中的一块数据后，要加载新数据覆盖它之前，需要同步，确保没有线程还在读旧数据。

思考题：以下情况是否需要 __syncthreads()？

  1. 线程将数据写入 shared[ty][tx]，随后又从 shared[ty][tx] 读取。
    • 答案：不需要。线程读取的是自己刚写入的值。
  2. 线程将数据写入 shared[ty][tx]，随后从 shared[tx][ty] 读取（交换了索引）。
    • 答案：需要。因为读取的位置可能是由其他线程写入的。
  3. 在情况2之后（已有一个屏障），线程再从 shared[TILE_WIDTH-1-tx][ty] 读取。
    • 答案：不需要。因为屏障后没有线程再修改 shared，所有读取操作都是安全的，顺序无关紧要。
  4. 在情况3之后，线程再执行 shared[ty][tx] = new_value。
    • 答案：需要。因为这次写入可能会覆盖其他线程在情况3中尚未读取的数据，构成读-写 hazard。

处理边界条件

上述分块内核假设矩阵宽度是分块宽度的整数倍。在实际应用中，我们需要处理任意大小的矩阵，包括非整数倍的情况和矩形矩阵。

核心策略：条件加载与归零
在将数据从全局内存加载到共享内存时，检查要加载的索引是否在矩阵有效范围内。

  • 如果在范围内，正常加载。
  • 如果超出范围（例如，在矩阵右边缘或下边缘的分块中），则向共享内存的对应位置写入 0.0f。

这样处理的好处是：

  1. 在计算阶段，即使线程访问了共享内存中这些“填充”的零值，乘加运算的结果也不受影响（任何数乘以零还是零）。
  2. 无需在计算循环内部添加复杂的边界判断，保持了内循环的高效性。

修改后的加载逻辑示例：

// 加载 M 的分块，带边界检查
int load_row_M = row;
int load_col_M = m * TILE_WIDTH + tx;
if (load_row_M < width && load_col_M < width) {
    Mds[ty][tx] = d_M[load_row_M * width + load_col_M];
} else {
    Mds[ty][tx] = 0.0f;
}

// 加载 N 的分块，带边界检查
int load_row_N = m * TILE_WIDTH + ty;
int load_col_N = col;
if (load_row_N < width && load_col_N < width) {
    Nds[ty][tx] = d_N[load_row_N * width + load_col_N];
} else {
    Nds[ty][tx] = 0.0f;
}

其他必要的修改：

  • 分块数量计算：循环次数应为 (width + TILE_WIDTH - 1) / TILE_WIDTH，向上取整以确保覆盖整个矩阵。
  • 输出保护：在将结果 p_value 写回全局内存 d_P 时，必须检查 row 和 col 是否在 P 的有效范围内，防止越界写入。

总结

本节课中我们一起学习了如何利用共享内存优化矩阵乘法：

  1. 识别瓶颈：简单矩阵乘法的性能受限于全局内存带宽。
  2. 引入分块：通过将数据分块加载到共享内存，利用数据在线程块内的复用性，大幅减少全局内存访问。
  3. 编写内核：我们逐步构建了分块矩阵乘法内核，包括共享内存分配、协作加载、屏障同步和计算循环。
  4. 理解同步：深入探讨了 __syncthreads() 屏障的作用、规则和使用场景，这是保证并行程序正确性的关键。
  5. 处理边界：学习了通过条件加载和归零技术来处理非整数倍分块和矩阵边界，使算法更具通用性。

通过应用分块技术，我们成功将矩阵乘法的计算重点从昂贵的全局内存转移到了快速的共享内存，从而显著提升了程序的潜在性能，使其更充分地利用GPU的强大计算能力。在接下来的实验中，你将有机会亲手实现这一优化。

005：动态随机存取内存与内存访问合并

概述

在本节课中，我们将学习动态随机存取内存的组织结构和工作原理，包括其突发传输模式和多存储体并行机制。我们还将探讨GPU编程中的一个关键概念——内存访问合并，并学习如何编写代码以实现高效的内存带宽利用。

DRAM基础：动态随机存取内存

上一节我们介绍了并行计算的基本概念，本节中我们来看看现代计算机中用于大容量存储的核心技术：动态随机存取内存。

动态随机存取内存是一种存储技术，其中每个比特（0或1）都存储在一个电容器中。电容器可以看作是两个平行的金属板，电子聚集在其中一块板上，从而在两板之间形成电场。这个电场可以保持电荷，从而存储信息。

每个存储单元由一个电容器和一个晶体管组成。晶体管作为一个电压控制的开关，连接电容器和用于读写的数据线。然而，电容器上的电荷会随着时间的推移而泄漏，通常在约50毫秒后就会丢失。因此，DRAM需要定期（例如每30-40毫秒）刷新所有比特，这就是“动态”一词的由来。

由于每个比特只需要一个晶体管和一个电容器，DRAM单元可以做得非常小，从而允许在单个芯片上集成数十亿个存储单元。

DRAM的组织结构：核心阵列与突发模式

了解了DRAM的基本单元后，我们来看看这些单元是如何组织起来工作的。

一个典型的DRAM存储体包含一个核心阵列，大约有1000条数据线，每条数据线连接约1000个存储单元，总计约100万个比特。地址通过行解码器输入，解码后驱动相应的选择线，激活一整行（约1000个比特）的存储单元。这些比特通过数据线进入读出放大器，放大微小的电压变化，并将其锁存在锁存器中。随后，通过多路复用器，根据列地址选择出所需的特定比特或一小组比特。

DRAM的核心阵列内部没有时钟信号，访问任何地址所需的时间是相同的。然而，在芯片级别，DRAM与CPU或GPU之间的接口是时钟同步的。

关键点：每次从核心阵列读取时，必须读取整行数据。读取操作会耗尽该行所有存储单元的电容器电荷，因此之后必须将整行数据写回。这意味着，即使你只需要该行中的几个比特，DRAM也会提供整行数据。

这种特性引出了突发传输模式的概念。与其只读取请求的几个比特，不如利用已经读取的整行数据，连续传输多个相邻的比特。这样，在接口等待下一次核心阵列访问的较长延迟期间，可以持续传输数据，从而更有效地利用内存接口带宽。

在现代系统中，DRAM几乎总是以突发模式运行，以最大化接口效率。如果应用程序不需要突发传输中的所有数据，可以简单地忽略它们。

提升带宽：多存储体并行

仅靠突发模式还不足以充分利用内存接口的潜力。为了进一步提升带宽，现代DRAM芯片采用了多存储体结构。

多个存储体可以并行工作。当一个存储体正在进行核心阵列访问（耗时较长）时，另一个存储体可以利用接口传输数据。通过交错不同存储体的数据传输，可以显著减少接口的空闲时间，从而成倍增加有效带宽。

在实际设计中，一个DRAM芯片可能包含16、32甚至上百个存储体，以进一步隐藏核心访问延迟，实现接近接口峰值带宽的持续数据传输。

GPU编程实践：内存访问合并

我们已经了解了DRAM如何高效地提供数据。现在，从GPU编程的角度来看，我们需要确保线程的内存访问模式能够充分利用这种突发传输特性，这就是内存访问合并。

在C/C++中，多维数组在内存中是按行线性存储的。例如，一个4x4的矩阵，在内存中的顺序是第一行的四个元素，接着是第二行的四个元素，依此类推。

考虑一个简单的矩阵乘法内核（非分块版本）：

float Pvalue = 0;
for (int k = 0; k < Width; ++k) {
    Pvalue += M[row * Width + k] * N[k * Width + col];
}
P[row * Width + col] = Pvalue;

让我们分析线程的内存访问模式：

  • 对于矩阵N的访问：同一线程块中，相邻的线程（threadIdx.x相邻）在循环的同一迭代中，会访问N矩阵中相邻的列元素。由于数组按行存储，这些列元素在内存地址上是连续的。因此，这些访问可以被合并——GPU可以一次性从DRAM获取一个包含这些连续地址数据的突发块，然后分发给各个线程，从而高效利用带宽。
  • 对于矩阵M的访问：同一线程块中，相邻的线程（threadIdx.x相邻）在循环的同一迭代中，会访问M矩阵中相同的行元素（因为row相同）。这看起来似乎很好，可以共享数据。但在非分块内核中，每个线程独立进行全局内存访问，实际上每个线程都会发起一次对该相同地址的加载请求。更糟糕的是，如果线程在threadIdx.y维度上相邻，它们会访问不同行的元素，这些元素在内存中相隔Width个元素，访问模式是分散的，无法合并，导致内存带宽利用率极低。

因此，在简单的矩阵乘法内核中，对M的访问模式是非合并的，这是其性能低下的一个重要原因。

解决方案：通过共享内存实现“转角”

在上一讲的分块矩阵乘法中，我们引入共享内存不仅是为了数据复用，也巧妙地解决了内存访问合并的问题。

在从全局内存加载数据块到共享内存时，我们让线程块中的所有线程协作加载整个数据块。具体操作如下：

  1. 合并访问全局内存：每个线程负责从全局内存中读取数据块的一个元素到共享内存。我们安排线程的读取索引，使得相邻线程读取的是全局内存中地址连续的元素。这样，对全局内存的访问是完美合并的，无论是加载M的数据块还是N的数据块。
  2. 在共享内存中重新组织数据：数据被加载到共享内存后，线程可以根据计算的需要，以任何模式（包括之前非合并的模式）从共享内存中高效地读取数据，因为共享内存位于芯片内部，延迟极低。

这种技术有时被称为“转角”，因为它逻辑上改变了数据的访问方向。我们通过在共享内存中暂存数据，将对全局内存的非合并访问，转换成了合并访问，从而显著提升了内存带宽的利用率。

判断合并访问的简单经验：在访问全局内存时，确保线程的访问地址偏移量与threadIdx.x呈简单的线性关系（例如 base_address + threadIdx.x），这样相邻线程访问相邻地址，有利于合并。应避免访问地址与threadIdx.y相关或带有较大乘数的复杂偏移。

常量内存与缓存

在卷积等算法中，所有线程都会重复读取相同的、不变的数据（如卷积核或滤波器）。GPU为此提供了常量内存。常量内存位于芯片上，并通过一个缓存机制进行加速。

缓存是一小块高速内存，用于存储最近从较慢内存（如全局内存、常量内存）中读取的数据副本。它利用了程序的两种局部性：

  1. 时间局部性：最近被访问的数据很可能在短期内再次被访问。
  2. 空间局部性：当某个内存位置被访问时，其附近的内存位置很可能在短期内也被访问。

当线程请求的数据已经在缓存中（称为“命中”），则可以直接从快速缓存中提供，速度远快于访问全局内存。对于只读的常量数据，GPU的常量缓存设计可以更简单高效，提供接近共享内存的访问速度（约5个周期）。

缓存与共享内存（暂存器）的区别：

  • 缓存：由硬件自动管理。硬件决定哪些数据保留在缓存中，程序员通常无法直接控制。
  • 共享内存：由程序员显式管理。程序员决定何时、何地、如何将数据放入共享内存以及如何使用它，提供了更灵活和可预测的性能优化手段。

在卷积示例中，将卷积核（掩膜）数据放入常量内存，可以利用常量缓存让所有线程快速共享这些只读数据。

总结

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

  1. DRAM的工作原理：基于电容存储，需要定期刷新，通过核心阵列、行/列解码器进行访问。
  2. 提升DRAM带宽的技术：突发传输模式利用单次行访问获取连续数据；多存储体并行通过交错多个存储体的操作来隐藏延迟，两者结合以实现高持续带宽。
  3. 内存访问合并：GPU编程的关键优化技术。线程对全局内存的访问应遵循“相邻线程访问相邻地址”的模式，以匹配DRAM的突发传输特性，从而最大化有效内存带宽。
  4. 利用共享内存实现合并：通过先将数据以合并模式从全局内存加载到共享内存，再在共享内存中进行任意模式的计算，可以解决某些算法中固有的非合并访问问题。
  5. 常量内存与缓存：对于所有线程只读的常量数据，使用常量内存可以利用GPU上的常量缓存，实现高速数据共享。

理解这些内存层次结构的特点和优化方法，对于编写高性能的GPU程序至关重要。

006：卷积与常量内存

在本节课中，我们将学习卷积计算模式，并探讨如何利用CUDA的常量内存来优化卷积操作。卷积是高性能计算、信号处理以及深度学习中的核心操作。我们将从基础概念开始，逐步深入到使用共享内存进行优化的分块卷积算法。

卷积基础回顾

上一节我们介绍了卷积的基本概念。卷积涉及两个数据集：输入数据集 N 和一个通常较小的掩码（或滤波器）M。输出数据集 P 的每个元素都是通过将掩码与输入数据的对应区域进行点积运算得到的。

对于一个一维卷积，计算输出元素 P[i] 的公式如下：

P[i] = Σ (M[j] * N[i + j - radius])， 其中 j = 0 到 mask_width-1

这里，radius = mask_width / 2（向下取整）。当掩码超出输入边界时，我们将其视为与零相乘。

常量内存的使用

在卷积计算中，掩码 M 对于所有输出线程都是只读且恒定的。这使得它成为使用常量内存的理想候选。常量内存通过缓存机制提供快速访问，特别适合被所有线程频繁读取的少量数据。

以下是使用常量内存的关键步骤：

  1. 声明常量内存：在全局作用域使用 __constant__ 限定符声明掩码。

    __constant__ float M_c[MASK_WIDTH][MASK_WIDTH];

  2. 从主机复制数据：使用 cudaMemcpyToSymbol 函数将主机上的掩码数据复制到设备的常量内存中。

    cudaMemcpyToSymbol(M_c, host_mask, MASK_WIDTH * MASK_WIDTH * sizeof(float));

  3. 在核函数中访问：由于常量内存变量具有文件作用域，核函数可以直接通过变量名 M_c 访问它，而无需作为参数传入。

使用常量内存可以减少对全局内存的访问压力，并利用缓存提升性能。

优化策略：数据复用与分块

基础的卷积核函数是内存受限的，因为每个输出元素都需要从全局内存加载整个掩码和对应的输入数据。为了提升性能，我们需要复用从全局内存加载的数据。

数据复用分析

考虑一个使用 5x5 掩码的二维卷积。对于输入数据中间的一个元素，它会被多少个不同的输出元素需要？答案是 25 个。因为掩码覆盖的 25 个位置在计算不同的输出时，都会用到这个中心输入元素。这种输入数据被多个输出复用的特性，是我们进行性能优化的基础。

分块卷积策略

为了利用数据复用，我们引入共享内存，采用分块计算策略。核心思想是将输入数据块加载到共享内存中，让块内的多个线程共享这些数据，从而减少对全局内存的访问。

在实现分块卷积时，一个关键决策是如何分配线程的工作。主要有两种并行化策略：

  1. 输出并行化：每个线程负责计算一个输出元素。这是最直观的方式。
  2. 输入并行化：每个线程负责将一个输入元素加载到共享内存中。然后，由块内的一部分线程（通常是前 tile_width 个线程）负责计算输出。

策略比较与选择

以下是两种策略的简要比较：

  • 输出并行化：
    • 优点：逻辑直观，所有计算线程都保持活跃。
    • 挑战：加载输入数据到共享内存时，由于需要加载额外的“光晕”区域，可能导致部分线程需要加载多个元素，引起分支分歧和潜在的合并访问问题。
  • 输入并行化：
    • 优点：加载阶段所有线程工作量均衡，易于实现高效的合并内存访问。
    • 挑战：只有一部分线程参与计算，会造成一定的计算资源浪费（部分线程闲置）。但随着掩码增大，闲置比例会降低。

对于初学者，输入并行化策略通常更易于实现和理解，也是后续实验推荐的方法。它通过清晰的职责分离（所有线程加载，部分线程计算）简化了边界处理和共享内存的索引计算。

输入并行化分块卷积示例

以下以二维卷积为例，概述输入并行化策略的关键步骤：

  1. 定义变量：

    • tile_width：输出块的大小（例如 16x16）。
    • block_size：线程块的大小。因为需要加载光晕区域，所以 block_size = tile_width + mask_width - 1（例如 5x5 掩码对应 20x20）。
    • 确保 block_size * block_size 不超过 GPU 的每块线程数限制（通常是 1024）。

  2. 计算网格和块尺寸：

    • 网格尺寸基于覆盖整个输出所需的 tile_width 大小块的数量。
    • 块尺寸设置为 block_size。

  3. 加载数据到共享内存：

    • 每个线程根据其线程ID (tx, ty) 计算对应的输入坐标。
    • 进行边界检查，如果坐标有效，则从全局内存加载输入值到共享内存数组 tile 中；否则，将 tile 中对应位置置为 0。
    • 使用 __syncthreads() 确保所有加载完成。

  4. 计算输出：

    • 只有满足 (tx < tile_width) && (ty < tile_width) 的线程参与计算。
    • 这些线程使用共享内存中的 tile 数据和常量内存中的掩码 M_c 进行卷积计算。
    • 计算完成后，再次检查输出坐标是否有效，然后将结果写回全局内存。

这种策略将复杂的边界处理和加载逻辑与计算逻辑分离，使得核函数结构更清晰。

总结

本节课我们一起学习了卷积计算及其在CUDA上的优化。

  • 我们回顾了卷积的基本原理和公式。
  • 介绍了常量内存的特性与使用方法，它适合存储像卷积掩码这样只读且恒定的数据。
  • 深入分析了卷积中的数据复用特性，这是性能优化的关键。
  • 探讨了基于共享内存的分块卷积策略，比较了输出并行化和输入并行化两种方法。
  • 详细阐述了输入并行化策略的实现步骤，这是推荐初学者使用的方案，它通过让所有线程加载数据、部分线程计算的方式，平衡了实现复杂度和性能。

理解这些概念和策略，将为你在后续实验中实现高效的三维分块卷积打下坚实的基础。

007：卷积性能分析与归约模式入门 🚀

概述

在本节课中，我们将深入分析卷积算法的性能，并引入一种新的并行计算模式——归约。我们将探讨如何通过分块策略提升卷积的内存重用率，并学习归约模式的基本概念、并行实现方法及其在GPU上的性能考量。

卷积性能分析回顾 🔍

上一节我们介绍了三种并行化卷积的方法：输出并行、带共享内存的输入并行以及纯输入并行。本节中，我们来看看如何量化这些方法的性能提升，特别是内存带宽的减少。

内存重用率分析

卷积算法的核心优化在于减少对全局内存的访问。通过将输入数据块加载到共享内存，多个输出线程可以重用这些数据。

一维卷积重用公式

对于一个分块大小为 T、掩码宽度为 M 的一维卷积：

  • 从全局内存加载到共享内存的元素数量：T + M - 1
  • 被替换的全局内存访问次数（即共享内存访问次数）：T * M
  • 平均重用率（带宽减少因子）为：
    重用率 = (T * M) / (T + M - 1)

二维卷积重用公式

将一维公式扩展到二维，假设输出块大小为 T x T，掩码大小为 M x M：

  • 加载到共享内存的元素数量：(T + M - 1)^2
  • 被替换的全局内存访问次数：T^2 * M^2
  • 平均重用率为：
    重用率 = (T^2 * M^2) / ((T + M - 1)^2)
    这本质上是一维重用率的平方。

边界效应

上述分析主要针对内部数据块。对于边界块，由于部分“光晕”区域位于输入数据之外，实际加载的数据量更少，因此重用率会略高。但在大数据集下，边界块占比很小，其影响可以忽略。

线程束分化与控制流

在输入并行方法中，并非所有线程都参与计算，这导致了线程束分化。例如，在一个 32x32 的线程块中处理 30x30 的输出块时：

  • 加载阶段：有30个线程束（warp）存在控制流分化（部分线程加载数据，部分线程写入0）。
  • 计算阶段：有31个线程束存在控制流分化（部分线程计算，部分线程不计算）。

虽然存在分化，但由于是简单的 if 条件（无 else 分支），性能损失相对较小。

现代GPU的算力与带宽挑战 ⚡

随着GPU发展，计算能力的提升速度远快于内存带宽，这使得数据重用变得至关重要。

案例分析

我们通过“字节/浮点操作”比率来分析：

  • 未分块卷积：每个输出元素需要加载一个输入元素（4字节）并执行两次浮点操作（乘和加），因此比率为 2 字节/FLOP。
  • GPU算力需求：要达到峰值算力利用率，所需的“重用因子”可以通过 (峰值算力 / 基于带宽的算力) 计算。

以下是不同年代GPU的需求变化：

  1. 2010年GPU：需要约 13倍 重用才能充分利用算力。
  2. 2020年GPU：需要约 52倍 重用。
  3. 2023年H100 GPU：得益于高带宽内存（HBM），需求降至约 26倍 重用，但仍高于2010年水平。

对卷积算法的启示

  • 对于一维卷积，即使使用较大的掩码（如宽度55），也难以满足现代GPU的高重用需求。
  • 对于二维卷积，由于重用潜力是掩码大小的平方，一个 9x9 的掩码就能提供81倍的潜在重用，更容易接近饱和GPU算力。
  • 线程协作是进一步提升重用率的有效技术，通过让每个线程处理多个输出元素，可以增加数据块大小，从而在共享内存容量允许范围内获得更高的重用。

并行计算模式：归约 🌳

现在，我们转向一种新的、广泛使用的并行计算模式——归约。

什么是归约？

归约是将一组输入值通过一个二元运算符合并为单个值的过程。

  • 运算符：必须是结合律和交换律的，例如加法、乘法、最大值、最小值。
  • 恒等值：运算符的恒等元（如加法为0，乘法为1）。
  • 应用场景：求总和、求极值、在MapReduce框架中汇总结果等。

并行归约算法

并行归约通常采用树形结构：

  1. 第一步：N 个元素两两合并，产生 N/2 个部分结果。
  2. 后续步骤：重复合并，每一步活动元素数量减半。
  3. 经过 log₂(N) 步后，得到最终结果。

该算法是工作高效的，总操作次数为 N-1，与串行算法相同。

在GPU上实现归约的挑战

尽管树形归约在理论上高效，但在GPU上直接实现会面临并行度递减的问题：

  • 初始阶段：需要大量线程（如 N/2），可能超过GPU硬件限制。
  • 后续阶段：活动线程数迅速减少，导致硬件利用率下降，最后几步几乎是串行的。
  • CUDA模型限制：线程数量在核函数启动时固定，无法动态减少。

分块归约策略

为了解决上述问题，我们采用分块策略：

  1. 加载到共享内存：每个线程块（包含 M 个线程）从全局内存加载 2M 个元素到共享内存。
  2. 块内归约：在共享内存中，使用树形归约将 2M 个元素合并为1个值。这需要 log₂(2M) 步。
  3. 写回全局内存：每个线程块将其结果写回全局内存。如果还有多个块，可以启动另一个归约核函数或将结果传回CPU处理。

块内归约代码示例（以加法为例）

以下是核心归约循环的伪代码，假设 partial_sum 是共享内存数组，已加载了 2*blockDim.x 个数据：

__shared__ float partial_sum[2*BLOCK_SIZE];
unsigned int t = threadIdx.x;
unsigned int index = 2 * t; // 每个线程负责的起始索引

for (unsigned int stride = 1; stride <= blockDim.x; stride *= 2) {
    // 只有满足条件的线程参与本步计算
    if (t % stride == 0) {
        partial_sum[index] += partial_sum[index + stride];
    }
    __syncthreads(); // 确保本步所有写入完成
}
// 最终结果存储在 partial_sum[0] 中
if (t == 0) {
    output[blockIdx.x] = partial_sum[0]; // 写回全局内存
}

关键点：

  • stride 变量控制合并的距离，每一步翻倍。
  • t % stride == 0 条件用于筛选每一步中活跃的线程。
  • __syncthreads() 屏障至关重要，用于确保前一步的所有写入对下一步的读取可见。
  • 最后，仅用线程0将结果写回全局内存，无需额外的 __syncthreads()，因为该线程参与了最后一步计算，数据对其立即可见。

总结 🎯

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

  1. 卷积性能量化：我们推导了一维和二维卷积在分块策略下的内存重用率公式，并讨论了边界效应和线程束分化的影响。
  2. 算力与带宽的平衡：通过案例分析，我们理解了现代GPU对数据重用日益增长的需求，以及二维卷积和线程协作技术在应对这一挑战中的优势。
  3. 归约模式入门：我们学习了归约的基本概念、并行树形算法及其工作高效性。
  4. GPU归约实现：我们探讨了在GPU上实现归约时面临的并行度递减挑战，并介绍了通过分块和共享内存来解决这一问题的策略，给出了核心代码示例。

归约是一种基础且强大的并行模式，在后续课程和许多实际应用中都会频繁出现。理解其原理和实现方法对编写高效的并行程序至关重要。

008：归约与扫描算法详解

在本节课中，我们将深入学习并行计算中两个至关重要的模式：归约和扫描。我们将从回顾归约算法开始，探讨其优化方法，然后深入讲解扫描算法的原理与实现，特别是两种主流的并行扫描算法：Kogge-Stone算法和Brent-Kung算法。

归约算法回顾

上一节我们介绍了归约模式的基本概念。归约是指将一组输入值通过一个二元运算符（如求和、求最大值）合并为单个值的过程。该运算符必须是结合律和交换律的。

在GPU上实现归约时，我们通常利用共享内存。对于一个包含M个线程的线程块，我们从全局内存读取2M个值到共享内存，并确保良好的内存合并访问。

归约的核心步骤是一个循环，其跨度从1开始，每次迭代翻倍，直到达到线程块维度。在每一步中，每个活跃线程将其对应的两个元素进行归约操作，并将结果写回共享内存。每次迭代后，活跃线程数减半。

以下是归约循环的核心代码逻辑：

for (int stride = 1; stride <= blockDim.x; stride *= 2) {
    __syncthreads(); // 确保上一步的写入完成
    if (threadIdx.x % (2 * stride) == 0) {
        shared_array[threadIdx.x] = shared_array[threadIdx.x] + shared_array[threadIdx.x + stride];
    }
}

__syncthreads() 在此至关重要，它确保所有线程在读取数据前，上一步的写入操作已经完成。缺少同步会导致数据竞争和错误结果。

归约的线程组织优化

我们之前讨论的归约方法存在资源利用率低的问题。在后续步骤中，大量线程处于空闲状态，甚至整个线程束（Warp）可能只有一个线程在工作。

本节我们来看一种更优的线程组织方法。其核心思想是紧凑存储中间结果，使得活跃线程始终是连续的，从而更好地利用线程束。

在优化版本中，每个线程不再处理固定间隔的两个元素，而是处理两个相距较远的元素（初始跨度为线程数的一半），并将结果紧凑地存储在共享数组的前半部分。这样，活跃线程在每一步都是连续的。

以下是优化后归约循环的逻辑变化：

for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
    __syncthreads();
    if (threadIdx.x < stride) {
        shared_array[threadIdx.x] += shared_array[threadIdx.x + stride];
    }
}

这种方法在前几步避免了线程束分化，所有活跃线程都位于连续的线程束中，提高了计算资源的利用率。

并行算法开销

即使并行算法是“工作高效”的（即执行的操作总数与串行算法相同），它仍然存在并行开销。这些开销源于并行化本身引入的额外工作。

以下是并行归约中主要的开销来源：

  • 线程索引计算：每个线程需要计算自己在网格和块中的位置。
  • 数据移动：将数据从全局内存加载到共享内存，这本身不参与计算。
  • 线程同步：使用 __syncthreads() 等屏障带来的开销。
  • 控制流：循环和条件判断（如 if 语句）在并行代码中更为复杂。

理解这些开销有助于我们设计和优化更高效的并行算法。

扫描算法介绍

现在，让我们转向扫描算法。扫描，也称为前缀和，是归约的扩展。归约产生一个最终结果，而扫描产生一个数组，其中每个元素是输入数组从开始到该位置所有元素的归约结果。

对于一个输入数组 [x0, x1, ..., xn-1] 和满足结合律的二元运算符 ⊕，其包含性扫描的输出数组 [y0, y1, ..., yn-1] 定义为：

  • y0 = x0
  • y1 = x0 ⊕ x1
  • y2 = x0 ⊕ x1 ⊕ x2
  • ...
  • yn-1 = x0 ⊕ x1 ⊕ ... ⊕ xn-1

扫描在并行计算中用途广泛，例如流压缩、排序、分配资源等。

Kogge-Stone扫描算法

Kogge-Stone算法是一种经典的并行前缀扫描算法，灵感来源于硬件加法器设计。它通过多次迭代计算前缀和。

该算法的核心步骤是：

  1. 将数据从全局内存加载到共享内存。
  2. 进行 log2(n) 次迭代，跨度 stride 从1开始，每次翻倍。
  3. 在每次迭代中，对于线程索引 j >= stride 的活跃线程，执行操作：dest[j] = src[j] ⊕ src[j - stride]。
  4. 迭代完成后，共享内存中即包含扫描结果。

直接实现会导致读写冲突，因为我们在读取 src 的同时可能正在向同一个数组写入。为了解决这个问题，我们引入双缓冲技术。

双缓冲技术

双缓冲通过维护两个共享内存数组（例如 buffer0 和 buffer1）来消除读写依赖。在每次迭代中，我们从源缓冲区读取，向目标缓冲区写入，然后交换两个缓冲区的指针角色。

以下是使用双缓冲的Kogge-Stone扫描循环伪代码：

src = buffer0;
dst = buffer1;
for (int stride = 1; stride < n; stride *= 2) {
    __syncthreads(); // 等待所有数据就绪（包括初始加载）
    if (threadIdx.x >= stride) {
        dst[threadIdx.x] = src[threadIdx.x] ⊕ src[threadIdx.x - stride];
    } else {
        // 非活跃线程需要复制其值到目标缓冲区，以保持数据完整
        dst[threadIdx.x] = src[threadIdx.x];
    }
    // 交换源和目标指针
    temp = src; src = dst; dst = temp;
}
__syncthreads();
// 最终结果在 `src` 指向的缓冲区中

双缓冲允许我们每迭代仅使用一次 __syncthreads()，提高了性能。

Brent-Kung扫描算法

Kogge-Stone算法虽然延迟低，但其操作数量为 O(n log n)，并非工作高效的（串行扫描仅需 O(n) 次操作）。Brent-Kung算法提供了另一种选择，它分两阶段进行，总操作数量更少，约为 2n，是工作高效的。

Brent-Kung算法的两个阶段：

  1. 归约阶段（上行）：构建一棵平衡二叉树，从叶子节点向根节点计算部分和。这一步与优化后的归约类似，最终根节点得到全局和。
  2. 分发阶段（下行）：从根节点开始，向下遍历树，将必要的部分和值加到相应的节点上，从而为每个位置生成最终的前缀和。

Brent-Kung算法用更多的步骤（2 log n - 2 步）换取了更少的总体操作，在计算资源受限或问题规模较大时可能更有优势。

总结

本节课我们一起深入学习了并行计算中的归约与扫描。

  • 我们回顾了基础的归约算法，并探讨了通过改变线程组织方式来优化资源利用率的技巧。
  • 我们分析了并行算法固有的开销来源。
  • 我们引入了扫描（前缀和）的概念及其广泛应用。
  • 我们详细讲解了Kogge-Stone扫描算法的原理，并利用双缓冲技术解决了实现中的同步问题。
  • 最后，我们简要介绍了工作效率更高的Brent-Kung扫描算法，它通过上行-下行两阶段过程来减少总操作数。

理解这些核心算法模式对于设计高效GPU程序至关重要。在接下来的实验中，你将有机会亲手实现这些算法。

009：并行扫描算法与原子操作

在本节课中，我们将要学习两种并行扫描算法（前缀和）的对比，并深入探讨原子操作的概念及其在并行编程中的应用。我们将从回顾上周的Coggystone算法开始，然后介绍Brent-Kung算法，分析两者在延迟和工作效率上的权衡。最后，我们将引入原子操作，解释其工作原理，并探讨其在直方图计算等应用中的使用。

回顾：并行扫描与Coggystone算法

上一节我们介绍了并行扫描（前缀和）的概念以及Coggystone算法。并行扫描的定义如下：给定一个包含 n 个元素的数组 x[0] 到 x[n-1]，返回一个前缀和数组，其中每个输出元素 y[i] 是输入数组前 i+1 个元素在某个二元结合运算符（如加法）下的累积结果。

公式：
y[i] = x[0] op x[1] op ... op x[i]

Coggystone算法在 log₂(n) 次并行迭代中完成计算，但其总工作量是 O(n log n)，而顺序算法仅需 O(n) 次操作。因此，Coggystone算法的工作效率较低，但其延迟较小，因为并行步骤少。

Brent-Kung算法：基于平衡二叉树

本节中，我们来看看另一种基于平衡二叉树思想的并行扫描算法——Brent-Kung算法。该算法旨在提高工作效率，使其接近顺序算法的 O(n)，但代价是增加了并行步骤（即同步点），从而可能增加延迟。

Brent-Kung算法分为两个阶段：

  1. 归约阶段（Fan-in）：自底向上构建一棵平衡二叉树，将部分和汇聚到树的“根”节点（实际上是数组的最后一个元素）。
  2. 分布阶段（Fan-out）：自顶向下遍历树，将累积的部分和分发出去，以计算最终的前缀和结果。

该算法在共享内存数组 T 上操作，其大小为线程块大小的两倍。

归约阶段内核代码分析

以下是归约阶段的核心循环代码逻辑：

for (int stride = 1; stride < 2 * blockDim.x; stride *= 2) {
    __syncthreads(); // 确保数据已加载到共享内存
    int index = (threadIdx.x + 1) * stride * 2 - 1;
    if (index < 2 * blockDim.x && (index - stride) >= 0) {
        T[index] += T[index - stride];
    }
}

代码解释：

  • stride 从1开始，每次迭代翻倍，直到达到 2 * blockDim.x。
  • index 的计算确保了线程操作不重叠，避免了读写冲突，因此无需像Coggystone算法那样将读写步骤分开。
  • 每个活动线程从 T[index - stride] 读取数据，加到 T[index] 上，并写回。

分布阶段内核代码分析

归约阶段完成后，我们得到了部分和。分布阶段则负责将这些部分和组合成最终的前缀和：

for (int stride = blockDim.x / 2; stride >= 1; stride /= 2) {
    __syncthreads(); // 技术上在第一次迭代可能不需要，但为安全起见保留
    int index = (threadIdx.x + 1) * stride * 2 - 1;
    if ((index + stride) < 2 * blockDim.x) {
        T[index + stride] += T[index];
    }
}

代码解释：

  • stride 从 blockDim.x / 2 开始，每次迭代减半，直到1。
  • 此阶段线程将 T[index] 的值加到 T[index + stride] 上，方向与归约阶段相反。
  • 经过此阶段，共享内存数组 T 中包含了完整的前缀和结果。

算法对比与权衡

以下是Coggystone算法与Brent-Kung算法的关键对比：

  • 工作量：Coggystone为 O(n log n)，Brent-Kung为 O(n)（更高效）。
  • 步骤数（延迟）：Coggystone约需 log₂(n) 步，Brent-Kung约需 2 * log₂(n) 步（更多同步，延迟可能更高）。
  • 线程利用率：对于处理相同数量的元素，Brent-Kung所需的线程数约为Coggystone的一半（因为每个Brent-Kung线程加载2个元素）。
  • 适用场景：由于GPU计算能力强，通常更青睐延迟低的Coggystone算法。只有当运算符本身非常复杂耗时，需要极力减少计算次数时，才会考虑使用Brent-Kung算法。

分层扫描处理大规模数据

无论是Coggystone还是Brent-Kung，单个线程块能处理的数据量是有限的（例如1024或2048个元素）。为了扫描数十亿元素的大数组，我们需要采用分层方法：

  1. 第一层内核：每个线程块扫描自己负责的一块输入数据，生成该块内的局部前缀和，并将该块所有元素归约后的块总和写入一个辅助数组。
  2. 第二层内核：对存储块总和的辅助数组执行扫描。这样，辅助数组的第 i 个元素就包含了前 i 个输入块的总和。
  3. 第三层内核：另一个内核启动，每个线程块（除了第一个）从辅助数组中读取对应的累积块总和，并将其加到本线程块之前计算好的局部前缀和的每个元素上。第一个块的结果已经是最终的全局前缀和，无需调整。

关键点：必须通过完成内核并启动新内核的方式，来确保一个内核对全局内存的写入对后续内核可见。试图在单个内核内通过线程块间同步来实现此目的是危险且容易导致死锁的。

包含扫描与排除扫描

之前讨论的都是包含扫描，即输出 y[i] 包含输入 x[i]。
排除扫描的定义是：输出 y[i] 包含的是前 i 个输入 x[0] 到 x[i-1] 的累积结果，y[0] 是运算符的单位元（如加法中的0）。

公式（排除扫描）：
y[0] = identity
y[i] = x[0] op x[1] op ... op x[i-1] (对于 i > 0)

两者算法本质相同，只是输出位置有偏移。在实践中，可以通过在包含扫描结果前补单位元，或直接调整算法中的写入索引来相互转换。

原子操作基础

现在，我们转向并行编程中的另一个核心概念——原子操作。当多个线程可能写入同一内存位置时，我们需要原子操作来保证结果的正确性。

原子操作将一个“读-修改-写”序列作为一个不可分割的单元执行。相对于其他针对同一地址的原子操作，它们保证不会出现执行交错，从而避免数据竞争。

为什么需要原子操作？

考虑两个线程都要对全局变量 X（初始为0）执行 X++ 操作。如果不使用原子操作，可能会发生以下交错执行：

  1. 线程A读取 X=0。
  2. 线程B读取 X=0。
  3. 线程A计算 0+1=1，写入 X=1。
  4. 线程B计算 0+1=1，写入 X=1。
    最终 X=1，而不是正确的 2。这就是数据竞争。

使用原子操作（如 atomicAdd）可以保证：无论线程A和B谁先执行，最终 X 都会等于 2。原子性保证了操作的顺序执行（不交错），但不规定线程A和B的执行先后顺序。

CUDA中的原子操作

CUDA提供了多种原子操作 intrinsic 函数，例如：

  • atomicAdd(int* address, int val): 原子性地将 val 加到 address 指向的值上，返回旧值。
  • atomicSub, atomicExch, atomicMin, atomicMax, atomicCAS (比较并交换) 等。

重要提示：原子操作是序列化的，对同一地址的多次原子操作会排队执行，这可能成为性能瓶颈，因此应谨慎使用。

使用原子操作实现同步的陷阱

理论上，可以用原子操作和全局变量实现跨线程块的全局同步（类似 __syncthreads 但针对整个网格）。但强烈不建议这样做。原因在于，如果启动的线程块数量超过GPU能同时容纳的数量，先执行的线程块会在同步点等待尚未被调度的线程块，而后者因为GPU资源已满无法启动，从而导致死锁。

直方图计算：原子操作的应用实例

直方图是一种常见的数据分析操作，它将大量数据点归类到有限的“桶”中，并统计每个桶的频数。

简单的并行直方图算法

一个直观的并行化方法是：为每个输入数据元素分配一个线程。每个线程根据其数据值，找到对应的直方图桶，然后使用原子操作将该桶的计数值加1。

以下是一个简单直方图内核的伪代码思路：

__global__ void histogram_simple(unsigned char* data, int data_size, unsigned int* histogram) {
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;

    for (int i = tid; i < data_size; i += stride) {
        unsigned char bin = data[i]; // 确定桶索引 (0-255)
        atomicAdd(&histogram[bin], 1); // 原子递增
    }
}

算法分析：

  • 优点：实现简单，并行度高。
  • 缺点：
    1. 原子操作串行化：所有线程都可能竞争更新少数几个桶（特别是高频值），导致严重的序列化瓶颈。
    2. 内存访问模式：对 histogram 数组的访问是完全随机的（取决于输入数据），无法合并，导致全局内存访问效率低下。
    3. 输入读取：如果像最初例子中那样给每个线程分配连续的数据块，则对输入 data 的读取也无法合并。上面改进后的伪代码使用了跨步循环 (i += stride)，使得同一warp的线程读取相邻内存地址，实现了对输入数据的合并访问。

性能瓶颈与优化方向

简单原子直方图的主要瓶颈在于原子操作的竞争。优化策略的核心是减少原子操作的冲突。

一种关键技术是私有化：

  1. 每个线程（或线程块）先使用自己的局部直方图（在寄存器或共享内存中）进行统计。
  2. 局部直方图只被单个线程（或线程块内的线程）访问，因此无需原子操作，速度极快。
  3. 在所有局部统计完成后，再将局部直方图的结果累加到全局直方图中。这一步仍然需要原子操作，但冲突次数大大减少，因为每个线程（或线程块）只向每个桶原子加一次（其局部计数值），而不是为每个数据元素都加一次。

我们将在下一节课中详细探讨私有化等直方图优化技术。

总结

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

  1. Brent-Kung扫描算法：一种基于平衡二叉树、工作效率更高（O(n)）但步骤更多（2*log₂(n)）的并行前缀和算法。我们分析了其归约和分布两个阶段的代码逻辑。
  2. 分层扫描：通过多内核协作，将局部扫描结果与块总和的扫描结果相结合，以处理远超单个线程块容量的大规模数据。
  3. 原子操作：解释了原子操作作为“读-修改-写”的不可分割单元，如何用于解决数据竞争，并强调了其序列化特性可能带来的性能问题。
  4. 直方图计算：介绍了使用原子操作的简单并行直方图算法，并指出了其存在的原子竞争和内存访问效率问题，为后续的优化讲座做了铺垫。

理解这些算法和概念的权衡，对于在GPU上设计高效、正确的并行程序至关重要。

010：原子操作与稀疏矩阵向量乘法

概述

在本节课中，我们将继续探讨原子操作，理解其硬件实现原理和性能影响，并将其应用于直方图计算。随后，我们将转向并行稀疏矩阵向量乘法，学习如何高效处理现实世界中大量存在的稀疏数据。

原子操作回顾与硬件实现

上一节我们介绍了原子操作的基本概念。本节中我们来看看原子操作在硬件层面是如何实现的，以及这为何会对性能产生影响。

原子操作，特别是我们用于直方图的读-修改-写操作，需要从内存中读取一个值，在处理器内部修改它，然后再写回内存。

当一个流多处理器执行原子操作时，它需要访问DRAM。来自其他SM的对同一全局地址的原子操作必须被协调，确保它们相对于彼此是序列化的。在整个读-修改-写周期内，其他SM或同一SM内的其他线程都无法访问该特定内存位置。

其性能影响在于，每个读-修改-写操作都包含两次完整的内存访问延迟。所有对同一内存位置的原子操作在时间上是串行化的。这意味着原子操作的吞吐量，即应用程序对同一位置连续执行原子操作的速率，将受到整个读-修改-写序列总延迟的限制。由于每次内存访问可能需要数百个周期，这严重限制了SM的有效利用率和内存带宽。

原子操作的性能优化：缓存与共享内存

为了缓解全局内存原子操作的性能瓶颈，后来的GPU架构引入了硬件改进，将原子操作移至更靠近处理器的缓存中执行。

其思想是将操作转移到最需要共享数据的处理器附近的缓存中。对于多个SM，这意味着芯片边缘的缓存。原子操作在芯片上的缓存中执行，延迟大幅降低至数十个周期，但访问同一位置的操作仍然是串行化的。

更好的方案是使用共享内存进行原子操作。共享内存距离SM仅几个周期，因此原子操作的吞吐量可以显著提高。然而，共享内存是线程块私有的，无法在线程块之间共享数据。

直方图计算的私有化技术

如果我们需要跨线程块的全局原子操作，但又想利用共享内存的高性能，可以通过算法改写来实现，这种方法称为私有化。

私有化依赖于我们并不关心操作顺序这一事实。对于直方图，每个桶的操作只是加一，我们只关心所有加一操作都发生，而不关心其顺序。

我们可以让每个线程块拥有直方图的一个私有副本。每个线程块使用共享内存原子操作高效地计算其数据部分。最后，所有线程块将其私有副本的结果通过（次数少得多的）全局内存原子操作累加到全局直方图中。

以下是实现步骤：

  1. 在共享内存中分配一个私有直方图副本。
  2. 使用至少256个线程（对应256个直方图桶）将私有副本初始化为零。
  3. 执行同步操作确保初始化完成。
  4. 每个线程处理其分配的数据部分，但使用共享内存原子操作更新私有直方图副本。
  5. 再次同步，确保所有线程完成计算。
  6. 最后，使用256个线程并行地将私有直方图每个桶的结果通过原子加操作累加到全局直方图的对应桶中。

私有化技术要求所执行的操作满足结合律和交换律。对于直方图加法，两者都满足。另一个约束是直方图必须足够小，以便其副本能够放入共享内存。

稀疏矩阵向量乘法简介

许多现实世界的系统和应用涉及稀疏数据，即矩阵中大部分元素为零。直接在包含大量零的稠密表示上进行计算会浪费计算资源和内存带宽。

并行稀疏矩阵向量乘法面临的关键挑战是其不规则性。与稠密矩阵乘法不同，它缺乏规整的数据访问模式和数据复用，难以从编译器优化中受益。

获得良好性能的关键在于：

  • 减少控制流分歧。
  • 平衡线程间的负载。
  • 通过优化数据布局实现良好的内存合并访问，以最大化内存带宽利用率。

稀疏矩阵存储格式

为了高效处理稀疏矩阵，人们开发了多种压缩存储格式。

CSR格式

压缩稀疏行格式是一种在CPU上流行的高效格式。

  • data 数组：按行优先顺序存储所有非零元素值。
  • col_index 数组：存储每个非零元素所在的列索引。
  • row_ptr 数组：存储每行非零元素在 data 和 col_index 数组中的起始位置。数组长度为行数加一，最后一个元素是非零元素总数。

在CSR格式上实现SpMV的核函数中，每个线程负责计算矩阵的一行与向量的点积。由于每行非零元素数量不同，会导致线程执行循环迭代次数不同（控制流分歧），并且对 data 和 col_index 数组的访问可能无法合并。

ELL格式

为了在GPU上获得更好的内存合并访问，提出了ELLPACK格式。

  1. 填充：找到矩阵中最长的行（非零元素最多），将所有其他行填充到此长度，填充值为零。
  2. 转置：将填充后的矩阵进行转置存储。

这样，在计算时，所有线程在每次迭代中访问 data 和 col_index 数组中相邻的位置，从而实现了内存访问的合并。由于所有行长度相同，不再需要 row_ptr 数组。缺点是当行长度差异很大时，填充会引入大量零元素，浪费存储和计算。

COO格式

坐标格式简单地为每个非零元素存储一个三元组：行索引、列索引和值。

  • 三个数组：row_index, col_index, data，长度均为非零元素总数。
  • 三元组的顺序可以任意排列，不影响结果。

COO格式易于并行化，可以为每个非零元素启动一个线程。每个线程读取一个三元组，计算 data[i] * x[col_index[i]]，然后需要将结果累加到 y[row_index[i]] 中。这会导致多个线程可能写入同一个 y 元素，因此需要使用原子操作，可能引发写冲突。

混合格式

结合ELL和COO的优点，形成混合格式。

  1. 选择一个适中的行长度进行ELL填充（例如，所有行的平均长度），将能放入ELL格式的非零元素进行填充和转置。
  2. 剩余无法放入ELL格式的非零元素，则使用COO格式存储。

这样，大部分计算通过高效的ELL格式完成，且合并访问良好；少部分剩余元素通过COO格式处理，减少了整体填充量。

JDS格式：通过排序实现负载均衡

另一种优化思路是通过对矩阵行进行排序来改善负载均衡。

  1. 排序：按每行非零元素数量降序排列矩阵行。非零元素最多的行排在最前面。
  2. 记录置换：存储一个 row_perm 数组，记录新行号对应的原始行号。
  3. 存储：对排序后的矩阵使用类似CSR的格式存储（有时称为JDS格式）。

这样，当将行分配给线程块时，一个线程块内的各行长度相近，减少了线程块内的负载不平衡。计算完成后，需要通过 row_perm 数组将结果写回 y 向量的正确位置。

总结

本节课我们一起学习了原子操作的硬件实现细节及其性能影响，并掌握了通过共享内存私有化技术来优化直方图计算。随后，我们转向了稀疏矩阵向量乘法，探讨了CSR、ELL、COO等多种稀疏矩阵存储格式及其在GPU上的实现策略，包括如何通过填充转置、混合格式以及行排序等技术来优化内存访问、减少控制流分歧和改善负载均衡，从而高效处理不规则的计算任务。

011：稀疏矩阵格式与机器学习入门

概述

在本节课中，我们将继续学习稀疏矩阵的存储格式，特别是JDS格式及其转置版本，并探讨这些格式在GPU上的性能影响。随后，我们将转向机器学习和深度学习的基础概念，了解其核心思想、应用领域以及多层感知机的基本结构。

稀疏矩阵向量乘法回顾

上一节我们介绍了压缩稀疏行格式、坐标格式和ELLPACK格式。本节中，我们来看看如何通过引入间接层和排序来进一步优化不规则数据的处理。

稀疏矩阵向量乘法的目标是计算 Y = A * X，其中矩阵A包含大量零元素。我们讨论过的格式旨在高效存储和访问这些非零元素。

压缩稀疏行格式

CSR格式包含三个数组：

  • data: 存储所有非零元素的值。
  • col_index: 存储每个非零元素对应的列索引。
  • row_ptr: 存储每行非零元素在data和col_index数组中的起始和结束位置（行数+1个元素）。

ELLPACK格式

ELLPACK格式通过填充零元素使所有行具有相同的长度（最大非零元数），然后进行转置，以实现线程间更好的内存合并访问。

坐标格式与混合格式

COO格式简单地存储每个非零元的行索引、列索引和值。我们可以将其与ELLPACK结合，选择一个小于实际最大行长度的阈值，将超出阈值的元素用COO格式存储，从而减少填充开销。

JDS格式及其转置

尽管使用了混合格式，我们仍然面临负载不均衡的问题，因为长行和短行可能被分到同一个线程块或线程束中。

JDS格式的核心思想

JDS格式通过引入一个间接层来解决负载不均衡问题。具体步骤如下：

  1. 对矩阵的行按照非零元素数量进行降序排序。
  2. 使用一个row_perm数组记录排序后行号与原行号的对应关系。
  3. 将排序后的行数据按CSR类似的方式存储。

这样，相邻线程处理的行具有相似的长度，提高了线程利用率，减少了线程束内部分线程空闲等待的时间。

JDS格式的内存访问问题

然而，简单的JDS格式重新引入了内存访问不连续的问题，因为每行的第一个元素在内存中并不相邻。

JDS转置格式

为了解决内存合并访问问题，我们可以对JDS格式的数据进行转置。与ELLPACK不同，JDS转置不需要填充零元素，因为我们已经通过排序解决了负载均衡问题。

以下是JDS转置格式的关键组件：

  • data_t: 转置后的非零值数组。
  • col_index_t: 转置后的列索引数组。
  • col_ptr: 列指针数组，指示转置后每一“列”（即原JDS排序后行数据的每一列）的起始位置。col_ptr[i+1] - col_ptr[i] 给出了在第 i 次循环迭代中活跃的线程数量。

JDS转置的代码逻辑

在核函数中，每个线程负责计算一个行（排序后的虚拟行）的点积。线程通过比较自己的虚拟行号与col_ptr数组的差值来决定在每次循环迭代中是否活跃。活跃的线程从data_t和col_index_t数组中连续地读取数据，实现了高效的内存合并访问。计算完成后，通过row_perm数组将结果写回输出向量Y的正确位置。

核心计算循环的伪代码逻辑如下：

float sum = 0.0f;
int sec = 0;
while (col_ptr[sec+1] - col_ptr[sec] > row) { // 如果该列活跃线程数大于当前行号，则该线程需参与计算
    int index = col_ptr[sec] + row;
    sum += data_t[index] * x[col_index_t[index]];
    sec++;
}
y[row_perm[row]] = sum; // 通过排列数组写回原行

不同稀疏矩阵格式的应用场景

以下是针对不同矩阵结构推荐使用的格式：

  • 随机矩阵/图：通常使用ELLPACK格式。由于度的分布相对均匀，填充带来的开销可以接受。
  • 行长度方差大的矩阵（如三角矩阵）：推荐使用JDS格式。通过排序有效平衡负载。
  • 带状矩阵：通常使用ELLPACK格式。各行非零元数量相近，转置后能获得良好的内存聚合访问。
  • 极度稀疏的矩阵：直接使用COO格式即可，数据量小，处理简单。
  • 具有规则非零结构的矩阵（如三对角矩阵）：可以使用更专用的格式（如对角线存储），将结构信息编码到数据结构中，无需存储索引。

此外，还有许多其他格式和混合格式，如压缩稀疏列、分块CSR等，可根据具体应用选择。

机器学习与深度学习引言

从稀疏矩阵计算过渡到新的领域，我们现在开始探讨机器学习和深度学习。这些技术利用计算能力从数据中学习模式和关系，解决了许多难以用传统算法形式化描述的问题。

机器学习定义

机器学习是一种构建应用程序的方法，其内部逻辑并非完全由人类预先定义。它通过从大量标注数据（输入-输出对）中学习，自动调整模型参数，从而对新的输入产生合理的输出。这个过程称为训练。

学习任务的类型

机器学习可以应用于多种任务：

  • 分类：将输入映射到预定义的类别（如图像识别、缺陷检测）。
  • 回归：预测连续值（如温度预测、曲线拟合）。
  • 转录：将非结构化数据转换为文本（如语音识别、OCR）。
  • 翻译：将一种语言的符号序列转换为另一种语言。
  • 异常检测：识别数据中的异常模式（如入侵检测）。

从特征工程到表示学习

传统机器学习依赖于特征工程，即由人类专家设计并提取对任务有用的数据特征（如图像的边缘、音频的频谱）。然而，对于许多复杂问题（如识别一辆汽车），定义合适的特征非常困难。

深度学习属于表示学习，它通过多层非线性处理单元（神经网络）自动学习数据的层次化特征表示。较低层学习简单特征（如边缘），较高层将这些特征组合成更复杂的抽象概念（如物体部件），最终用于决策。

多层感知机

一个基础的深度学习模型是多层感知机，它是神经网络的基本形式。

  • 感知机：一个简单的线性分类器，计算输入向量 x 和权重向量 w 的点积，加上偏置 b，然后通过一个激活函数（如符号函数）产生输出。公式为：output = sign(w·x + b)。它在多维空间中定义一个超平面。
  • 线性不可分问题：单个感知机无法解决非线性可分问题，如异或逻辑。
  • 多层感知机：通过堆叠多个感知机层可以解决非线性问题。一个典型的三层MLP包括：
    1. 输入层：接收原始数据。
    2. 隐藏层：包含多个感知机，学习数据的中间表示。其计算为：h = g(W1 * x + b1)，其中 g 是非线性激活函数（如ReLU）。
    3. 输出层：生成最终结果。对于分类任务，通常使用 softmax 函数将输出转换为概率分布。公式为：y = softmax(W2 * h + b2)。

理论上，一个足够大的单隐藏层MLP可以近似任何函数。但在实践中，更深的网络通常能更快、更有效地学习复杂的特征层次。

模型训练

模型参数（权重和偏置，统称为 θ）需要通过训练来学习。给定一个训练样本 (x, y)，过程如下：

  1. 前向传播：用当前参数 θ 和输入 x 计算网络输出 ŷ。
  2. 损失计算：使用损失函数 L(ŷ, y)（如均方误差、交叉熵）衡量预测输出 ŷ 与真实标签 y 之间的差异。
  3. 反向传播：计算损失函数相对于每个参数的梯度 ∂L/∂θ。
  4. 参数更新：使用优化算法（如梯度下降）沿梯度反方向更新参数：θ_new = θ_old - α * ∂L/∂θ，其中 α 是学习率。

通过在整个训练集上迭代这个过程，模型逐渐调整其参数以最小化总体损失。

总结

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

  1. JDS格式及其转置：通过按行非零元数排序解决负载不均衡，并通过转置数据解决内存访问不连续的问题，从而在GPU上实现高效的稀疏矩阵向量乘法。
  2. 机器学习基础：了解了机器学习通过数据驱动的方式解决难以形式化描述的问题，其核心是从示例中学习映射关系。
  3. 深度学习与表示学习：认识了深度学习作为表示学习的一种，能够自动学习数据的层次化特征，免除了复杂的手工特征工程。
  4. 多层感知机：学习了MLP的基本结构，包括输入层、隐藏层和输出层，以及它如何通过非线性组合解决复杂分类问题。我们还简要了解了模型训练的前向传播和反向传播过程。

下一讲，我们将更深入地探讨神经网络训练的具体细节，包括损失函数、梯度下降以及更复杂的网络结构。

012：多层感知机与卷积神经网络

在本节课中，我们将学习多层感知机的前向与反向传播过程，并了解卷积神经网络的基本概念及其优势。

概述

上一节我们介绍了神经网络的基本概念。本节中，我们将深入探讨多层感知机的具体工作流程，包括前向传播进行推理，以及反向传播进行训练。随后，我们将引入卷积神经网络，并解释其为何在图像处理等任务中比传统的全连接网络更具优势。

多层感知机的前向与反向传播

多层感知机是一种经典的前馈神经网络。其核心在于通过一系列层（输入层、隐藏层、输出层）对输入数据进行变换，最终得到预测结果。

前向传播：推理过程

给定网络参数（权重 W 和偏置 b，合称为 θ）以及输入 X，前向传播的目标是产生一个预测标签 y（假设是分类任务）。

对于一个三层网络（输入层、隐藏层、输出层）：

  1. 输入层到隐藏层：计算 FC1 = W1 · X + b1，其中 W1 是权重矩阵，b1 是偏置向量。然后，将结果 FC1 通过一个激活函数（如Sigmoid），得到隐藏层输出 S1。
  2. 隐藏层到输出层：计算 FC2 = W2 · S1 + b2。同样，将 FC2 通过激活函数，得到输出 S2。
  3. 最终分类：对于分类任务，我们通常不会直接使用 S2 中的最大值（argmax），因为它在训练时导数不友好。相反，我们会使用 Softmax 函数将 S2 转换为一个概率分布 K，其中每个元素 K_i 代表属于第 i 类的概率。

前向传播公式总结：
Y = Softmax( W2 · Sigmoid( W1 · X + b1 ) + b2 )

反向传播：训练过程

训练的目标是找到一组参数 θ，使得网络在所有训练数据上的总损失最小。我们通过定义一个损失函数 E 来衡量预测值与真实标签之间的差距，例如均方误差：E = (y - t)^2，其中 t 是真实标签。

为了最小化 E，我们需要计算损失函数相对于每个权重和偏置的梯度（即导数 ∂E/∂W, ∂E/∂b），然后使用梯度下降法更新参数。

由于网络是层层嵌套的，我们使用链式法则从输出层向输入层反向计算梯度：

  1. 计算损失 E 对网络最终输出 Y 的梯度：∂E/∂Y。
  2. 计算 Y 对上一层输出 S2 的梯度：∂Y/∂S2，这取决于 Softmax 的导数。
  3. 将两者相乘得到 ∂E/∂S2 = (∂E/∂Y) * (∂Y/∂S2)。
  4. 继续反向计算 ∂E/∂FC2、∂E/∂S1、∂E/∂FC1，最终得到 ∂E/∂W1 和 ∂E/∂b1。

反向传播核心：利用链式法则，将误差从输出层逐层反向传递，并利用前向传播时存储的中间值（如 S1, S2）高效计算梯度。

激活函数与Softmax

在多层感知机中，激活函数为网络引入了非线性，使其能够学习更复杂的模式。以下是常见的激活函数：

  • Sigmoid (Logistic) 函数：σ(x) = 1 / (1 + e^{-x})。它将输入映射到 (0, 1) 区间，其导数可以用自身表示：σ'(x) = σ(x) * (1 - σ(x))，计算方便。
  • 修正线性单元 (ReLU)：ReLU(x) = max(0, x)。计算速度快，但负值区域的导数为0。变体如 Leaky ReLU (f(x) = max(αx, x)， α≈0.01) 解决了这个问题。
  • Softmax 函数：用于多分类任务的输出层，将一组实数转换为概率分布。对于类别 C，给定向量 Z，第 i 类的概率为：
    K_i = e^{Z_i} / Σ_{j=0}^{C-1} e^{Z_j}
    Softmax 的输出所有元素之和为1，每个元素在0到1之间，可被解释为概率。

梯度下降与训练

我们使用梯度下降法来更新网络参数。基本步骤如下：

  1. 前向传播：输入一个（或一批）训练样本 X，计算网络输出 K。
  2. 计算损失：根据真实标签 t 和输出 K 计算损失 E。
  3. 反向传播：计算损失 E 关于所有参数 θ 的梯度 ∇E。
  4. 参数更新：沿梯度反方向更新参数：θ_new = θ_old - ε * ∇E，其中 ε 是学习率。

为了提高效率并避免陷入局部最优，实践中常采用：

  • 随机梯度下降 (SGD)：每次随机选取一个样本计算梯度并更新。
  • 小批量梯度下降：每次选取一小批样本计算平均梯度并更新，这是最常用的方法。
  • 带动量的梯度下降：更新方向不仅考虑当前梯度，还累积一部分历史梯度，有助于加速收敛并稳定训练过程。

过拟合

在训练过程中，如果模型在训练集上表现越来越好，但在未见过数据（验证集）上表现变差，则发生了过拟合。这意味着模型过度学习了训练数据中的噪声和细节，而非一般性规律。

防止过拟合的策略包括：

  • 早停：当验证集损失不再下降时停止训练。
  • 使用验证集：用未参与训练的数据监控模型泛化能力。
  • 理解本质：我们最终目标是逼近一个能处理无限可能输入的函数，而非完美拟合有限的训练样本。

从多层感知机到卷积神经网络

尽管多层感知机功能强大，但在处理像图像这样的高维数据时存在明显缺陷。例如，一个28x28的手写数字图像有784个像素。若使用一个仅含10个神经元的隐藏层，仅第一层就需要 784 * 10 + 10 = 7850 个参数。对于更大图像，参数数量会爆炸式增长，导致计算和存储成本极高，且训练困难。

卷积神经网络通过两种关键思想解决了这些问题：

  1. 局部连接：每个神经元只与输入图像的一小块局部区域（如5x5）相连，而非全部像素。这基于“图像中有意义的特征通常出现在局部”的观察。
  2. 权值共享：同一个卷积核（滤波器）被滑动应用于整个输入图像的不同位置。这意味着无论特征出现在图像的哪个角落，都使用相同的权重来检测它。这大大减少了需要学习的参数数量，并使网络对输入的位置变化具有平移不变性。

卷积操作公式（以2D图像为例）：
Y[i, j] = Σ_{m=0}^{K1-1} Σ_{n=0}^{K2-1} W[m, n] * X[i+m, j+n] + b
其中 W 是 K1 x K2 的卷积核，X 是输入图像，Y 是输出特征图，b 是偏置。

经典CNN结构：LeNet-5

LeNet-5是一个用于手写数字识别的早期成功CNN模型，其结构清晰地展示了CNN的组成：

  1. 卷积层 (C1)：使用6个5x5的卷积核，从32x32输入生成6个28x28的特征图。
  2. 池化层 (S2)：对每个特征图进行2x2下采样（如取最大值），得到6个14x14的特征图。池化用于降维和提供一定程度的平移不变性。
  3. 卷积层 (C3)：使用16个5x5的卷积核，与S2层的特征图进行连接（非全连接），生成16个10x10的特征图。
  4. 池化层 (S4)：再次进行2x2下采样，得到16个5x5的特征图。
  5. 卷积层 (C5)：使用120个5x5的卷积核，生成120个1x1的特征图（可视为一个120维向量）。
  6. 全连接层 (F6)：将C5的120维向量通过全连接层映射到84维。
  7. 输出层：最后通过全连接层映射到10个输出（对应0-9数字），并使用高斯连接等技巧（现代网络已较少使用）。

这种“卷积+池化”交替，最后接全连接层的结构，成为后来许多CNN模型的基础范式。

为什么CNN更有效？

  1. 参数效率高：通过局部连接和权值共享，显著减少了参数量。
  2. 平移不变性：相同的滤波器扫描整个图像，使网络不关心特征出现的位置。
  3. 层次化特征提取：浅层卷积核学习边缘、颜色等低级特征；深层卷积核组合低级特征，形成更复杂的高级特征（如形状、物体部件）。
  4. 更适合高维空间数据：如图像、音频、视频，其中数据在空间或时间上具有强相关性。

自2012年AlexNet在ImageNet竞赛中凭借GPU加速的深度CNN大幅领先传统计算机视觉方法后，CNN迅速成为计算机视觉领域的主流方法。

总结

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

  1. 多层感知机的工作原理，包括利用前向传播进行预测，以及利用反向传播和链式法则计算梯度进行训练。
  2. 常见的激活函数（Sigmoid, ReLU）和用于多分类的Softmax函数。
  3. 使用梯度下降及其变体（小批量、带动量）来优化网络参数。
  4. 机器学习中的常见问题——过拟合及其应对策略。
  5. 卷积神经网络的核心思想：局部连接和权值共享，以及它们如何解决全连接网络在处理图像数据时的瓶颈。
  6. 通过LeNet-5的结构，了解了CNN的典型组成：卷积层、池化层和全连接层的堆叠。
  7. CNN在参数效率、平移不变性和层次化特征学习方面的优势，这使其在图像识别等任务上取得了巨大成功。

下一节，我们将更具体地探讨CNN中各层的维度计算、实现细节以及在GPU上进行优化的策略。

013：卷积神经网络层实现 🧠

在本节课中，我们将要学习卷积神经网络中卷积层和池化层的实现细节，包括其计算原理、顺序代码结构、并行化策略，以及如何通过矩阵乘法来优化卷积运算。

概述

上一节我们介绍了CNN的基本概念和结构。本节中，我们来看看如何具体实现CNN中的核心计算层——卷积层和池化层。我们将从顺序代码开始，逐步分析其并行化潜力，并探讨如何利用GPU架构进行高效计算。

卷积层计算原理

卷积层的输入称为X，包含B个不同的图像。每个图像有C个通道（例如，RGB图像C=3）。每个通道有H x W个像素。因此，输入数据集是一个四维张量：[B, C, H, W]。

卷积层的权重集（滤波器）用于提取特征。假设我们有M个特征图（滤波器），所有滤波器大小相同，为K x K。每个特征图对应C个通道，每个通道有一个独立的K x K滤波器。因此，权重张量的维度是 [M, C, K, K]。

卷积操作包括逐点乘法和求和。首先，在每个通道上进行卷积运算，然后将所有通道的结果求和，从而“折叠”通道维度。在边界处，我们忽略滤波器会超出输入范围的输出像素，因此输出尺寸会变小。

输出Y的尺寸计算如下：

  • H_out = H - K + 1
  • W_out = W - K + 1

输出Y也是一个四维张量：[B, M, H_out, W_out]。

顺序代码实现

理解并行化的第一步是分析顺序（CPU）版本的代码。以下是卷积层前向传播的顺序伪代码：

// 假设数组已适当线性化，此处使用多维索引便于理解
for (int b = 0; b < B; ++b) { // 遍历批次中的图像
    for (int m = 0; m < M; ++m) { // 遍历输出特征图
        for (int h = 0; h < H_out; ++h) { // 遍历输出高度
            for (int w = 0; w < W_out; ++w) { // 遍历输出宽度
                float sum = 0.0f;
                for (int c = 0; c < C; ++c) { // 遍历输入通道
                    for (int p = 0; p < K; ++p) { // 遍历滤波器高度
                        for (int q = 0; q < K; ++q) { // 遍历滤波器宽度
                            // 计算输入像素和滤波器权重的乘积并累加
                            sum += X[b][c][h + p][w + q] * W[m][c][p][q];
                        }
                    }
                }
                Y[b][m][h][w] = sum; // 存储结果
            }
        }
    }
}

这段代码清晰地展示了七层嵌套循环。我们的并行化目标就是将最外层的几个循环（图像、特征图、输出像素）映射到GPU的线程和线程块上。

并行化策略分析

在GPU上实现卷积时，我们需要考虑如何将计算任务分配给大量线程。以下是主要的并行化维度：

  1. 跨图像并行：批次B中的不同图像可以完全独立处理。
  2. 跨特征图并行：不同的输出特征图M可以并行计算。
  3. 跨输出像素并行：每个输出特征图内的像素(H_out, W_out)可以并行计算。
  4. 跨输入通道并行：对输入通道C的求和操作理论上也可以并行，但这需要线程间同步（例如原子操作）来合并结果。

选择哪种并行化策略取决于具体问题的参数（如图像大小、滤波器大小、特征图数量）和目标GPU的架构。通常，输出像素数量巨大，因此跨输出像素并行是主要手段。当特征图数量较少时，仅靠像素并行可能不足以充分利用GPU，此时可能需要结合其他维度。

基础GPU内核设计

一个基础的GPU卷积内核可能将线程块分配给输出特征图中的一个图块（Tile）。假设每个线程块计算一个 TILE_WIDTH x TILE_WIDTH 大小的输出像素块。

以下是线程块和网格的映射方式示例：

  • 网格X维度：映射到输出特征图索引 M。
  • 网格Y维度：线性化映射到单个特征图内所有图块的位置（需要将二维的图块索引转换为一维）。
  • 网格Z维度：映射到批次中的图像索引 B。

主机端代码可能如下所示：

#define TILE_WIDTH 16
dim3 blockDim(TILE_WIDTH, TILE_WIDTH, 1);
int W_grid = ceil((float)W_out / TILE_WIDTH);
int H_grid = ceil((float)H_out / TILE_WIDTH);
int Y_grid = W_grid * H_grid; // 单个特征图所需的线程块数量
dim3 gridDim(M, Y_grid, B); // 网格维度
convolutionForwardKernel<<<gridDim, blockDim>>>(...);

在内核中，每个线程需要根据其线程ID和块ID计算出它负责的精确输出位置 (b, m, h, w)，然后执行内层循环（遍历通道和滤波器）来计算该像素的值。

这种基础实现的缺点是输入数据重用性差。同一个输入像素会被多个计算不同输出特征图或相邻输出像素的线程多次读取，导致全局内存带宽成为瓶颈。

通过矩阵乘法优化卷积

为了更高效地利用GPU的矩阵计算单元（如Tensor Cores）和优化数据复用，可以将卷积运算转换为矩阵乘法。

核心思想是“展开”输入数据和滤波器：

  1. 滤波器矩阵：将每个 M x C x K x K 的滤波器展开成矩阵 W‘ 的一行，该矩阵大小为 [M, C*K*K]。
  2. 输入数据矩阵：对于每个输出位置，将对应的输入感受野（C x K x K）展开成矩阵 X‘ 的一列。如果有 H_out * W_out 个输出位置，则 X‘ 的大小为 [C*K*K, H_out*W_out]。

那么，卷积运算 Y = conv(X, W) 就等价于矩阵乘法 Y‘ = W‘ * X‘，其中结果矩阵 Y‘ 的大小为 [M, H_out*W_out]，可以重新变形为 [M, H_out, W_out]。

优势：

  • 可以利用高度优化的矩阵乘法库（如cuBLAS）或手写的高效矩阵乘法内核。
  • 在矩阵乘法中，通过分块计算可以显著提高输入数据和滤波器权重的数据复用率，减少对全局内存的访问。

注意：在实际实现中，我们通常不会在物理内存中真正创建展开后的矩阵 X‘ 和 W‘，因为这会增加内存开销和数据搬运成本。相反，我们使用“隐式展开”，即在矩阵乘法内核中，通过索引计算直接访问原始张量 X 和 W 中的元素，模拟展开后的访问模式。

池化（子采样）层实现

池化层通常跟在卷积层之后，用于降低特征图的空间维度，减少参数和计算量，同时提供一定的平移不变性。

最常见的池化操作是最大池化或平均池化。它在一个 N x N 的窗口内进行操作，并用该窗口的统计量（最大值或平均值）代替该区域。输出尺寸为：

  • H_pool = floor(H_in / N)
  • W_pool = floor(W_in / N)

顺序伪代码如下（以平均池化为例）：

for (int b = 0; b < B; ++b) {
    for (int m = 0; m < M; ++m) {
        for (int x = 0; x < H_pool; ++x) {
            for (int y = 0; y < W_pool; ++y) {
                float sum = 0.0f;
                for (int p = 0; p < N; ++p) {
                    for (int q = 0; q < N; ++q) {
                        sum += inputY[b][m][N*x + p][N*y + q]; // inputY 是卷积层的输出
                    }
                }
                float average = sum / (N * N);
                outputS[b][m][x][y] = sigmoid(average + bias); // 可能通过一个激活函数
            }
        }
    }
}

内核融合优化：池化层计算简单，数据复用率低。一个重要的优化是“内核融合”，即将池化层的计算直接合并到前一个卷积层的内核中。这样，卷积层的输出不必写回全局内存，再由另一个内核读回，而是直接在芯片上（寄存器或共享内存）进行池化操作，最后只将池化后的较小结果写回全局内存。这大大减少了数据移动的开销。

总结

本节课中我们一起学习了卷积神经网络核心层的实现。

  1. 我们首先分析了卷积层的计算原理和顺序代码结构。
  2. 接着，我们探讨了在GPU上并行化卷积运算的各种策略，包括跨图像、特征图和输出像素进行并行。
  3. 然后，我们介绍了一种强大的优化手段：将卷积运算转换为矩阵乘法，从而利用GPU高效的数据复用模式和专用硬件。
  4. 最后，我们简要介绍了池化层的实现，并提到了通过内核融合来减少数据移动的重要优化技术。

理解这些底层实现细节对于在GPU上高效运行深度学习模型至关重要。在实际项目中，需要根据具体的层参数和硬件特性，灵活选择和组合这些策略以达到最佳性能。

014：子图同构与动态负载均衡 🚀

在本节课中，我们将学习一个重要的图挖掘问题——子图同构，并探讨如何利用CUDA的动态负载均衡技术来高效解决该问题。我们将从问题定义出发，逐步深入到算法实现、性能挑战以及优化策略。

概述 📋

子图同构是图论中的一个基本问题，旨在从大型数据图中找出所有与给定查询图结构相同的子图。这个问题在生物信息学、社交网络分析、自然语言处理等领域有广泛应用。然而，由于其计算复杂度呈指数级增长，且需要处理海量数据，因此对算法的性能和可扩展性提出了极高要求。

问题定义与挑战 ⚙️

子图同构问题涉及两个图：数据图 ( G_d ) 和查询图 ( G_q )。目标是找出数据图中所有与查询图同构的子图。同构意味着两个图在结构上完全相同，即节点和边一一对应。

计算复杂度

子图同构是一个NP完全问题，计算复杂度极高。例如，在一个包含6500万节点和18亿边的图中，仅查找大小为6的团（完全子图）就可能产生数千亿个实例。随着查询图规模的增大，实例数量呈指数级增长，这对内存和计算资源提出了巨大挑战。

现有方法的局限性

传统的解决方法主要分为深度优先搜索（DFS）和广度优先搜索（BFS）。CPU上通常采用DFS，因为它易于实现且内存占用较低，但并行度有限。GPU上更适合BFS，因为它可以充分利用大量线程并行处理，但内存需求巨大，难以扩展。

算法基础 🧠

搜索算法

子图同构的搜索算法基于逐步匹配节点。首先，根据查询图的最小度数过滤数据图中的节点，排除不可能匹配的节点。这一过程称为“剥离”。然后，将查询图转换为有向无环图（DAG），并按照拓扑顺序遍历数据图，逐步匹配节点。

关键操作

搜索过程中的关键操作包括：

  1. 掩码操作：根据度数或对称性剪枝，排除不可能匹配的节点。
  2. 交集操作：生成新候选节点时，需要计算邻居节点的交集。这是算法中最耗时的操作，占总计算时间的99%以上。

GPU上的优化策略 ⚡

现有工作：Parsic

Parsic是首个在GPU上实现DFS的子图同构解决方案。它通过二进制编码和快速交集操作提升了性能，并采用了两种并行策略：

  • 节点每块：每个线程块处理一个父节点，块内的线程负责遍历其子图。
  • 边每块：进一步细分工作，每个线程块处理一条边，以增加并行度。

负载均衡问题

Parsic采用静态分配策略，将线程块固定分配给数据图中的节点。然而，数据图中节点的度数分布不均，导致某些线程块工作负载过重，而其他线程块空闲，从而引起负载失衡。

动态负载均衡 🔄

线程块调度

在CUDA中，线程块是独立调度的。每个线程块可以独立执行，且顺序不受其他块影响。这种设计提供了透明的可扩展性，但同时也带来了负载均衡的挑战。如果线程块的工作负载差异较大，会导致某些流多处理器（SM）空闲，而其他SM过载，从而降低整体性能。

动态负载均衡的实现

为了解决负载失衡问题，我们引入了动态负载均衡机制。具体步骤如下：

  1. 启动最大并发线程块：通过cudaOccupancyMaxActiveBlocksPerMultiprocessor API确定可同时运行的最大线程块数，并以此启动内核。
  2. 工作队列管理：使用全局原子操作分配节点ID给空闲线程块。当线程块完成当前任务后，将其加入工作队列。
  3. 工作卸载：繁忙的线程块可以将其DFS栈中的部分工作卸载给空闲线程块，通过全局内存和原子操作实现通信和同步。

性能提升

通过动态负载均衡，我们显著减少了负载失衡。例如，在某些测试案例中，最忙SM与平均SM的工作负载差异从6倍降低到1.06倍。这使得整体性能提升了7倍以上，并且在多GPU环境下也表现出良好的可扩展性。

总结 🎯

本节课我们一起学习了子图同构问题的定义、挑战及其在GPU上的优化策略。通过引入动态负载均衡，我们有效解决了静态分配导致的负载失衡问题，显著提升了算法性能。关键要点包括：

  • 子图同构是NP完全问题，计算复杂度高，需要高效且可扩展的解决方案。
  • GPU上的DFS实现面临负载均衡挑战，动态负载均衡通过工作卸载和全局通信机制有效解决了这一问题。
  • 线程块独立调度的特性使得动态负载均衡成为可能，但需要仔细设计通信和同步机制。

希望本节课的内容能帮助你更好地理解如何利用CUDA的高级特性优化复杂计算问题。

015：嘉宾讲座 - 关系图神经网络的高效编程与编译框架

概述

在本节课中，我们将学习关系图神经网络的高效编程与编译框架。我们将探讨图神经网络的基本概念、面临的性能挑战，以及一个名为Hector的系统如何通过创新的中间表示和代码生成技术来解决这些挑战。

图神经网络背景与目标 🎯

图神经网络的目标场景是大规模图表示学习。以下是两个典型例子：

例子1：链接预测
对于电子商务公司，他们可能希望追踪每位买家购买的产品，并向消费者推荐潜在商品。问题在于：某位消费者是否可能成为某特定商品的潜在买家？如果我们将此建模为一个网络，其中买家和商品都是节点，购买历史是链接，那么这就可以建模为链接预测问题。如果买家对商品感兴趣，他们之间就应该存在一条链接。

例子2：节点分类
考虑一个引文网络。我们可以将计算机科学领域的出版物建模为一个网络，其中文章是节点，并可以给它们贴上不同主题的标签，例如系统会议、计算机视觉或自然语言处理。我们有文章引用的历史记录。当出现一个新节点时，我们想了解它属于哪个领域，这可以建模为节点分类问题。

图表示学习的好处在于，它可以学习图中嵌入的高级结构化信息。图表示学习所做的是编码这些信息，并将其转化为低维稠密向量。我们有一个模型，它接收图并进行学习，然后模型会输出表示，即每个节点的低维稠密向量。一旦我们有了这些低维稠密向量，就可以将它们输入到一个分类器函数中，以判断例如某个节点是否是一篇系统领域的论文。

关系图神经网络简介

近年来，关系图神经网络被广泛采用。RGNN是GNN的一种变体。GNN是一种深度学习模型，根据输入图数据的稀疏矩阵进行操作。它以图作为输入，输出节点的结构。RGNN考虑了不同的类型，并为每种边类型和节点类型分配权重。

我们以RGAT层的执行为例。RGAT是关系图神经网络的一种变体。每个节点的特征向量结果，其输入是每个节点的嵌入向量。嵌入向量是从描述节点的现有信息中产生的向量，例如，它可能是从产品描述中学习到的向量。经过这个RGAT层后，每个节点的输出特征现在包含了与其有边相连的节点的信息。

图神经网络有两个阶段。让我们从后往前看：

  • 第二阶段是节点聚合。例如，对于节点1，边01和21上的消息以加权方式聚合。这里的消息和权重都是在第一阶段生成的。权重是注意力，由模型产生，用于强调或减弱消息。
  • 第一阶段是消息生成阶段。对于节点1，我们需要计算两条边上的注意力和消息。为了生成消息，RGAT模型对源节点特征应用一个线性层。注意力是通过一个以源节点特征、目标节点特征和边类型为参数的函数计算出来的。这是一个接收源节点嵌入、目标节点嵌入以及根据边类型的多个权重的过程。

与节点1的例子类似，这个RGAT层将对所有节点应用相同的过程。

RGNN执行面临的三大挑战

RGNN执行主要面临三个挑战：

1. 固有的内存密集型特性
我们展示了现有性能最佳编译器Graphiler的推理时间分解。我们在两个数据集上通过两个模型执行一个生成节点特征的模型。这两个模型是HGT和RGAT，两个数据集是MUTAG和Freebase 15K。PyTorch Geometric的图变换器计算更密集，MUTAG是一个较小的图。如图所示，分析器在索引和复制上花费了显著的时间。这些是为了准备高度专业化的内核（如矩阵乘法）而进行的大量数据移动。许多计算不是矩阵乘法，它们被归类为其他计算，并涉及向量点积和矩阵向量乘法等例子，它们通常是内存密集型的。

2. PyTorch编程接口缺乏张量视图导致不必要的数据复制
让我们以消息生成为例，同样的问题也存在于RGAT模型的边注意力计算中。为了生成消息，我们需要在每条边上进行矩阵乘法。典型的做法是使用PyTorch提供的批处理矩阵乘法。在PyTorch中，批处理矩阵乘法本质上是通用矩阵乘法的三维扩展。所有三个张量（两个输入张量和输出张量）都具有相同的前导维度。第一个维度始终是样本数量。对于第一个维度中的每个切片，我们从第一个输入得到一个二维矩阵，从第二个输入得到另一个二维矩阵，我们对它们应用矩阵乘法，然后得到一个输出。通过对样本数量次重复此操作，我们可以得到一个三维输出。基本上，它的作用是在第一个维度上有多个样本切片，然后在每个切片中执行通用矩阵乘法。

因此，在每条边上的矩阵乘法就是批处理矩阵乘法中的一个样本。在进行批处理矩阵乘法之前，权重张量必须被具体化，这导致了冗余的浪费问题。基本上，在PyTorch中，你可以访问张量，当你进行批处理矩阵乘法时，你会写类似 C = A @ B 的语句，这里的B是一个具体的张量。它不是一系列指针，你必须将这些权重移动到这个新张量中并具体化，这导致了非常大的内存消耗。

3. 由于内核与数据访问方案耦合导致的高优化成本
对于异构图，我们可以将每种边类型的邻接矩阵存储在单独的稀疏矩阵中，或者也可以存储在一个稀疏矩阵中。如果我们要为这些选项的每种组合维护内核，我们将需要维护和优化数量成倍增长的内核。

Hector框架概述

为了系统地解决这些挑战，我们提出了Hector中间表示和代码生成框架。

  • Hector捕获了RGNN模型的关键特性，并提供了减少内存访问的机会。
  • 通过中间表示调度和具体化来避免图复制。
  • Hector捕获模型语义，并生成具有灵活数据访问方案的代码，以消除冗余数据复制，避免高优化成本。
  • Hector将模型、数据布局和运算符特定调度彼此解耦。

接下来，我将通过一个例子说明Hector的工作流程和两项优化。

Hector工作流程与软件架构

这是Hector的工作流程和软件架构。为了方便起见，我们在编程接口提供了类似DGL的API。代码生成器将处理被装饰的函数。基本上，在Python中，这是应用装饰器的语法糖。当我们提供一个语法时，我们提供一个装饰器，被装饰的函数将由我们的代码生成器处理。

然后，代码生成将处理被装饰的函数。API调用首先被转换到中间表示级别。中间表示扩展了Python语言以表达模型语义。然后代码被降低到运算符内级别中间表示。在这个级别，它指定了CUDA内核将如何生成。这里它指定了边类型线性运算符应该被降低到一个GEMM内核，并指定了加载和存储方案。

代码生成器然后基于模板生成CUDA内核，编译并在PyTorch运行时注册。然后，在GNN执行期间，PyTorch运行时将使用这个CUDA内核进行模型执行。

中间表示级别优化：线性运算符重排序

让我们看看中间表示级别中间表示是如何工作的。这是一个例子。中间表示级别中间表示将模型语义定义为一序列运算符，并以调整张量布局的方式呈现。注意，我们这里使用了一个不同于引文网络的图。图中的节点是论文、作者和机构。边表示关系，例如论文引用另一篇论文，以及作者属于某个机构。在我们的模型中，在边消息生成期间，每条边上都有一个矩阵乘法。为了在中间表示级别表达这一点，我们做的是创建一个边上的迭代，并指定矩阵乘法。用Python语法表示，我们从边源节点获取特征，切片权重张量，然后将结果显示在边元素中。

我们在这里提出的一项优化是线性运算符重排序。这发生在中间表示级别，目的是减少计算复杂性和中间数据集大小。乘法次数是 N*P*Q + N*P*Q。第一个维度是N，然后我们有K和Q作为矩阵乘法的两个维度。通过重新排序这三个矩阵之间的乘积计算，我们可以将第二项中的N因子减少到K。基本上，这里K是隐藏维度，N是边的数量。我们将处理大规模图，所以N可能非常大，可能达到数百万。但这里的隐藏维度非常有限，因为这是节点的低维稠密向量表示，所以隐藏维度通常很低，大约几百个元素。因此，中间数据大小也缩小了。

运算符内级别优化：消除冗余权重复制与复合具体化

接下来，我们将展示运算符内级别中间表示如何帮助消除冗余权重复制并启用复合具体化。

这里我展示的是RGAT层中的边消息生成。这可以通过一个GEMM内核有效地完成。如果我们想使用一个GEMM内核，它包含三个步骤：

  1. 首先，根据源节点索引数组收集源节点特征。
  2. 其次，它在这里进行GEMM计算。基本上，这里的每个虚线矩形可以是一系列块。我们将这些矩阵乘法分配给几组块来处理。在一个简单的GEMM内核中，我们会进行分块，然后将这些分块乘法分配给每个块。但在这里，我们做的是三种不同类型的矩阵乘法，它们有不同的权重，输入以类似的方式放入内核。所以发生的情况是，我们有三个块组，每个组处理一种类型的边。例如，对于引用边，我们将专门使用该边类型的权重。虚线框显示了三种不同的矩阵乘法组。当我们在这里做矩阵乘法时，我们需要做一些指针操作。基本上，输入是在节点域中。第一个维度是节点索引。然后我们在这些边上执行这些矩阵乘法，所以我们需要从节点维度转换到边维度。我们查看存储每条边源节点索引的数组，并使用统计信息进行指针方案，将它们引导到相应的块中。
  3. 第三阶段，我们将这些结果分散到存储边中间数据的张量缓冲区中。这是一个分散方案。因为我们遵循将相同类型的边的消息存储在一起的方案，所以你可以看到基本上有三个段。前两个元素彼此相邻存储，因此对于引用地址，第二部分是第三到第五个切片，这是用于撰写地址，最后两个用于雇佣地址。

当我们生成具有自定义加载/存储方案的内核时，我们只保留权重的一个副本，并在矩阵乘法期间从中加载正确的矩阵，从而消除了由于PyTorch缺乏张量视图而导致的冗余权重复制。相比之下，在PyTorch中，它会进行批处理矩阵乘法，我们想要得到这个，我们需要输入两个张量，这两个张量都是三维的，并且我们需要具体化存储每条边权重的第二个张量。这意味着我们需要为每条边读取相应的权重并进行复制，这造成了非常大的权重复制开销。

我们在这里做的另一项优化是复合具体化。这消除了冗余计算和内存占用。我们注意到这里的一些计算是重复的。边消息仅依赖于源节点和边类型。例如，如果一个作者写了多篇论文，那么这些边上的消息是相同的。为了利用这一点，输出张量将具有不同的布局。原本输出张量的每一行都有边特征，但现在每一行都有一个对应于唯一节点-边类型对的特征。这种变化在降低到运算符内级别时反映在布局选择中。输入矩阵X的分散和输出矩阵Y的分散会相应改变。当涉及到后续内核时，布局的变化将得到反映，并且在降低过程中，内核的收集方案将相应改变。

运算符内级别中间表示指定了CUDA内核将如何生成。Hector原型基于两个模板：一个是GEMM模板，另一个是遍历模板。GEMM模板涉及灵活的加载和存储方案，允许配置使用不同的收集和分散方案。遍历模板是一个通用的边稀疏内核，我们尽可能将其运算符降低到GEMM以获得最佳性能。

性能评估

我们在单GPU机器上评估了Hector系统。异构图由DGL库和OGB库提供。DGL代表深度图库，是亚马逊的一个GNN执行框架。OGB代表开放图基准，它为算法开发者和图处理系统开发者存储公开可用的图。

将最优化后的Hector代码与现有系统进行比较，Hector的性能始终更优。Hector也更节省内存，因为它不会在大图上触发任何内存不足错误，而现有系统会触发内存不足错误。我们展示了两个部分：第一部分是在三个模型上的推理，第二部分是在这三个模型上的训练，我们在从大到小的不同图上进行。例如，Wiki KG2是从维基百科创建的大型图，MUTAG较小。我们在不同模型和不同图上展示这些执行，并与基线进行比较，我们看到性能更好，并且我们没有触发任何内存不足错误。

为了展示Hector如何缓解内存密集型问题，我们回到推理时间分解图，并将分析器的性能与Hector进行比较。为了清晰起见，我们放大了图的部分。如图所示，Hector几乎没有索引和复制操作，因为它避免了数据移动。它确实有少量从主机到设备的数据传输，以传递指定块分配的元数据。与批处理矩阵乘法相比，批处理矩阵乘法需要从设备到设备复制以创建大的中间张量，因此会有大量的索引和复制操作，这在内存占用方面效率低下，并且引入了实际的内存密集型操作。

我们还研究了这些优化如何影响执行。我们展示了一个表格，比较了不同模型、图和优化组合的训练。我们这里有两个优化：C代表复合具体化，R代表线性运算符重排序，然后我们可以一起使用它们。我们可以意识到的一点是，对于不同的训练目标组合（无论是训练还是推理、模型和图），实际上最优的优化组合是不同的。这是因为对于图神经网络，执行是输入依赖的。例如，当一个图的度较大时，意味着每个节点的平均边数较多。当度较大时，我们将有更密集的边操作。此外，当我们进行矩阵乘法，从节点域收集数据以产生边消息数据时，我们将在那里进行更多的计算。当输入不同时，计算的特征也不同。我们在这里展示了一个分解示例。A和Freebase 15k是两个图。我们展示了使用两种优化、仅使用一种优化以及不使用任何优化的情况。蓝色部分表示GEMM内核的执行时间，橙色部分是遍历内核的时间，灰色部分是PyTorch操作的时间。

总结

本节课中，我们一起学习了关系图神经网络的高效编程与编译框架。我们首先了解了图神经网络的基本概念和应用场景。然后，我们深入探讨了RGNN执行面临的三大挑战：固有的内存密集型特性、PyTorch接口导致的数据复制以及内核与数据方案耦合的高优化成本。

为了系统性地解决这些问题，我们介绍了Hector框架。Hector通过其创新的中间表示和代码生成器，实现了模型、数据布局和调度的解耦。我们重点讲解了两项关键优化：在中间表示级别通过线性运算符重排序减少计算和内存占用；在运算符内级别通过生成自定义GEMM内核消除冗余权重复制，并利用复合具体化进一步减少冗余计算。

性能评估表明，Hector在多种模型和图数据上均优于现有系统，并且具有更高的内存效率。这些优化效果会根据具体的图和模型特征而变化，为未来的自动优化选择提供了研究方向。

通过本讲，我们看到了如何通过编译器和系统层面的创新，来显著提升图神经网络这类内存密集型应用的执行效率。

016：CUDA流与高级并行编程模型

概述

在本节课中，我们将学习CUDA流的概念，它如何通过任务并行性来提升GPU程序的性能。我们还将探讨如何通过流来重叠数据传输与内核计算。最后，我们将简要介绍CUDA之外的其他并行编程模型，如OpenCL、HIP和OpenACC，了解它们的设计理念和适用场景。

CUDA流：任务并行性的实现

上一节我们讨论了GPU的内存层次结构和优化策略。本节中，我们来看看如何利用CUDA流实现任务级别的并行性，以进一步提升程序性能。

CUDA流是一个操作队列，其中可以包含内核启动、内存拷贝等操作。默认情况下，所有操作都在一个默认流中顺序执行。通过创建多个流，我们可以让不同流中的操作并行执行，从而充分利用GPU的计算资源和DMA引擎。

流的基本概念与工作原理

一个流代表一个按FIFO顺序执行的操作序列。主机线程将操作（如内存拷贝、内核启动）放入流队列，然后GPU驱动异步地从队列中取出并执行这些操作。操作在单个流内严格按顺序执行，但不同流中的操作可以并行执行。

利用流重叠计算与数据传输

以下是通过将数据分段并使用多个流来重叠数据传输与内核计算的理想化流程：

  1. 阶段一：将第一段数据从主机内存传输到设备内存（H2D）。此阶段无法与其他操作重叠。
  2. 阶段二：启动内核计算第一段数据，同时将第二段数据传输到设备。此时，计算与数据传输开始重叠。
  3. 阶段三：进入稳定流水线状态。此时可以并行执行三个操作：将第三段数据传入设备（H2D）、计算第二段数据、将第一段结果传回主机（D2H）。

为了实现这种三路并行，我们至少需要三个流。

代码示例：使用双缓冲流

以下是使用两个流实现双缓冲的基本代码结构。我们为不同的数据段分配独立的内存空间，并在两个流之间交替执行操作。

cudaStream_t stream0, stream1;
cudaStreamCreate(&stream0);
cudaStreamCreate(&stream1);

// 假设已将数据分段，并为每段分配了主机和设备内存 (h_a_seg0, d_a_seg0 等)

// 对偶数段使用 stream0
cudaMemcpyAsync(d_a_seg0, h_a_seg0, size, cudaMemcpyHostToDevice, stream0);
cudaMemcpyAsync(d_b_seg0, h_b_seg0, size, cudaMemcpyHostToDevice, stream0);
vectorAddKernel<<<gridDim, blockDim, 0, stream0>>>(d_a_seg0, d_b_seg0, d_c_seg0, N);
cudaMemcpyAsync(h_c_seg0, d_c_seg0, size, cudaMemcpyDeviceToHost, stream0);

// 对奇数段使用 stream1
cudaMemcpyAsync(d_a_seg1, h_a_seg1, size, cudaMemcpyHostToDevice, stream1);
cudaMemcpyAsync(d_b_seg1, h_b_seg1, size, cudaMemcpyHostToDevice, stream1);
vectorAddKernel<<<gridDim, blockDim, 0, stream1>>>(d_a_seg1, d_b_seg1, d_c_seg1, N);
cudaMemcpyAsync(h_c_seg1, d_c_seg1, size, cudaMemcpyDeviceToHost, stream1);

// 等待所有流完成
cudaStreamSynchronize(stream0);
cudaStreamSynchronize(stream1);

关键点：

  • cudaMemcpyAsync 和内核启动的 <<<...>>> 语法中，最后一个参数指定了流。
  • 内核启动的第三个参数是动态共享内存大小，这里设为0。
  • 需要使用 cudaStreamSynchronize 来等待特定流中的操作完成。

流的硬件演进与注意事项

从软件队列到硬件Hyper-Q

早期GPU的流支持完全由驱动软件实现，这可能导致队头阻塞问题。例如，一个流中的拷贝操作可能会阻塞另一个流中已就绪的拷贝操作。

现代GPU（如Kepler架构及以后）引入了 Hyper-Q 等硬件特性，为每个引擎（如拷贝引擎、计算引擎）提供了多个硬件工作队列。这允许不同流中的操作真正并发执行，消除了虚假的流间依赖，最高支持32路并发。

内核并发执行的考量

当多个流中的内核同时执行时，它们会在GPU上分时共享计算资源。这可能不是你想要的，因为一个内核可能会被另一个内核拖慢。你需要根据应用需求，决定是让内核串行执行以尽快完成第一个任务，还是允许并发以获得整体吞吐量。

如何选择分段大小？

分段大小（即每个流处理的数据量）的选择是一个权衡：

  • 分段过小：内核启动、DMA传输的固定开销占比过大，GPU可能因线程块过少而无法被充分利用。
  • 分段过大：流水线启动和排空阶段的非重叠开销占比过大，无法有效隐藏数据传输延迟。

最佳分段大小取决于具体的内核特性、GPU硬件（SM数量、内存带宽）和数据规模。通常需要通过实验（性能剖析）或自动调优来确定。

张量核心：面向AI的专用硬件

上一节我们介绍了通用的任务并行性。本节中，我们来看看为特定计算模式（如矩阵乘法）设计的专用硬件单元——张量核心。

张量核心是GPU中为加速机器学习（尤其是低精度矩阵乘累加运算）而引入的专用硬件。它们支持16位（half）或8位整数等格式，能显著提升矩阵运算吞吐量。

编程模型：Warp级矩阵乘累加

CUDA提供了 wmma（Warp Matrix Multiply Accumulate）API 来使用张量核心。编程模型围绕“片段”进行操作，一个片段代表一个子矩阵块。

以下是使用张量核心进行矩阵乘法的简化代码框架：

#include <cuda_fp16.h>
#include <cuda_runtime.h>

// 声明用于存储矩阵块的片段
wmma::fragment<wmma::matrix_a, M, N, K, half, wmma::row_major> a_frag;
wmma::fragment<wmma::matrix_b, M, N, K, half, wmma::col_major> b_frag;
wmma::fragment<wmma::accumulator, M, N, K, float> c_frag;

// 1. 从全局内存加载矩阵块到片段
wmma::load_matrix_sync(a_frag, d_a, stride_a);
wmma::load_matrix_sync(b_frag, d_b, stride_b);

// 2. 将累加器片段初始化为零
wmma::fill_fragment(c_frag, 0.0f);

// 3. 使用张量核心执行乘累加运算 D = A * B + C
wmma::mma_sync(c_frag, a_frag, b_frag, c_frag);

// 4. 将结果片段存储回全局内存
wmma::store_matrix_sync(d_c, c_frag, stride_c, wmma::mem_row_major);

关键点：

  • 操作是warp同步的，需要整个warp（32个线程）共同参与。
  • 数据在“片段”对象中管理，对程序员部分透明。
  • 使用张量核心可以带来数量级的性能提升，尤其是在AI和HPC领域。

其他并行编程模型简介

OpenCL：开放标准跨平台加速

OpenCL是一个开放的、跨平台的并行计算框架，支持CPU、GPU、DSP、FPGA等多种设备。其编程模型与CUDA非常相似。

核心概念对应关系：

  • 工作项 -> CUDA线程
  • 工作组 -> CUDA线程块
  • 本地内存 -> CUDA共享内存
  • 私有内存 -> CUDA寄存器

OpenCL代码结构更显式，需要更多设置代码来获取平台、设备上下文等信息。

HIP：CUDA代码的可移植层

HIP是AMD推出的C++运行时API和内核语言，旨在将CUDA代码便捷地移植到AMD GPU上。其API与CUDA高度相似，hipify工具可以自动化大部分转换工作。

代码示例：向量加法

// HIP 代码与 CUDA 代码极其相似
hipMalloc(&d_a, size);
hipMemcpy(d_a, h_a, size, hipMemcpyHostToDevice);

hipLaunchKernelGGL(vectorAddKernel, gridDim, blockDim, 0, 0, d_a, d_b, d_c, N);

hipMemcpy(h_c, d_c, size, hipMemcpyDeviceToHost);
hipDeviceSynchronize();

OpenACC：基于编译指示的并行编程

OpenACC采用“编译器指令”（pragma）的方式，让程序员通过注解来提示编译器何处可以并行化，由编译器负责生成并行代码（包括数据移动和内核启动）。它追求以最小的代码改动获得可观的加速比。

代码示例：矩阵乘法（添加OpenACC指令）

#pragma acc data copyin(M[0:Mh*Mw], N[0:Nh*Nw]) copyout(P[0:Mh*Nw])
{
    #pragma acc parallel loop
    for(int i = 0; i < Mh; i++) {
        #pragma acc loop
        for(int j = 0; j < Nw; j++) {
            float sum = 0.0f;
            for(int k = 0; k < Mw; k++) {
                sum += M[i*Mw + k] * N[k*Nw + j];
            }
            P[i*Nw + j] = sum;
        }
    }
}

特点与权衡：

  • 优点：代码改动小，可维护性高，具备一定的可移植性（可针对CPU/GPU编译）。
  • 缺点：性能严重依赖编译器的优化能力；对复杂、不规则的计算模式优化效果有限；程序员对底层控制较弱。

总结

本节课中我们一起学习了CUDA流的高级用法，理解了如何通过任务并行性来重叠数据传输与计算，从而更充分地利用GPU硬件。我们还探讨了选择合适分段大小的考量因素。接着，我们介绍了专为AI计算设计的张量核心及其编程接口。最后，我们概览了CUDA之外的几种并行编程模型（OpenCL、HIP、OpenACC），了解了它们不同的设计哲学和适用场景，认识到在编程便利性、性能控制和可移植性之间总是存在权衡。掌握这些概念有助于你为不同的应用场景选择合适的工具和优化策略。

017：替代编程模型与高性能计算展望 🚀

在本节课中，我们将学习两种主要的替代并行编程模型：OpenACC（基于编译制导）和MPI（基于消息传递）。我们将了解它们的基本概念、执行模型，并通过示例代码理解其工作原理。最后，我们将探讨并行编程的历史、现状以及未来面临的挑战。

OpenACC 编译制导模型 🧠

上一节我们介绍了CUDA编程模型。本节中，我们来看看一种更高级的抽象模型——OpenACC。OpenACC允许程序员通过向顺序代码中添加特殊的编译指令（Pragmas）来指示编译器如何自动并行化代码，以及如何管理GPU内存的拷贝。

OpenACC 执行模型

OpenACC的执行模型与CUDA的线程块和线程概念类似，但使用了不同的术语：

  • Gangs： 类似于CUDA中的线程块。
  • Workers： 类似于CUDA中线程块内的线程。

并行性表达

以下是OpenACC中表达不同层级并行性的方式：

1. 冗余执行

#pragma acc parallel num_gangs(1024)
{
    // 这段代码将被1024个gangs冗余执行
    a = b + c;
}

这段代码声明了1024个gangs，每个gang都会执行相同的代码块。这适用于需要多次独立运行相同计算的情况，例如蒙特卡洛模拟。

2. 循环并行化（Gang级）

#pragma acc parallel num_gangs(1024)
{
    #pragma acc loop gang
    for (int i = 0; i < 2048; i++) {
        // 循环的2048次迭代被分配到1024个gangs上
        // 每个gang处理大约2次迭代
        c[i] = a[i] + b[i];
    }
}

这里，外层的parallel指令创建了gangs，内层的loop gang指令将循环迭代映射到这些gangs上。

3. 嵌套并行化（Gang和Worker级）

#pragma acc parallel num_gangs(1024) num_workers(32)
{
    #pragma acc loop gang
    for (int i = 0; i < 2048; i++) {
        #pragma acc loop worker
        for (int j = 0; j < 512; j++) {
            // 外层i循环在gangs间并行
            // 内层j循环在每个gang的workers间并行
            c[i][j] = a[i][j] + b[i][j];
        }
    }
}

这个例子创建了1024个gangs，每个gang有32个workers。总并行度为 1024 * 32 = 32768。

规约操作

在并行循环中进行规约操作（如求和）时，需要显式告知编译器：

float sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
for (int i = 0; i < N; i++) {
    sum += a[i] * b[i]; // 并行规约求和
}

reduction(+:sum) 子句告诉编译器这是一个对变量sum的加法规约。编译器会为每个并行执行单元创建私有副本，最后合并结果。需要注意的是，浮点数的加法不满足结合律，并行规约可能导致与顺序执行略有不同的结果，这是被允许的。

消息传递接口（MPI）模型 📨

上一节我们看了基于共享内存抽象的OpenACC。本节中，我们来看看面向分布式内存系统的编程模型——MPI。MPI是高性能计算领域使用最广泛的编程模型之一，用于在由网络连接的多个计算节点（每个节点有自己的内存）上编写并行程序。

MPI 基本概念

在MPI模型中：

  • 每个进程运行在独立的地址空间中。
  • 进程间通过发送和接收消息来通信和协调。
  • 没有全局共享内存。

MPI 程序结构

一个典型的MPI程序遵循以下模式：

1. 初始化与信息获取

#include <mpi.h>

int main(int argc, char** argv) {
    int rank, size;

    MPI_Init(&argc, &argv); // 初始化MPI环境
    MPI_Comm_rank(MPI_COMM_WORLD, &rank); // 获取当前进程的ID（rank）
    MPI_Comm_size(MPI_COMM_WORLD, &size); // 获取进程总数

    // ... 并行计算代码 ...

    MPI_Finalize(); // 结束MPI环境，同步所有进程
    return 0;
}
  • MPI_Init 必须在所有其他MPI调用之前执行。
  • MPI_Comm_rank 返回当前进程在通信域（如MPI_COMM_WORLD）中的唯一标识（从0开始）。
  • MPI_Comm_size 返回通信域中的进程总数。
  • MPI_Finalize 清理MPI环境，并确保所有进程都完成计算后才退出。

2. 点对点通信：发送与接收

MPI的核心操作是MPI_Send和MPI_Recv。

  • 阻塞发送:

    MPI_Send(void* buf, int count, MPI_Datatype datatype,
         int dest, int tag, MPI_Comm comm);
    • buf: 发送数据缓冲区的起始地址。
    • count: 发送的元素数量。
    • datatype: 发送数据的MPI数据类型（如MPI_INT, MPI_FLOAT）。
    • dest: 目标进程的rank。
    • tag: 消息标签，用于区分不同类型的信息。
    • comm: 通信域。

  • 阻塞接收:

    MPI_Recv(void* buf, int count, MPI_Datatype datatype,
         int source, int tag, MPI_Comm comm,
         MPI_Status* status);
    • source: 发送进程的rank，可使用MPI_ANY_SOURCE接收来自任何进程的消息。
    • status: 返回接收操作的状态信息（如实际接收的数据量、发送者rank等）。
    • 其他参数与MPI_Send对应。

阻塞通信意味着MPI_Send和MPI_Recv调用要等到通信操作（至少是本地部分）安全完成后才返回。 这通常能避免额外的内存拷贝，但需要仔细设计通信顺序以防止死锁。

MPI 向量加法示例

以下是一个简化的MPI向量加法示例，采用“主从”模式：

  • Rank 0 进程作为服务器：负责生成数据、分发数据、收集结果。
  • 其他进程作为计算节点：接收数据、进行计算、返回结果。

服务器进程（Rank 0）伪代码逻辑：

  1. 分配并初始化完整的输入向量A和B。
  2. 将向量A和B的片段依次发送给每个计算节点（rank 1 到 rank size-1）。
  3. 调用MPI_Barrier等待所有计算节点接收完数据。
  4. 依次从每个计算节点接收结果片段，并组合成最终结果向量。
  5. 释放资源并结束。

计算节点进程（Rank > 0）伪代码逻辑：

  1. 分配本地内存以存储接收到的A和B片段。
  2. 按顺序从服务器（rank 0）接收A片段和B片段。
  3. 执行本地向量加法（此处可替换为CUDA内核调用，实现MPI+CUDA混合编程）。
  4. 调用MPI_Barrier与所有进程同步。
  5. 将结果发送回服务器进程。

MPI+CUDA混合编程是现代超算应用的常见模式：使用MPI在节点间协调和通信，使用CUDA在节点内的GPU上进行加速计算。

并行编程的挑战与未来展望 🔮

本节课我们一起学习了OpenACC和MPI。最后，我们来探讨并行编程领域仍然存在的挑战和未来的发展方向。

并行编程的历史与现状

  • 过去：并行计算曾是利基领域，主要由研究生和国家实验室的研究人员为了获得极致性能而使用。
  • 转折点：21世纪初，由于功耗墙限制，单核性能提升放缓，多核/众核处理器成为主流。GPU因其强大的并行计算能力和高内存带宽，从图形处理领域进入通用计算领域。
  • 现在：GPU无处不在（从手机到超算），CUDA等编程模型降低了并行编程门槛，数百万开发者可以编写GPU程序。并行加速带来的性能提升直接转化为成本（如数据中心功耗）和体验（如手机续航）的优势。

当前面临的主要挑战

尽管取得了巨大进展，并行编程仍然面临诸多挑战：

1. 软件工程实践
并行性是一种控制抽象，而传统的软件工程强调数据抽象（如信息隐藏、封装）。这两者之间存在内在冲突。一个经典问题是“同步异常”：很难在不破坏封装的前提下，让派生类安全地扩展现有同步对象（如线程安全队列）的行为。

2. 表达依赖与顺序
在复杂并行任务中，明确表达“任务A必须在任务B之前完成”这类依赖关系仍然很困难。CUDA主要通过线程块内的同步和内核间的顺序启动来表达依赖。虽然引入了动态并行（内核启动内核），但开销较大，且表达能力仍有局限。

3. 确定性与可重现性
并行程序，特别是存在数据竞争的程序，其执行结果可能因线程调度、硬件差异等因素而不同。这给调试带来了极大困难（例如，添加或删除一个printf语句可能就会改变bug的行为）。并行调试器本身也面临挑战：当一个线程命中断点时，其他线程、其他流处理器、子内核等应处于什么状态？

4. 性能可移植性
为一种GPU架构优化的代码，在另一种GPU或未来架构上可能无法高效运行。开发者需要关注众多可能成为瓶颈的因素（寄存器数量、共享内存大小、线程块配置等），并针对不同平台进行调整。

性能评估指标

在评估并行程序时，需要选择合适的指标：

  • 加速比: Speedup = T_sequential / T_parallel
    • 比较并行版本与优化后的最佳顺序版本的运行时间。
    • 固定问题规模。
  • 效率: Efficiency = Speedup / P （P为处理器数量）
    • 衡量对并行资源的利用程度。
  • 可扩展性: 随着处理器数量增加，加速比是否线性增长或效率是否保持稳定。
    • 受限于问题中的固有顺序部分（阿姆达尔定律）。
  • 扩展加速比: 随着处理器数量增加，同时线性增加问题规模，看效率是否能维持。这对于吞吐量导向或受内存容量限制的应用更为合适。

未来展望

并行编程的未来充满机遇：

  • 硬件演进： 随着制程工艺进步放缓，通过增加核心数量、发展晶圆级计算、近内存计算等新架构来提升性能变得更为重要。
  • 编程模型： 需要更高级、更易于表达依赖和同步的抽象，同时不牺牲性能。领域特定语言（DSL）和智能编译器是研究方向。
  • 应用驱动： 人工智能、大数据分析、科学计算等数据密集型应用是推动并行计算发展的核心动力。性能直接关系到产品竞争力和运营成本。

总结：并行编程已经从尖端技术转变为必备技能。尽管OpenACC、MPI、CUDA等工具提供了强大的能力，但在易用性、可维护性、可移植性和确定性方面仍有很长的路要走。这为未来的研究者（也许就是在座的各位！）留下了广阔的空间去创新和解决这些挑战。

本节课中我们一起学习了 OpenACC编译制导模型的基本用法、MPI消息传递模型的核心概念和编程模式，并深入探讨了并行编程领域的历史、现状以及未来面临的核心挑战与机遇。希望这些知识能为你在高性能计算领域的进一步探索打下坚实的基础。

posted @ 2026-03-29 09:31  布客飞龙II  阅读(40)  评论(0)    收藏  举报