【CUDA】Kernel Atomic Stream
【CUDA】Kernel & Atomic & Stream
Kernel
线程同步
在 CUDA 编程中，线程同步是极为关键的环节。
cudaDeviceSynchronize()
cudaDeviceSynchronize(); 函数的作用在于确保针对一个问题的所有内核操作都完成后，才安全地开始下一个操作。可以将其看作是一个屏障。它通常在 int main() {} 或者其他非 __global__ 函数中被调用。由于 GPU 线程是异步执行的，它们的执行顺序并不确定。例如在一些复杂的计算流程中，如果后续的计算依赖于前面内核操作的结果，那么使用这个函数就能避免在前面结果还未就绪时就开始后续操作，从而防止错误的产生。
__syncthreads()
__syncthreads(); 则是用于在内核中设置线程执行的屏障。当多个线程可能会对相同的内存位置进行操作时，就需要使用它。比如，在某些情况下，一个线程可能还在对内存中的某个区域进行处理，而另一个线程已经完成了相关任务并可能会修改该内存区域。如果没有这个同步操作，就可能出现数值不稳定和错误。例如在处理一个复杂的数据结构，多个线程同时对其不同部分进行读写操作时，就需要在合适的点使用 __syncthreads() 来保证数据的正确性。
// 模拟位移动作的核函数，第二个位移动作，依赖第一个位移动作结果
__global__ void shiftBitsSecondStep(int* data) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < DATA_LENGTH) {
// 这里使用__syncthreads进行同步，确保第一个位移动作完成
__syncthreads();
data[idx] <<= 2;
}
}
__syncwarps()
__syncwarps(); 函数专门用于同步一个线程束（warp）内的所有线程。线程束是 GPU 执行模型中的一个重要概念，通过同步线程束内的线程，可以更好地协调它们的操作，确保在某些特定的计算场景下，线程束内的操作按照预期顺序进行。
// CUDA核函数，使用__syncwarps进行线程束内同步
__global__ void warpSyncExample(int* array) {
int tid = threadIdx.x + blockIdx.x * blockDim.x;
if (tid < ARRAY_SIZE) {
// 每个线程先进行一次自增操作
array[tid]++;
// 使用__syncwarps进行线程束内同步，确保本线程束内所有线程都完成了自增操作
__syncwarps();
// 再进行一次自增操作，这里可以看到同步后的效果，两次自增是连贯的按顺序的（在每个线程束内）
array[tid]++;
}
}
以向量加法和位移动作的并行计算为例，如果要对数组 a = [1, 2, 3, 4] 和 b = [5, 6, 7, 8] 进行向量加法并将结果存储在 c 中，然后再给 c 中的每个元素加 1，按照数学运算顺序（PEDMAS），我们需要确保所有的加法操作完成后再进行加 1 的操作。如果不进行线程同步，就可能出现某个元素在还未完成加法时就被加 1，从而得到错误的输出向量。同样，在并行化位移动作时，如果一个位移动作依赖于前一个位移动作的结果，那么就必须使用同步操作来确保前一个位移动作完成后再进行下一个。
线程安全性
当一段代码被称为 “线程安全” 时，意味着它可以被多个线程同时运行而不会导致竞态条件或其他意外行为。在 CUDA 中，竞态条件是指一个线程在另一个线程完成之前就开始下一个任务。为了防止这种情况，我们使用 cudaDeviceSynchronize() 函数。就好像一群线程在赛跑，有些线程由于各种原因先到达终点，我们需要手动告诉这些 “获胜” 的线程在终点等待落后的线程，这样才能保证整个过程的正确性。如果涉及到使用不同的 CPU 线程调用多个 GPU 内核，可以参考相关的资料进一步深入了解。
SIMD/SIMT（单指令，多线程）
与 CPU 的 SIMD（单指令多数据）类似，GPU 具有单指令多线程（SIMT）特性。在传统的顺序执行 for 循环时，每个线程可以运行 for 循环的单个迭代，这样看起来就只花费一个迭代的时间。而且，当迭代数量增加时，如果有足够的核心，它可以线性地扩展执行速度（当然，如果没有足够的核心来并行处理所有独立的 for 循环迭代，速度提升会受限）。与 CPU 相比，GPU 的 SIMT 具有一些优势，例如指令按顺序发出，不需要分支预测，并且其控制相对简单，这为容纳更多的核心提供了空间。
线程层次结构
根据 thread - hierarchy 的说明，每个块中的线程数量是有限制的。因为一个块中的所有线程都期望驻留在同一个流式多处理器核心上，并且必须共享该核心有限的内存资源。https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#thread-hierarchy 告诉我们在当前的 GPU 上，一个线程块最多可以包含 1024 个线程，每个线程束包含 32 个线程，每个块理论上最多有 32 个线程束。了解这个线程层次结构对于合理设计 CUDA 程序中的线程布局和资源分配非常重要。
数学内建函数
CUDA 提供了仅设备可用的硬件指令用于基本的数学运算，这些被称为数学内建函数。虽然可以使用主机设计的操作，如 log()（主机），但使用 logf()（设备）会运行得更快。这些数学内建函数允许在设备 / GPU 上进行高效的数学运算。可以向 nvcc 编译器传递 -use_fast_math 参数来转换为这些仅设备可用的操作，不过会有几乎难以察觉的精度误差。另外，使用 --fmad = true 可以进行融合乘加操作，进一步提高某些数学计算的效率。
通过对 CUDA 这些线程同步、安全性以及相关特性的深入理解，我们能够在编写 GPU 程序时更加得心应手，充分发挥 GPU 的强大计算能力，为处理复杂的计算任务提供有力的支持。无论是深度学习中的大规模矩阵运算，还是其他科学计算领域的密集型计算，都能从这些特性的合理运用中受益。
接下来是两部分代码示例，分别是向量加法和矩阵乘法的实现以及性能分析。
vector_add
下面程序实现了一个简单的向量加法（Vector Addition）并通过 CPU 和 GPU 的 1D（一维）和 3D（三维）版本进行性能对比。
- CPU版本 (
VectorAddCpu): 使用传统的 for 循环在 CPU 上执行向量加法。 - GPU 1D版本 (
VectorAddGpu1D): 在 1D 网格和线程块结构下执行向量加法。 - GPU 3D版本 (
VecotrAddGpu3D): 在 3D 网格和线程块结构下执行向量加法，模拟三维网格的加法运算。
#include <cuda_runtime.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <iostream>
const int N = 10000000;
// 一维网格的block大小
const int BLOCK_SIZE_1D = 1024;
// 三维网格的block大小
const int BLOCK_SIZE_3D_X = 16;
const int BLOCK_SIZE_3D_Y = 8;
const int BLOCK_SIZE_3D_Z = 8;
void VectorAddCpu(float* a, float* b, float* c, int n) {
for (int i = 0; i < n; ++i) {
c[i] = a[i] + b[i];
}
}
__global__ void VectorAddGpu1D(float* a, float* b, float* c, int n) {
// 注意我们的grid是1维的，block也是1维度的
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) {
c[i] = a[i] + b[i];
}
}
__global__ void VecotrAddGpu3D(float* a, float* b, float* c, int nx, int ny,
int nz) {
// 计算当前线程在三维网格中的坐标
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
int k = blockIdx.z * blockDim.z + threadIdx.z;
if (i < nx && i < ny && k < nz) {
// 三维数组在一维化后的索引
int idx = i + j * nx + k * nx * ny;
if (idx < nx * ny * nz) {
c[idx] = a[idx] + b[idx];
}
}
}
void InitVector(float* vec, int n) {
for (int i = 0; i < n; ++i) {
vec[i] = (float)rand() / RAND_MAX;
}
}
double GetTime() {
struct timespec ts;
// CLOCK_MONOTONIC表示一个单调递增的时钟，从某个固定时间点开始（通常是系统启动时），不受系统时间的修改影响。
clock_gettime(CLOCK_MONOTONIC, &ts);
return ts.tv_sec + ts.tv_nsec * 1e-9;
}
int main() {
float *h_a, *h_b, *h_c_cpu, *h_c_gpu_1d, *h_c_gpu_3d;
float *d_a, *d_b, *d_c_1d, *d_c_3d;
size_t size = N * sizeof(float);
h_a = (float*)malloc(size);
h_b = (float*)malloc(size);
h_c_cpu = (float*)malloc(size);
h_c_gpu_1d = (float*)malloc(size);
h_c_gpu_3d = (float*)malloc(size);
srand(time(NULL));
InitVector(h_a, N);
InitVector(h_b, N);
cudaMalloc(&d_a, size);
cudaMalloc(&d_b, size);
cudaMalloc(&d_c_1d, size);
cudaMalloc(&d_c_3d, size);
cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);
// 等同于N/BLOCK_SIZE 上取整
int num_blocks = (N + BLOCK_SIZE_1D - 1) / BLOCK_SIZE_1D;
// 设定三维数组的各个维度
int nx = 100, ny = 100, nz = 1000; // 10000000
dim3 block_size_3d(BLOCK_SIZE_3D_X, BLOCK_SIZE_3D_Y, BLOCK_SIZE_3D_Z);
// 计算每个轴方向上，至少需要多少个线程块
dim3 num_blocks_3d((nx + block_size_3d.x - 1) / block_size_3d.x,
(ny + block_size_3d.y - 1) / block_size_3d.y,
(nz + block_size_3d.z - 1) / block_size_3d.z);
printf("Performing warm-up runs...\n");
for (int i = 0; i < 3; ++i) {
VectorAddCpu(h_a, h_b, h_c_cpu, N);
VectorAddGpu1D<<<num_blocks, BLOCK_SIZE_1D>>>(d_a, d_b, d_c_1d, N);
VecotrAddGpu3D<<<num_blocks_3d, block_size_3d>>>(d_a, d_b, d_c_3d, nx,
ny, nz);
cudaDeviceSynchronize();
}
printf("Benchmarking CPU implementation...\n");
double cpu_total_time = 0.0;
for (int i = 0; i < 20; ++i) {
double start_time = GetTime();
VectorAddCpu(h_a, h_b, h_c_cpu, N);
double end_time = GetTime();
cpu_total_time += end_time - start_time;
}
double cpu_avg_time = cpu_total_time / 20.0;
printf("Benchmarking GPU 1D implementation...\n");
double gpu_1d_total_time = 0.0;
for (int i = 0; i < 100; ++i) {
cudaMemset(d_c_1d, 0, size);
double start_time = GetTime();
VectorAddGpu1D<<<num_blocks, BLOCK_SIZE_1D>>>(d_a, d_b, d_c_1d, N);
cudaDeviceSynchronize();
double end_time = GetTime();
gpu_1d_total_time += end_time - start_time;
}
double gpu_1d_avg_time = gpu_1d_total_time / 100.0;
cudaMemcpy(h_c_gpu_1d, d_c_1d, size, cudaMemcpyDeviceToHost);
bool correct_1d = true;
for (int i = 0; i < N; ++i) {
if (fabs(h_c_cpu[i] - h_c_gpu_1d[i]) > 1e-4) {
correct_1d = false;
std::cout << i << " cpu: " << h_c_cpu[i] << " != " << h_c_gpu_1d[i]
<< std::endl;
break;
}
}
printf("1D Results are %s\n", correct_1d ? "correct" : "incorrect");
printf("Benchmarking GPU 3D implementation...\n");
double gpu_3d_total_time = 0.0;
for (int i = 0; i < 100; ++i) {
cudaMemset(d_c_3d, 0, size);
double start_time = GetTime();
VecotrAddGpu3D<<<num_blocks_3d, block_size_3d>>>(d_a, d_b, d_c_3d, nx,
ny, nz);
cudaDeviceSynchronize();
double end_time = GetTime();
gpu_3d_total_time += end_time - start_time;
}
double gpu_3d_avg_time = gpu_3d_total_time / 100.0;
cudaMemcpy(h_c_gpu_3d, d_c_3d, size, cudaMemcpyDeviceToHost);
bool correct_3d = true;
for (int i = 0; i < N; ++i) {
if (fabs(h_c_cpu[i] - h_c_gpu_3d[i]) > 1e-4) {
correct_3d = false;
std::cout << i << " cpu: " << h_c_cpu[i] << " != " << h_c_gpu_3d[i]
<< std::endl;
break;
}
}
printf("3D Results are %s\n", correct_3d ? "correct" : "incorrect");
printf("CPU averge time: %f milliseconds\n", cpu_avg_time * 1000);
printf("GPU 1D averge time: %f milliseconds\n", gpu_1d_avg_time * 1000);
printf("GPU 3D averge time: %f milliseconds\n", gpu_3d_avg_time * 1000);
printf("Speedup (CPU vs GPU 1D): %fx\n", cpu_avg_time / gpu_1d_avg_time);
printf("Speedup (CPU vs GPU 3D): %fx\n", cpu_avg_time / gpu_3d_avg_time);
printf("Speedup (GPU 1D vs GPU 3D): %fx\n",
gpu_1d_avg_time / gpu_3d_avg_time);
free(h_a);
free(h_b);
free(h_c_cpu);
free(h_c_gpu_1d);
free(h_c_gpu_3d);
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c_1d);
cudaFree(d_c_3d);
return 0;
}
输出如下：（设备：4060 8G）
Performing warm-up runs...
Benchmarking CPU implementation...
Benchmarking GPU 1D implementation...
1D Results are correct
Benchmarking GPU 3D implementation...
3D Results are correct
CPU averge time: 15.989620 milliseconds
GPU 1D averge time: 0.748264 milliseconds
GPU 3D averge time: 0.749384 milliseconds
Speedup (CPU vs GPU 1D): 21.368963x
Speedup (CPU vs GPU 3D): 21.337017x
Speedup (GPU 1D vs GPU 3D): 0.998505x
可以看到相对于 CPU，GPU 在执行该向量加法计算时有着大约 21 倍的加速。无论是 1D 还是 3D 配置，GPU 的表现都非常接近，表明 1D 配置对于这种任务已经足够有效。虽然使用了 3D 网格和线程块配置，性能与 1D 配置相比甚至有所下降，表明对于这种简单的向量加法任务，3D 配置并不比 1D 配置更高效。在具有明确三维布局的任务中，使用 3D 配置有助于提高计算的并行度和数据访问效率。
matmul
下面这段代码主要实现了矩阵乘法运算，并对比了在 CPU 和 GPU 上执行该运算的性能差异（通过计算平均执行时间及相应的加速比来体现）。代码中先在 CPU 上用常规的三层嵌套循环方式实现矩阵乘法，同时利用 CUDA 编写了对应的 GPU 并行版本的矩阵乘法函数，之后通过多次运行两种实现方式来进行性能测试与对比。
#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
// 3x2 @ 2x4 = 3x4
const int M = 256; // Number of rows in A and C
const int K = 512; // Number os columns in A and row in B
const int N = 256; // Number of columns in B and C
const int BLOCK_SIZE = 32;
void MatmulCPU(float* A, float* B, float* C, int m, int k, int n) {
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
float sum = 0;
for (int p = 0; p < k; ++p) {
sum += A[i * k + p] * B[p * n + j];
}
C[i * m + j] = sum;
}
}
}
__global__ void MatmulGPU(float* A, float* B, float* C, int m, int k, int n) {
// x轴方向是水平方向对应于二维数组中的列，y是竖直方向对应于行
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < m && col < n) {
float sum = 0;
for (int p = 0; p < k; ++p) {
sum += A[row * k + p] * B[p * n + col];
}
C[row * n + col] = sum;
}
}
void InitMatrix(float* mat, int rows, int cols) {
for (int i = 0; i < rows * cols; ++i) {
mat[i] = (float)rand() / RAND_MAX;
}
}
double GetTime() {
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC, &ts);
return ts.tv_sec + ts.tv_nsec * 1e-9;
}
int main() {
float *h_A, *h_B, *h_C_cpu, *h_C_gpu;
float *d_A, *d_B, *d_C;
int size_A = M * K * sizeof(float);
int size_B = K * N * sizeof(float);
int size_C = M * N * sizeof(float);
h_A = (float*)malloc(size_A);
h_B = (float*)malloc(size_B);
h_C_cpu = (float*)malloc(size_C);
h_C_gpu = (float*)malloc(size_C);
srand(time(NULL));
InitMatrix(h_A, M, K);
InitMatrix(h_B, K, N);
cudaMalloc(&d_A, size_A);
cudaMalloc(&d_B, size_B);
cudaMalloc(&d_C, size_C);
cudaMemcpy(d_A, h_A, size_A, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, size_B, cudaMemcpyHostToDevice);
// Define grid and block dimensions
dim3 blockDim(BLOCK_SIZE, BLOCK_SIZE);
dim3 gridDim((N + BLOCK_SIZE - 1) / BLOCK_SIZE,
(M + BLOCK_SIZE - 1) / BLOCK_SIZE);
// Warm-up runs
printf("Performing warm-up runs...\n");
for (int i = 0; i < 3; i++) {
MatmulCPU(h_A, h_B, h_C_cpu, M, K, N);
MatmulGPU<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, K, N);
cudaDeviceSynchronize();
}
// Benchmark CPU implementation
printf("Benchmarking CPU implementation...\n");
double cpu_total_time = 0.0;
for (int i = 0; i < 20; i++) {
double start_time = GetTime();
MatmulCPU(h_A, h_B, h_C_cpu, M, K, N);
double end_time = GetTime();
cpu_total_time += end_time - start_time;
}
double cpu_avg_time = cpu_total_time / 20.0;
// Benchmark GPU implementation
printf("Benchmarking GPU implementation...\n");
double gpu_total_time = 0.0;
for (int i = 0; i < 20; i++) {
double start_time = GetTime();
MatmulGPU<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, K, N);
cudaDeviceSynchronize();
double end_time = GetTime();
gpu_total_time += end_time - start_time;
}
double gpu_avg_time = gpu_total_time / 20.0;
// Print results
printf("CPU average time: %f microseconds\n", (cpu_avg_time * 1e6f));
printf("GPU average time: %f microseconds\n", (gpu_avg_time * 1e6f));
printf("Speedup: %fx\n", cpu_avg_time / gpu_avg_time);
// Free memory
free(h_A);
free(h_B);
free(h_C_cpu);
free(h_C_gpu);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
return 0;
}
输出如下：
Performing warm-up runs...
Benchmarking CPU implementation...
Benchmarking GPU implementation...
CPU average time: 103018.952200 microseconds
GPU average time: 108.085850 microseconds
Speedup: 953.121541x
可以发现GPU对矩阵乘法的加速是非常可观的，而深度学习底层中又充斥着大量矩阵乘法，所有GPU对深度学习的重要性可见一斑。
Atomic
在 CUDA 编程的世界里，原子操作（Atomic Operations）是确保并行计算正确性和内存安全的关键概念。今天，让我们深入探讨一下原子操作的奥秘。
一、原子操作的基本概念
原子操作的“原子”一词来源于物理学中的不可分割概念，意味着一个操作在内存位置上是完整且不可中断地被一个线程执行，在该操作完成之前，其他线程无法访问或修改相同的内存位置，这有效防止了竞态条件（Race Conditions）的出现。然而，由于在原子操作期间对单个内存块的操作进行了限制，其执行速度会略有损失，但它以牺牲一定速度为代价确保了内存操作的安全性，这种安全性是由硬件层面来保证的。
二、整数原子操作
CUDA 提供了丰富的整数原子操作函数：
- atomicAdd(int address, int val)*：原子性地将值
val加到address所指向的内存位置的值上，并返回该内存位置的旧值。例如，在多个线程对一个共享的计数器变量进行累加操作时，使用atomicAdd可以确保每个加法操作不会被其他线程干扰，从而得到正确的累加结果。 - atomicSub(int address, int val)*：与
atomicAdd类似，不过是执行减法操作，原子性地从address指向的内存位置的值中减去val，并返回旧值。 - atomicExch(int address, int val)*：原子性地将
address处的值与val进行交换，并返回address处的旧值。这在一些需要原子性地更新并获取旧值的场景中非常有用，比如在实现某些特定的状态更新机制时。 - atomicMax(int address, int val)*：原子性地将
address指向的内存位置的值设置为当前值与val中的最大值。 - atomicMin(int address, int val)*：与
atomicMax相对，将address处的值设置为当前值与val中的最小值。 - atomicAnd(int address, int val)*：原子性地对
address处的值和val执行按位与操作。 - atomicOr(int address, int val)*：原子性地对
address处的值和val执行按位或操作。 - atomicXor(int address, int val)*：原子性地对
address处的值和val执行按位异或操作。 - atomicCAS(int address, int compare, int val)*：原子性地比较
address处的值与compare值，如果相等，则将其替换为val，并返回原始值。这个操作常用于实现无锁数据结构中的一些复杂逻辑，如原子性的条件更新操作。
三、浮点原子操作
对于浮点类型，也有相应的原子操作。从 CUDA 2.0 开始，提供了 atomicAdd(float* address, float val) 函数，用于原子性地将 val 加到 address 所指向的浮点型内存位置的值上，并返回旧值。而对于双精度浮点变量的原子操作，从 CUDA Compute Capability 6.0 开始支持，使用 atomicAdd(double* address, double val) 函数。
四、原子操作的实现原理
现代 GPU 配备了特殊的硬件指令来高效地执行这些原子操作，在硬件层面常使用比较与交换（Compare-and-Swap，CAS）等技术。我们可以将原子操作类比为一种非常快速的、硬件级别的互斥锁（mutex）操作。例如，一个简单的原子加法操作可以想象成以下伪代码的执行过程：
lock(memory_location)
old_value = *memory_location
*memory_location = old_value + increment
unlock(memory_location)
return old_value
在 CUDA 中，我们甚至可以自己实现一个软件层面的原子加法操作，如以下代码示例：
__device__ int softwareAtomicAdd(int* address, int increment) {
__shared__ int lock;
int old;
if (threadIdx.x == 0) lock = 0;
__syncthreads();
while (atomicCAS(&lock, 0, 1)!= 0); // Acquire lock
old = *address;
*address = old + increment;
__threadfence(); // Ensure the write is visible to other threads
atomicExch(&lock, 0); // Release lock
return old;
}
这个自定义的 softwareAtomicAdd 函数通过使用 atomicCAS 来获取锁，然后执行加法操作，最后释放锁并返回旧值，虽然它在效率上可能不如硬件原生的原子操作，但有助于我们理解原子操作背后的互斥和同步逻辑。
五、原子操作与互斥（Mutex）
互斥（Mutex）在多线程编程中用于防止多个线程同时访问共享资源。在 CUDA 中，我们可以构建自己的互斥结构并在核函数中使用原子操作来实现互斥功能。例如：
// 我们的互斥结构
struct Mutex {
int *lock;
};
// 初始化互斥锁
__host__ void initMutex(Mutex *m) {
cudaMalloc((void**)&m->lock, sizeof(int));
int initial = 0;
cudaMemcpy(m->lock, &initial, sizeof(int), cudaMemcpyHostToDevice);
}
// 获取互斥锁
__device__ void lock(Mutex *m) {
while (atomicCAS(m->lock, 0, 1)!= 0) {
// 自旋等待
}
}
// 释放互斥锁
__device__ void unlock(Mutex *m) {
atomicExch(m->lock, 0);
}
// 核函数示例，展示互斥锁的使用
__global__ void mutexKernel(int *counter, Mutex *m) {
lock(m);
// 临界区
int old = *counter;
*counter = old + 1;
unlock(m);
}
在上述代码中，我们定义了 Mutex 结构，通过 initMutex 函数初始化互斥锁，在 mutexKernel 核函数中，使用 lock 和 unlock 函数来保护对共享计数器 counter 的操作，确保在同一时间只有一个线程能够进入临界区（Critical Section）对计数器进行修改，从而避免了数据竞争。
原子操作在 CUDA 编程中扮演着至关重要的角色，无论是在简单的共享资源计数场景，还是在复杂的无锁数据结构实现中，都离不开原子操作的支持。理解原子操作的原理和用法，能够帮助我们更好地开发高效、正确的 CUDA 并行程序，充分发挥 GPU 的强大计算能力，在深度学习、科学计算等众多领域实现更出色的性能表现。
代码示例
#include <cuda_runtime.h>
#include <stdio.h>
const int NUM_THREADS = 1000;
const int NUM_BLOCKS = 1000;
__global__ void IncrementCounterNonAtomic(int* counter) {
int old = *counter;
int new_value = old + 1;
*counter = new_value;
}
__global__ void IncrementCounterAtomic(int* counter) {
int a = atomicAdd(counter, 1);
}
double GetTime() {
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC, &ts);
return ts.tv_sec + ts.tv_nsec * 1e-9;
}
int main() {
int h_counterNonAtomic = 0;
int h_counterAtomic = 0;
int *d_counterNonAtomic, *d_counterAtomic;
// cudaError_t cudaMalloc(void** devPtr, size_t size);
// void**
// devPtr:指向设备内存指针的指针。这个指针会被赋值为设备内存的起始地址。
cudaMalloc((void**)&d_counterNonAtomic, sizeof(int));
cudaMalloc((void**)&d_counterAtomic, sizeof(int));
cudaMemcpy(d_counterNonAtomic, &h_counterNonAtomic, sizeof(int),
cudaMemcpyHostToDevice);
cudaMemcpy(d_counterAtomic, &h_counterAtomic, sizeof(int),
cudaMemcpyHostToDevice);
double start_time = GetTime();
IncrementCounterNonAtomic<<<NUM_BLOCKS, NUM_THREADS>>>(d_counterNonAtomic);
double end_time = GetTime();
double non_atomic_time = end_time - start_time;
start_time = GetTime();
IncrementCounterAtomic<<<NUM_BLOCKS, NUM_THREADS>>>(d_counterAtomic);
end_time = GetTime();
double atomic_time = end_time - start_time;
cudaMemcpy(&h_counterNonAtomic, d_counterNonAtomic, sizeof(int),
cudaMemcpyDeviceToHost);
cudaMemcpy(&h_counterAtomic, d_counterAtomic, sizeof(int),
cudaMemcpyDeviceToHost);
printf("Non-atomic counter value: %d\t time = %f\n", h_counterNonAtomic,
non_atomic_time);
printf("Atomic counter value: %d\t time = %f\n", h_counterAtomic,
atomic_time);
cudaFree(d_counterAtomic);
cudaFree(d_counterNonAtomic);
return 0;
}
Stream
在 CUDA 编程的世界里，CUDA Streams（流）是实现高效并行计算的重要概念。今天，就让我们深入探究 CUDA Streams 的奥秘，了解它如何助力我们更好地利用 GPU 资源。
一、CUDA Streams 的直观理解
想象一下“河流中的溪流”，CUDA Streams 中的操作就如同溪流中的水流，只能沿着时间的正向流动，就像一条时间轴。例如，先进行数据复制（时间步 1），接着执行一些计算（时间步 2），然后再将数据复制回来（时间步 3），这就是 Streams 的基本思想。在 CUDA 中，我们可以同时拥有多个流，并且每个流都有自己独立的时间线。这使得我们能够重叠操作，从而更充分地利用 GPU 的强大计算能力。
以训练大规模语言模型为例，若只是单纯地花费大量时间将所有的标记数据加载进 GPU 又加载出来，效率会非常低下。而 Streams 允许我们在数据移动的同时持续进行计算，它引入了一种名为“预取（prefetching）”的软件抽象概念，也就是在数据实际被需要之前就提前进行数据的移动，这样就巧妙地隐藏了数据移动所带来的延迟。
二、代码中的 Streams 体现
- 默认流与流的标识：默认流即 stream 0，也被称为 null 流。例如，当我们像这样启动一个内核：
myKernel<<<gridSize, blockSize>>>(args);，实际上它使用的就是 null 流（0），这等价于myKernel<<<gridSize, blockSize, 0, 0>>>(args);。在 CUDA 内核的执行配置中，<<<gridDim, blockDim, Ns, S>>>这种形式用于指定相关参数，其中S（cudaStream_t）就是用来指定相关联的流，它是一个可选参数，默认值为 0。 - 创建不同优先级的流：我们可以创建具有不同优先级的流，如 stream 1 和 stream 2。这意味着在运行时它们会按照特定的顺序执行，从而让我们对内核的并发执行有了更多的控制权。以下是创建具有不同优先级流的示例代码：
// 创建具有不同优先级的流
int leastPriority, greatestPriority;
CHECK_CUDA_ERROR(cudaDeviceGetStreamPriorityRange(&leastPriority, &greatestPriority));
CHECK_CUDA_ERROR(cudaStreamCreateWithPriority(&stream1, cudaStreamNonBlocking, leastPriority));
CHECK_CUDA_ERROR(cudaStreamCreateWithPriority(&stream2, cudaStreamNonBlocking, greatestPriority));
三、相关重要概念
- Pinned Memory（固定内存）：可以将其理解为被锁定在特定位置且不能被操作系统随意移动的内存。当我们想要将数据移动到 GPU 并在其上进行计算时，它就发挥了重要作用。因为如果操作系统移动了数据，GPU 在寻找数据时就会出错，导致段错误（segfault）。以下是分配固定内存的代码示例：
// 分配固定内存
float* h_data;
cudaMallocHost((void**)&h_data, size);
- Events（事件）：
- 测量内核执行时间：通过在核函数启动前后放置事件，可以精确地测量内核的执行时间。例如：
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, stream);
kernel<<<grid, block, 0, stream>>>(args);
cudaEventRecord(stop, stream);
cudaEventSynchronize(stop);
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
- 流间同步：事件还可以用于创建不同流之间的依赖关系，确保一个操作在另一个操作完成之后才开始。
- 重叠计算与数据传输：能够标记数据传输的完成，从而通知计算可以在该数据上开始进行。
- Callbacks（回调函数）：借助回调函数，我们可以构建一个计算管道。当 GPU 上的一个操作完成时，它可以触发 CPU 上的另一个操作开始，而 CPU 上的这个操作可能又会为 GPU 安排更多的工作任务。例如：
void CUDART_CB MyCallback(cudaStream_t stream, cudaError_t status, void *userData) {
printf("GPU operation completed\n");
// 触发下一批工作
}
kernel<<<grid, block, 0, stream>>>(args);
cudaStreamAddCallback(stream, MyCallback, nullptr, 0);
CUDA Streams 为我们在 GPU 编程中提供了强大的工具，无论是优化计算资源的利用，还是处理复杂的多任务并行场景，它都有着不可替代的作用。通过深入理解 Streams 及其相关概念，并熟练运用到实际的 CUDA 编程中，我们能够让 GPU 发挥出更强大的性能，在深度学习、科学计算等众多领域取得更好的计算效果。
stream_basics
#include <cuda_runtime.h>
#include <stdio.h>
// #val：CUDA 函数的字符串形式（通过宏的字符串化操作 # 获得）
#define CHECK_CUDA_ERROR(val) Check((val), #val, __FILE__, __LINE__)
template <typename T>
void Check(T err, const char* func, const char* const file, const int line) {
if (err != cudaSuccess) {
fprintf(stderr, "CUDA error at %s:%d code=%d(%s) \"%s\" \n", file, line,
static_cast<unsigned int>(err), cudaGetErrorString(err), func);
exit(EXIT_FAILURE);
}
}
__global__ void VectorAdd(const float* A, const float* B, float* C,
int num_elements) {
int i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < num_elements) {
C[i] = A[i] + B[i];
}
}
int main(void) {
int num_elements = 50000;
size_t size = num_elements * sizeof(float);
float *h_A, *h_B, *h_C;
float *d_A, *d_B, *d_C;
cudaStream_t stream1, stream2;
h_A = (float*)malloc(size);
h_B = (float*)malloc(size);
h_C = (float*)malloc(size);
for (int i = 0; i < num_elements; ++i) {
h_A[i] = rand() / (float)RAND_MAX;
h_B[i] = rand() / (float)RAND_MAX;
}
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_A, size));
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_B, size));
CHECK_CUDA_ERROR(cudaMalloc((void**)&d_C, size));
// 创建流
CHECK_CUDA_ERROR(cudaStreamCreate(&stream1));
CHECK_CUDA_ERROR(cudaStreamCreate(&stream2));
// 异步地将Host数据复制到Device
CHECK_CUDA_ERROR(
cudaMemcpyAsync(d_A, h_A, size, cudaMemcpyHostToDevice, stream1));
CHECK_CUDA_ERROR(
cudaMemcpyAsync(d_B, h_B, size, cudaMemcpyHostToDevice, stream2));
// 启动核
int threads_per_block = 256;
int block_per_grid =
(num_elements + threads_per_block - 1) / threads_per_block;
VectorAdd<<<block_per_grid, threads_per_block, 0, stream1>>>(d_A, d_B, d_C,
num_elements);
// 异步地将数据复制回Host
CHECK_CUDA_ERROR(
cudaMemcpyAsync(h_C, d_C, size, cudaMemcpyDeviceToHost, stream1));
// 同步流
CHECK_CUDA_ERROR(cudaStreamSynchronize(stream1));
CHECK_CUDA_ERROR(cudaStreamSynchronize(stream2));
for (int i = 0; i < num_elements; ++i) {
if (fabs(h_A[i] + h_B[i] - h_C[i]) > 1e-5) {
fprintf(stderr, "Result verification failed at element %d!\n", i);
exit(EXIT_FAILURE);
}
}
printf("Test PASSED\n");
CHECK_CUDA_ERROR(cudaFree(d_A));
CHECK_CUDA_ERROR(cudaFree(d_B));
CHECK_CUDA_ERROR(cudaFree(d_C));
CHECK_CUDA_ERROR(cudaStreamDestroy(stream1));
CHECK_CUDA_ERROR(cudaStreamDestroy(stream2));
free(h_A);
free(h_B);
free(h_C);
return 0;
}
stream_advanced
#include <cuda_runtime.h>
#include <stdio.h>
#include <iostream>
#define CHECK_CUDA_ERROR(val) Check((val), #val, __FILE__, __LINE__)
template <typename T>
void Check(T err, const char* const func, const char* const file,
const int line) {
if (err != cudaSuccess) {
fprintf(stderr, "CUDA error at %s:%d code=%d(%s) \"%s\" \n", file, line,
static_cast<unsigned int>(err), cudaGetErrorString(err), func);
exit(EXIT_FAILURE);
}
}
__global__ void Kernel1(float* data, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
data[idx] *= 2.0f;
}
}
__global__ void Kernel2(float* data, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
data[idx] += 1.0f;
}
}
// CUDART_CB 是一个宏，用于标记回调函数的声明，以确保它与 CUDA
// 运行时的调用约定兼容。
void CUDART_CB MyStreamCallback(cudaStream_t stream, cudaError_t status,
void* userData) {
printf("Stream callback: Operation completed\n");
}
int main(void) {
const int N = 1000000;
size_t size = N * sizeof(float);
float *h_data, *d_data;
cudaStream_t stream1, stream2;
cudaEvent_t event;
std::cout << event << std::endl;
// cudaMallocHost在主机（Host）上分配固定内存，不会被OS换出到磁盘
CHECK_CUDA_ERROR(cudaMallocHost(&h_data, size));
CHECK_CUDA_ERROR(cudaMalloc(&d_data, size));
for (int i = 0; i < N; ++i) {
h_data[i] = static_cast<float>(i);
}
// 创建不同优先级的流
int leaset_proority, greatest_priority;
// 获取当前设备支持的流优先级范围
CHECK_CUDA_ERROR(
cudaDeviceGetStreamPriorityRange(&leaset_proority, &greatest_priority));
// 创建一个流（stream1），具有指定优先级（最低优先级）。
// cudaStreamNonBlocking标志表示该流是非阻塞流。
CHECK_CUDA_ERROR(cudaStreamCreateWithPriority(
&stream1, cudaStreamNonBlocking, leaset_proority));
CHECK_CUDA_ERROR(cudaStreamCreateWithPriority(
&stream2, cudaStreamNonBlocking, greatest_priority));
// 创建事件
CHECK_CUDA_ERROR(cudaEventCreate(&event));
// 异步操作和事件记录
CHECK_CUDA_ERROR(
cudaMemcpyAsync(d_data, h_data, size, cudaMemcpyHostToDevice, stream1));
Kernel1<<<(N + 255) / 256, 256, 0, stream1>>>(d_data, N);
// 记录事件，在 stream1 中记录事件，表示 stream1 的某些操作完成。
CHECK_CUDA_ERROR(cudaEventRecord(event, stream1));
// 让 stream2 等待 event 触发后才开始执行其后续操作。
CHECK_CUDA_ERROR(cudaStreamWaitEvent(stream2, event, 0));
Kernel2<<<(N + 255) / 256, 256, 0, stream2>>>(d_data, N);
// 为 stream2 添加一个回调函数 MyStreamCallback。
// 当 stream2 的所有任务完成后，回调函数会被主机线程调用。
CHECK_CUDA_ERROR(cudaStreamAddCallback(stream2, MyStreamCallback, NULL, 0));
CHECK_CUDA_ERROR(
cudaMemcpyAsync(h_data, d_data, size, cudaMemcpyDeviceToHost, stream2));
CHECK_CUDA_ERROR(cudaStreamSynchronize(stream1));
CHECK_CUDA_ERROR(cudaStreamSynchronize(stream2));
for (int i = 0; i < N; ++i) {
float expected = (static_cast<float>(i) * 2.0f) + 1.0f;
if (fabs(h_data[i] - expected) > 1e-5) {
fprintf(stderr, "Result verification failed at element %d!\n", i);
exit(EXIT_FAILURE);
}
}
printf("Test PASSED\n");
CHECK_CUDA_ERROR(cudaFreeHost(h_data));
CHECK_CUDA_ERROR(cudaFree(d_data));
CHECK_CUDA_ERROR(cudaStreamDestroy(stream1));
CHECK_CUDA_ERROR(cudaStreamDestroy(stream2));
CHECK_CUDA_ERROR(cudaEventDestroy(event));
}
参考：https://github.com/Infatoshi/cuda-course/tree/master/05_Writing_your_First_Kernels