向量加法

agicy2026/3/9大约 13 分钟

向量加法

题目描述

编写一个 GPU 程序，对两个包含 32 位浮点数的向量执行逐元素加法。程序应接收两个长度相等的输入向量，并生成一个包含它们之和的输出向量。

实现要求

不允许使用外部库。
solve 函数签名必须保持不变。
最终结果必须存储在向量 $\mathbf{C}$ 中。

示例

示例 1

Input:  A = [1.0, 2.0, 3.0, 4.0]
        B = [5.0, 6.0, 7.0, 8.0]
Output: C = [6.0, 8.0, 10.0, 12.0]

示例 2

Input:  A = [1.5, 1.5, 1.5]
        B = [2.3, 2.3, 2.3]
Output: C = [3.8, 3.8, 3.8]

约束条件

输入向量 $\mathbf{A}$ 和 $\mathbf{B}$ 长度相同。
$1 \le N \le 100,000,000$ 。
性能测试将在 $N = 25,000,000$ 的规模下进行。

解题思路

向量加法是一个典型的数据并行（Data Parallelism）问题。由于每个元素的加法运算 $C[i] = A[i] + B[i]$ 都是独立的，我们可以利用 GPU 的大规模并行线程来同时计算多个元素。

模式分析

本题属于经典的 Map（映射）模式。在这种模式中，每个输入元素独立地通过相同的函数（此处为加法操作）映射到一个输出元素。由于计算过程没有数据依赖（如前缀和或归约），各线程之间完全不需要同步，这使得它成为最易于并行的算法结构。

从数据访问角度看，这是一个 Gather（聚集） 类型的操作：每个线程读取两个输入源 $A[i]$ 和 $B[i]$ ，并写入一个目标 $C[i]$ 。由于输入和输出索引一一对应且对齐，它也符合 Output-Centered（输出中心） 的设计思想，即以输出位置为锚点，分配线程去计算该位置的值。

效率分析

算术强度（Arithmetic Intensity）

向量加法的计算非常轻量，每个线程执行 1 次浮点加法操作，但需要进行 3 次 4 字节的全局内存访问（读取 $A[i]$ ，读取 $B[i]$ ，写入 $C[i]$ ）。

计算量：1 FLOP
访存量： $3 \times 4 = 12$ Bytes
算术强度： $1 / 12 \approx 0.083$ FLOPs/Byte

受限类型

由于算术强度极低，远低于任何现代 GPU 的峰值计算/带宽比（通常在 10~100+ FLOPs/Byte），因此这是一个典型的 Memory Bandwidth Bound（显存带宽受限） 内核。

程序的性能瓶颈完全取决于 GPU 的显存带宽，而非计算核心的算力。即使使用算力更强的 GPU，如果显存带宽没有提升，该程序的运行速度也不会有显著改善。优化方向应集中在确保内存合并访问（Coalesced Access）以最大化带宽利用率。

调度策略

采用 一维线性网格（1D Grid） 调度是最自然的选择。我们将 $N$ 个任务均匀分配给多个线程块（Thread Block）。

Grid 维度：根据向量总长度 $N$ 和每个 Block 的线程数动态计算。
Block 维度：通常选择 128、256 或 512 等 32 的倍数，以匹配 GPU 的 Warp 调度机制，最大化 Occupancy（占用率）。
线程索引：使用 $idx = blockIdx.x \times blockDim.x + threadIdx.x$ 计算全局唯一索引，直接对应数据下标。

这种简单的调度方式配合 GPU 的硬件调度器，能够自动处理成千上万个 Block 的执行顺序，无需用户手动管理任务队列。

核心策略

线程映射：将每个待计算的元素索引 $i$ 映射到一个唯一的 GPU 线程。
全局索引计算：利用 CUDA 内置变量计算当前线程对应的全局索引。
- blockIdx.x：当前线程块的索引。
- blockDim.x：每个线程块包含的线程数。
- threadIdx.x：当前线程在块内的索引。
- 全局索引公式：idx = blockIdx.x * blockDim.x + threadIdx.x。
边界检查：确保计算出的索引 idx 小于向量长度 $N$ ，防止越界访问。

进阶优化：网格跨步循环 (Grid-Stride Loop)

虽然直接的一对一映射（每个线程处理一个元素）对于 $N$ 小于最大线程总数的情况是可行的，但在实际应用中， $N$ 往往远大于 GPU 能同时调度的线程数。此外，为了提高程序的健壮性和灵活性，推荐使用 网格跨步循环。

在内存带宽受限（Memory Bound）的场景下，通过限制 Block 的总数量并采用 Grid-Stride Loop，可以带来显著的资源优势：

减少上下文开销：GPU 需要为每个 Block 分配寄存器和共享内存资源。过多的 Block（尤其是当 $N$ 非常大时）会增加调度器的负担。限制 Grid 大小（例如设为 GPU SM 数量的 8~16 倍）可以减少这些管理开销。
避免资源争抢：在某些极端情况下，如果 Kernel 占用了过多的寄存器或共享内存，过大的 Grid 可能导致后续 Block 无法及时被调度，甚至在复杂的生产环境中与其他 Kernel 争抢资源。
提升指令级并行：每个线程串行处理多个元素，有助于掩盖指令流水线的延迟，增加单个线程的有效工作量。
原理：每个线程不仅处理索引 idx 的元素，还会在处理完当前元素后，跳过整个网格的大小（Grid Size），去处理下一个需要计算的元素 idx + gridDim.x * blockDim.x，直到覆盖所有数据。
优势：
- 解耦：代码逻辑与具体的 Grid 尺寸解耦，无论启动多少个 Block，程序都能正确运行。
- 重用：增加了每个线程的工作量，有助于掩盖指令延迟和初始化开销。
- 调试：即使配置为 1 个 Block 1 个 Thread，程序也能串行正确执行，便于调试。

代码实现

你的代码是完全正确的，这是非常标准的 Element-Wise Mapping（逐元素映射）解法。

解法一：逐元素映射 (Element-Wise Mapping)

这是最直观的 GPU 编程模式：每个线程负责计算输出向量中的一个元素。

#include <cuda_runtime.h>

/**
 * @brief CUDA Kernel for vector addition
 * 
 * 每个线程计算一个元素：C[i] = A[i] + B[i]
 */
__global__ void vector_add(const float* A, const float* B, float* C, int N) {
    // 计算全局索引
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    
    // 边界检查：防止线程索引超出向量长度
    if (i < N) {
        C[i] = A[i] + B[i];
    }
}

/**
 * @brief Host function to launch the kernel
 */
extern "C" void solve(const float* A, const float* B, float* C, int N) {
    // 线程块大小：通常选择 128, 256 或 512
    int threadsPerBlock = 256;
    
    // 网格大小：向上取整，确保有足够的线程覆盖所有 N 个元素
    // 公式：(N + threadsPerBlock - 1) / threadsPerBlock
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    // 启动 Kernel
    vector_add<<<blocksPerGrid, threadsPerBlock>>>(A, B, C, N);
    
    // 同步：确保 GPU 计算完成后再返回（虽然某些 OJ 可能隐式处理，但显式同步是好习惯）
    cudaDeviceSynchronize();
}

代码解析

索引计算：
int i = blockIdx.x * blockDim.x + threadIdx.x;
这是 CUDA 中最基础的线性索引公式，能够将二维的线程层次结构（Grid -> Block -> Thread）映射到一维的线性数据空间。
边界检查：
if (i < N)
由于 blocksPerGrid 是向上取整计算的，启动的总线程数 blocksPerGrid * threadsPerBlock 往往会略大于 N。如果不加这个检查，多余的线程会访问数组越界内存，导致未定义行为或程序崩溃。
Kernel 配置：
<<<blocksPerGrid, threadsPerBlock>>>
这种配置保证了即便是 $N = 100,000,000$ 这样的大数组，也能一次性启动足够的线程来处理。在现代 GPU 上，这种配置是非常高效的。

进阶优化：网格跨步循环 (Grid-Stride Loop)

虽然你的代码对于本题已经足够完美，但在工业界代码中，我们常采用 网格跨步循环。它的好处是解耦了数据规模 N 和网格大小 blocksPerGrid。即使限制了 blocksPerGrid 的数量（例如为了控制并发量或适应特定硬件限制），代码依然能正确运行。

__global__ void vector_add_grid_stride(const float* A, const float* B, float* C, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;

    // 线程在处理完当前元素后，跳跃 stride 个位置处理下一个元素
    for (int i = idx; i < N; i += stride) {
        C[i] = A[i] + B[i];
    }
}

/**
 * @brief 进阶优化：向量化内存访问 (Vectorized Memory Access)
 * 
 * 使用 float4 类型一次性加载/存储 4 个浮点数。
 * 这可以显著减少内存指令的数量，提高显存带宽利用率。
 * 
 * 注意：输入数组的地址必须是对齐的（通常 cudaMalloc 分配的地址满足要求）。
 */
__global__ void vector_add_float4(const float* A, const float* B, float* C, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;

    // 1. 处理主要的对齐部分
    // 将指针强制转换为 float4*，每个线程处理 4 个元素
    // 注意：这里假设 N 足够大且内存对齐，实际工程中需要更严谨的边界处理
    int N4 = N / 4;
    const float4* A4 = reinterpret_cast<const float4*>(A);
    const float4* B4 = reinterpret_cast<const float4*>(B);
    float4* C4 = reinterpret_cast<float4*>(C);

    for (int i = idx; i < N4; i += stride) {
        float4 a = A4[i]; // 一次加载 4 个 float (128 bits)
        float4 b = B4[i];
        float4 c;
        
        c.x = a.x + b.x;
        c.y = a.y + b.y;
        c.z = a.z + b.z;
        c.w = a.w + b.w;
        
        C4[i] = c; // 一次存储 4 个 float
    }

    // 2. 处理剩余的元素 (Peeling / Remainder Loop)
    // 从 4 的倍数位置开始，逐个处理剩余元素
    int remainder_start = N4 * 4;
    for (int i = remainder_start + idx; i < N; i += stride) {
        C[i] = A[i] + B[i];
    }
}

进阶优化二：向量化内存访问 (Vectorized Memory Access)

除了 Grid-Stride Loop，另一个极其重要的优化手段是 向量化内存访问。

原理：GPU 的显存控制器通常以 32 字节或 128 字节为粒度进行事务传输。如果每个线程只读取一个 float (4 字节)，虽然 Warp 内的线程合并访问可以利用带宽，但使用更宽的数据类型（如 float4，16 字节）可以进一步减少内存指令的总数，降低指令发射开销，并提高总线利用率。
实现：使用 CUDA 内置的 float4 类型，将指针强制转换后进行读写。
注意：必须确保内存地址是 16 字节对齐 的（cudaMalloc 分配的内存默认是对齐的），并且要正确处理无法被 4 整除的尾部数据。

进阶优化三：异步拷贝与双缓冲 (Async Copy & Double Buffering)

在 Ampere (SM80) 及之后的架构（如 A100, H100, RTX 30/40 系列）中，可以使用 异步数据拷贝 (LDGSTS) 指令，配合 流水线 (Pipeline) 机制，将从全局内存（Global Memory）到共享内存（Shared Memory）的数据加载与计算完全重叠。

虽然向量加法这种简单的 Element-wise 操作通常不需要 Shared Memory，但为了展示“最强大”的访存优化技术，我们可以构建一个基于 cp.async 的流水线模型。

#include <cuda_runtime.h>
#include <cuda/pipeline> // 需要 CUDA 11+

/**
 * @brief 使用 Ampere 架构异步拷贝的向量加法
 * 
 * 利用 Shared Memory 作为缓冲区，实现 Global -> Shared 的异步加载
 * 以及 Shared -> Register 的计算流水线。
 * 
 * 注意：这种方法对于简单的向量加法是"杀鸡用牛刀"，
 * 主要用于展示在计算密集型任务（如 GEMM, Convolution）中的核心优化范式。
 */
__global__ void vector_add_async(const float* A, const float* B, float* C, int N) {
    // 共享内存缓冲区：双缓冲 (Double Buffering)
    // 两个阶段：0 (加载下一块), 1 (计算当前块)
    // 每个线程处理 4 个 float，Block 处理 blockSize * 4 个 float
    extern __shared__ float4 shared_mem[]; 
    float4* smem_A = shared_mem; // A 的缓冲区
    float4* smem_B = smem_A + blockDim.x * 2; // B 的缓冲区紧随 A 之后
    
    int tid = threadIdx.x;
    int idx = blockIdx.x * blockDim.x + tid;
    int stride = gridDim.x * blockDim.x;
    
    // 向量化指针
    const float4* A4 = reinterpret_cast<const float4*>(A);
    const float4* B4 = reinterpret_cast<const float4*>(B);
    float4* C4 = reinterpret_cast<float4*>(C);
    int N4 = N / 4;

    // 创建流水线对象
    cuda::pipeline<cuda::thread_scope_thread> pipe = cuda::make_pipeline();

    // -----------------------------------------------------------
    // 阶段 1: 预取 (Prefetch) 第一块数据到 Shared Memory (Buffer 0)
    // -----------------------------------------------------------
    int i = idx;
    if (i < N4) {
        // 异步拷贝 Global -> Shared
        // 使用 pipe.producer_acquire() 获取“生产”权限（此处简化）
        pipe.producer_acquire();
        // cp.async: 从 A4[i] 拷贝 16 字节到 smem_A[tid] (Buffer 0)
        cuda::memcpy_async(&smem_A[tid], &A4[i], sizeof(float4), pipe);
        cuda::memcpy_async(&smem_B[tid], &B4[i], sizeof(float4), pipe);
        pipe.producer_commit(); // 提交本次批次的异步拷贝请求
    }

    // -----------------------------------------------------------
    // 阶段 2: 主循环流水线
    // -----------------------------------------------------------
    // 我们使用简单的单缓冲逻辑演示核心 API，真实双缓冲需要轮转 buffer索引
    // 这里为了逻辑清晰，展示 "Wait -> Compute -> Store -> Fetch Next" 模式
    
    for (; i < N4; i += stride) {
        // 等待当前批次拷贝完成
        pipe.consumer_wait(); 
        
        // 计算：从 Shared Memory 读取到寄存器进行加法
        float4 a = smem_A[tid];
        float4 b = smem_B[tid];
        float4 c;
        c.x = a.x + b.x; c.y = a.y + b.y; 
        c.z = a.z + b.z; c.w = a.w + b.w;
        
        // 写回结果到 Global Memory
        C4[i] = c;
        
        // 释放消费权限
        pipe.consumer_release();

        // 预取下一块数据 (如果还有)
        int next_i = i + stride;
        if (next_i < N4) {
            pipe.producer_acquire();
            cuda::memcpy_async(&smem_A[tid], &A4[next_i], sizeof(float4), pipe);
            cuda::memcpy_async(&smem_B[tid], &B4[next_i], sizeof(float4), pipe);
            pipe.producer_commit();
        }
    }
}

原理：
- cp.async (PTX: cp.async.ca.shared.global)：绕过寄存器，直接由硬件 DMA 引擎将数据从 Global Memory 搬运到 Shared Memory。
- 计算与传输重叠：在数据搬运的同时，SM 可以执行其他独立的指令（如上一块数据的计算）。
- 流水线 (Pipeline)：通过多级缓冲（Multi-stage buffering），让 GPU 永远处于“既在算，又在搬”的饱和状态。

注：对于纯向量加法，直接从 Global Memory 读到寄存器（LDG 指令）通常已经足够快且能通过 L1/L2 Cache 缓存，cp.async 带来的额外 Shared Memory 读写开销可能反而得不偿失。但在矩阵乘法 (GEMM) 或卷积等重度复用数据的场景下，这是目前最强大的优化范式。

进阶优化四：零拷贝 (Zero-Copy) 与统一内存 (Unified Memory)

如果 $N$ 极大（甚至超过了显存容量），或者数据本身就在 Host 端且只用一次，那么传统的 cudaMemcpy (Host->Device) -> Kernel -> cudaMemcpy (Device->Host) 流程可能不是最优的。

我们可以利用 NVIDIA GPU 的 Unified Memory (UM) 或 Zero-Copy 技术，让 GPU 直接通过 PCIe/NVLink 总线访问 Host 内存。

适用场景：
- 数据量大到无法全部放入显存。
- 数据只被访问一次（Streaming access），此时传输开销与计算开销无法通过多次复用被摊薄。
- 异构计算平台（如 NVIDIA Jetson, Grace Hopper），CPU 和 GPU 共享物理内存。

/**
 * @brief 使用 Zero-Copy 的 Host 代码示例
 * 
 * 注意：Kernel 代码本身不需要改变，变的是内存分配和调用方式。
 */
void solve_zero_copy(const float* A_h, const float* B_h, float* C_h, int N) {
    float *A_pinned, *B_pinned, *C_pinned;
    
    // 1. 分配 Pinned Memory (页锁定内存)
    // 这种内存物理地址固定，GPU 可以直接通过 PCIe DMA 访问
    cudaHostAlloc(&A_pinned, N * sizeof(float), cudaHostAllocMapped);
    cudaHostAlloc(&B_pinned, N * sizeof(float), cudaHostAllocMapped);
    cudaHostAlloc(&C_pinned, N * sizeof(float), cudaHostAllocMapped);
    
    // 复制数据到 Pinned Memory (这一步是 CPU 内存拷贝，非常快)
    memcpy(A_pinned, A_h, N * sizeof(float));
    memcpy(B_pinned, B_h, N * sizeof(float));

    // 2. 获取 Device 指针
    float *A_d, *B_d, *C_d;
    cudaHostGetDevicePointer(&A_d, A_pinned, 0);
    cudaHostGetDevicePointer(&B_d, B_pinned, 0);
    cudaHostGetDevicePointer(&C_d, C_pinned, 0);

    // 3. 启动 Kernel
    // GPU 线程执行时，会直接通过 PCIe 读取 Host 端的 A_pinned 和 B_pinned
    // 并将结果写回 C_pinned
    // 这一步彻底省去了显式的 cudaMemcpy H2D 和 D2H
    int blockSize = 256;
    int numBlocks = (N + blockSize - 1) / blockSize;
    vector_add<<<numBlocks, blockSize>>>(A_d, B_d, C_d, N);
    
    cudaDeviceSynchronize();
    
    // 结果已经直接在 C_pinned (Host 端) 可用了
    memcpy(C_h, C_pinned, N * sizeof(float));
    
    cudaFreeHost(A_pinned);
    cudaFreeHost(B_pinned);
    cudaFreeHost(C_pinned);
}

优势：
- 流水线隐式重叠：GPU 读取 A/B 时，PCIe 传输自动与计算重叠。
- 突破显存限制：可以处理远超显存大小的数据集。
劣势：
- PCIe 瓶颈：PCIe 带宽（如 Gen4 x16 ~32GB/s）远低于显存带宽（如 A100 ~2000GB/s）。如果计算不够密集，GPU 可能会长时间等待 PCIe 数据，导致性能严重下降。

对于本题来说，你的写法更加简洁明了，完全符合要求。

复杂度分析

时间复杂度： $O(N/P)$ ，其中 $N$ 是向量长度， $P$ 是 GPU 上的并行处理器数量。在理想情况下，计算被完美并行化。
空间复杂度： $O(1)$ ，除了输入输出数组外，只需要极常数的额外空间用于存储索引变量。