CUDA编程基础：向量加法、矩阵乘法、错误处理与性能分析工具

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

CUDA（Compute Unified Device Architecture）是NVIDIA推出的通用并行计算平台。其核心思想是将GPU视为一个高度并行的协处理器，通过异构计算模型，让CPU（Host）负责逻辑控制和串行任务，GPU（Device）负责大规模数据并行计算。理解CUDA编程，首先要建立三个基本概念：线程层次结构、内存层次结构和核函数（Kernel）。

线程层次从细到粗分为：Thread（线程） → Block（线程块） → Grid（网格）。一个Grid包含多个Block，每个Block包含多个Thread。这种层次结构直接映射到GPU的硬件架构：Thread映射到CUDA Core，Block映射到SM（Streaming Multiprocessor），Grid映射到整个GPU设备。内存层次则包括：每个Thread私有的寄存器和局部内存，每个Block共享的共享内存（Shared Memory），以及所有Thread可见的全局内存（Global Memory）、常量内存和纹理内存。

上图展示了CUDA编程模型的核心架构。核函数通过 kernel<<<Grid, Block>>>(args) 语法启动，其中Grid和Block可以是二维或三维的，便于处理图像等二维数据。每个Thread通过内置变量 threadIdx、blockIdx、blockDim 和 gridDim 来定位自己在整个计算空间中的唯一索引，从而决定处理哪部分数据。

方法流程：向量加法与矩阵乘法的CUDA实现

向量加法是最简单的CUDA入门案例。假设有两个长度为N的浮点数组A和B，计算C = A + B。CPU串行版本使用一个for循环，而CUDA版本则让每个Thread负责计算一个元素。核心代码如下：

// 核函数定义
__global__ void vectorAdd(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];
    }
}

// 主函数中的调用
int main() {
    int N = 1 << 20;  // 1,048,576 个元素
    size_t bytes = N * sizeof(float);

    // 1. 分配主机内存并初始化
    float *h_A = (float*)malloc(bytes);
    float *h_B = (float*)malloc(bytes);
    float *h_C = (float*)malloc(bytes);
    // ... 初始化 h_A 和 h_B ...

    // 2. 分配设备内存
    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A, bytes);
    cudaMalloc(&d_B, bytes);
    cudaMalloc(&d_C, bytes);

    // 3. 将数据从主机拷贝到设备
    cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice);

    // 4. 配置核函数启动参数
    int threadsPerBlock = 256;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    // 5. 启动核函数
    vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);

    // 6. 将结果拷贝回主机
    cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);

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

矩阵乘法则更具挑战性。对于 C = A × B（A为M×K，B为K×N），每个Thread计算C矩阵的一个元素。朴素实现中，每个Thread需要从全局内存读取A的一行和B的一列，计算K次乘加。这种实现存在大量重复的全局内存访问，性能较差。优化版本使用共享内存分块（Tiling）技术：将矩阵划分为多个子块，每个Block负责一个子块的计算，先将子块数据加载到共享内存，再从共享内存中读取数据进行计算，大幅减少全局内存访问次数。

// 使用共享内存优化的矩阵乘法核函数（分块大小为 TILE_SIZE）
#define TILE_SIZE 16

__global__ void matMulTiled(const float *A, const float *B, float *C,
                            int M, int K, int N) {
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];

    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    float sum = 0.0f;

    for (int tile = 0; tile < (K + TILE_SIZE - 1) / TILE_SIZE; ++tile) {
        // 协作加载A和B的子块到共享内存
        if (row < M && tile * TILE_SIZE + threadIdx.x < K)
            As[threadIdx.y][threadIdx.x] = A[row * K + tile * TILE_SIZE + threadIdx.x];
        else
            As[threadIdx.y][threadIdx.x] = 0.0f;

        if (col < N && tile * TILE_SIZE + threadIdx.y < K)
            Bs[threadIdx.y][threadIdx.x] = B[(tile * TILE_SIZE + threadIdx.y) * N + col];
        else
            Bs[threadIdx.y][threadIdx.x] = 0.0f;

        __syncthreads();  // 确保子块加载完成

        for (int k = 0; k < TILE_SIZE; ++k)
            sum += As[threadIdx.y][k] * Bs[k][threadIdx.x];

        __syncthreads();  // 确保所有线程都使用完共享内存后再加载下一块
    }

    if (row < M && col < N)
        C[row * N + col] = sum;
}

错误处理：CUDA编程的必备防线

CUDA API函数和核函数执行都可能失败，但核函数是异步执行的，错误不会立即返回。因此，必须建立严格的错误检查机制。常用的错误检查模式如下：

// 宏定义错误检查
#define CHECK_CUDA_ERROR(call) do { \
    cudaError_t err = (call); \
    if (err != cudaSuccess) { \
        fprintf(stderr, "CUDA Error at %s:%d - %s\n", \
                __FILE__, __LINE__, cudaGetErrorString(err)); \
        exit(EXIT_FAILURE); \
    } \
} while(0)

// 使用示例
CHECK_CUDA_ERROR(cudaMalloc(&d_A, bytes));
CHECK_CUDA_ERROR(cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice));

// 核函数启动后必须检查错误
vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);
CHECK_CUDA_ERROR(cudaGetLastError());  // 捕获核函数启动错误
CHECK_CUDA_ERROR(cudaDeviceSynchronize());  // 等待核函数执行完毕并捕获执行错误

常见错误陷阱：忘记检查 cudaDeviceSynchronize() 的返回值，导致核函数内部错误（如越界访问）被静默忽略；在循环中频繁调用 cudaMalloc 而不复用内存，造成性能下降和内存碎片；使用 cudaMemcpy 时方向参数写反（HostToDevice 与 DeviceToHost 混淆）。

性能分析工具：nvprof与Nsight Systems

性能分析是CUDA开发中不可或缺的环节。NVIDIA提供了两代性能分析工具：nvprof（命令行工具，已弃用但仍在教学和简单场景中使用）和Nsight Systems（新一代图形化工具）。

工具	定位	主要功能	适用场景
nvprof	命令行性能分析器	核函数耗时、内存传输耗时、占用率统计	快速定位瓶颈、脚本化分析
Nsight Systems	系统级可视化分析器	时间线视图、CPU/GPU并发分析、API调用追踪	整体性能调优、发现CPU-GPU同步问题
Nsight Compute	核函数级分析器	SM利用率、内存吞吐、指令级分析、瓶颈诊断	深入优化单个核函数

使用nvprof分析程序的基本命令如下：

# 分析所有CUDA活动
nvprof ./vector_add

# 只分析核函数耗时，按耗时排序
nvprof --print-gpu-trace ./vector_add

# 输出到文件
nvprof --log-file profile.txt ./vector_add

# Nsight Systems 命令行版本
nsys profile --stats=true -o my_profile ./vector_add

典型的nvprof输出会显示每个核函数的名称、网格/块配置、占用率、全局内存读写量、耗时等关键指标。通过分析这些数据，可以判断程序是计算密集型（核函数耗时长，占用率高）还是内存密集型（内存吞吐接近峰值，计算单元空闲），从而采取不同的优化策略。

案例复盘：向量加法的性能分析实战

假设我们运行上述向量加法程序，使用nvprof分析得到以下关键数据：

== Profiling result ==
Time(%)  Time(us)  Calls  Name
  78.2%   1250.3    1     vectorAdd
  21.8%    348.7    2     [CUDA memcpy HtoD]
  0.0%      1.2     1     [CUDA memcpy DtoH]

从这个结果可以看出：核函数耗时占比高达78%，而内存拷贝（HtoD）也占用了21.8%。对于向量加法这种简单计算，核函数耗时应该远小于内存传输时间。如果核函数耗时异常高，可能的原因包括：线程配置不合理（Block大小不是32的倍数，导致线程束浪费）、全局内存访问未合并（相邻线程访问不相邻的内存地址）、核函数启动开销过大（Grid和Block维度设置不当）。

优化建议：对于向量加法，Block大小通常设为128或256（32的倍数），Grid大小按公式 (N + blockSize - 1) / blockSize 计算。如果N很大，可以考虑使用网格跨步循环（Grid-Stride Loop），让每个线程处理多个元素，减少线程创建开销。此外，使用 cudaMallocManaged（统一内存）可以简化编程，但可能带来性能损失，需权衡使用。

在Nsight Systems的时间线视图中，可以直观看到CPU和GPU的活动重叠情况。理想状态下，CPU在等待GPU计算时应该执行其他任务（如准备下一批数据），而不是阻塞在 cudaMemcpy 或 cudaDeviceSynchronize 上。通过流（Stream）技术，可以实现数据传输与计算的重叠，进一步提升整体吞吐量。

最后，性能优化是一个迭代过程：先通过分析工具定位瓶颈，再针对性地修改代码，然后重新分析验证效果。切忌盲目优化，一定要用数据说话。掌握nvprof和Nsight工具的使用，是CUDA开发者从“能跑”走向“跑得快”的关键一步。