Nvidia-CUDA-并行编程笔记-全-
Nvidia CUDA 并行编程笔记（全）
001：NVIDIA CUDA中的HelloWorld程序写法 🚀
在本节课中，我们将学习如何编写和运行CUDA程序。我们将从最简单的“Hello World”程序开始，逐步探索其不同变体，并了解如何执行多线程程序以及如何在GPU上分配内存。
概述
本节课是CUDA编程的入门。我们将首先了解一个基本的CUDA程序结构，然后分析其执行流程。接着，我们会探讨如何通过多线程来并行执行任务，并学习如何管理CPU与GPU之间的异步执行和数据传输。课程的核心目标是理解CUDA内核（Kernel）的启动机制以及线程层次结构。
CPU上的Hello World程序
首先，我们来看一个用C语言编写的标准“Hello World”程序。这个程序与CUDA无关，完全在CPU上运行。
#include <stdio.h>
int main() {
printf("Hello World\n");
return 0;
}
使用GCC编译器编译并运行此程序，会在屏幕上输出“Hello World”。
gcc hello_world.c -o hello_world
./hello_world
第一个CUDA Hello World程序
现在，我们来看一个在GPU上运行的“Hello World”程序。它与CPU版本有三个关键区别。
#include <stdio.h>
#include <cuda_runtime.h> // 区别1：包含CUDA头文件
// 区别2：使用 __global__ 关键字声明GPU内核函数
__global__ void d_kernel() {
printf("Hello World\n");
}
int main() {
// 区别3：使用特殊语法启动内核，<<<1, 1>>> 表示用1个线程块，每个块1个线程
d_kernel<<<1, 1>>>();
cudaDeviceSynchronize(); // 等待GPU内核执行完成
return 0;
}
关键概念解析
- 头文件：
#include <cuda_runtime.h>是编写CUDA代码所必需的，它包含了CUDA相关的函数和类型定义。 - 内核函数：
__global__关键字告诉编译器，d_kernel是一个GPU内核函数，它将在GPU上执行，而不是在CPU上。 - 内核启动：
d_kernel<<<1, 1>>>()是启动内核的语法。<<<1, 1>>>中的两个参数定义了线程的配置：- 第一个
1表示线程块（Block）的数量。 - 第二个
1表示每个线程块中的线程（Thread）数量。 - 因此，
<<<1, 1>>>表示启动一个包含1个线程的线程块，即总共只有1个线程执行该内核。
- 第一个
注意：如果省略 cudaDeviceSynchronize()，程序可能不会输出任何内容。这是因为CPU在启动内核后不会等待GPU完成，而是继续执行并可能提前结束程序。cudaDeviceSynchronize() 的作用是让CPU等待所有GPU线程执行完毕。
多个内核的调用顺序
上一节我们介绍了单个内核的启动。本节中我们来看看当程序中有多个内核调用时，它们的执行顺序是怎样的。
__global__ void d_kernel() {
printf("Hello World\n");
}
int main() {
d_kernel<<<1, 1>>>();
d_kernel<<<1, 1>>>();
d_kernel<<<1, 1>>>();
cudaDeviceSynchronize();
printf("on CPU\n");
return 0;
}
执行流程分析
- CPU按顺序将三个内核调用放入一个默认的队列（称为“流”，Stream）中。
- GPU从这个队列中按先进先出（FIFO）的顺序取出并执行内核。
cudaDeviceSynchronize()确保CPU等待队列中所有内核执行完毕。- 最后，CPU执行
printf(“on CPU\n”)。
因此，输出结果是确定性的：三行“Hello World”，然后是一行“on CPU”。内核的执行是串行的（一个接一个），但这是GPU内部的串行，CPU在启动它们后就去等待了。
CPU与GPU的异步执行
我们调整一下代码顺序，观察CPU和GPU异步执行的影响。
int main() {
d_kernel<<<1, 1>>>();
d_kernel<<<1, 1>>>();
printf("on CPU\n");
cudaDeviceSynchronize(); // 注意：同步调用放到了打印语句之后
return 0;
}
可能的输出结果
由于CPU和GPU是异步执行的，输出顺序变得不确定。以下是几种可能的情况：
on CPU->Hello World->Hello WorldHello World->on CPU->Hello WorldHello World->Hello World->on CPU
关键点：
- CPU启动内核后立即继续执行后面的
printf语句。 - 两个内核在默认流中仍按顺序执行（
Hello World总是先于第二个Hello World打印）。 cudaDeviceSynchronize()保证了两个内核最终都会执行完毕，但无法控制它相对于CPUprintf的执行时机。
混合CPU与GPU打印任务
让我们看一个更复杂的例子，其中混合了CPU和GPU的打印任务。
int main() {
d_kernel<<<1, 1>>>();
printf("CPU1\n");
d_kernel<<<1, 1>>>();
printf("CPU2\n");
d_kernel<<<1, 1>>>();
printf("CPU3\n");
cudaDeviceSynchronize();
printf("on CPU\n");
return 0;
}
并行性分析
在这个例子中，以下操作可以并行执行：
- GPU执行三个
d_kernel打印“Hello World”。 - CPU执行三个
printf打印“CPU1”、“CPU2”、“CPU3”。
因此，你可能会看到“CPU1/2/3”与“Hello World”交错打印的输出。然而，由于 cudaDeviceSynchronize() 的存在，最后的 ”on CPU” 一定会在所有GPU内核执行完毕后才打印。
一个常见的输出顺序是：CPU1 CPU2 CPU3 -> 三个 Hello World -> on CPU。这是因为CPU执行速度通常很快，在GPU完成大量工作前就可能完成了自己的打印任务。
使用多线程并行打印
前面的例子每个内核只使用一个线程。CUDA的强大之处在于大规模并行。以下是如何用多个线程打印“Hello World”。
__global__ void d_kernel() {
printf("Hello World\n");
}
int main() {
// 启动1个线程块，该块中包含32个线程
d_kernel<<<1, 32>>>();
cudaDeviceSynchronize();
return 0;
}
编译并运行此程序，将会看到“Hello World”被打印了32次，每个线程执行一次 printf 语句。
线程配置公式：<<<NumBlocks, ThreadsPerBlock>>>
NumBlocks：线程网格（Grid）中线程块的数量。ThreadsPerBlock：每个线程块中的线程数量。- 总线程数 =
NumBlocks * ThreadsPerBlock
并行计算示例：打印平方数
上一节我们介绍了如何用多线程执行相同任务。本节中我们来看看如何让每个线程处理不同的数据。假设我们想并行计算并打印0到99的平方。
思路分析
目标是让线程i计算并打印 i * i。我们需要一种方法让每个线程知道自己独特的索引（ID）。
__global__ void square_kernel() {
int tid = threadIdx.x; // 获取当前线程在线程块内的索引（0到ThreadsPerBlock-1）
printf(“Square of %d is %d\n”, tid, tid * tid);
}
int main() {
int n = 100;
// 启动1个线程块，包含100个线程
square_kernel<<<1, n>>>();
cudaDeviceSynchronize();
printf(“on CPU\n”);
return 0;
}
代码解释：
threadIdx.x是一个内置变量，表示线程在其所属线程块中的x维度索引。- 由于我们只使用了一个线程块（
<<<1, n>>>），threadIdx.x的范围是0到n-1。 - 每个线程根据自己独有的
tid计算平方并打印。
注意：由于所有线程是并行执行的，打印输出的顺序可能是混乱的（即非从0到99的顺序）。这是并行编程中常见的现象，如果需要顺序输出，则需要额外的同步或排序操作。
核心概念回顾
在深入代码之前，我们先简要回顾一些核心的计算机系统概念，这有助于理解CUDA的线程模型。
- 进程（Process）：正在执行的程序实例。它拥有独立的地址空间，包含代码、数据、堆栈等资源。例如，打开的Word应用程序就是一个进程。
- 线程（Thread）：进程内的一个执行流。它是“轻量级进程”，与同进程的其他线程共享代码和全局数据，但拥有独立的栈和寄存器。例如，Word中的拼写检查功能可以是一个独立的线程。
- 操作系统（OS）：管理硬件资源和为用户程序提供服务的软件。例如，Linux， Windows。
- 硬件组件：
- 缓存（Cache）：靠近处理器的小容量高速内存，用于加速数据访问。
- 主存（RAM）：容量更大的内存，用于存储正在运行的程序和数据。
- 核心（Core）：执行线程的物理处理单元。一个核心同一时刻只能执行一个线程。线程可能会在不同的核心之间调度（“跳跃”）。
在CUDA中，我们创建成千上万个线程，它们被组织成线程块和网格，在GPU的众多核心上并行执行。
总结
本节课我们一起学习了CUDA并行编程的基础。我们从最简单的“Hello World”程序出发，理解了CUDA内核函数（用 __global__ 修饰）的概念以及如何用 <<<…>>> 语法启动内核。我们探讨了CPU与GPU的异步执行特性，并学习了使用 cudaDeviceSynchronize() 进行同步。通过示例，我们看到了如何配置多线程（<<<NumBlocks, ThreadsPerBlock>>>）来并行执行任务，以及如何利用 threadIdx.x 让每个线程处理不同的数据。这些概念是构建更复杂、高性能CUDA应用程序的基石。
002：CUDA程序流程与cudaMemcpy操作 🚀
概述
在本节课中，我们将学习CUDA程序的基本流程，特别是CPU与GPU之间如何通过cudaMemcpy函数进行数据传输。我们将通过具体的代码示例，理解如何为GPU分配内存、启动内核以及将结果复制回CPU。
回顾与引入
上一节我们介绍了如何编写和启动一个简单的CUDA内核，并使用多个线程并行打印数字的平方。本节中，我们来看看如何将计算结果存储到数组中，并将其从GPU传回CPU。
我们将从一个简单的C语言串行代码开始，目标是将其并行化。
// 串行C代码示例
for (int i = 0; i < 100; i++) {
a[i] = i * i;
}
以下是实现相同功能的CUDA并行代码。
#include <stdio.h>
__global__ void square(int *a) {
int tid = threadIdx.x;
a[tid] = tid * tid;
}
int main() {
int n = 100;
int a[n];
int *d_a;
cudaMalloc((void**)&d_a, n * sizeof(int));
square<<<1, n>>>(d_a);
cudaMemcpy(a, d_a, n * sizeof(int), cudaMemcpyDeviceToHost);
for (int i = 0; i < n; i++) {
printf("%d ", a[i]);
}
cudaFree(d_a);
return 0;
}
代码解析
以下是上述代码的关键步骤解析：
- 主机（CPU）内存分配：在CPU上声明数组
a。 - 设备（GPU）内存分配：使用
cudaMalloc在GPU上分配内存，d_a是一个位于CPU上的指针，但它指向GPU上的内存地址。 - 内核启动：使用
<<<1, n>>>语法启动square内核，其中n是线程数。每个线程计算其threadIdx.x的平方，并写入GPU数组d_a的对应位置。 - 数据回传：使用
cudaMemcpy将计算结果从GPU内存 (d_a) 复制回CPU内存 (a)。第四个参数cudaMemcpyDeviceToHost指明了传输方向。 - 内存释放：使用
cudaFree释放GPU上分配的内存。
核心概念与注意事项
在深入更多示例前，我们需要理解几个核心概念：
- CPU与GPU内存分离：CPU和GPU拥有独立的内存空间（显存）。在离散GPU上，CPU无法直接访问GPU内存，反之亦然。因此，必须显式地进行内存分配和数据拷贝。
cudaMemcpy是阻塞操作：cudaMemcpy调用会阻塞CPU线程，直到所有先前启动的内核执行完毕且数据传输完成。因此，在这个例子中，我们不需要显式调用cudaDeviceSynchronize()。- PCIe总线瓶颈：CPU与GPU之间的数据传输通过PCI Express总线进行，其带宽有限。频繁或大量数据传输可能成为程序性能瓶颈。因此，应确保在GPU上进行足够的计算，以抵消数据传输的开销。
CPU/GPU内存访问示例
为了加深对内存分离的理解，我们来看几个例子。
示例1：错误的全局变量访问
以下代码尝试在内核中访问一个在CPU上定义的全局变量字符串指针，这将导致编译错误。
char *msg = "Hello World"; // 在CPU内存中
__global__ void hello() {
printf("%s\n", msg); // 错误！GPU内核无法直接访问CPU内存
}
int main() {
hello<<<1, 32>>>();
cudaDeviceSynchronize();
return 0;
}
关键点：在离散GPU上，定义在主机（CPU）上的变量不能被设备（GPU）代码直接访问。
示例2：使用预处理器宏
如果使用预处理器宏 #define，代码则可以正常工作。
#define MSG "Hello World" // 预处理器指令，在编译前进行文本替换
__global__ void hello() {
printf("%s\n", MSG); // 正确！编译前MSG已被替换为字符串字面量
}
int main() {
hello<<<1, 32>>>();
cudaDeviceSynchronize();
return 0;
}
关键点：#define 是预处理器指令，它在编译之前将 MSG 直接替换为 "Hello World" 字符串字面量。对于内核来说，它只是在打印一个常量字符串，而非访问CPU内存地址。
典型的CUDA程序流程 📊
一个标准的CUDA程序通常遵循以下步骤，我们将其与所需的CUDA API对应起来：
- 从文件系统加载数据到CPU内存：使用标准C/C++文件操作（如
fread）。 - 将数据从CPU内存传输到GPU内存：
- 使用
cudaMalloc在GPU上分配内存。 - 使用
cudaMemcpy(..., cudaMemcpyHostToDevice)将数据从主机复制到设备。
- 使用
- 在GPU上执行内核：使用
kernel_name<<<grid_size, block_size>>>(parameters)启动内核。 - 将结果从GPU内存传输回CPU内存：使用
cudaMemcpy(..., cudaMemcpyDeviceToHost)。 - 在CPU上使用或输出结果：使用
printf打印或fwrite写回文件。
重要提示：程序员有责任维护CPU和GPU上数据的一致性。通常，建议使用不同的变量名来区分它们（例如，h_a 表示主机数组，d_a 表示设备数组），以避免错误传递指针。
实践练习：字符串操作与边界问题
让我们通过一个操作字符串的练习，来学习如何正确处理GPU内存和线程索引。
任务描述
我们有一个字符串 "GdkknVnqkc"，每个字符在字母表中都是前一个字符（例如，‘H’的前一个字符是‘G’）。我们想用GPU线程并行地将每个字符递增，最终得到 "HelloWorld"。
有问题的初始代码
以下代码存在一个常见的边界错误。
__global__ void decode(char *a, int len) {
int tid = threadIdx.x;
if (tid < len) {
a[tid]++; // 所有线程都对字符进行加一操作
}
}
int main() {
char cpu_a[] = "GdkknVnqkc";
int len = strlen(cpu_a);
char *gpu_a;
cudaMalloc((void**)&gpu_a, (len + 1) * sizeof(char)); // 分配len+1字节，包含'\0'
cudaMemcpy(gpu_a, cpu_a, (len + 1) * sizeof(char), cudaMemcpyHostToDevice);
decode<<<1, len + 1>>>(gpu_a, len + 1); // 启动了 len+1 个线程
cudaMemcpy(cpu_a, gpu_a, (len + 1) * sizeof(char), cudaMemcpyDeviceToHost);
printf("%s\n", cpu_a);
cudaFree(gpu_a);
return 0;
}
问题分析：字符串 "GdkknVnqkc" 长度为10，但C风格字符串以 \0 结尾，所以 cpu_a 实际占用11字节。代码分配了11字节，并启动了11个线程。第11个线程（tid=10）操作的是字符串终止符 \0，将其递增为 \1，导致字符串失去终止符。后续 printf 可能打印乱码或导致段错误。
正确的代码
正确的做法是只对有效字符进行操作，不改变终止符。
decode<<<1, len>>>(gpu_a, len); // 只启动 len 个线程，不处理终止符
关键点：在并行处理数组（尤其是字符串）时，必须仔细管理线程索引与数据边界，确保不会意外修改或访问越界。
线程组织与限制 ⚙️
最后，我们简要探讨一下CUDA的线程组织。目前我们只使用了 <<<1, n>>> 中的第二个参数，它定义了一个线程块（Thread Block）中的线程数量。
- 线程块大小限制：每个线程块的最大线程数通常是1024（取决于GPU架构）。这意味着，使用
<<<1, n>>>这种格式，n不能超过1024。 - 如何启动更多线程：为了启动超过1024个线程（例如，处理8000个元素的数组），我们需要使用第一个参数来指定网格（Grid）中线程块的数量。这是下节课的重点。
例如，以下调用试图启动8000个线程，但可能会失败，因为单个线程块无法容纳8000个线程。
kernel<<<1, 8000>>>(); // 可能失败，超出线程块最大限制
我们将学习如何使用网格和线程块来灵活地组织大量线程。
总结
本节课中我们一起学习了：
- CUDA程序基本流程：包括主机与设备内存分配、数据传输和内核执行。
cudaMemcpy函数：用于在CPU和GPU之间复制数据，方向由cudaMemcpyHostToDevice和cudaMemcpyDeviceToHost控制。- CPU与GPU内存分离：这是CUDA编程的核心概念，必须为需在GPU上处理的数据在设备端创建副本。
- 数据边界管理：在线程化代码中，必须谨慎处理数组索引，防止越界访问。
- 线程组织初步：了解了线程块的概念及其线程数量限制，为学习网格化线程组织打下基础。
掌握这些基础是编写正确、高效CUDA程序的关键。下一节，我们将深入探讨CUDA的线程层次结构（网格和线程块），以支持大规模并行计算。
003：CUDA中的线程组织与线程块 🧵
概述
在本节课中，我们将深入学习CUDA中的线程组织模型。我们将探讨线程网格、线程块和线程的层次结构，理解如何通过多维索引来唯一标识每个线程，并通过实例学习如何为一维和二维数据结构编写并行初始化代码。
线程组织的层次结构
在上一讲中，我们学习了典型的CUDA程序流程。本节中，我们来看看CUDA如何组织海量的线程来执行并行计算。
当一个内核被启动时，它会以一个线程网格的形式执行。这个网格是一个三维数据结构，可以将其想象成一个魔方。整个魔方就是网格，而构成魔方的每一个小立方体就是一个线程块。
- 网格：是内核启动的最高层次组织，由线程块构成。
- 线程块：是网格的组成部分，内部包含线程。
- 线程：是实际执行计算的最小单元，位于线程块内部。
网格和线程块都可以是一维、二维或三维的，这为处理不同维度的数据（如图像、矩阵、体积数据）提供了灵活性。
标识线程与线程块
为了在并行计算中准确定位每个线程，我们需要一套坐标系统。
网格维度与块索引
网格的尺寸（即每个维度上有多少个线程块）可以通过内置变量访问：
gridDim.x：网格在x方向的尺寸（线程块数量）。gridDim.y：网格在y方向的尺寸。gridDim.z：网格在z方向的尺寸。
网格中每个线程块的位置由块索引唯一确定：
blockIdx.x：当前线程块在网格x方向上的索引。blockIdx.y：当前线程块在网格y方向上的索引。blockIdx.z：当前线程块在网格z方向上的索引。
索引范围：如果gridDim.x = 10，那么blockIdx.x的取值范围是 0 到 9。
线程块维度与线程索引
每个线程块的尺寸（即每个维度上有多少个线程）可以通过内置变量访问：
blockDim.x：线程块在x方向的尺寸（线程数量）。blockDim.y：线程块在y方向的尺寸。blockDim.z：线程块在z方向的尺寸。
线程块中每个线程的位置由线程索引唯一确定：
threadIdx.x：当前线程在线程块x方向上的索引。threadIdx.y：当前线程在线程块y方向上的索引。threadIdx.z：当前线程在线程块z方向上的索引。
索引范围：如果blockDim.x = 5，那么threadIdx.x的取值范围是 0 到 4。
内核启动配置示例
以下代码展示了如何配置并启动一个多维的内核：
dim3 grid(2, 3, 4); // 定义一个2x3x4的网格（共24个线程块）
dim3 block(5, 6, 7); // 定义一个5x6x7的线程块（共210个线程）
// 启动内核，共创建 24 * 210 = 5040 个线程
myKernel<<<grid, block>>>();
cudaDeviceSynchronize();
在内核函数myKernel内部，我们可以使用上述变量。例如，以下条件只会在一个特定的线程（位于第一个线程块的第一个线程）中为真：
if (blockIdx.x == 0 && blockIdx.y == 0 && blockIdx.z == 0 &&
threadIdx.x == 0 && threadIdx.y == 0 && threadIdx.z == 0) {
// 只有这个线程会执行这里的代码
}
应用实例：计算全局唯一线程ID
在实际编程中，我们经常需要为每个线程计算一个在整个网格范围内唯一的ID，以便处理像数组这样的线性内存。
情况一：单线程块，二维线程组织
假设我们有一个大小为 N * M 的一维数组，我们使用一个二维线程块（N x M 个线程）来初始化它。
内核启动配置：
#define N 5
#define M 6
dim3 block(N, M); // 5x6的线程块
myKernel<<<1, block>>>(matrix); // 网格中只有1个线程块
在内核中，计算唯一ID的公式为：
int i = threadIdx.x * blockDim.y + threadIdx.y;
matrix[i] = i; // 为每个位置赋予唯一值
公式解释：threadIdx.x 是行索引，blockDim.y（固定为6）是总列数，threadIdx.y 是列索引。该公式将二维索引映射到一维数组。
情况二：多线程块，一维线程组织
同样初始化一个 N * M 的数组，现在我们使用多个一维线程块。每个线程块包含 M 个线程，总共启动 N 个线程块。
内核启动配置：
#define N 5
#define M 6
myKernel<<<N, M>>>(matrix); // N个块，每个块M个线程
在内核中，计算唯一ID的公式变为：
int i = blockIdx.x * blockDim.x + threadIdx.x;
matrix[i] = i;
公式解释：blockIdx.x 是线程块索引，blockDim.x（固定为6）是每个块的线程数，threadIdx.x 是块内线程索引。该公式结合了块索引和块内线程索引来生成全局唯一ID。
这两种方法都能成功地将数组初始化为 [0, 1, 2, ..., 29]，展示了组织线程的灵活性。
总结
本节课我们一起学习了CUDA中线程组织的核心概念：
- CUDA执行模型是分层的：网格(Grid) -> 线程块(Block) -> 线程(Thread)。
- 网格和线程块都可以配置为一维、二维或三维结构，以适应不同维度的数据。
- 使用
blockIdx、threadIdx、gridDim、blockDim这些内置变量可以唯一标识每个线程并获取执行配置信息。 - 关键在于计算全局唯一的线程ID，其公式取决于网格和线程块的维度配置。我们通过两个初始化数组的实例演示了如何推导和应用这个公式。
理解线程组织是编写高效CUDA并行程序的基础，它直接决定了任务如何被分解并映射到GPU的数千个核心上执行。
004：GPU计算层次与线程发散
在本节课中，我们将学习如何为大规模数据启动CUDA内核，深入理解GPU的计算层次结构，并探讨线程发散问题及其对性能的影响。
概述
上一节我们介绍了线程在GPU上的组织方式，包括网格、线程块以及线程ID的计算。本节中，我们将看看如何为大规模数据启动内核，理解GPU的底层执行单元——线程束，并分析条件分支如何导致线程发散，从而影响性能。
为大规模数据启动内核
当数据规模巨大时，我们需要计算合适的线程块数量来启动内核。以下是一个示例，展示了如何根据用户指定的数据大小动态计算所需的线程块数量。
#include <stdio.h>
#include <cuda.h>
#define BLOCK_SIZE 1024
__global__ void initVectorKernel(unsigned int *vector, unsigned int vector_size) {
unsigned int id = blockIdx.x * blockDim.x + threadIdx.x;
if (id < vector_size) {
vector[id] = id;
}
}
int main(int argc, char *argv[]) {
unsigned int n = atoi(argv[1]);
unsigned int *d_vector, *h_vector;
cudaMalloc(&d_vector, n * sizeof(unsigned int));
h_vector = (unsigned int*)malloc(n * sizeof(unsigned int));
unsigned int num_blocks = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
printf("Number of blocks: %u\n", num_blocks);
initVectorKernel<<<num_blocks, BLOCK_SIZE>>>(d_vector, n);
cudaMemcpy(h_vector, d_vector, n * sizeof(unsigned int), cudaMemcpyDeviceToHost);
for (int i = 0; i < n; i++) {
printf("%u ", h_vector[i]);
}
printf("\n");
free(h_vector);
cudaFree(d_vector);
return 0;
}
以下是代码中需要解决的两个关键问题及其解决方案：
- 整数除法问题：使用
(n + BLOCK_SIZE - 1) / BLOCK_SIZE代替n / BLOCK_SIZE，确保当n不能被BLOCK_SIZE整除时，能分配足够的线程块。 - 内存越界访问：在内核函数中添加条件判断
if (id < vector_size)，确保只有有效的线程ID才会执行内存写入操作，防止访问非法内存地址。
矩阵平方应用：从CPU到GPU的并行化
我们以一个矩阵平方运算为例，展示如何逐步将其从CPU版本并行化到GPU上，并分析性能变化。
CPU版本
CPU版本的矩阵平方使用三层嵌套循环，时间复杂度为 O(n³)。
void matrixSquareCPU(float *matrix, float *result, int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
float sum = 0.0f;
for (int k = 0; k < n; k++) {
sum += matrix[i * n + k] * matrix[k * n + j];
}
result[i * n + j] = sum;
}
}
}
对于一个64x64的矩阵，CPU执行时间约为 1.52毫秒。
GPU版本1：仅并行化外层循环
在第一个GPU版本中，我们仅并行化了最外层的 i 循环。每个线程负责计算结果矩阵的一整行。
__global__ void matrixSquareGPUv1(float *matrix, float *result, int n) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) {
for (int j = 0; j < n; j++) {
float sum = 0.0f;
for (int k = 0; k < n; k++) {
sum += matrix[i * n + k] * matrix[k * n + j];
}
result[i * n + j] = sum;
}
}
}
// 内核启动: <<<1, n>>>
这个版本存在两个问题：
- 每个线程内部仍有大量串行工作（
j和k循环）。 - 并行度仅为
n（例如64），远未充分利用GPU的数千个核心。
因此，其执行时间（6.39毫秒）甚至比CPU版本更慢。
GPU版本2：并行化两个外层循环
为了提升并行度，我们同时并行化 i 和 j 循环。每个线程负责计算结果矩阵中的一个单独元素。
__global__ void matrixSquareGPUv2(float *matrix, float *result, int n) {
int id = blockIdx.x * blockDim.x + threadIdx.x;
int i = id / n;
int j = id % n;
if (i < n && j < n) {
float sum = 0.0f;
for (int k = 0; k < n; k++) {
sum += matrix[i * n + k] * matrix[k * n + j];
}
result[i * n + j] = sum;
}
}
// 内核启动: <<<(n*n + BLOCK_SIZE -1)/BLOCK_SIZE, BLOCK_SIZE>>>
这个版本取得了显著改进：
- 并行度：提升到
n * n（例如4096）。 - 串行工作：每个线程仅剩一个
k循环。
其执行时间大幅降低至 0.1毫秒，优于CPU版本。
注意：最内层的 k 循环涉及对同一内存位置（sum）的累加，如果直接并行化会导致数据竞争，需要同步机制，这将在后续课程中讨论。
GPU计算层次结构
理解GPU的硬件执行模型对编写高效CUDA程序至关重要。其计算层次从高到低如下：
- 线程：最小的执行单元，类似于一个“工人”。
- 线程束：由 32个线程 组成的基本执行单元。线程束内的所有线程以 SIMD 方式执行，即单指令，多数据。它们步调一致地执行相同的指令，但操作不同的数据。
- 线程块：由多个线程束组成。一个线程块最多可包含1024个线程。
- 块内同步：线程块内的线程可以通过
__syncthreads()函数进行同步。 - 共享内存：线程块拥有一个高速的、块内线程可共享的片上内存（Shared Memory），能显著减少全局内存访问延迟。
- 块内同步：线程块内的线程可以通过
- 流式多处理器：GPU上的一个处理单元，可同时管理多个线程块（通常1-8个）。一个SM包含大量核心（例如128个）。
- GPU设备：由多个SM组成，可同时运行数万甚至数十万个线程。
关键点：锁步执行仅发生在线程束内部。不同线程束之间的执行是异步的。
线程发散
线程发散是GPU编程中一个重要的性能考量因素。它发生在一个线程束内的线程需要执行不同的指令路径时。
发散的原因与影响
由于线程束以SIMD方式执行，当遇到条件分支（如 if-else）时，如果线程束内有些线程满足条件，有些不满足，GPU会先执行满足条件分支的线程，其他线程等待（执行空操作）；然后再执行另一分支的线程。这导致了串行化执行，降低了效率。
__global__ void divergentKernel(int *data) {
int id = threadIdx.x;
if (id % 2 == 0) {
data[id] *= 2; // 偶数线程执行
} else {
data[id] += 1; // 奇数线程执行
}
}
在上面的内核中，一个线程束内的奇偶线程将分两批执行，造成发散。
何时会发生发散？
导致发散的关键不是条件语句本身，而是条件在同一个线程束内是否对所有线程评估出相同的结果。
- 不会发散的例子：
if (threadIdx.x < 32)。对于一个线程束（0-31），条件对所有线程都为真，因此没有发散。 - 会发散的例子：循环次数依赖于线程ID，
for(int i=0; i<threadIdx.x; i++)。线程束内不同线程的循环次数不同，最后一个线程（迭代31次）将决定整个线程束的执行时间，导致其他线程等待，引发发散。
避免发散的一个技巧
在某些特定情况下，可以重写代码来消除条件分支。
原始代码（可能发散）：
if (x == y) {
x = z;
} else if (x == z) {
x = y;
}
优化后代码（无分支）：
x = (x == y) * z + (x == z) * y;
// 或者更通用的： x = z * (x == y) + y * (x == z);
优化后的代码通过逻辑运算和算术运算替代了条件判断，使得线程束内所有线程执行相同的指令流，从而避免了发散。但这种方法并非在所有情况下都适用。
总结
本节课中我们一起学习了：
- 如何为大规模数据计算并启动合适数量的线程块，并注意处理整数除法和内存越界问题。
- 通过矩阵平方的例子，实践了将CPU代码并行化到GPU的逐步优化过程，认识到提高并行度和减少线程内串行工作的重要性。
- 深入理解了GPU的计算层次结构，特别是线程束的SIMD执行模式。
- 分析了线程发散问题的成因、影响以及避免方法。核心在于意识到：线程束内执行路径的分化是性能杀手，编写代码时应尽量让同一个线程束内的线程执行相同的指令流。
掌握这些概念是编写高效CUDA程序的基础。下一节课，我们将探讨GPU的内存模型。
005：CUDA中的线程发散（续）与内存模型概述
在本节课中，我们将深入学习CUDA编程中的线程发散问题，并通过具体示例探讨如何优化代码以减少发散。随后，我们将从计算部分转向内存部分，初步了解GPU的内存层次结构。
回顾与概述
上一节我们介绍了GPU的计算层次结构、典型的矩阵平方应用并行化方法，以及处理大规模数据时内核启动的常见问题。我们还深入探讨了线程束（Warp）的概念、执行方式（锁步、SIMD），并初步了解了线程发散（Thread Divergence）问题。
本节我们将继续深入线程发散问题，通过更多示例学习如何识别和缓解发散。最后，我们将开始新的主题——GPU内存，了解其基本组织方式。
线程发散深入分析
示例一：多分支条件语句
首先，我们回顾上节课结束时的一个简单示例。以下内核启动了一个包含32个线程（即一个线程束）的核函数：
__global__ void D_kernel(float* vector, int vector_size) {
int id = threadIdx.x + blockIdx.x * blockDim.x;
if (id < vector_size) {
if (id % 2 == 1) {
// 语句 S1 (奇数ID线程执行)
} else {
// 语句 S2 (偶数ID线程执行)
}
// 语句 S3 (所有线程执行)
}
}
关键观察点：对于这个 if-else 结构，无论线程ID是奇数还是偶数，线程束执行它都需要 2步。第一步，一部分线程执行 S1，另一部分闲置；第二步，另一部分线程执行 S2，之前执行的线程闲置。虽然存在发散，但执行步数是固定的。
练习：复杂条件分支
现在考虑一个更复杂的条件分支：
if (x == y) {
a = z + w;
} else if (x == z) {
a = y + w;
} else {
a = y + z;
}
问题：
- 线程束执行此
if-else-if结构最多需要多少步？ - 如何消除此处的线程发散？
分析与解答：
- 执行步数：在最坏情况下，线程束中的线程可能分别进入三个不同的条件分支。因此，线程束需要按顺序让进入每个分支的线程组依次执行，总共需要 3步。
- 消除发散：我们可以利用数学等价替换，将条件分支转换为一个无分支的表达式，从而让所有线程执行相同的操作：
a = (x == y) * (z + w) + (x == z) * (y + w) + (x == w) * (y + z);
但更优雅的解法是利用这个公式：a = (y * (x == z) + z * (x == y) + w * (x == w)) / ((x == y) + (x == z) + (x == w));
当x等于y、z或w时，此公式能正确计算出对应的a值，且所有线程执行相同的计算，只需1步，彻底消除了发散。
示例二：包含循环的昂贵分支
考虑以下内核，它同样由32个线程（一个线程束）执行：
__global__ void T_kernel(float* vector, int vector_size) {
int id = threadIdx.x + blockIdx.x * blockDim.x;
if (id <= 0) { // 仅线程ID为0满足条件
vector[0] = 0;
for(int i = 1; i <= 100; i++) {
vector[0] += i; // 计算1到100的和
}
} else {
vector[id] = 1;
}
}
问题：线程束执行这个 if-else 结构需要多少步？
分析与解答：
- 只有
id=0的线程会进入if分支，执行vector[0] = 0;（1步）和100次循环（100步），共101步。 - 其余31个线程 (
id=1到31) 在if阶段闲置。 - 然后，这31个线程进入
else分支，执行vector[id] = 1;（1步），而id=0的线程在此阶段闲置。 - 因此，线程束总共需要 101 + 1 = 102步 来执行此结构。性能瓶颈在于
id=0线程的昂贵循环。
优化：我们可以将循环计算（求和）替换为闭合公式（等差数列求和公式）：
if (id <= 0) { vector[0] = 100 * 101 / 2; }
这样，if 分支的耗时从101步降为1步。优化后，线程束执行整个 if-else 结构只需要 2步（if分支1步，else分支1步）。
进一步优化以消除发散：虽然步数减为2，但线程依然发散（一部分执行if，另一部分执行else）。我们可以利用整数运算技巧，将其合并为单一步骤：
vector[id] = 1 + ( ( -(id >> 31) ) & (100*101/2 - 1) );
这个表达式利用了整数右移和位操作，使得当 id=0 时，结果为 100*101/2；当 id>0 时，结果为 1。所有线程执行相同的指令序列，仅需1步，完全消除了发散。
示例三：Switch-Case 语句
考虑以下使用 switch-case 的内核：
__global__ void D_kernel(float* vector, int vector_size) {
int id = threadIdx.x + blockIdx.x * blockDim.x;
switch(id) {
case 0: vector[0] = 0; break;
case 1: vector[1] = vector[0]; break;
case 2: vector[2] = vector[id-2]; break;
case 3: vector[3] = vector[6]; break;
case 4: vector[4] = vector[5] + 4; break;
case 5: vector[5] = vector[4] - 4; break;
case 6: vector[6] = vector[3]; break;
case 7: vector[7] = vector[7] * 2; break;
case 8: vector[8] = vector[4] + id; break;
case 9: vector[9] = 9; break;
}
}
问题：线程束执行此 switch 结构需要多少步？
分析与解答：
switch 语句在逻辑上等价于一系列的 if-else if 判断。线程束需要让进入每个 case 的线程组依次执行。因为有10个 case，所以在最坏情况下需要 10步。
优化：通过观察，我们可以合并一些执行相同或等价操作的 case，从而减少步数：
case 0和case 2：由于线程锁步执行，当id=2执行vector[2] = vector[0];时，vector[0]已被id=0的线程初始化为0。因此两者效果可合并（将vector对应位置赋0）。case 3和case 6：都执行vector[3] = vector[6];和vector[6] = vector[3];（具体值取决于执行顺序，但操作对称）。case 4和case 8：case 8执行vector[8] = vector[4] + 8;，这与case 4的vector[4] = vector[5] + 4;在数据流上相关，在某些条件下可考虑合并逻辑。
通过合并，可以将执行步数从10步减少到7步或更少。
示例四：屏障同步与死锁
考虑以下包含全局屏障的内核：
__global__ void example_kernel() {
int id = threadIdx.x + blockIdx.x * blockDim.x;
if (id < 16) {
printf("Inside if\n");
__syncthreads(); // 全局屏障
} else {
printf("Inside else\n");
}
}
问题：这个程序的输出是什么？会发生什么？
分析与解答：
- 假设内核启动了一个包含32个线程的线程束。
- 前16个线程 (
id=0到15) 进入if分支，打印 “Inside if”，然后到达__syncthreads()屏障并等待。 - 后16个线程 (
id=16到31) 进入else分支，打印 “Inside else”，但永远不会到达__syncthreads()屏障。 - 因此，先到达屏障的16个线程将无限期等待另外16个线程到达屏障，而这是永远不会发生的。这就导致了 死锁（Deadlock）。
关键点：线程束内的线程如果因为条件分支而无法全部到达同步点，就会导致死锁。这是编写CUDA代码时需要特别注意的问题。
GPU内存模型概述
前面我们主要关注了GPU的计算方面。现在，让我们开始了解GPU的另一个核心方面——内存。
GPU拥有复杂的内存层次结构，设计目标是隐藏全局内存的高延迟并提高数据访问带宽。以下是其主要的层次：
线程 (Thread) -> 寄存器 (Registers) [最快，每个线程私有]
|
V
线程块 (Block) -> 共享内存 (Shared Memory) [较快，块内线程共享]
|
V
全局内存 (Global Memory) / 设备内存 (Device Memory) [较慢，所有线程可访问]
|
V
纹理内存 (Texture Memory) 和 常量内存 (Constant Memory) [特殊用途，只读]
各级内存详解
- 寄存器 (Registers)
- 位置与速度：位于芯片上，是速度最快的内存。
- 作用域：每个线程私有。编译器会尽可能将局部变量分配到寄存器中。
- 容量：数量有限。如果寄存器不足，变量会“溢出”到更慢的本地或全局内存。
-
共享内存 (Shared Memory)
- 位置与速度：位于芯片上，速度仅次于寄存器，但远快于全局内存。
- 作用域：同一个线程块内的所有线程共享。是块内线程通信和协作的主要桥梁。
- 用途：常用于存储需要被块内线程频繁访问的数据片段，例如矩阵乘法中的瓦片（Tile）。
-
全局内存 (Global Memory)
- 位置与速度：即GPU的板载DRAM（设备内存）。容量大（数GB），但延迟非常高（数百个时钟周期）。
- 作用域：对所有网格中的所有线程可见。是主机（CPU）与设备（GPU）之间传输数据的主要区域，也是不同线程块之间通信的媒介（尽管不直接同步）。
- 带宽：现代GPU的全局内存带宽很高（可达数百GB/s），但高延迟仍是主要性能瓶颈。
-
纹理内存 (Texture Memory) 和 常量内存 (Constant Memory)
- 特性：它们是全局内存中经过特殊优化、具有独立缓存的空间。
- 只读性：内核通常只能读取它们的内容。
- 纹理内存：针对具有空间局部性的访问模式（如图像处理中的邻近像素访问）进行了优化。
- 常量内存：适合存储所有线程都需要读取的常量数据。当所有线程读取相同地址时，性能极高（广播机制）。
- 容量：很小（常量内存通常64KB，纹理内存缓存也有限）。
性能启示
编程时应遵循以下原则以提升性能：
- 优先使用寄存器：尽量使用局部变量。
- 合理利用共享内存：将需要重复访问的数据从全局内存缓存到共享内存。
- 优化全局内存访问：通过合并访问（Coalesced Access）来最大化内存带宽利用率（后续课程详述）。
- 善用常量与纹理内存：将只读的常量或具有空间局部性的数据放置于此。
总结
本节课我们深入探讨了CUDA中的线程发散问题。通过多个示例，我们学习了：
- 如何分析条件分支、循环和
switch语句导致的线程束执行步数。 - 如何运用数学等价变换和位操作等技巧，将发散代码重构为无分支代码，从而显著提升性能，甚至完全消除发散。
- 线程发散在同步点（如
__syncthreads()）可能引发死锁，需要谨慎处理。
最后，我们开启了GPU内存主题的学习，概述了GPU的内存层次结构，包括寄存器、共享内存、全局内存、纹理内存和常量内存的特性与用途。理解这个层次结构对于编写高性能CUDA程序至关重要。
在接下来的课程中，我们将更深入地研究内存访问模式、共享内存的使用以及如何实现高效的内存访问。
006：GPU内存类型、延迟与局部性 🧠
概述
在本节课中，我们将深入学习GPU的内存层次结构，理解带宽与延迟的概念，并探讨如何利用局部性原理来优化CUDA程序的性能。
GPU内存层次结构回顾
上一节我们深入探讨了线程发散问题。本节中，我们来看看GPU的内存组织。
GPU的内存层次结构与计算层次结构非常相似，了解这种层次结构有助于我们优化应用程序性能，例如利用共享内存来降低延迟。
以下是GPU的典型内存层次结构：
- 寄存器：每个线程私有，速度最快，容量有限。
- 共享内存：每个线程块共享，可编程，速度快。
- L1缓存/共享内存：每个流多处理器（SM）配置，部分可编程（作为共享内存）。
- L2缓存：所有SM共享，用于实现快速原子操作等。
- 全局内存：所有线程和网格共享，主机与设备通信的主要通道，容量大但延迟高。
- 常量内存：只读，所有线程共享，容量小（64KB），带宽高。
- 纹理内存：只读，所有线程共享，容量小（约2KB），针对2D空间局部性优化。
核心概念公式：
- 全局内存延迟：
400 - 800个时钟周期 - 共享内存延迟：
20 - 30个时钟周期 - 全局内存带宽：
~200 GB/s - 共享内存带宽：
~1 TB/s
带宽与延迟
带宽
带宽衡量的是单位时间内能够传输的数据量，而不是数据传输的速度。在GPU中，高带宽意味着可以同时进行大量的数据读写操作。
类比：飞机运输少量人员速度快（低延迟），火车运输大量人员速度慢但运量大（高带宽）。GPU类似于火车，虽然单次数据访问可能较慢，但可以并行处理大量数据访问，从而获得高吞吐量。
以下是提高带宽的几种通用技术：
- 数据重用：在设备上复用已存在的数据，减少主机与设备间的数据传输。
- 数据压缩：在传输前压缩数据，减少传输量。条件是：
压缩时间 + 传输压缩数据时间 < 传输原始数据时间。 - 重计算：在某些情况下，与其将结果从设备传回主机，不如在主机上重新计算，以节省回传的带宽。
延迟
延迟是指完成一次操作（如内存读写）所需的时间。理想情况下延迟应为零，但内存访问延迟不可避免。
在CPU上，主要通过以下方式隐藏内存访问延迟：
- 缓存层次结构（L1, L2, L3, L4缓存）。
- 多线程/多进程：当一个线程等待IO时，调度其他线程执行。
- 写回缓存：先写入高速缓存，CPU无需等待数据写入慢速内存。
- 指令级并行：执行不依赖于正在等待IO的指令的下一条指令。
- 流水线：利用硬件不同单元并行工作。
在GPU上，由于线程数量极多，每个线程可用的缓存容量很小，因此通过缓存隐藏延迟的效果有限。GPU主要依靠大规模多线程来隐藏延迟。
延迟隐藏原理：当一个线程块因内存访问（IO）而暂停时，流多处理器（SM）会立即切换到另一个就绪的线程块执行计算。通过让大量线程块交替进行计算和IO，可以保持SM始终繁忙，从而隐藏IO延迟，甚至获得超线性加速。
超线性加速示例：
- 顺序执行（无延迟隐藏）：4个线程块，每个需5个时间单位，共需
4 * 5 = 20个单位时间。 - 利用多线程隐藏延迟执行：4个线程块交错执行，仅需约8个单位时间即可全部完成。
- 加速比：
20 / 8 ≈ 2.5，超过了线程块数量（4）的线性比例（1），故为超线性加速。
局部性原理
局部性原理对CPU和GPU的性能都至关重要，缓存就是为利用局部性而设计的。GPU线程块内的线程会利用L1和L2缓存。
GPU的L1缓存（64KB）是可配置的，程序员可以指定一部分作为可编程的共享内存（Scratchpad Memory），另一部分作为硬件管理的L1缓存。程序员应积极利用共享内存来存储被频繁访问的数据片段，以利用局部性。
在GPU上，空间局部性尤为重要，它与我们后续将学到的内存合并访问技术紧密相关。
局部性的类型
在CPU架构中，主要利用两种局部性：
-
时间局部性：如果某个内存位置被访问，那么它很可能在不久的将来再次被访问。
- 代码示例：
a[i] = 10; // 第一次访问（写） a[i] = a[i] + 5; // 再次访问（读+写） b = a[i]; // 第三次访问（读）
- 代码示例：
-
空间局部性：如果某个内存位置被访问，那么它附近的内存位置很可能也会被访问。
- 代码示例：
for (int i = 0; i < N; i++) { sum += a[i]; // 访问 a[0], a[1], a[2]... 这些地址在内存中是连续的 }
- 代码示例：
缓存的设计就是为了利用这两种局部性。当访问某个字（word）时，其所在的整个内存块会被载入缓存。如果后续访问同一位置（时间局部性）或其附近位置（空间局部性），就可以快速从缓存中获取数据。
案例分析：矩阵乘法的循环顺序优化
以下是两个完成矩阵乘法 C = A * B 的代码片段，假设矩阵以行主序方式存储。
代码片段 1:
for (int i = 0; i < N; i++) {
for (int j = 0; j < P; j++) {
for (int k = 0; k < M; k++) {
C[i][j] += A[i][k] * B[k][j]; // B 被以列主序访问
}
}
}
代码片段 2:
for (int i = 0; i < N; i++) {
for (int k = 0; k < M; k++) {
for (int j = 0; j < P; j++) {
C[i][j] += A[i][k] * B[k][j]; // A, B, C 均被以行主序访问
}
}
}
问题：哪个代码片段性能更好？
分析与答案：
- 在行主序存储下，按行顺序访问数组能充分利用空间局部性，因为相邻元素被一起载入缓存。
- 代码片段1：
C[i][j]和A[i][k]是行主序访问，但B[k][j]是列主序访问（内层循环k变化最快）。这导致对B的访问无法利用缓存，性能较差。 - 代码片段2：
C[i][j]、A[i][k]和B[k][j]全部是行主序访问（对于B，外层k慢，内层j快，相当于访问B的行）。这能充分利用空间局部性，性能更好。
实测结果通常显示，代码片段2的执行时间显著短于代码片段1（例如，4.7秒 vs 9.5秒），因为它所有矩阵的访问模式都与内存布局一致，缓存命中率更高。
总结
本节课我们一起学习了GPU详细的内存层次结构及其特点，理解了带宽与延迟的区别以及GPU如何利用大规模多线程隐藏延迟。我们还回顾了局部性原理（时间局部性与空间局部性），并通过矩阵乘法的例子，看到了访问模式如何显著影响程序性能。理解这些概念是进行后续GPU内存优化（如使用共享内存、实现合并访问）的基础。
007：内存合并与AoS对SoA 🚀
概述
在本节课中，我们将要学习GPU内存访问中一个至关重要的性能优化概念——内存合并。我们将探讨什么是内存合并，以及如何通过选择合适的数据结构布局（数组结构体 与 结构体数组）来最大化内存合并的效率，从而显著提升GPU程序的性能。
上一节回顾
在上一节中，我们介绍了内存层次结构、带宽以及空间局部性和时间局部性的概念。我们通过CPU上的例子看到，遵循内存的存储顺序（如行主序）进行访问可以带来显著的性能提升。
本节中，我们将把焦点转移到GPU上，看看除了局部性之外，另一个对GPU内存性能影响巨大的因素——内存合并。
什么是内存合并？ 🤔
内存合并是GPU硬件层面的一项优化技术。其核心定义如下：
当一个线程束（32个线程）访问一个连续的内存块（例如，一个32字大小的块）时，这些零散的内存访问请求会被硬件“合并”成一个单一的、更宽的内存事务。
这意味着，原本可能需要32次独立内存访问的操作，现在可能只需要1次或几次合并后的访问即可完成。这极大地减少了内存延迟，提升了内存带宽的利用率。
合并访问示例
考虑一个经典的初始化数组的内核：
__global__ void initArray(int *A) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
A[tid] = 0; // 合并的内存访问
}
在这个例子中，线程0访问A[0]，线程1访问A[1]，...，线程31访问A[31]。由于这些元素在内存中是连续存储的，并且位于同一个内存块内，因此这32次访问会被硬件合并成一次高效的内存事务。
非合并访问的代价
如果没有内存合并，每个线程的加载/存储指令都需要独立的内存周期。在最坏情况下，一个线程束可能需要32个内存周期，这将导致严重的性能瓶颈，使得GPU的计算能力无法充分发挥。
内存合并的内部实现 🛠️
内存合并是由GPU硬件自动完成的，程序员无法直接控制其开关。其工作方式与线程束的执行模型紧密相关：
- 线程束以锁步方式执行指令。
- 当线程束执行一条内存指令（加载或存储）时，内存控制器会检查所有线程请求的内存地址。
- 如果这些地址落在一个对齐的、连续的内存块（例如32字节、64字节或128字节，具体取决于GPU架构和访问类型）内，硬件就会将这些请求合并为一个内存事务。
关键要点：内存合并是硬件的固有特性。程序员的责任是通过组织数据和访问模式来促成合并访问的发生。
内存访问模式分析 📊
根据线程束内线程访问内存的方式，我们可以将访问模式分为三类：
以下是三种典型的内存访问模式图示：
- 完全合并访问：线程束中所有线程访问连续且位于同一内存块内的地址。这是最理想的情况。
![]()
-
部分合并访问：线程束中只有部分线程的访问地址是连续的，其他线程的访问地址分散在不同的内存块中。这会导致多个内存事务。
![]()
-
完全非合并访问：线程束中每个线程访问的地址都位于完全不同的、不连续的内存块中。这会导致最差的性能，需要最多（32次）的内存事务。
![]()
代码示例：步长访问与随机访问 💻
让我们通过具体代码来理解不同的访问模式。
1. 完全合并访问
// 经典示例：连续访问
A[tid] = value;
2. 完全非合并访问（通过大跨度步长实现）
// 假设 chunk_size = 32 (或更大)
int start = tid * chunk_size;
int end = start + chunk_size;
for (int i = start; i < end; i++) {
A[i] = value; // 每个线程访问完全独立的内存块
}
当chunk_size为32时，线程0访问元素0-31，线程1访问元素32-63...，彼此毫无重叠，导致完全非合并。
3. 部分合并访问
// 假设 chunk_size = 2
int start = tid * chunk_size;
int end = start + chunk_size;
for (int i = start; i < end; i++) {
A[i] = value;
}
当chunk_size为2时，前16个线程（tid 0-15）访问的元素（如0-1, 2-3, ..., 30-31）可能位于同一个或两个内存块内，从而产生部分合并。这通常需要2次内存事务。
4. 随机访问
// 索引来自另一个输入数组，访问模式不可预测
int index = input[tid];
A[index] = value;
这种模式在图形算法（如BFS、SSSP）中很常见，其合并效率完全取决于input数组中的数据分布，通常是低效的。
数据结构布局：AoS vs SoA 🏗️
数据在内存中的组织方式（数据结构布局）对内存合并有决定性影响。主要分为两种：
数组结构体
数组结构体 是CPU编程中常见的形式。它将一个对象的所有属性打包在一起，然后创建该对象的数组。
// AoS (Array of Structures)
struct Student {
int roll_no;
double marks;
char name[50];
};
Student students[1000]; // 数组的每个元素都是一个完整的Student结构体
- CPU优势：利于空间局部性。当CPU线程访问
students[i].roll_no时，很可能紧接着访问students[i].marks，后者可能已被缓存。 - GPU劣势：当GPU线程束访问不同学生的
roll_no时，访问地址的跨度是sizeof(Student)，这会导致跨步访问，严重阻碍内存合并。
结构体数组
结构体数组 是GPU优化中推荐的形式。它将所有对象的同一属性分别存储在独立的数组中。
// SoA (Structure of Arrays)
struct StudentData {
int roll_no[1000];
double marks[1000];
char name[1000][50];
};
StudentData students; // 一个结构体，其成员是多个数组
- GPU优势：利于内存合并。当GPU线程束访问不同学生的
roll_no时（即访问roll_no[0],roll_no[1], ...），这些地址是连续的，可以实现完美的合并访问。 - CPU劣势：可能破坏空间局部性，因为访问一个学生的所有属性需要从三个不同的数组位置读取。
GPU上的AoS与SoA实现对比 ⚖️
以下是两种布局在GPU内核中访问方式的对比：
AoS 访问方式（性能较差）
// 内核函数
__global__ void processAoS(Student *d_students) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
// 跨步访问，不利于合并
d_students[tid].roll_no = tid;
d_students[tid].marks = tid * 10.0;
}
// 内存分配与拷贝
Student *d_students;
cudaMalloc(&d_students, N * sizeof(Student));
cudaMemcpy(d_students, h_students, N * sizeof(Student), cudaMemcpyHostToDevice);
SoA 访问方式（性能更优）
// 内核函数
__global__ void processSoA(int *roll_no, double *marks) {
int tid = blockIdx.x * blockDim.x + threadIdx.x;
// 连续访问，利于合并
roll_no[tid] = tid;
marks[tid] = tid * 10.0;
}
// 内存分配与拷贝（需要为每个数组成员单独分配）
int *d_roll_no;
double *d_marks;
cudaMalloc(&d_roll_no, N * sizeof(int));
cudaMalloc(&d_marks, N * sizeof(double));
cudaMemcpy(d_roll_no, h_roll_no, N * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(d_marks, h_marks, N * sizeof(double), cudaMemcpyHostToDevice);
性能对比：在实际测试中，使用SoA布局的内核执行时间可能仅为使用AoS布局内核的1/3，这巨大的性能提升主要归功于高效的内存合并。
总结
本节课中我们一起学习了GPU并行编程中两个核心的优化概念：
- 内存合并：GPU硬件将线程束内连续的、对齐的内存访问合并为单个宽内存事务的机制，是获得高内存带宽的关键。
- AoS与SoA：数据结构布局的选择直接影响内存合并的效率。
- AoS 在CPU上利用局部性有优势，但在GPU上会导致跨步访问，破坏合并。
- SoA 是GPU编程的推荐模式，它能确保对同一属性的访问是连续的，从而实现完美的内存合并，带来显著的性能提升。
作为程序员，在设计GPU程序时，应当时刻考虑数据的访问模式，并优先采用SoA布局来最大化内存合并的效益。
008：GPU中的共享内存 🧠
概述
在本节课中，我们将要学习CUDA编程中一个非常强大的特性：共享内存。我们将了解什么是共享内存、如何声明和使用共享内存变量，以及为什么它在某些场景下至关重要。我们还将通过示例来探讨使用共享内存时可能遇到的同步问题。
共享内存简介
上一节我们介绍了内存合并（Memory Coalescing）的概念，它对于提升GPU内存访问效率至关重要。本节中，我们来看看另一种可以显著提升性能的内存：共享内存。
共享内存是GPU上L1缓存的一部分，但它与CPU缓存有一个关键区别：共享内存是可由程序员显式控制的。在CPU上，数据如何进入L1、L2、L3缓存完全由硬件管理，程序员无法干预。而在GPU上，程序员可以指定哪些数据应该存放在共享内存中，这为我们优化程序提供了极大的灵活性。
共享内存主要有两个关键特性：
- 按线程块分配：每个线程块（Thread Block）都拥有自己独立的一块共享内存。即使多个线程块被调度到同一个流式多处理器（SM）上执行，它们各自的共享内存也是隔离的。
- 线程块内共享：在同一个线程块内的所有线程，可以看到并修改共享内存中的同一份数据副本。这使得共享内存成为线程块内线程间通信和协作的理想场所。
共享内存通常用于以下两种情况：
- 数据重用：当线程块内需要反复访问一小块数据时，将其加载到共享内存中可以极大降低访问延迟。
- 线程协作：当需要在线程块内的线程之间进行同步或协调时（例如，实现一个屏障或归约操作），共享内存可以作为共享的工作空间。
如何声明共享内存
在CUDA内核（Kernel）中声明共享内存变量非常简单，只需在变量声明前加上 __shared__ 限定符即可。
以下是声明共享内存数组和变量的语法：
// 在内核函数中声明一个共享内存数组
__shared__ float shared_array[1024];
// 在内核函数中声明一个共享内存变量
__shared__ unsigned int shared_variable;
重要说明：
- 这些声明必须写在CUDA内核函数内部。
- 如果声明时不加
__shared__，例如float local_array[1024];，那么这个数组将是线程局部的（Thread-local），每个线程都有自己的副本，线程间无法直接访问彼此的数据。 - 加上
__shared__后，该变量在线程块内是共享的。例如，一个包含1024个线程的线程块，所有线程都操作shared_array的同一份物理存储。
示例与分析：理解共享内存的并发访问
为了深入理解共享内存的共享特性以及可能引发的并发问题，我们来看一个具体的代码示例。
假设我们启动一个内核，只使用一个线程块，该线程块包含1024个线程。内核代码如下：
__global__ void myKernel() {
__shared__ int s;
int i = threadIdx.x;
if (i == 0) {
s = 0; // 语句A：线程0初始化s为0
}
if (i == 1) {
s += 1; // 语句B：线程1给s加1
}
if (i == 100) {
s += 2; // 语句C：线程100给s加2
}
if (i == 0) {
printf("%d\n", s); // 语句D：线程0打印s的值
}
}
我们需要思考：这个程序可能的输出结果是什么？
关键概念回顾：线程束（Warp）与锁步执行
在分析之前，必须回忆GPU的执行模型：线程以32个为一组（称为一个Warp）进行调度。同一个Warp内的线程是“锁步”（Lock-step）执行的，即它们同时执行同一条指令（尽管可能因为分支而产生部分线程不活跃的情况）。
在我们的例子中：
- 线程0和线程1属于同一个Warp（Warp 0）。
- 线程100属于另一个Warp（Warp 3）。
并发执行与可能的输出
由于不同Warp之间的执行顺序是不确定的，语句A、B、C、D可能以多种方式交织执行，导致最终 s 的值不确定。以下是几种可能的执行顺序及其结果：
- 顺序执行 (A -> B -> C -> D)：
s最终为 0 + 1 + 2 = 3。输出 3。 - 线程100先执行 (C -> A -> B -> D)：线程100读取未初始化的
s（垃圾值），加2后写入。随后A将其覆盖为0，B加1。最终s为 1。输出 1。 - A和C并发执行产生竞争：线程0执行
s=0和线程100执行s+=2同时发生。s+=2操作在底层可能被分解为“读-改-写”三步。如果线程100读取了垃圾值，计算了垃圾值+2，但在写入前，线程0写入了0，线程1写入了1，最后线程100的写入覆盖了结果。最终s可能为 垃圾值+2。 - 另一种竞争情况：类似情况3，但线程1的
s+=1在线程100最终写入后才执行，最终s可能为 垃圾值+3。 - B和C并发执行：线程1的
s+=1和线程100的s+=2并发执行，且s已被初始化为0。如果线程100的写入后发生，s最终为 2。输出 2。
结论：这个简单的程序由于存在数据竞争（Data Race），可能产生多种输出（如1, 2, 3, 垃圾值+2, 垃圾值+3），其行为是非确定性的。
这个例子清晰地展示了当多个线程不加协调地访问共享内存时会导致的问题。要获得确定性的、正确的结果，我们必须引入同步机制。
共享内存的典型使用模式
让我们回到一个更实际的例子，看看共享内存如何被应用，以及同步为何必不可少。
问题：有一个1024x1024的矩阵 M。我们想对每个元素 M[i][j] 进行如下操作：M[i][j] = M[i][j] + M[i][j+1]（对最后一列不操作）。我们计划将每一行分配给一个线程块（共1024个线程块），每个线程块内的1024个线程处理该行的1024个元素。
为了优化性能，我们决定将一行数据从全局内存加载到共享内存中进行计算。实现步骤可分为三步：
- 声明并加载数据：在线程块内声明一个共享内存数组，然后将全局内存中对应行的数据复制进来。
- 执行计算：所有线程使用共享内存中的数据执行
shared_data[j] = shared_data[j] + shared_data[j+1]。 - 写回结果：将计算后的共享内存数据写回全局内存的对应行。
这个流程中存在两个关键的同步问题：
问题一：步骤间的屏障
- 在步骤1和步骤2之间，必须确保线程块内所有线程都已完成数据加载，才能开始计算。否则，某些线程可能还在加载数据，而另一些线程已经开始使用未完全加载的数据进行计算。
- 同理，在步骤2和步骤3之间，必须确保所有线程都已完成计算，才能将结果写回全局内存。
问题二：计算步骤中的数据竞争
- 在步骤2中，线程
j需要读取shared_data[j]和shared_data[j+1]。而线程j+1同时也在读取shared_data[j+1]和shared_data[j+2]。如果操作是就地更新（in-place update），并且没有同步，就可能发生我们前面例子中类似的竞争条件，导致结果错误。
为了解决这些问题，CUDA提供了 __syncthreads() 函数。它是一个线程块级别的屏障，调用该函数会阻塞线程块内的所有线程，直到每个线程都到达这个调用点，然后所有线程才能继续执行后续指令。
一个正确的、使用了同步的伪代码结构如下：
__global__ void matrixKernel(float* M) {
__shared__ float row_data[1024];
int tid = threadIdx.x;
int row_start = blockIdx.x * 1024;
// 步骤1：协作加载数据
row_data[tid] = M[row_start + tid];
__syncthreads(); // 等待所有线程加载完毕
// 步骤2：执行计算（注意处理边界）
if (tid < 1023) { // 最后一列不计算
float temp = row_data[tid] + row_data[tid + 1];
__syncthreads(); // 确保所有读取操作完成后再写入
row_data[tid] = temp;
}
__syncthreads(); // 等待所有计算完成
// 步骤3：写回结果
M[row_start + tid] = row_data[tid];
}
注意：实际处理此类相邻元素依赖的计算时，模式可能更复杂（例如使用双缓冲），但此示例展示了同步的基本用法。
总结
本节课中我们一起学习了CUDA中共享内存的核心知识：
- 共享内存的本质：它是GPU上可由程序员控制的高速、片上内存，以线程块为单位进行分配和管理。
- 声明方法：使用
__shared__限定符在内核中声明变量或数组。 - 核心价值：用于数据重用以降低延迟，以及实现线程块内的线程间协作与通信。
- 关键挑战：并发访问共享内存会引入数据竞争，导致非确定性和错误的结果。
- 解决方案：必须使用
__syncthreads()等同步原语来协调线程块内线程的执行顺序，确保在关键操作点（如数据加载后、计算后）所有线程达成一致。
理解并正确使用共享内存与同步，是编写高效、正确CUDA并行程序的关键一步。在下一节课中，我们将更详细地探讨 __syncthreads() 以及其他同步机制。
009：动态共享内存与存储体冲突 🚀
概述
在本节课中，我们将要学习CUDA编程中的线程同步机制、动态共享内存的分配与管理，并初步了解共享内存的存储体冲突问题。上一节我们介绍了共享内存的基本概念和使用方法，本节中我们来看看如何确保线程间的正确协作，以及如何更灵活地使用共享内存。
线程同步：__syncthreads() 🔄
在上一节中，我们看到了由于线程的异步执行，可能导致数据竞争和不确定的输出。为了确保线程块内所有线程在执行到某个点之前都已完成特定操作，我们需要使用同步机制。
CUDA提供了 __syncthreads() 函数，它是一个线程块内的屏障。其工作原理如下：
- 所有属于同一线程块的线程在执行到
__syncthreads()时都会暂停。 - 线程会在此等待，直到同一线程块内的所有其他线程都到达这个同步点。
- 一旦所有线程都到达，它们才会一起继续执行之后的代码。
公式/代码描述：
// 线程块内所有线程执行至此
__syncthreads(); // 同步点：所有线程在此等待彼此
// 所有线程同步后，从此处继续执行
同步示例分析
回顾上节课末尾可能产生非确定输出的代码。通过添加 __syncthreads()，我们可以强制线程执行顺序，从而得到确定性的结果（输出始终为3）。
关键点：
__syncthreads()只对同一线程块内的线程有效。不同线程块的线程之间不会通过此函数进行同步。- 如果将
__syncthreads()放在条件语句（如if）内部，而并非块内所有线程都能满足该条件，那么未能进入该条件的线程将永远无法到达同步点，导致死锁。 - 对于同一线程束（Warp） 内的线程，由于SIMD（单指令多数据）的锁步执行特性，它们天然地按步骤执行，在某些简单情况下可能不需要显式同步。但涉及不同线程束间的协作时，
__syncthreads()是必需的。
配置L1缓存与共享内存 ⚙️
GPU的L1缓存/共享内存总大小通常为64KB。程序员可以根据内核的需求，调整用于L1缓存和共享内存的空间分配比例。
这通过 cudaDeviceSetCacheConfig 函数实现：
代码描述：
cudaDeviceSetCacheConfig(cudaFuncCachePrefer option);
参数 option 说明：
cudaFuncCachePreferNone: 默认。编译器自动根据内核代码决定分配策略。cudaFuncCachePreferShared: 偏好共享内存。分配 48KB共享内存 + 16KB L1缓存。cudaFuncCachePreferL1: 偏好L1缓存。分配 16KB共享内存 + 48KB L1缓存。cudaFuncCachePreferEqual: 均等分配。分配 32KB共享内存 + 32KB L1缓存。
使用示例：
cudaDeviceSetCacheConfig(cudaFuncCachePreferL1); // 为该内核设置L1偏好配置
myKernel<<<grid, block>>>(...); // 启动内核
注意： 如果内核声明的共享内存数组大小超过了配置的共享内存容量，编译器会报错。
动态共享内存 🧠
之前我们使用的共享内存，其大小在编译时就是确定的（例如 __shared__ int s_data[256];）。然而，有时我们需要根据运行时的输入数据来决定共享内存的大小，这时就需要使用动态共享内存。
分配动态共享内存
动态共享内存的分配分为两个步骤：
- 内核启动时指定大小：在内核调用（
<<< >>>）的第三个执行配置参数中，指定需要分配的共享内存字节数。 - 内核内声明：在内核中，使用
extern __shared__关键字声明一个未指定大小的共享内存数组。
代码描述：
// 主机端代码
int size = 1024; // 运行时决定的大小
myKernel<<<grid, block, size * sizeof(float)>>>(...);
// 设备端（内核）代码
__global__ void myKernel(...) {
extern __shared__ float s_data[]; // 动态共享内存数组
// 使用 s_data，其大小为启动内核时指定的 size * sizeof(float)
}
使用动态共享内存
动态共享内存通常作为一个“大块”内存进行分配。如果需要将其用作多个不同大小的数组，可以在内核内通过指针计算来手动划分。
示例：将一块动态共享内存分为两部分
__global__ void myKernel(int n) {
extern __shared__ int s_data[];
int *s_part1 = s_data; // 第一部分起始地址
int *s_part2 = &s_data[n/2]; // 第二部分起始地址（偏移 n/2）
// 现在可以独立使用 s_part1 和 s_part2
}
共享内存的存储体冲突 ⚡️（初步了解）
共享内存为了实现高带宽并行访问，被组织成多个存储体。当前GPU架构通常有32个存储体（与线程束大小对应）。
- 理想情况：一个线程束中的32个线程分别访问32个不同存储体中的地址。这些访问可以同时进行，实现高带宽。
- 存储体冲突：如果同一个线程束内的多个线程访问同一个存储体中的不同地址，这些访问将被迫串行化，从而显著降低访问效率。
类比：
想象一个有32个柜台的银行（32个存储体）。如果32个人（一个线程束）每人去一个不同的柜台办理业务，所有业务可以同时处理。但如果有多个人都要去同一个柜台，他们就必须排队，导致整体速度变慢。
简单示例：
假设共享内存按4字节（例如int）进行存储体映射。
- 无冲突访问：线程0访问地址0，线程1访问地址128，线程2访问地址256...（地址间隔较大，可能落入不同存储体）。
- 有冲突访问：线程0访问地址0，线程1访问地址4，线程2访问地址8...（这些地址可能位于同一存储体，导致冲突）。
存储体冲突是优化共享内存访问性能时需要重点考虑的问题。我们将在后续课程中详细讨论其原理和避免方法。
总结 🎯
本节课中我们一起学习了：
- 线程同步：使用
__syncthreads()函数来同步线程块内的线程，确保数据操作的顺序性和正确性。 - 内存配置：如何使用
cudaDeviceSetCacheConfig来配置L1缓存和共享内存的分配比例，以适应不同内核的需求。 - 动态共享内存：如何在内核启动时动态指定共享内存的大小，并在内核中使用
extern __shared__进行声明和划分。 - 存储体冲突：初步了解了共享内存的存储体组织结构，以及多线程访问同一存储体时可能引发的性能下降问题。
掌握这些知识，你将能够编写出协作更安全、内存使用更灵活、并为后续性能优化打下基础的CUDA内核。下一节我们将更深入地探讨存储体冲突及矩阵计算等高级优化技术。
010：存储体冲突、纹理内存、常量内存及计算能力 🚀
概述
在本节课中，我们将继续深入学习GPU内存体系结构。我们将详细探讨共享内存中的存储体冲突问题，并通过示例理解其成因和影响。接着，我们将介绍两种特殊的只读内存：纹理内存和常量内存，了解它们的特性、声明方式和使用场景。最后，我们将简要介绍CUDA计算能力的概念及其与GPU硬件家族的关系。
存储体冲突详解
上一节我们介绍了共享内存的基本概念和__syncthreads()同步原语。本节中，我们来看看共享内存内部的组织结构——存储体，以及由此引发的“存储体冲突”问题。
共享内存被组织成32个存储体。对共享内存的访问由线程束（32个线程）执行。理解以下关键点至关重要：
- 对同一存储体的访问是顺序执行的。
- 对同一存储体中不同字的访问会导致顺序化，从而引发存储体冲突。
- 对同一存储体中同一字的访问不会导致顺序化，因此不会引发存储体冲突。
- 连续的字存储在相邻的存储体中。
以下是理解存储体访问模式的几个核心示例：
无存储体冲突的访问模式
当线程束中的每个线程访问连续且位于不同存储体中的字时，所有访问可以并行完成，没有冲突。
__global__ void noBankConflict(float* output) {
__shared__ float s_data[1024];
int tid = threadIdx.x;
// 每个线程访问不同存储体中的元素
s_data[tid] = tid * 1.0f; // 乘以1仅为示意，无实际影响
// ... 其他操作
}
在这个例子中，thread0访问s_data[0]（存储体0），thread1访问s_data[1]（存储体1），以此类推。所有访问并行发生。
导致存储体冲突的访问模式
当线程束中的多个线程访问同一存储体中的不同字时，这些访问会被顺序化，导致性能下降。
__global__ void bankConflict(float* output) {
__shared__ float s_data[1024];
int tid = threadIdx.x;
// 导致冲突的访问模式：多个线程访问同一存储体
s_data[tid * 32] = tid * 1.0f;
// ... 其他操作
}
在这个例子中，thread0访问s_data[0]（存储体0），thread1访问s_data[32]（存储体0），thread2访问s_data[64]（存储体0）。由于它们访问的是存储体0中的不同字，这些请求会被顺序处理，形成存储体冲突。这被称为“多路存储体冲突”（例如，本例是32路冲突）。
存储体冲突可视化
为了更直观地理解，我们可以通过以下图示来分析几种访问模式：
- 场景A：每个线程访问不同存储体中的唯一字。这是无冲突的理想情况。
- 场景B：两个线程（如thread0和thread16）访问同一存储体中的不同字。这会导致2路存储体冲突。
- 场景C：访问模式看似杂乱，但每个线程访问的存储体仍然是唯一的。这同样是无冲突的。
- 场景D & E：多个线程访问同一存储体中的同一个字。根据规则3，这属于特殊情况，访问可以广播给所有请求线程，因此不会引发存储体冲突。
优化优先级：如果一个程序同时存在未合并的全局内存访问和存储体冲突，应优先优化全局内存访问，使其合并。因为全局内存访问延迟（400-800周期）远高于共享内存。在全局内存访问优化好后，再着手解决存储体冲突问题。
纹理内存
纹理内存是一种全局的、只读的内存空间，针对具有空间局部性的访问模式（如图像处理中的2D邻域访问）进行了优化。
声明与绑定
使用纹理内存需要以下步骤：
- 声明纹理引用：使用
texture类型声明一个变量。
其中，texture<float, 2, cudaReadModeElementType> texRef;float是数据类型，2表示是2D纹理，cudaReadModeElementType表示以元素形式读取。
- 绑定到设备内存：将纹理引用绑定到已分配的CUDA数组。
cudaArray* cuArray; // ... 分配并初始化cuArray (例如使用cudaMallocArray, cudaMemcpyToArray) cudaBindTextureToArray(texRef, cuArray);
在核函数中访问
在核函数内部，使用tex2D()函数来读取纹理内存。
float value = tex2D(texRef, x, y);
应用示例：图像旋转
以下是一个简化的核函数，演示如何使用纹理内存实现图像旋转（或仿射变换）：
__global__ void transformKernel(float* output, int width, int height, float theta) {
unsigned int x = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int y = blockIdx.y * blockDim.y + threadIdx.y;
if (x < width && y < height) {
// 计算旋转后的坐标 (u, v)
float u = x * cosf(theta) - y * sinf(theta);
float v = x * sinf(theta) + y * cosf(theta);
// 从纹理内存中读取旋转后坐标处的像素值
// 注意：tex2D会自动处理边界和插值（如果设置）
float pixel = tex2D(texRef, u + 0.5f, v + 0.5f);
// 写入输出数组（全局内存）
output[y * width + x] = pixel;
}
}
在这个例子中，输入图像被绑定到纹理内存texRef。每个线程计算其对应输出像素(x, y)在输入图像中对应的源位置(u, v)，然后通过tex2D高效地读取该位置的值。纹理内存硬件会对(u, v)附近的像素进行缓存，非常适合这种访问模式。
重要提示：纹理内存是只读的。计算结果需要写入到全局内存（如output数组）中。
常量内存
常量内存是另一种全局只读内存，大小固定（通常为64KB）。它用于存储需要被所有线程频繁访问的常量数据。常量内存通过缓存实现高速访问。
声明与初始化
- 声明常量变量：使用
__constant__限定符在全局作用域声明。__constant__ int constData[256];
- 从主机初始化：常量内存只能从主机（CPU）代码进行初始化，使用
cudaMemcpyToSymbol函数。int h_data = 10; cudaMemcpyToSymbol(constData, &h_data, sizeof(int));
在核函数中访问
在核函数中，可以像读取普通全局变量一样读取常量内存。
__global__ void myKernel(int* result) {
int tid = threadIdx.x + blockIdx.x * blockDim.x;
result[tid] = constData[0] * tid; // 所有线程读取同一个常量值
}
简单示例
以下代码展示了常量内存的完整使用流程：
// 1. 声明常量内存变量
__constant__ unsigned int myConst;
__global__ void constantKernel(unsigned int* data) {
// 3. 在核函数中访问常量内存
int tid = threadIdx.x;
data[tid] = myConst; // 所有线程都将自己的输出设置为常量值
}
int main() {
unsigned int h_val = 10;
unsigned int* d_data;
cudaMalloc(&d_data, 32 * sizeof(unsigned int));
// 2. 从主机复制数据到常量内存
cudaMemcpyToSymbol(myConst, &h_val, sizeof(unsigned int));
constantKernel<<<1, 32>>>(d_data);
// ... 同步并打印d_data，所有32个元素值都为10
cudaFree(d_data);
return 0;
}
CUDA计算能力
计算能力（Compute Capability）是一个由主版本号和次版本号组成的版本号（例如8.6），它代表了GPU硬件的架构和功能集。
- 硬件特性标识：不同的计算能力版本支持不同的硬件特性（如双精度原子操作、动态并行、张量核心等）。
- 编译指定：在使用
nvcc编译时，可以通过-arch=sm_xx标志指定目标计算能力（例如-arch=sm_86），以确保编译器使用该版本支持的指令和优化。 - 与CUDA版本区别：CUDA版本（如CUDA 11.0）是软件工具包的版本。计算能力是GPU硬件的属性。高版本的CUDA工具包支持更多计算能力的设备。
主要计算能力版本与特性
| 计算能力 | 架构家族 | 新增关键特性示例 |
|---|---|---|
| 1.x | Tesla | 基础CUDA支持 |
| 2.x | Fermi | 线程束内投票、共享内存原子操作 |
| 3.x | Kepler | 统一内存、动态并行性（从3.5开始） |
| 5.x, 6.x | Maxwell, Pascal | 改进的功耗效率、统一内存增强 |
| 7.x | Volta, Turing | 张量核心（AI加速）、独立的线程调度 |
| 8.x | Ampere | 硬件支持的异步拷贝、第三代张量核心 |
| 9.x | Hopper (最新) | 第四代张量核心、Transformer引擎 |
注意：计算能力4.x被跳过。计算能力10.x（Lovelace）和11.x（Blackwell）对应更新的架构，部分仍在开发或发布初期。
选择编译目标时，需要根据程序使用的特性和目标部署平台的GPU来确定合适的计算能力版本。
总结
本节课中我们一起学习了CUDA内存体系的几个高级主题。
我们深入分析了共享内存的存储体冲突，明白了其产生条件是线程束内多个线程访问同一存储体中的不同字，并通过代码示例学会了如何识别和避免这类冲突。
接着，我们探讨了两种特殊的只读内存：纹理内存和常量内存。纹理内存针对空间局部性访问优化，非常适合图像处理等算法；常量内存则通过缓存为所有线程提供对常量数据的高速访问。我们掌握了它们的声明、初始化和在核函数中访问的方法。
最后，我们介绍了CUDA计算能力的概念，理解了它是标识GPU硬件功能和架构的版本号，并知道了如何在编译时指定目标计算能力。
至此，我们已经完成了对CUDA内存子系统主要组成部分的学习。从全局内存、共享内存、寄存器到纹理和常量内存，这些知识将帮助我们更好地优化CUDA程序，充分发挥GPU的并行计算潜力。
011：数据竞争与同步 🚦
概述
在本节课中，我们将要学习并行编程中的核心概念——数据竞争与同步。我们将从数据竞争的定义和必要条件入手，通过多个代码示例深入理解其产生的原因和潜在影响。接着，我们将探讨如何通过不同的同步机制来避免数据竞争，确保并行程序的正确性。课程内容将涵盖互斥、原子操作等基本概念，并通过实际案例（如成绩统计、最短路径计算）来巩固理解。
什么是同步？
在并行编程的视角下，当多个线程协作完成特定任务时，同步指的是这些线程之间的协调。当计算任务要求线程之间必须相互协调，以确保任务能够正确完成且不产生数据竞争时，就需要同步。
上一节我们介绍了GPU内存的各个方面，本节中我们来看看线程间的协调与同步。
何时不需要同步？
并非所有并行操作都需要同步。一个典型的例子是初始化数组。例如，将数组所有元素设置为零。每个线程独立地将其负责的数组元素赋值为零，线程之间无需通信或协调即可完成任务。这是一个“令人尴尬的并行”操作。
然而，在今天的课程中，我们将看到多个需要线程间协调的场景。
同步的成本
在深入同步机制之前，需要理解一个关键点：同步操作的成本非常高。其开销远高于内存访问和计算操作。
以下是各类操作的大致开销排序（从低到高）：
- 计算：最快，通常在一个GPU周期内完成。
- 内存访问：较慢，访问寄存器、共享内存、全局内存的延迟依次增加（从几个周期到数百个周期）。
- 同步：最慢，开销最大。
因此，在设计CUDA程序时，优化顺序应为：
- 首先，尝试设计无需同步的算法。
- 其次，优化内存访问模式（如合并访问、避免存储体冲突）。
- 最后，再考虑优化计算本身（如避免线程发散）。
课程路线图
在接下来的几讲中，我们将系统学习同步相关的所有内容。以下是我们的学习路线：
以下是未来五讲的核心主题：
- 数据竞争：理解其定义和产生条件。
- 互斥：确保临界区一次只被一个线程执行。
- 实现同步的三种主要方式：
- 原子操作
- 锁
- 屏障（例如，我们已在内存章节见过的
__syncthreads()，它用于同步线程块内的线程）
- 两种重要的并行计算模式：
- 归约
- 前缀和
本节课，我们将重点学习数据竞争，并初步了解两个线程间的互斥实现。
课堂练习：识别数据竞争 (1)
让我们通过一个具体的编程问题来理解数据竞争。
任务描述：
我们有一个结构体 Point，包含属性 x 和 y。声明一个包含 n 个 Point 的数组。每个线程处理四个连续的元素（例如，线程0处理元素0-3，线程1处理元素4-7，依此类推）。每个线程需要：
- 计算其负责的四个
x值的平均值。 - 检查这四个点的
y值：- 如果任何一个
y值大于步骤1计算出的x平均值，则将该四个点的所有y值替换为该平均值。 - 否则，将这四个
y值累加到一个全局变量global_sum中。
- 如果任何一个
- 程序最后需要打印出被设置为平均值的元素数量。
问题：请尝试构思代码实现，并指出代码中可能存在的潜在问题（特别是数据竞争）。
分析与解答：
- 独立计算部分：每个线程计算自己四个元素的
x平均值、检查y值、以及可能的替换操作，这些都是独立的，不会引发数据竞争。 - 数据竞争风险点1 - 全局累加：当多个线程同时执行
global_sum = global_sum + y_value时，会发生数据竞争。这条语句在底层会被分解为“读取global_sum-> 计算新值 -> 写回global_sum”三个步骤。如果两个线程几乎同时读取了旧的global_sum，分别加上自己的y_value后写回，那么后写入的结果会覆盖前一个，导致部分加法丢失。- 示例：设
global_sum初值为4。线程A想加2，线程B想加3。正确结果应为9。但如果它们同时读取到4，线程A计算4+2=6并写回，线程B计算4+3=7并写回（覆盖了6），最终结果可能是6或7，而不是9。
- 示例：设
- 数据竞争风险点2 - 全局计数：同样，当多个线程需要递增一个全局计数器
count（用于统计有多少个元素被设置为平均值）时，count++操作也会引发类似的数据竞争。
这个练习清晰地展示了当多个线程并发读写共享变量时，数据竞争如何导致程序结果错误。
数据竞争的必要条件
从上面的例子，我们可以总结出数据竞争发生的四个必要条件，必须全部满足才会发生：
- 存在多个线程。
- 这些线程访问（读或写）同一个共享内存位置。
- 至少有一个线程执行的是写操作。
- 这些访问是并发进行的（没有强制的时间顺序）。
如果缺少其中任何一个条件，数据竞争就不会发生。
避免数据竞争的方法
既然数据竞争需要四个条件同时成立，那么避免数据竞争的方法就是破坏其中至少一个条件。以下是四种对应策略：
以下是避免数据竞争的四种主要思路：
- 顺序执行：只使用一个线程，破坏“多个线程”的条件。这牺牲了并行性。
- 数据私有化：为每个线程创建共享变量的私有副本，让线程操作自己的数据，破坏“共享内存位置”的条件。最后可能需要合并这些私有结果。
- 使用屏障分离读写：通过屏障同步点，确保所有读操作完成后，再进行写操作（或反之），破坏“并发访问”的条件。例如，先让所有线程并行读取数据，然后同步(
__syncthreads)，再让一个线程执行写操作。 - 实现互斥：确保对于临界区（访问共享资源的代码段），同一时刻只有一个线程能够进入，破坏“并发访问”的条件。这通常通过锁或原子操作实现。
课堂练习：分析互斥代码 (2)
考虑以下两个线程（T1和T2）的代码片段，其中 flag 是一个初始值未知的共享变量，S1 和 S2 代表两个临界区（例如，更新共享变量的操作）。
// 线程 T1
while(flag); // 循环等待，直到 flag 为 0
S1; // 临界区
flag = 1;
// 线程 T2
while(!flag); // 循环等待，直到 flag 为 1
S2; // 临界区
flag = 0;
问题：这段代码是否存在数据竞争？它保证了什么属性？
分析与解答：
- 是否存在数据竞争？ 从高层逻辑分析，这段代码的设计目的是实现互斥。通过分析
flag初始值为0或1的各种执行交错情况，可以发现S1和S2永远不会同时执行。虽然对flag的读写是并发的，但简单的内存赋值（如flag=1）在硬件层面通常是原子的单条指令，其并发执行不会导致数据竞争引发的不确定性错误。因此，这段代码本身没有会导致程序错误的数据竞争。 - 它保证了什么？
- 互斥：它确实保证了
S1和S2这两个临界区不会同时被执行。 - 执行顺序：进一步分析可以发现，该代码还隐含地保证了
S2总是先于S1执行。
- 互斥：它确实保证了
- 潜在问题：这段代码存在一个严重缺陷——可能发生死锁或线程饥饿。例如，若
flag初始值为1，T2会直接执行S2并将flag设为0。如果T2结束后T1才开始运行，它会将flag设为1，然后陷入while(flag);的无限循环，因为再也没有线程会将flag改回0。
这个练习说明，即使避免了数据竞争，不正确的同步设计也可能导致死锁等其他并发问题。
课堂练习：成绩统计与数据竞争 (3)
任务描述：有80个学生，学号即数组索引。给定一个存储学生分数的数组 marks[]。评分标准为：S(>=90), A(80-89), B(70-79), C(60-69), D(50-59), E(40-49), U(<40)。需要统计每个等级的学生人数。
问题：如何用CUDA内核实现？会存在数据竞争吗？如何避免？
分析与解答：
- 实现思路1（按学生并行）：启动80个线程（或更多），每个线程处理一个学生的分数。线程根据分数判断等级，然后递增对应等级的全局计数器（例如，
count_S++,count_A++）。- 数据竞争：存在。如果多个学生属于同一等级，多个线程就会同时尝试递增同一个计数器，导致数据竞争。
- 避免方法：需要使用原子操作（如
atomicAdd）来保护对每个等级计数器的递增操作。
- 实现思路2（按等级并行）：启动7个线程（对应7个等级）。每个线程遍历整个80人的分数数组，统计属于自己等级的人数。
- 数据竞争：不存在。每个线程独立操作自己的私有计数器，最后写入结果。线程间没有共享的写入位置。
- 缺点：每个线程都要遍历整个数组，计算量是思路1的7倍，当学生数量极大时效率低下。
- 综合策略：可以结合两种思路。例如，使用多个线程块，每个块内先按思路2进行局部归约（避免块内数据竞争），然后再使用原子操作或更高级的归约方法在块间进行全局累加。
案例分析：并行最短路径算法
任务描述：使用GPU计算从源城市（如图中的A）到所有其他城市的最短路径。假设使用类似Bellman-Ford的算法，在每轮迭代中，每个顶点尝试通过其邻居来松弛到其他顶点的距离。
核心代码逻辑：
__global__ void shortestPathKernel(Graph graph, int* distance) {
int tid = blockIdx.x * blockDim.x + threadIdx.x; // 线程ID对应顶点ID
if (tid >= graph.numVertices) return;
int my_dist = distance[tid];
for (each neighbor ‘n’ of vertex ‘tid’) {
int weight = getEdgeWeight(tid, n); // 获取边权
int alt_dist = my_dist + weight;
if (alt_dist < distance[n]) { // 如果找到更短路径
distance[n] = alt_dist; // 更新邻居的距离
}
}
}
问题：找出上述内核代码中至少两处可能发生数据竞争的地方，并分析其影响。
分析与解答：
- 竞争点1（影响正确性）：
if (alt_dist < distance[n])中的读取distance[n]与distance[n] = alt_dist;中的写入distance[n]可能并发发生。- 场景：顶点A和顶点D同时试图更新顶点C的距离。C的初始距离为无穷大。A计算的
alt_dist为3，D计算的为7。两者都读取到distance[C] = 无穷大，并都判断alt_dist < 无穷大成立。如果它们随后并发地执行写操作，最终distance[C]的值可能是3或7，取决于谁最后写入。而正确的最短路径是3。这导致了错误的结果。
- 场景：顶点A和顶点D同时试图更新顶点C的距离。C的初始距离为无穷大。A计算的
- 竞争点2（不影响最终正确性，但影响效率）：
int my_dist = distance[tid];中的读取与另一个线程正在执行的distance[tid] = ...写入可能并发发生。- 场景：顶点A正在更新
distance[B]为7。同时，顶点B正在读取自己的distance[B]以计算其邻居的距离。B可能读到旧的无穷大值，从而基于错误信息进行松弛。然而，在Bellman-Ford的多轮迭代中，下一轮迭代B会读到更新后的正确值7。因此，这不会影响算法的最终正确性，但可能会增加达到收敛所需的迭代轮数，影响性能。
- 场景：顶点A正在更新
这个案例说明，数据竞争不仅可能导致结果错误，也可能导致性能下降。在并行图算法中，同步至关重要。
总结
本节课中我们一起学习了并行编程中数据竞争与同步的基础知识。
- 数据竞争 发生在多个线程并发读写共享内存且至少有一个写操作时，它会导致程序行为不确定和结果错误。
- 我们通过多个实例（累加、计数、成绩统计、最短路径）深入分析了数据竞争的产生场景和后果。
- 数据竞争有四个必要条件，破坏任一条件即可避免竞争。
- 避免数据竞争的主要方法包括顺序化、数据私有化、使用屏障和实现互斥。
- 我们初步探讨了通过标志变量实现简单互斥的方法，并指出了其可能存在的死锁缺陷。
- 最后，我们以并行最短路径算法为例，分析了实际应用中复杂的数据竞争问题及其对正确性和性能的影响。
理解数据竞争是编写正确、高效CUDA程序的基础。在接下来的课程中，我们将学习更健壮、更高效的同步机制，如原子操作和锁。
012：互斥锁 🔒
概述
在本节课中，我们将学习如何使用控制流和数据流（即不使用原子操作或硬件屏障等特殊指令）来实现两个线程间的互斥锁。我们将分析几种不同的实现方法，探讨它们各自在互斥性、进展性和死锁方面的优缺点，并最终介绍经典的Peterson算法。
上节回顾
在上一节中，我们介绍了同步、互斥锁和临界区的概念。我们通过数据库的例子了解了数据竞争，分析了导致数据竞争的四个必要条件，并探讨了如何通过消除这些条件来避免数据竞争。课程最后，我们分析了一个用于计算城市间最短距离的GPU内核代码，并指出了其中因多个线程同时读写共享数组而可能引发的数据竞争问题。
同步的三种方式
实现同步通常有三种主要方式：
- 控制流 + 数据流：仅使用程序逻辑（如循环、条件判断）和普通变量进行同步。
- 原子操作：使用硬件支持的原子指令。
- 屏障：使用硬件支持的同步屏障。
其中，原子操作和屏障需要硬件支持，而控制流+数据流的方式则完全在软件层面实现。
任务：实现两线程互斥锁
我们的目标是：仅使用控制流和数据流（不使用原子操作、锁或屏障），为两个线程实现互斥锁机制。
以下是几种尝试方案及其分析。
方案一：简单标志法
我们首先尝试一个简单的方案，使用一个共享的布尔标志 flag。
// 线程 T1 的代码
while(!flag); // 等待
S1; // 临界区
flag = false;
// 线程 T2 的代码
while(flag); // 等待
S2; // 临界区
flag = true;
假设：flag 初始值为 false。
分析：
- 互斥性：❌ 不保证。如果
flag初始值为true，两个线程的while循环条件会同时为假，导致它们同时进入临界区。 - 进展性：❌ 不保证。执行顺序被固定：必须先执行
S2，然后才能执行S1。如果T2不执行，T1将永远等待。 - 死锁：✅ 不会发生。因为
flag的值非真即假，不会导致两个线程都卡在while循环中。
结论：此方案存在严重缺陷，无法保证互斥性。
方案二：对称等待法
我们尝试修改方案一，让两个线程对称地等待对方。
// 线程 T1 的代码
while(!flag);
S1;
flag = false;
// 线程 T2 的代码
while(flag);
S2;
flag = true;
假设：flag 初始值为 false。flag 被声明为 volatile，确保读写直接作用于主存，避免缓存一致性问题。
分析：
- 互斥性：✅ 保证。
flag在同一时刻只能有一个值（true或false），因此只有一个线程能通过while循环。 - 进展性：❌ 不保证。如果
flag初始为true且T1永远不想进入临界区，那么想进入临界区的T2将永远在while(flag);处等待。 - 死锁：❌ 可能发生。考虑以下执行顺序：
T1执行flag = false;（假设之前已执行完S1）。T2执行flag = true;。T1和T2几乎同时到达各自的while循环。
此时，flag的最终值为true（后写入者生效）。T1看到while(!true)为假，进入S1。T2看到while(true)为真，陷入无限等待。这虽然不是一个典型的“双方互相等待”的死锁，但导致了T2饥饿，从效果上看类似于死锁。更精确地说，在并发执行时，可能导致一个线程永远等待。
结论：此方案保证了互斥性，但无法保证进展性，且在某些执行顺序下会导致线程饥饿/死锁。
方案三：数组声明法 (Lock Version 1)
我们尝试实现通用的 lock() 和 unlock() 函数，使用一个标志数组 flag[2]。
volatile bool flag[2] = {false, false};
void lock(int tid) {
int me = tid; // 0 或 1
int other = 1 - me;
flag[me] = true; // 声明自己希望进入临界区
while(flag[other] == true); // 如果对方也希望进入，则等待
}
void unlock(int tid) {
flag[tid] = false; // 声明自己离开临界区
}
分析：
- 互斥性：✅ 保证。只有当对方线程的标志为
false时，当前线程才能进入临界区。 - 进展性：✅ 保证。即使一个线程不参与（其
flag为false），另一个线程也能顺利进入临界区。 - 死锁：❌ 可能发生。考虑以下执行顺序：
T0执行flag[0] = true;T1执行flag[1] = true;- 现在
flag[0]和flag[1]都为true。 T0执行while(flag[1] == true);陷入等待。T1执行while(flag[0] == true);陷入等待。
双方都在等待对方将标志置为false，形成死锁。
结论：此方案保证了互斥性和进展性，但存在死锁风险。
方案四：牺牲变量法 (Lock Version 2)
我们引入一个共享的 victim 变量来打破对称性，避免死锁。
volatile int victim;
void lock(int tid) {
int me = tid;
victim = me; // 主动“牺牲”自己，让对方优先
while(victim == me); // 如果自己仍然是“牺牲者”，则等待
}
void unlock(int tid) {
// 解锁操作简单，无需修改 victim
// victim 的值会在下一次 lock 竞争时被覆盖
}
分析：
- 互斥性：✅ 保证。
victim只能为 0 或 1。最后执行victim = me;的线程会成为真正的“牺牲者”，在其while循环中等待。 - 进展性：❌ 不保证。如果只有一个线程（如
T0）调用lock()，它会执行victim = 0;，然后发现victim == 0为真，从而在while循环中无限等待，即使没有竞争者。 - 死锁：✅ 不会发生。
victim的最终值会确定一个唯一的等待线程，另一个线程总能进入临界区。
结论：此方案保证了互斥性且无死锁，但无法保证进展性（单个线程无法独立工作）。
最终方案：Peterson算法 (Lock Version 3)
Peterson算法巧妙地结合了声明意向 (flag[]) 和礼貌谦让 (victim) 两种机制。
volatile bool flag[2] = {false, false};
volatile int victim;
void lock(int tid) {
int me = tid;
int other = 1 - me;
flag[me] = true; // 步骤1：声明自己希望进入
victim = me; // 步骤2：主动让对方优先
// 步骤3：等待条件：当对方想进入，并且自己是被礼让者时，继续等待
while(flag[other] == true && victim == me);
}
void unlock(int tid) {
flag[tid] = false; // 离开临界区，取消声明
}
分析：
- 互斥性：✅ 保证。进入临界区的条件是：
flag[other] == false（对方不想进）或victim == other（对方是牺牲者）。两个线程不可能同时使对方的条件为假。 - 进展性：✅ 保证。如果只有一个线程想进入（
flag[other] == false），它将立即通过while循环。即使两个线程都想进入，victim变量也能确保其中之一（最后设置victim的线程）会等待，另一个则进入。 - 死锁：✅ 不会发生。
victim变量确保了至少有一个线程能通过while循环。不会出现双方无限等待的情况。
核心思想：
flag[me] = true：大声宣布“我想进去”。victim = me：礼貌地说“您先请”。while(...)：等待的条件是“您想进，并且您还让我先请”，那我就继续等。否则（您不想进，或者该我进了），我就进去。
总结
本节课我们一起学习了如何仅用控制流和数据流实现两线程互斥锁。我们经历了多个方案的迭代：
- 简单标志法：互斥性和进展性均无法保证。
- 对称等待法：保证互斥性，但无进展性，可能饥饿。
- 数组声明法 (V1)：保证互斥性和进展性，但存在死锁风险。
- 牺牲变量法 (V2)：保证互斥性且无死锁，但无进展性。
- Peterson算法 (V3)：完美地同时保证了互斥性、进展性和无死锁。
Peterson算法是一个经典的软件互斥解决方案，它优雅地展示了如何通过共享变量和程序逻辑来解决复杂的同步问题。需要注意的是，它只适用于两个线程的场景。在下一节课中，我们将探讨如何将互斥锁的概念扩展到多个线程。
013：互斥锁（续）及CUDA中的原子操作 🔒⚛️
概述
在本节课中，我们将要学习如何为多个线程（n个线程）实现互斥锁，并深入了解CUDA中的原子操作。我们将从适用于多线程的“面包店算法”开始，然后探讨原子操作的概念、其与互斥锁的区别，以及几种常见的原子指令。
面包店算法 🍞
上一节我们介绍了如何为两个线程实现互斥锁。本节中我们来看看如何为n个线程实现互斥锁，这需要使用“面包店算法”。
面包店算法由Leslie Lamport提出，它模拟了面包店或服务台取号排队的机制。每个线程在进入临界区前会获取一个“号码”（标签），然后根据号码大小顺序进入临界区，从而保证互斥和先进先出（FIFO）的顺序。
以下是该算法的核心伪代码：
bool flag[n]; // 初始化为 false
int label[n]; // 初始化为 0
void lock(int me) {
flag[me] = true; // 表示线程 me 希望进入临界区
label[me] = 1 + max(label[0], ..., label[n-1]); // 获取当前最大标签值加1
while ( (存在 k != me 且 flag[k] == true) &&
( (label[k], k) < (label[me], me) ) ) {
// 等待：如果存在其他希望进入且“号码”更小的线程，则循环等待
}
}
void unlock(int me) {
flag[me] = false; // 表示线程 me 已离开临界区
}
核心概念解释：
flag[me]: 表示线程me是否希望进入临界区。label[me]: 线程me的“排队号码”。- 字典序比较
(label[k], k) < (label[me], me): 先比较标签label，如果标签相同，则比较线程IDk和me。这确保了当多个线程获得相同标签时，线程ID较小的线程优先。
算法要点与讨论
- 标签生成非原子性：
label[me] = 1 + max(...)这一步不是原子操作。如果多个线程同时计算最大值并获取标签，它们可能得到相同的标签值。 - 字典序比较的必要性：如果只比较标签 (
label[k] < label[me])，当多个线程获得相同标签时，它们可能同时进入临界区，破坏互斥性。因此需要引入线程ID作为决胜局。 - 内存可见性：
flag和label数组应声明为volatile，或在支持缓存一致性的系统中确保写入能及时对其他线程可见，否则可能因缓存导致线程读取到旧值。 - 在GPU上的局限性：面包店算法在GPU的同一个Warp内执行会导致死锁。因为Warp内的32个线程是锁步（SIMD）执行的。如果所有线程都执行
lock函数，它们会同时设置flag、计算相同的label，然后大部分线程会在while循环中等待。而唯一满足条件（如线程ID最小）的线程也因为Warp锁步执行而无法前进，从而导致整个Warp死锁。
从互斥锁到原子操作 🔄
上一节我们探讨了基于软件的互斥锁实现。本节中我们来看看GPU上更常用、更底层的同步原语——原子操作。
为什么需要原子操作？
考虑一个简单的例子：两个线程（来自不同的线程块）同时对一个全局内存变量x进行自增操作 x++。
// 假设 x 指向全局内存，初始值为0
// 内核启动配置: <<<2, 1>>> 表示2个线程块，每个块1个线程
__global__ void kernel(int *x) {
x[0]++; // 这不是原子操作！
}
x++ 在底层可能被编译成三条机器指令：
- LOAD: 将
x[0]的值从内存加载到寄存器R。 - ADD: 将寄存器
R的值加1。 - STORE: 将寄存器
R的值存回x[0]所在的内存。
由于两个线程的执行可能交错进行，最终x的值可能是1或2，而不是预期的2。这就是数据竞争。
原子操作可以解决这个问题。它将“读取-修改-写入”这一系列操作打包成一个不可分割的、连续的整体操作。在执行原子操作期间，其他线程无法访问目标内存位置，从而保证了结果的正确性。
原子操作 vs. 互斥锁
- 互斥锁：是一种更高级别的软件抽象，用于保护一段代码（临界区）。线程在进入前获取锁，退出时释放锁。在GPU上，由于线程数量极多，锁的开销很大，容易导致性能下降和死锁。
- 原子操作：是硬件支持的底层指令，直接作用于单个内存位置（如加、减、比较交换等）。它更轻量，适用于简单的数据更新场景。在GPU编程中，原子操作比互斥锁更常用。
原子操作的特性与应用场景
- 硬件支持：原子操作需要硬件支持，CUDA GPU提供了丰富的原子指令。
- 内存范围：原子操作可用于全局内存和共享内存。
- 常见指令：
atomicAdd(address, val): 原子加法。atomicCAS(address, compare, val): 原子比较并交换（可实现锁）。atomicMin/atomicMax(address, val): 原子求最小/最大值。
- 应用场景：统计直方图、求和、求极值、实现无锁数据结构等。
现实世界类比：银行ATM取款 💳
假设一个联合账户有1000元存款，双胞胎兄弟同时从两台ATM机各取1000元。ATM机的逻辑是：
- 检查余额
>=取款金额。 - 如果满足，则扣款并吐出现金。
如果步骤1和2不是原子操作，可能发生以下交错执行：
- 两台ATM都检查余额（1000 >= 1000，通过）。
- 两台ATM都扣款（1000 - 1000 = 0）。
- 两台ATM都吐出现金1000元。
结果：账户余额为0，但用户总共取走了2000元，银行损失1000元。
原子操作可以确保“检查余额”和“扣款”作为一个不可分割的单元执行，从而避免此问题。
CUDA原子加法示例 ➕
为了解决之前x++的数据竞争问题，我们可以使用CUDA的atomicAdd函数。
__global__ void safe_kernel(int *x) {
atomicAdd(&x[0], 1); // 原子地将 x[0] 的值增加1
}
函数原型：int atomicAdd(int* address, int val);
address: 要进行原子操作的目标内存地址。val: 要加上的值。- 返回值：返回该地址在操作前的旧值。
使用atomicAdd后，无论启动多少个线程块和线程，也无论执行顺序如何，最终x[0]的值都将是所有线程执行次数的总和（例如，启动配置为<<<K1, K2>>>，则最终值为 K1 * K2）。
性能注意：原子操作会导致对同一内存地址的访问串行化，即多个线程必须排队依次执行该操作。这会成为性能瓶颈，应尽量避免在热点循环中频繁使用原子操作。
总结
本节课中我们一起学习了：
- 面包店算法：一种为n个线程实现互斥锁的算法，通过取号和字典序比较来保证互斥与顺序。但它在GPU的Warp内执行会导致死锁。
- 原子操作：硬件支持的、不可分割的单一内存操作。它解决了多线程下的数据竞争问题，是GPU上重要的同步机制。
- 原子操作与互斥锁的区别：原子操作更底层、更轻量，适用于简单的数据更新；互斥锁更高级，用于保护代码段，但在GPU上开销大。
atomicAdd指令：CUDA中用于原子加法的函数，它能确保计数的正确性，但会引入串行化开销。
下一节课，我们将深入探讨更强大的atomicCAS（比较并交换）指令，并学习如何使用它来实现自旋锁。
014：CUDA中的多种原子函数与屏障 🚀
在本节课中，我们将学习如何使用CUDA中的原子函数来解决数据竞争问题，并深入探讨屏障（Barrier）这一重要的同步机制。我们将从回顾单源最短路径算法中的数据竞争问题开始，然后学习如何使用原子操作（如 atomicMin 和 atomicCAS）来确保代码的正确性。接着，我们将了解如何利用 atomicCAS 实现锁（Lock）和单线程执行区（Single Section）。最后，我们将正式定义屏障，并介绍CUDA中不同级别的屏障同步。
回顾：单源最短路径算法中的数据竞争 🔄
在之前的课程中，我们学习了单源最短路径算法。我们注意到，在并行执行时，代码中存在数据竞争问题，这可能导致计算出错误的距离值。
问题出现在以下代码段中：
if (alt_dist < dist[neighbor]) {
dist[neighbor] = alt_dist;
}
当多个线程同时读取和更新同一个邻居节点（例如节点C）的距离时，可能会发生数据竞争。例如，线程A和线程B可能都读取到节点C的初始距离为无穷大，然后分别尝试将其更新为3和6。最后写入的值（可能是6）将覆盖掉更短的正确距离（3），从而导致错误。
上一节我们介绍了数据竞争的问题，本节中我们来看看如何使用原子操作来解决它。
使用原子操作解决数据竞争 ⚛️
为了消除数据竞争，我们需要确保“读取-比较-写入”这一系列操作是原子性的，即不可分割的。CUDA提供了 atomicMin 函数，它正是用于此目的。
我们可以将上述易出错的代码替换为：
atomicMin(&dist[neighbor], alt_dist);
atomicMin 函数会原子性地比较 dist[neighbor] 的当前值与 alt_dist 的值，并将 dist[neighbor] 更新为两者中的较小值。这样，无论线程执行顺序如何，每个节点的距离都会被正确地更新为所有可能路径中的最小值。
深入原子比较与交换（atomicCAS） 🔧
atomicCAS（Compare And Swap）是一个功能强大的原子操作，其语法如下：
int atomicCAS(int* address, int compare, int val);
它的工作流程是：
- 原子性地读取
address指针处的值。 - 将该值与
compare参数进行比较。 - 如果相等，则将
address处的值设置为val。 - 无论是否交换，函数都返回
address处原来的旧值。
atomicCAS 有两个典型应用：
- 实现锁（Lock），用于保护临界区。
- 实现单线程执行区（Single Section），确保只有任意一个线程执行特定代码块。
使用 atomicCAS 实现锁 🔒
简单的 atomicCAS 调用本身并不能实现锁。我们需要将其与循环结合，构建一个自旋锁。
以下是使用 atomicCAS 在GPU上实现一个能正确处理线程束（Warp）内同步的锁的示例：
// 假设 lock_var 初始化为 0
__shared__ int lock_var;
// 尝试获取锁
int old;
do {
old = atomicCAS(&lock_var, 0, 1); // 尝试将锁从0设为1
} while (old != 0); // 如果old不是0，说明锁已被占用，继续循环
// 临界区代码
// ...
// 释放锁
lock_var = 0;
重要说明：这个基础实现在CPU或不同Warp的线程间可以工作，但在同一个Warp的线程间会导致死锁。因为Warp内的线程是同步执行的，一个线程获得锁后，其他线程会在循环中等待，而获得锁的线程也在等待循环中的其他线程，形成死锁。
为了解决Warp内死锁问题，需要结合条件判断和线程掩码，确保同一时刻只有一个Warp线程进入临界区，而其他线程被正确屏蔽。这通常涉及更复杂的逻辑。
使用 atomicCAS 实现单线程执行区 🎯
单线程执行区要求只有任意一个线程执行某段代码，其他线程无需等待，直接跳过。
使用 atomicCAS 可以轻松实现：
// 假设 single_flag 初始化为 0
__shared__ int single_flag;
int old = atomicCAS(&single_flag, 0, 1); // 尝试将标志从0设为1
if (old == 0) {
// 只有一个线程会成功进入此区域
// 单线程执行区代码
// ...
}
// 其他线程直接跳过，继续执行后续代码
在这个实现中，single_flag 在首次被成功后就不再是0，因此后续线程的 atomicCAS 操作都会失败，从而不会进入if块。线程执行完单线程区后，也无需重置 single_flag。
理解屏障（Barrier） 🚧
屏障是一种同步点，要求所有参与线程都必须到达该点后，才能有任何线程继续执行后续代码。
类比：想象一群朋友约定去踢足球。大家从各自家中（异步）出发，到一个集合点汇合（屏障）。只有所有人都到达集合点后，大家才一起前往足球场。
在CUDA中，存在不同级别的屏障：
- 隐式全局屏障：内核（Kernel）的结束点本身就是一个隐式屏障，主机代码会等待所有GPU线程完成。
- 线程块内屏障：使用
__syncthreads()函数。它同步同一个线程块内的所有线程。 - 全局屏障：从CUDA 9开始，支持
__syncgrid()（或类似机制，具体API需查阅文档），可以同步网格（Grid）中所有线程块内的线程。 - 线程束内屏障：对于Warp内的线程，由于它们以锁步方式执行指令，因此本身就存在隐式的执行屏障，通常不需要显式同步。
总结 📚
本节课我们一起学习了以下核心内容：
- 原子函数：我们回顾了使用
atomicMin解决单源最短路径算法中的数据竞争问题。 atomicCAS：深入学习了其原理，并探讨了如何用它来实现锁（解决临界区互斥访问）和单线程执行区。特别要注意在GPU上实现锁时需考虑Warp内线程的死锁问题。- 屏障：我们正式定义了屏障作为同步点的概念，并介绍了CUDA中不同层次的同步机制，包括隐式内核屏障、
__syncthreads()线程块内屏障以及更新的全局屏障支持。
通过掌握原子操作和屏障，你可以编写出正确、高效且同步良好的CUDA并行程序，避免数据竞争和协调线程间的执行顺序。
015：屏障、归约与CUDA中的前缀和 🚀
概述
在本节课中，我们将要学习CUDA编程中的三个核心概念：屏障（Barrier）、归约（Reduction）和前缀和（Prefix Sum）。我们将了解屏障如何实现线程间的同步，如何利用归约高效地计算数组的总和，以及如何实现和应用前缀和算法。
屏障（Barrier）回顾与深入
上一节我们介绍了原子操作，本节我们继续探讨屏障这一同步机制。
屏障是一个程序点，要求所有线程都必须到达该点后，任何线程才能继续执行后续指令。在CUDA中，内核（Kernel）的结束处对所有GPU线程来说是一个隐式的全局屏障。从CUDA 9开始，grid.sync()函数可以作为跨GPU所有线程的全局屏障。而__syncthreads()函数则充当线程块（Thread Block）内所有线程的屏障。
屏障示例分析
以下是使用__syncthreads()的一个示例代码：
__global__ void kernel(unsigned int* vector, unsigned int vector_size) {
unsigned int id = threadIdx.x;
vector[id] = id; // 每个线程根据其ID初始化数组元素
__syncthreads(); // 线程块内屏障
if (id < vector_size - 1 && vector[id + 1] != id + 1) {
printf("__syncthreads doesn't work\n");
}
}
以下是代码逻辑的逐步分析：
- 每个线程根据其
threadIdx.x（即ID）为数组vector的对应位置赋值。 - 调用
__syncthreads()，确保线程块内所有线程都完成了初始化操作。 - 每个线程（除了最后一个）检查其下一个元素的值是否等于
id + 1。
潜在问题：如果数组大小vector_size超过了一个线程块的最大线程数（例如1024），则需要启动多个线程块。__syncthreads()仅同步同一个线程块内的线程。因此，可能出现以下情况：
- 第一个线程块（线程ID 0-1023）的线程完成了初始化，并通过了屏障。
- 第二个线程块（线程ID 1024-2047）的线程仍在执行初始化。
- 此时，第一个线程块中ID为1023的线程执行检查：它试图访问
vector[1024]（即id+1），但这个位置的值尚未被第二个线程块的线程初始化。这可能导致检查条件成立，从而打印错误信息。
解决方案：如果使用grid.sync()替代__syncthreads()，则可以确保所有GPU线程（跨所有线程块）都完成初始化后才进行检查，从而避免此问题。
屏障的数据同步作用
__syncthreads()不仅提供控制流同步（所有线程到达后才继续），还提供数据同步。它确保在屏障点之前，线程块内所有线程对内存的写入操作，对屏障点之后所有线程的读取操作是可见的。
这是通过内存栅栏（Memory Fence） 操作实现的。内存栅栏强制将线程的写入刷新到主内存，并确保后续读取从主内存获取最新值，从而保证内存操作的顺序性和可见性。
__threadfence_block(): 在__syncthreads()内部隐式调用，确保线程块内的数据同步。__threadfence(): 确保GPU内所有线程的数据同步。
归约（Reduction）
上一节我们介绍了屏障，本节中我们来看看归约。归约是一种将一组数据（例如一个数组）聚合为单个值或少量值的操作，例如计算总和、最大值、最小值等。
为什么需要归约？
考虑计算一个包含N个元素的数组的总和。一种简单的方法是使用原子操作：
__global__ void atomic_sum(int* array, int* sum, int n) {
int id = threadIdx.x;
if (id < n) {
atomicAdd(sum, array[id]); // 原子加法
}
}
这种方法的问题是严重的顺序性：尽管有多个线程，但atomicAdd操作是串行的，同一时刻只有一个线程能成功更新sum。这无法充分利用GPU的并行能力，性能甚至可能差于CPU串行代码。
归约的基本思想
归约通过树形结构并行地合并数据。其核心思想是：在每一步，将数据两两分组进行合并，使数据量减半，经过 log₂N 步后得到最终结果。
操作前提：归约操作（如加法、求最大值）必须是可结合（Associative） 的。即 (a ⊕ b) ⊕ c = a ⊕ (b ⊕ c)。加法、乘法、求最大值、最小值、按位与/或/异或等操作满足结合律。而减法、除法则不满足。
归约算法实现
以下是一种常见的归约实现方式（假设数组大小为2的幂，且只使用一个线程块）：
__global__ void reduction_sum(int* input, int* output, int n) {
// 假设input数组已在GPU上，且n为2的幂
// 使用共享内存或全局内存进行原地操作
for (int offset = n / 2; offset > 0; offset >>= 1) {
int id = threadIdx.x;
if (id < offset) {
input[id] += input[id + offset]; // 将相隔offset的元素相加
}
__syncthreads(); // 每一步都需要同步
}
if (threadIdx.x == 0) {
*output = input[0]; // 最终结果在input[0]
}
}
算法步骤：
- 初始时，有N个元素，使用 N/2 个线程。
- 每个线程将位置
id的元素与位置id + offset的元素相加，结果存回id位置。offset初始为N/2。 - 调用
__syncthreads()确保本步所有加法完成。 offset减半（offset >>= 1），重复步骤2-3，直到offset为0。- 此时，
input[0]中存储的就是所有元素的总和。
图解：对于数组 [2, 3, 5, 6, 1, 9, 4, 7]
- 第1步 (
offset=4):[2+1, 3+9, 5+4, 6+7, 1, 9, 4, 7]->[3, 12, 9, 13, ...] - 第2步 (
offset=2):[3+9, 12+13, 9, 13, ...]->[12, 25, ...] - 第3步 (
offset=1):[12+25, ...]->[37, ...]
结果37即为总和。
优势：
- 并行度高：每一步都使用大量线程。
- 步数少：仅需 log₂N 步，而非串行的 N 步。
- 无需原子操作：通过巧妙的索引计算和同步避免数据竞争。
注意事项：
- 需要处理数组大小非2的幂的情况（例如，在数组末尾填充0）。
- 当使用多个线程块时，需要先在各块内进行归约，然后再对块的结果进行归约，并可能使用
grid.sync()或启动多个内核。
前缀和（Prefix Sum）
在学习了归约之后，我们来看一个与之相关但功能不同的算法——前缀和，也称为扫描（Scan）。
什么是前缀和？
对于一个输入数组 [x₀, x₁, x₂, ..., xₙ₋₁]，其前缀和输出数组 [y₀, y₁, y₂, ..., yₙ₋₁] 定义为：
y₀ = x₀y₁ = x₀ + x₁y₂ = x₀ + x₁ + x₂- ...
yₙ₋₁ = x₀ + x₁ + ... + xₙ₋₁
即，输出数组的每个元素是输入数组从起始到当前位置的累积和。
前缀和的应用
一个典型应用是工作项分配。假设有多个线程，每个线程需要向一个全局工作列表（数组）中推送不同数量的数据项。我们想知道每个线程推送的数据在全局数组中的起始和结束索引。
例如：
- 线程0推送4个项。
- 线程1推送3个项。
- 线程2推送9个项。
- 线程3推送5个项。
我们可以构建一个计数数组 counts = [4, 3, 9, 5]。计算其前缀和（从0开始）得到 prefix = [0, 4, 7, 16]。
- 线程0的项位于索引
prefix[0]到prefix[1]-1，即[0, 3]。 - 线程1的项位于索引
prefix[1]到prefix[2]-1，即[4, 6]。 - 以此类推。前缀和提供了每个线程数据在全局数组中的“偏移量”。
前缀和算法实现（Blelloch Scan）
以下是一种高效的前缀和并行算法实现（上行扫描）：
__global__ void prefix_sum(int* data, int n) {
// 假设 data 是全局数组，n 为2的幂，且只使用一个线程块
// 上行（Reduce）阶段：构建求和树
for (int offset = 1; offset < n; offset <<= 1) {
int id = threadIdx.x;
int temp = 0;
if (id >= offset) {
temp = data[id - offset]; // 读取前一个偏移量的值
}
__syncthreads(); // 分离读和写，避免数据竞争
if (id >= offset) {
data[id] += temp; // 累加到当前位置
}
__syncthreads(); // 同步本次迭代
}
// 此时 data 数组的最后一个元素 data[n-1] 是总和
// 下行（Downsweep）阶段（此处省略）：将部分和分散，得到最终前缀和
// ... (需要额外的步骤将 exclusive scan 转换为 inclusive scan，或反之)
}
算法逻辑（以上行阶段为例）：
offset从1开始，每次迭代翻倍（1, 2, 4, ...）。- 对于每个
offset，每个线程（除了前offset个）读取当前位置向前offset步的值。 - 使用
__syncthreads()确保所有读取操作在写入开始前完成，这是避免数据竞争的关键。 - 将读取的值累加到当前元素。
- 再次同步，确保本步所有计算完成，再进行下一步。
图解（简化上行阶段）：
初始: [1, 2, 3, 4, 5, 6, 7, 8]
offset=1: [1, 1+2, 2+3, 3+4, 4+5, 5+6, 6+7, 7+8] -> [1, 3, 5, 7, 9, 11, 13, 15]
offset=2: [1, 3, 1+5, 3+7, 5+9, 7+11, 9+13, 11+15] -> [1, 3, 6, 10, 14, 18, 22, 26]
offset=4: [1, 3, 6, 10, 1+14, 3+18, 6+22, 10+26] -> [1, 3, 6, 10, 15, 21, 28, 36]
最终data[7]=36是总和。要得到完整的前缀和数组[1,3,6,10,15,21,28,36]，还需要一个下行阶段。
关键点：
- 数据竞争：在读取
data[id - offset]和写入data[id]时，必须用屏障分离，防止其他线程同时修改正在读取的数据。 - 屏障位置：代码中使用了两个
__syncthreads()。第一个确保“读”在“写”之前完成；第二个确保本步所有计算完成，属于迭代间的控制同步。 - 复杂度：同样为 O(log N) 步。
其他应用
前缀和不仅用于上述偏移计算，还广泛应用于：
- 流压缩（Stream Compaction）：过滤数组中的有效元素。
- 基数排序（Radix Sort）的关键步骤。
- 构建数据结构，如芬威克树（Fenwick Tree）。
- 计算直方图（Histogram）的累积分布。
总结
本节课中我们一起学习了CUDA编程中三个重要的高级主题：
- 屏障：深入理解了
__syncthreads()和grid.sync()的同步范围，并通过例子认识到屏障选择不当可能导致逻辑错误。同时，明白了屏障兼具控制同步和数据同步（通过内存栅栏）的功能。 - 归约：掌握了将大规模数据并行聚合为单一值的高效方法。理解了归约的树形结构思想，学习了其并行实现代码，并认识到其相对于原子操作的巨大性能优势。
- 前缀和：学习了前缀和的概念及其在并行计算中的关键应用（如工作分配）。分析了一种并行前缀和算法的实现，特别注意了其中使用屏障来避免数据竞争的技巧。
这些算法是许多并行计算应用的基础构件，理解它们对于编写高效、正确的CUDA程序至关重要。
016：CUDA中的函数 🚀
概述
在本节课中，我们将要学习CUDA中不同类型的函数。到目前为止，我们只见过一种使用 __global__ 关键字定义的函数，我们称之为内核。但CUDA还提供了其他几种定义函数的方式。我们将详细探讨这些不同类型的函数，并通过多个示例来理解它们的设计理念和具体用法。如果时间允许，我们还将简要介绍CUDA的并行算法库——Thrust。
函数类型介绍
上一节我们介绍了CUDA编程的基本概念，本节中我们来看看CUDA中定义函数的几种不同方式。
在CUDA中，主要有三种方式声明函数，每种方式决定了函数在哪里被调用以及在哪里执行。
__global__函数（内核）- 使用
__global__关键字声明。 - 从主机（CPU）调用。
- 在设备（GPU）上执行。
- 必须返回
void类型。
- 使用
__device__函数（设备函数）- 使用
__device__关键字声明。 - 从设备（GPU）调用。
- 在设备（GPU）上执行。
- 通常由内核函数调用。
- 使用
__host__函数（主机函数）- 使用
__host__关键字声明。 - 从主机（CPU）调用。
- 在主机（CPU）上执行。
- 这与普通的C++函数行为一致。
__host__关键字的主要用途是与__device__结合，定义可在主机和设备上执行的函数。
- 使用
一个程序可以包含多种类型的函数，并且同一种函数可以被多次调用。
函数调用规则与示例
理解了基本类型后，我们通过一个具体的代码示例来深入理解这些函数之间的调用关系。
以下是示例代码的第一部分，展示了四种不同方式声明的函数：
// 一个既可在主机运行也可在设备运行的函数
__host__ __device__ void dh_func() {
printf("I can run on both CPU and GPU\n");
}
// 一个设备函数
__device__ void d_func() {
dh_func(); // 设备函数可以调用主机-设备函数
}
// 一个内核函数
__global__ void d_kernel() {
dh_func(); // 内核可以调用主机-设备函数
d_func(); // 内核可以调用设备函数
}
// 一个主机函数
__host__ void host_func() {
dh_func(); // 主机函数可以调用主机-设备函数
}
以下是代码的第二部分，即主函数：
int main() {
// 从主机（CPU）调用内核
d_kernel<<<1, 1>>>();
cudaDeviceSynchronize();
// 从主机调用主机函数
host_func();
// 从主机调用主机-设备函数
dh_func();
// 以下调用会导致编译错误，因为设备函数只能从设备调用
// d_func();
return 0;
}
基于以上代码，我们可以总结出函数调用的核心规则：
__device__函数：只能被__global__内核或其他__device__函数调用。__global__内核：只能从主机（CPU）调用。__host__ __device__函数：可以被主机或设备代码调用，但其函数体内的代码必须是主机和设备都能执行的“交集”。
为了更清晰地展示，以下是可能的函数调用关系图：
main->d_kernel(有效)main->host_func(有效)main->dh_func(有效)host_func->dh_func(有效)d_kernel->dh_func(有效)d_kernel->d_func(有效)d_func->dh_func(有效)
重要限制：__host__ __device__ 函数内部不能调用纯 __device__ 函数。因为 __host__ __device__ 函数可能从CPU被调用，而 __device__ 函数只能在GPU上执行，这会导致运行时错误。同理，其内部也不能使用 threadIdx.x 或 __syncthreads() 等GPU特有的变量或函数。
统一内存与函数
上一节我们探讨了函数调用的基本规则，本节中我们来看看函数如何与一种特殊的内存——统一内存进行交互。
统一内存（使用 cudaMallocManaged 分配）提供了一种在CPU和GPU之间自动同步的数据视图。这对于 __host__ __device__ 函数特别有用，因为它们可能从两端访问同一数据。
请看以下示例：
__host__ __device__ void func(int *counter) {
(*counter)++;
printf("Counter: %d\n", *counter);
}
__global__ void print_k(int *counter) {
func(counter); // 在GPU上调用func
}
int main() {
int *counter;
// 分配统一内存
cudaMallocManaged(&counter, sizeof(int));
*counter = 0; // 在CPU上初始化
printf("Main: %d\n", *counter); // 输出: 0
print_k<<<1, 1>>>(counter); // 启动内核
cudaDeviceSynchronize(); // 等待GPU完成
func(counter); // 在CPU上再次调用func
printf("Main: %d\n", *counter); // 输出: 2
cudaFree(counter);
return 0;
}
程序输出：
Main: 0
Counter: 1
Counter: 2
Main: 2
解释：
- 计数器在CPU上初始化为0。
- 内核
print_k在GPU上执行func，将计数器增加到1并打印。 - 由于是统一内存，GPU的修改对CPU可见。
- CPU再次调用
func，此时看到的值是1，将其增加到2并打印。 - 最终在
main中打印的值为2。
这个例子展示了统一内存如何简化CPU和GPU之间的数据共享。如果使用普通的 cudaMalloc（仅在GPU分配）或 malloc（仅在CPU分配），并尝试从另一端访问，则会导致段错误或读取垃圾值。
实践：数组增量示例
理论需要结合实践，本节我们通过一个具体的编程任务来巩固所学知识：编写一个既能被CPU串行执行，也能被GPU并行执行的数组元素自增函数。
目标：实现一个函数，对数组的所有元素执行 arr[i]++ 操作。
初始方案（CPU串行，GPU单线程串行）：
__host__ __device__ void func(int *arr, int n) {
for (int i = 0; i < n; ++i) {
arr[i]++;
}
}
int main() {
int n = 1024;
int *h_arr = new int[n];
int *d_arr;
cudaMalloc(&d_arr, n * sizeof(int));
// 初始化数组...
// CPU端串行执行
func(h_arr, n);
// GPU端串行执行（仅使用1个线程）
func<<<1, 1>>>(d_arr, n);
// ...
}
此方案在GPU上也是串行的，效率低下。
优化方案（GPU并行）：
为了利用GPU的并行能力，我们修改调用方式，让每个GPU线程处理一个数组元素。
__host__ __device__ void func(int *elem) {
(*elem)++; // 只操作单个元素
}
__global__ void d_parallel(int *arr, int n) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < n) {
func(&arr[idx]); // 每个线程调用func处理一个元素
}
}
int main() {
// ... 分配和初始化内存
// CPU端仍串行
for (int i = 0; i < n; ++i) func(&h_arr[i]);
// GPU端并行
int threadsPerBlock = 256;
int blocks = (n + threadsPerBlock - 1) / threadsPerBlock;
d_parallel<<<blocks, threadsPerBlock>>>(d_arr, n);
// ...
}
现在GPU执行是并行的，但CPU端仍需显式循环。
最终方案（统一的函数接口）：
我们可以让 func 函数自身根据调用者决定行为，提供更统一的接口。
__host__ __device__ void func(int *arr, int n) {
// 当从GPU并行调用时，n=1，每个线程增加一个元素
// 当从CPU串行调用时，n=数组长度，循环增加所有元素
for (int i = 0; i < n; ++i) {
arr[i]++;
}
}
__global__ void d_parallel_unified(int *arr, int n) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < n) {
// 每个线程告诉func只增加它负责的那个元素
func(&arr[idx], 1);
}
}
int main() {
// ... 分配和初始化内存
// CPU端调用
func(h_arr, n); // n是数组长度，内部循环
// GPU端调用
d_parallel_unified<<<blocks, threadsPerBlock>>>(d_arr, n);
// ...
}
这个设计让 func 函数更加灵活，能适应不同的执行环境。
扩展知识：Thrust库简介 ⚡
在课程的最后，我们简要介绍一下CUDA的Thrust库。Thrust是一个类似于C++标准模板库的并行算法库，它提供了丰富的数据结构（如向量）和算法（如排序、归约、前缀和）。
Thrust的特点：
- 高层抽象：程序员无需关心内核启动、线程管理等底层细节。
- 可移植性：相同的Thrust代码可以在CPU（使用STL后端）或GPU（使用CUDA后端）上运行，编译器根据设备选择实现。
- 高性能：库内部经过优化，能有效利用GPU的并行计算能力。
简单示例（归约求和）：
使用Thrust进行数组求和比手动编写归约内核简单得多：
#include <thrust/reduce.h>
#include <thrust/device_vector.h>
// ...
thrust::device_vector<int> d_vec(n);
// ... 填充数据
int sum = thrust::reduce(d_vec.begin(), d_vec.end(), 0, thrust::plus<int>());
Thrust库包含大量此类实用算法，能极大提升CUDA程序的开发效率。
总结
本节课中我们一起学习了CUDA中不同类型的函数及其调用规则。我们深入探讨了 __global__、__device__ 和 __host__ 函数的定义、调用限制以及它们与统一内存的协作。通过数组增量等示例，我们实践了如何编写既能用于CPU串行执行也能用于GPU并行执行的灵活函数。最后，我们简要了解了能简化并行编程的Thrust库。掌握这些函数类型是构建高效、复杂CUDA应用程序的基础。
017：CUDA中的支持功能
在本节课中，我们将要学习CUDA编程中的一些支持功能，包括如何使用Makefile自动化编译和运行流程，如何使用CUDA-GDB调试程序，以及如何使用性能分析工具（如nvprof）来评估和优化CUDA程序的性能。
Makefile基础
上一节我们介绍了CUDA中的函数类型。本节中，我们来看看如何利用Makefile来管理CUDA项目的编译和运行。
Makefile是一组包含目标和变量的命令集合，用于自动化编译、链接和清理项目文件。它特别适用于管理大型项目中文件之间的依赖关系。
以下是Makefile的基本结构：
target: prerequisites
command
- target：可以是生成的文件名，也可以是要执行的操作名称。
- prerequisites：是执行该目标所依赖的文件或其他目标。
- command：是生成目标或执行操作所需的一系列shell命令。注意：命令前必须使用Tab字符缩进，而不是空格。
一个简单的Makefile示例
以下是一个打印“Hello”的简单Makefile规则：
hello:
echo Hello
当在终端运行 make hello 命令时，它会执行 echo Hello，输出“Hello”。
用于CUDA程序的Makefile示例
让我们为一个打印“Hello”32次的简单CUDA程序编写Makefile。我们希望它能自动完成编译、运行和清理。
hello:
nvcc hello.cu -o hello
run: hello
./hello
clean: run
rm hello
这个Makefile包含三个规则：
hello：编译hello.cu文件，生成可执行文件hello。run：依赖于hello目标。在hello编译完成后，运行./hello。clean：依赖于run目标。在程序运行后，删除生成的可执行文件hello。
执行 make clean 将按顺序触发这三个目标。
在Makefile中使用变量
我们可以使用变量来使Makefile更清晰、更易于维护。
FILE1 = hello
FILE2 = run
FILE3 = clean
$(FILE1):
nvcc hello.cu -o hello
$(FILE2): $(FILE1)
./hello
$(FILE3): $(FILE2)
rm hello
变量通过 = 赋值，并通过 $(变量名) 引用。
使用 all 目标管理多个任务
all 目标常用于指定默认要构建的一系列目标。
all: hello1 hello2 hello3
hello1:
nvcc hello1.cu -o hello1
./hello1
hello2:
nvcc hello2.cu -o hello2
./hello2
hello3:
nvcc hello3.cu -o hello3
./hello3
clean:
rm hello1 hello2 hello3
运行 make all 将按顺序编译并运行 hello1、hello2 和 hello3，最后可以运行 make clean 进行清理。
Makefile中的条件与循环
Makefile支持条件判断和循环，以增加灵活性。
条件判断示例：
TARGET_CPU = x86_64
ifeq ($(TARGET_CPU), x86_64)
ARCH_FLAG = -m64
else
ARCH_FLAG = -m32
endif
循环示例：
run_tests:
for number in 1 2 3 4 5; do \
./my_program input$$number.txt output$$number.txt; \
done
综合示例：自动化测试多个图数据
假设我们有一个CUDA程序 vertex_removal.cu，需要对 test_bench.zip 中的多个图文件（t1.mtx 到 t10.mtx）进行测试。
all: vertex_removal
unzip -o examples/test_cases/test_bench.zip -d .
number=1; \
while [ $$number -le 10 ]; do \
if [ $$number -ne 5 ]; then \
./vertex_removal t$$number.mtx vr_graph_$$number.txt; \
fi; \
number=$$((number + 1)); \
done
vertex_removal:
nvcc vertex_removal.cu -o vertex_removal
clean:
rm -f vertex_removal *.txt
这个Makefile自动化了以下流程：
- 编译
vertex_removal.cu。 - 解压测试数据。
- 使用循环对除
t5.mtx外的所有测试图运行程序，并将结果输出到对应文件。 clean目标用于清理生成的可执行文件和输出文件。
使用CUDA-GDB调试
调试并行程序比调试串行程序更具挑战性，因为线程执行顺序的不确定性可能导致不同的结果。CUDA-GDB是GNU调试器(GDB)的扩展，专门用于调试运行在GPU硬件上的CUDA程序。
调试示例：空指针解引用
考虑以下存在错误的CUDA内核：
__global__ void myKernel(int *x) {
*x = 0; // 潜在的空指针解引用
}
int main() {
int *d_x; // 仅在主机端声明指针，未在设备上分配内存
myKernel<<<2, 10>>>(d_x);
cudaDeviceSynchronize();
return 0;
}
此代码中，设备指针 d_x 未被分配设备内存，导致内核中尝试解引用非法地址。
使用 cudaGetLastError 捕获错误
一种调试方法是使用 cudaGetLastError 函数。
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA error: %s\\n", cudaGetErrorString(err));
}
运行有问题的程序后，此代码会输出类似 CUDA error: an illegal memory access was encountered 的信息，提示存在非法内存访问。
使用CUDA-GDB进行交互式调试
更强大的方法是使用CUDA-GDB。
- 编译调试版本：使用
-G和-g标志编译CUDA代码，生成调试信息。nvcc -G -g my_program.cu -o my_program_debug - 启动调试：
cuda-gdb my_program_debug - 设置断点并运行：在
(cuda-gdb)提示符下：(cuda-gdb) break main (cuda-gdb) run - 分析错误：程序崩溃后，CUDA-GDB会显示错误位置（如
my_program.cu:6）、引发错误的线程信息（网格、块、线程ID）以及指针的值（可能为0x0），从而帮助定位空指针解引用问题。
常用的CUDA-GDB命令
info cuda kernels：列出当前所有CUDA内核的信息。info cuda threads：显示当前内核中线程的状态。cuda thread block thread：切换调试焦点到指定的线程块和线程。info cuda lanes：显示特定线程（通道）的详细信息。break filename.cu:23 if threadIdx.x == 1：在特定文件和行号设置条件断点。
性能分析：nvprof
性能分析是优化CUDA程序的关键步骤。nvprof 是NVIDIA提供的命令行性能分析工具，属于非侵入式分析，无需修改源代码。
使用nvprof进行性能分析
假设我们有一个程序调用了三个内核（kernel1、kernel2、kernel3）各100次。
-
直接分析：
nvprof ./my_cuda_programnvprof会输出每个内核的执行时间、调用次数以及占总时间的百分比。==PROF== Profiling result: Time(%) Time Calls Kernel 39.2% 1.23ms 100 kernel2 33.1% 1.04ms 100 kernel3 26.5% 0.83ms 100 kernel1 95.0% -- 300 CUDA API (cudaLaunchKernel)从结果可知，
kernel2耗时最多，是首要的优化目标。同时，CUDA内核启动开销也占据了显著时间。 -
分析特定代码区域：可以使用
nvprof的API在代码中标记分析范围。#include <nvToolsExt.h> nvtxRangePushA("MyCodeSection"); // ... 要分析的代码 ... nvtxRangePop();然后运行：
nvprof --analysis-metrics -o profile.nvvp ./my_cuda_program生成的
profile.nvvp文件可以用更可视化的工具（如NVIDIA Nsight Systems）进行深入分析。
nvprof 可以提供丰富的指标，包括内存吞吐量、缓存命中率、分支分化等，帮助开发者全面了解程序在GPU上的行为。
本节课中我们一起学习了CUDA编程中的三项重要支持功能：使用Makefile自动化项目构建流程，使用CUDA-GDB调试GPU内核代码，以及使用nvprof工具进行程序性能分析。掌握这些工具能显著提升CUDA程序的开发效率和运行性能。
018：动态并行、多GPU与线程束投票 🚀
概述
在本节课中，我们将学习CUDA编程中的三个高级主题：动态并行、多GPU编程以及线程束投票。动态并行允许从GPU内核内部启动新的内核，多GPU编程则涉及如何利用系统中的多个GPU进行协同计算，而线程束投票则提供了一种在同一个线程束内进行高效通信和决策的机制。
动态并行 🔄
上一节我们介绍了CUDA程序的构建、调试和性能分析。本节中，我们来看看动态并行。
动态并行在需要嵌套并行的情况下非常有用。例如，从CPU启动一个包含N个线程的内核，然后该内核的每个线程再启动另一个包含M个线程的内核。这种模式可以高效地模拟嵌套循环等场景，特别是当内层循环的迭代次数在运行时才能确定时。
动态并行的适用场景
以下是动态并行典型的适用场景：
- 分层数据结构：例如树结构，每个节点的子节点数量可能不同。遍历所有节点（外层循环），并对每个节点的可变数量的子节点进行处理（内层循环）。
- 递归算法：递归调用可以建模为递归树，树中每个节点的分支数量（递归调用次数）可能不同。
- 非均匀工作负载：任务可以自然地划分为独立的批次，但每个批次的工作量并不均匀。
注意：并非所有嵌套并行都需要动态并行。如果内外层循环的迭代次数在编译时已知，可以使用静态并行，即直接启动一个包含 N * M 个线程的单一内核，而不是使用动态并行。
动态并行示例
以下是一个简单的动态并行示例代码：
__global__ void child(int father) {
printf("Parent %d, Child %d\n", father, threadIdx.x);
}
__global__ void parent() {
printf("Parent %d\n", threadIdx.x);
child<<<1, 5>>>(threadIdx.x);
}
int main() {
parent<<<1, 3>>>();
cudaDeviceSynchronize();
return 0;
}
编译与运行：编译动态并行代码需要指定计算能力（至少3.5）并启用可重定位设备代码模式。
nvcc -arch=sm_35 -rdc=true dynamic_example.cu -o dynamic_example
程序输出分析：parent 内核由3个线程执行，每个线程打印自己的ID，然后启动一个包含5个线程的 child 内核。child 内核打印其父线程ID和自己的线程ID。由于内核启动是异步的，且多个 parent 线程并发执行，child 内核的输出顺序可能与 parent 线程的启动顺序交错，但属于同一个 parent 线程启动的5个 child 线程的输出顺序是固定的。
动态并行的计算与内存特性
- 计算特性：
- 父内核与子内核关联着各自的网格。
- 一个父内核可以启动多个子内核网格。
- 父内核和子内核可能异步执行。父内核线程在启动子内核后不会等待其完成，而是继续执行后续语句。
- 父网格在其所有子网格完成之前不会被视为完成。
- 内存特性：
- 父内核和子内核共享相同的全局内存和常量内存。
- 父内核和子内核拥有不同的本地内存和共享内存。不能将父内核的共享内存指针传递给子内核直接访问。
- 在子内核启动之前，父内核完成的所有全局内存操作对子内核可见。
- 子内核完成的所有全局内存操作，只有在父内核通过
cudaDeviceSynchronize()等方式同步等待子内核完成后，才对父内核可见。
多GPU编程 🖥️🖥️
上一节我们介绍了动态并行。本节中，我们来看看如何利用多个GPU进行并行计算。
使用多个GPU的主要优势在于可以获得更多的计算核心和更大的总内存容量，从而能够加速计算或处理更大规模的数据。
多GPU系统架构
典型的系统包含一个CPU和多个GPU。GPU之间可能存在两种连接方式：
- 点对点访问：某些GPU对（例如，通过PCIe桥接器或NVLink直接连接）可以直接访问彼此的内存，无需经过CPU，数据传输延迟较低。
- 通过CPU中转：对于没有直接连接的GPU对，数据交换需要通过CPU内存中转（GPU A -> CPU内存 -> GPU B），这会带来额外的开销。
多GPU编程基础
在CUDA中，系统中的GPU被编号为0到N-1。核心API是 cudaSetDevice(int device_id)，用于设置当前操作的GPU设备。在此调用之后，所有的内核启动、内存分配和数据传输操作都将在该设备上进行。
以下代码展示了基本的双GPU使用模式：
cudaSetDevice(0); // 切换到 GPU 0
kernel1<<<...>>>(...); // 在 GPU 0 上启动 kernel1
cudaMemcpyAsync(..., cudaMemcpyDeviceToHost); // 从 GPU 0 异步拷贝数据
cudaSetDevice(1); // 切换到 GPU 1
kernel2<<<...>>>(...); // 在 GPU 1 上启动 kernel2
cudaMemcpyAsync(..., cudaMemcpyDeviceToHost); // 从 GPU 1 异步拷贝数据
关键点：为了确保 kernel1 和 kernel2 能够真正在GPU 0和GPU 1上并发执行，必须使用异步内存拷贝函数 cudaMemcpyAsync，而不是阻塞式的 cudaMemcpy。否则，CPU会在第一个 cudaMemcpy 处等待，导致内核串行执行。
多GPU相关API
以下是用于多GPU编程的一些重要API：
cudaGetDeviceCount(int* count)：获取系统中可用的GPU数量。cudaDeviceCanAccessPeer(int* canAccessPeer, int device, int peerDevice)：检查device能否直接访问peerDevice的内存（点对点访问）。cudaDeviceEnablePeerAccess(int peerDevice, unsigned int flags)：启用从当前设备到peerDevice的点对点访问。
注意：点对点访问在硬件层面是对称的，但软件API需要显式地双向启用。一个设备最多只能与8个其他设备建立点对点连接。
枚举设备属性示例
以下代码展示了如何获取并打印系统中每个GPU设备的属性：
int device_count;
cudaGetDeviceCount(&device_count);
for (int i = 0; i < device_count; ++i) {
cudaDeviceProp prop;
cudaGetDeviceProperties(&prop, i);
printf("Device %d: %s\n", i, prop.name);
printf(" Compute Capability: %d.%d\n", prop.major, prop.minor);
// ... 打印其他属性
}
线程束投票 ✋
上一节我们介绍了多GPU编程。本节中，我们来看看线程束投票，这是一种在同一个线程束（Warp，通常32个线程）内进行高效协作的机制。
线程束投票函数对于需要根据某个条件在 warp 内进行集体决策的场景非常有用，例如，当需要根据数据依赖性为整个 warp 分配任务时。
线程束投票函数
CUDA 提供了以下内置的线程束投票函数：
__all(int predicate)：如果 warp 内所有活动线程的predicate值非零，则对所有线程返回1（真），否则返回0（假）。__any(int predicate)：如果 warp 内任意活动线程的predicate值非零，则对所有线程返回1，否则返回0。__ballot(int predicate)：返回一个32位的无符号整数掩码（mask），其中第n位对应于 warp 内第n个线程的predicate值（非零则为1，否则为0）。该掩码对所有线程返回相同的值。
__ballot 是 __all 和 __any 的泛化形式。
投票函数示例
考虑以下内核，它使用 __all 检查 warp 内所有线程的 threadIdx.x 是否都小于100：
__global__ void kernel() {
int tid = threadIdx.x;
int pred = (tid < 100);
int result = __all(pred);
if (tid % 32 == 0) { // 每个warp选一个代表线程打印
printf("Warp %d: __all(tid<100) = %d\n", tid/32, result);
}
}
// 启动配置：<<<1, 128>>> (4个warps)
对于前3个warps（线程0-95），所有线程都满足条件，result 为1。对于第4个warp（线程96-127），只有部分线程满足条件，result 为0。输出可能为：1, 1, 1, 0。
如果将 __all 替换为 __any，则第4个warp中也有线程满足条件，因此所有warps的 result 都为1。
__ballot 示例：
unsigned int mask = __ballot(threadIdx.x % 2 == 0); // 检查是否为偶数线程
对于每个warp，mask 是一个32位数，偶数位为1，奇数位为0。其十六进制表示为 0x55555555（二进制0101...）。
使用投票优化原子操作
线程束投票的一个典型应用是优化高争用的原子操作。考虑以下需要条件原子递增的代码：
// 原始代码：高争用
if (condition) {
atomicAdd(&counter, value);
}
如果warp中有多个线程满足 condition，它们会各自调用 atomicAdd，导致严重的资源争用。
优化版本使用 __ballot 和 __popc（计算掩码中1的位数）来减少原子操作调用次数：
// 优化代码：减少争用
unsigned int mask = __ballot(condition);
if ((threadIdx.x % 32) == 0) { // 每个warp选一个leader线程
int num_active = __popc(mask); // 计算满足条件的线程数
if (num_active > 0) {
atomicAdd(&counter, num_active * value); // Leader线程执行一次原子加
}
}
在这个优化版本中，无论warp中有多少线程满足条件，每个warp只执行一次原子操作，由选出的leader线程完成，从而显著降低了争用。
注意：投票函数的结果受线程活动状态影响。如果投票函数被放在条件语句（如 if(tid % 2 == 0)）内部，那么只有满足条件的“活动”线程会参与投票，不满足条件的线程被掩蔽（masked），不参与投票计算。
总结
本节课中我们一起学习了CUDA编程的三个高级主题。
- 动态并行：允许内核从设备端启动新的内核，适用于处理运行时才能确定工作量的嵌套并行任务，但需要注意其启动开销和内存访问规则。
- 多GPU编程：通过
cudaSetDevice管理多个GPU设备，利用异步操作实现真正的并发执行，并可以使用点对点访问优化GPU间的数据传输。 - 线程束投票：提供了
__all、__any和__ballot等函数，用于在warp内进行高效的集体决策和通信，常用于优化原子操作等高性能计算模式。
掌握这些高级特性，能够帮助你编写出更灵活、更高效、能够充分利用硬件资源的CUDA程序。
019：CUDA中的多维数组 🚀
在本节课中，我们将学习如何在CUDA中声明和使用多维数组。我们将通过三个具体的应用实例来深入理解这一概念：矩阵平方、二维最小值算法和模板计算。通过本课，你将学会如何更直观地处理多维数据，避免复杂的索引计算，并编写更易读的GPU代码。
概述 📋
在之前的课程中，我们使用一维数组来表示二维矩阵，并通过公式 i * matrix_size + j 来访问元素。本节课，我们将学习如何直接在CUDA中声明和访问二维（甚至更高维）数组，从而使代码更简洁、更易维护。
矩阵平方应用 🔲
上一节我们回顾了使用一维数组进行矩阵平方计算的方法。本节中，我们来看看如何使用真正的二维数组来实现相同的功能。
在之前的实现中，我们使用一维数组并通过公式计算索引来模拟二维访问。现在，我们将直接在GPU上声明二维数组。
CPU版本（使用一维数组）
void cpu_square(float* matrix, float* result, int matrix_size) {
for (int i = 0; i < matrix_size; i++) {
for (int j = 0; j < matrix_size; j++) {
for (int k = 0; k < matrix_size; k++) {
result[i * matrix_size + j] += matrix[i * matrix_size + k] * matrix[k * matrix_size + j];
}
}
}
}
GPU版本（使用二维数组）
以下是使用CUDA二维数组实现矩阵平方的关键步骤：
-
在GPU上声明二维数组：
__device__ float d_A[N][N]; __device__ float d_B[N][N]; -
主机端代码配置：
dim3 blocksPerGrid(1, 1); dim3 threadsPerBlock(N, N); // 使用 cudaMemcpyToSymbol 将数据从主机复制到设备 cudaMemcpyToSymbol(d_A, A, size); cudaMemcpyToSymbol(d_B, B, size); // 启动核函数 square_kernel<<<blocksPerGrid, threadsPerBlock>>>(N); // 将结果复制回主机 cudaMemcpyFromSymbol(C, d_B, size);
- 核函数实现：
__global__ void square_kernel(int n) { int i = threadIdx.y; int j = threadIdx.x; if (i < n && j < n) { float sum = 0.0f; for (int k = 0; k < n; k++) { sum += d_A[i][k] * d_A[k][j]; } d_B[i][j] = sum; } }
核心优势： 我们不再需要复杂的索引计算 i * n + j，可以直接使用 d_A[i][j] 这种直观的语法访问元素，代码可读性大大增强。
二维最小值算法 📉
上一节我们介绍了矩阵操作，本节中我们来看看一个图像处理中常见的算法——二维最小值算法。
该算法扫描输入矩阵中每个2x2的窗口，找出窗口内的最小值，并放入输出矩阵的对应位置。这类似于卷积操作。
算法定义
对于输出矩阵 O 中的元素 O[i][j]，其值为输入矩阵 A 中以下四个元素的最小值：
A[i-1][j], A[i-1][j+1], A[i][j], A[i][j+1]
边界情况（第一行和最后一列）保持不变。
CPU实现
void compute(float grid[N][N], float prev_grid[N][N]) {
for (int i = 1; i < N; i++) { // 忽略第一行
for (int j = 0; j < N-1; j++) { // 忽略最后一列
float min_val = fminf(prev_grid[i-1][j], prev_grid[i-1][j+1]);
min_val = fminf(min_val, prev_grid[i][j]);
min_val = fminf(min_val, prev_grid[i][j+1]);
grid[i][j] = min_val;
}
}
}
GPU实现
以下是GPU核函数的关键部分：
__global__ void min2d_kernel(int n) {
int i = blockIdx.y * blockDim.y + threadIdx.y;
int j = blockIdx.x * blockDim.x + threadIdx.x;
// 处理边界，避免第一行和最后一列
if (i > 0 && j < n-1) {
float min_val = fminf(d_prev_grid[i-1][j], d_prev_grid[i-1][j+1]);
min_val = fminf(min_val, d_prev_grid[i][j]);
min_val = fminf(min_val, d_prev_grid[i][j+1]);
d_grid[i][j] = min_val;
} else if (i == 0 || j == n-1) {
// 边界直接复制
d_grid[i][j] = d_prev_grid[i][j];
}
}
实现要点： 每个GPU线程负责计算一个输出元素。我们使用二维线程块来自然映射到矩阵的(i, j)坐标，并通过条件判断优雅地处理了边界情况。
模板计算 🧮
上一节我们学习了空间上的邻域操作，本节中我们来看看一个涉及多次迭代的模板计算。
模板计算中，每个单元格的新值由其上方两个邻居的当前值决定，公式为：
A[i][j] += A[i-1][j] + A[i-1][j+1]
该计算需要迭代多次（例如10次），且每次迭代后，新值成为下一次迭代的输入。
基础版本（使用额外空间）
最简单的实现是使用两个独立的数组，一个用于读取（上一轮结果），一个用于写入（本轮结果）。
GPU核函数示例：
__global__ void stencil_kernel(int n) {
int i = blockIdx.y * blockDim.y + threadIdx.y;
int j = blockIdx.x * blockDim.x + threadIdx.x;
if (i > 0 && j < n-1) {
d_grid_new[i][j] = d_grid_old[i][j] + d_grid_old[i-1][j] + d_grid_old[i-1][j+1];
}
}
// 每次迭代后，需要将 d_grid_new 的数据复制到 d_grid_old
进阶挑战：原地计算
现在，我们增加一个限制：不能使用额外的数组，必须在原输入矩阵上就地更新数据。
这就带来了数据竞争问题：当线程正在读取A[i-1][j]和A[i-1][j+1]来计算A[i][j]的新值时，另一个负责计算A[i-1][j]的线程可能正在写入它。
解决方案：按行顺序执行
为了消除数据竞争，我们可以将并行度限制在单行内。即，每次只计算一行：
- 从最后一行开始计算（因为它没有下方的依赖）。
- 等待该行所有线程计算完成后，再开始计算上一行。
- 依此类推，直到第一行。
主机端代码逻辑：
for (int iter = 0; iter < 10; iter++) {
for (int row = N-1; row >= 1; row--) { // 从底向上
// 1. 将当前行数据复制到设备
// 2. 启动核函数，只计算这一行
dim3 threadsPerBlock(N, 1); // 一行N个线程
stencil_inplace_kernel<<<1, threadsPerBlock>>>(N, row);
// 3. 将计算后的行复制回主机
}
}
核函数实现：
__global__ void stencil_inplace_kernel(int n, int current_row) {
int j = threadIdx.x; // 列索引
int i = current_row; // 行索引由主机传入
if (j < n-1) { // 忽略最后一列
d_grid[i][j] += d_grid[i-1][j] + d_grid[i-1][j+1];
}
}
性能权衡： 这种方法保证了正确性，但牺牲了并行度（从并行处理整个矩阵变为每次只并行处理一行）。另一种可能的方案是使用原子操作，但需要精心设计以避免将并行计算完全序列化。这可以作为课后思考题。
总结 🎯
本节课我们一起学习了CUDA中多维数组的核心用法：
- 声明与访问： 使用
__device__ float arr[N][M]在GPU上声明二维数组，并可直接用arr[i][j]语法访问，无需手动计算线性索引。 - 三个应用实例：
- 矩阵平方：展示了用二维数组简化代码逻辑。
- 二维最小值算法：演示了如何将图像处理类算法映射到二维线程网格，并处理边界条件。
- 模板计算：深入探讨了迭代计算中的依赖问题，并给出了通过限制执行顺序来实现原地计算的方案。
- 核心思想： 利用CUDA的二维/三维线程组织来自然匹配多维数据结构，使算法更清晰，并注意处理边界条件和迭代计算中的数据竞争问题。
通过掌握这些技巧，你可以更高效、更直观地在GPU上处理科学计算、图像处理等领域的多维数据问题。
浙公网安备 33010602011771号