第四节 第一个CUDA程序：向量加法的完整实践

从这一节开始，我们正式进入CUDA C的编码世界。向量加法（Vector Addition）是并行计算领域的“Hello World”，它虽然简单，却涵盖了CUDA编程模型中最核心的四个要素：内核函数（Kernel）、线程层次结构（Thread Hierarchy）、内存模型（Memory Model）以及主机与设备之间的数据搬运。通过亲手编写并运行这个程序，你将直观地理解<<< >>>启动语法背后的执行逻辑，掌握线程索引与数据下标之间的映射关系，并建立起一套完整的编译、运行与调试流程。本节将按照“原理剖析 → 代码实现 → 编译运行 → 调试排障”的路径，带你完成这一关键里程碑。

上图展示了向量加法程序的核心执行模型。主机端负责准备数据并通过PCIe总线拷贝到设备端全局内存，然后以网格（Grid）和线程块（Block）的形式组织成千上万个线程并行执行内核函数。每个线程通过内置变量blockIdx.x、threadIdx.x和blockDim.x计算出自己负责的数据元素下标，从而实现对数组元素的并行处理。理解这张图，就理解了CUDA并行编程的“骨架”。

4.1 内核启动语法与线程层次结构

4.1.1 核心原理：从串行到并行的思维转变

在传统C语言中，一个简单的for循环即可完成向量加法：

for (int i = 0; i < N; i++) {
    C[i] = A[i] + B[i];
}

这段代码在单核CPU上按顺序执行N次加法。CUDA的并行思路是：将循环体“向量化”为内核函数，让GPU的数千个线程同时执行这个循环体。每个线程只处理一个（或少数几个）数组元素，所有线程并行推进，从而将时间复杂度从O(N)降为O(1)（理想情况下）。

内核函数使用__global__修饰符声明，其返回值必须为void。启动内核时使用尖括号语法<<<gridDim, blockDim>>>，其中：

• gridDim：网格维度，即线程块（Block）的数量，可以是整数或dim3类型（支持二维/三维网格）。
• blockDim：线程块维度，即每个Block内的线程数量，同样支持dim3类型。

例如vecAdd<<<4, 256>>>(d_A, d_B, d_C, N);表示启动4个Block，每个Block包含256个线程，总共1024个线程参与计算。硬件上，GPU会将Block调度到流多处理器（SM）上执行，每个Block内的线程进一步划分为32个线程一组的线程束（Warp）进行调度。

4.1.2 线程索引的计算方法

线程在网格中的全局唯一ID（即数组下标）由以下公式计算：

int i = blockIdx.x * blockDim.x + threadIdx.x;

这个公式是CUDA编程中使用频率最高的表达式。其含义是：当前线程所在Block的编号（blockIdx.x）乘以每个Block的线程数（blockDim.x），再加上线程在Block内的编号（threadIdx.x）。对于一维场景，这个公式直观且高效。对于二维或三维场景，需要分别计算x、y、z三个维度的索引，并通常将y维度乘以x维度的宽度来线性化。

值得注意的是，网格和Block的维度上限：每个Grid的x/y/z维度最大为2^31-1、65535、65535；每个Block的线程总数不能超过1024，且x维度不能超过1024，y/z维度不能超过1024。这些限制由硬件架构决定，编程时必须遵守。

4.1.3 案例拆解：一个典型的向量加法内核

下面是一个完整的向量加法内核函数，注意其中包含边界检查——这是生产级代码的必备要素，因为当N不能被blockDim.x整除时，最后一个Block的部分线程会越界：

__global__ void vecAdd(float *A, float *B, float *C, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {  // 边界检查，防止越界访问
        C[i] = A[i] + B[i];
    }
}

假设N=1000，我们设置blockDim.x=256，则需要的Block数量为ceil(1000/256)=4。前3个Block（共768个线程）处理下标0~767，第4个Block的256个线程中只有前232个（下标768~999）执行有效计算，其余24个线程因不满足i < N而直接退出。这种“部分线程闲置”的现象在并行计算中不可避免，但通过合理选择Block大小可以最小化浪费。

4.2 完整程序编写与编译运行流程

4.2.1 程序骨架：主机代码与设备代码的协作

一个完整的CUDA程序包含主机代码（运行在CPU上）和设备代码（运行在GPU上）。主机代码负责内存管理、数据搬运和内核启动。以下是向量加法的完整实现：

#include <stdio.h>
#include <cuda_runtime.h>

#define N 1000
#define BLOCK_SIZE 256

// 内核函数（设备代码）
__global__ void vecAdd(float *A, float *B, float *C, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {
        C[i] = A[i] + B[i];
    }
}

int main() {
    // 1. 分配主机内存并初始化数据
    float *h_A = (float*)malloc(N * sizeof(float));
    float *h_B = (float*)malloc(N * sizeof(float));
    float *h_C = (float*)malloc(N * sizeof(float));
    for (int i = 0; i < N; i++) {
        h_A[i] = i * 1.0f;
        h_B[i] = i * 2.0f;
    }

    // 2. 分配设备全局内存
    float *d_A, *d_B, *d_C;
    cudaMalloc((void**)&d_A, N * sizeof(float));
    cudaMalloc((void**)&d_B, N * sizeof(float));
    cudaMalloc((void**)&d_C, N * sizeof(float));

    // 3. 将输入数据从主机拷贝到设备
    cudaMemcpy(d_A, h_A, N * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, N * sizeof(float), cudaMemcpyHostToDevice);

    // 4. 计算网格和Block维度，启动内核
    int gridSize = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
    vecAdd<<<gridSize, BLOCK_SIZE>>>(d_A, d_B, d_C, N);

    // 5. 将结果从设备拷贝回主机
    cudaMemcpy(h_C, d_C, N * sizeof(float), cudaMemcpyDeviceToHost);

    // 6. 验证结果
    for (int i = 0; i < N; i++) {
        if (h_C[i] != h_A[i] + h_B[i]) {
            printf("Error at index %d: %f + %f != %f\n", i, h_A[i], h_B[i], h_C[i]);
            break;
        }
    }
    printf("Vector addition completed successfully!\n");

    // 7. 释放内存
    free(h_A); free(h_B); free(h_C);
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    return 0;
}

4.2.2 编译与运行：nvcc编译器使用详解

CUDA程序使用nvcc编译器进行编译。nvcc是NVIDIA提供的专用编译器驱动，它会将设备代码（__global__函数）编译为PTX（并行线程执行）中间代码，再进一步编译为适用于具体GPU架构的SASS（流处理器汇编）代码。编译命令如下：

nvcc -arch=sm_80 -o vecAdd vecAdd.cu

参数说明：

• -arch=sm_80：指定目标GPU架构。sm_80对应Ampere架构（如A100），sm_86对应RTX 30系列，sm_90对应Hopper架构（H100）。如果不指定，nvcc会使用默认值，但可能无法发挥新硬件的全部特性。
• -o vecAdd：指定输出可执行文件名。
• vecAdd.cu：源文件，.cu扩展名是CUDA源文件的标准后缀。

运行程序：

./vecAdd

如果一切正常，终端将输出“Vector addition completed successfully!”。

4.2.3 调试流程：从错误信息到问题定位

CUDA编程中，错误排查是必备技能。常见的错误类型包括：

推荐在每次CUDA API调用后检查返回值：

cudaError_t err = cudaMemcpy(d_A, h_A, N * sizeof(float), cudaMemcpyHostToDevice);
if (err != cudaSuccess) {
    printf("cudaMemcpy failed: %s\n", cudaGetErrorString(err));
    exit(1);
}

此外，cuda-memcheck工具（现已更名为compute-sanitizer）可以检测内核中的越界访问和内存泄漏：

compute-sanitizer ./vecAdd

4.3 性能观察与优化初探

4.3.1 使用CUDA事件测量内核执行时间

向量加法虽然简单，但测量其执行时间有助于建立性能基准。CUDA事件（Event）是官方推荐的计时工具，它基于GPU硬件时钟，精度远高于CPU端的clock()函数：

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start);
vecAdd<<<gridSize, BLOCK_SIZE>>>(d_A, d_B, d_C, N);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
printf("Kernel execution time: %f ms\n", milliseconds);

对于N=1000的向量加法，内核执行时间通常在微秒级别，甚至低于计时精度。为了获得有意义的测量结果，建议将N增大到100万以上，并重复执行内核多次取平均值。

4.3.2 内存带宽与计算强度的初步分析

向量加法是典型的内存密集型（Memory-Bound）任务：每个线程执行一次浮点加法（2 FLOPs），但需要读取两个浮点数（8字节）并写入一个浮点数（4字节），总共12字节的内存流量。其算术强度仅为2 FLOPs / 12 Bytes ≈ 0.17 FLOPs/B，远低于现代GPU的平衡点（通常为10~20 FLOPs/B）。因此，向量加法的性能瓶颈在于内存带宽，而非计算能力。

以NVIDIA A100为例，其HBM2e显存带宽约为2TB/s。理论上，处理100万个元素的向量加法需要传输100万 × 12字节 = 12MB数据，耗时约12MB / 2TB/s = 6微秒。实际测量值通常会高于理论值，因为还存在内核启动开销、内存访问延迟等因素。