1.4 第一个实用 Kernel：向量加法

发表于 更新于 分类于 ， ，

通过一个完整的向量加法案例，走通 CUDA 编程的全流程：内存分配、数据传输、Kernel 编写、错误检查和性能测量，并对比 CPU 与 GPU 的执行效率。

📑 目录

1. 从 CPU 到 GPU：一个完整的开发流程

将计算任务从 CPU 搬到 GPU 就像把工作从一个能干的员工交给一支千人团队——你需要准备材料、分配任务、等待完成、收回成果。对应到 CUDA 编程，完整流程是：


graph LR
A["1. 分配 GPU 内存"] --> B["2. CPU→GPU 传数据"]
B --> C["3. 启动 Kernel"]
C --> D["4. GPU→CPU 取结果"]
D --> E["5. 释放 GPU 内存"]

每一步都有对应的 CUDA API：

步骤 API 说明
分配 GPU 内存 cudaMalloc() 类似 CPU 的 malloc()
Host→Device 传输 cudaMemcpy(..., HostToDevice) CPU 数据拷贝到 GPU
启动 Kernel kernel<<<grid, block>>>(args) GPU 并行执行
Device→Host 传输 cudaMemcpy(..., DeviceToHost) 结果拷回 CPU
释放 GPU 内存 cudaFree() 类似 CPU 的 free()

2. 向量加法 Kernel 实现

向量加法是最简单的并行问题——每个元素的计算完全独立，天然适合 GPU。给定两个长度为 N 的向量 A 和 B，计算：

$$
C[i] = A[i] + B[i], \quad i = 0, 1, \ldots, N-1
$$

2.1 Kernel 函数

1
2
3
4
5
6
7
8
9
__global__ void vectorAdd(const float* A, const float* B, float* C, int N) {
    // 计算全局线程索引
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    // 边界检查：防止越界访问
    if (idx < N) {
        C[idx] = A[idx] + B[idx];
    }
}

代码解析：

  • const float*：A 和 B 是只读输入
  • int idx：每个线程计算一个唯一的全局索引
  • if (idx < N)：最后一个 Block 可能有多余线程，必须检查边界
  • C[idx] = A[idx] + B[idx]：每个线程独立完成一个元素的加法

2.2 为什么需要边界检查

假设 N = 1000，Block 大小 = 256：

  • 需要 Block 数 = $\lceil 1000/256 \rceil = 4$
  • 总线程数 = 4 × 256 = 1024
  • 多出 24 个线程的索引为 1000~1023，超出数组范围

没有 if (idx < N) 检查，这些线程会访问非法内存地址，导致未定义行为或程序崩溃。

3. 主机端完整代码

3.1 内存管理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

int main() {
    int N = 1 << 20;  // 约 100 万个元素
    size_t size = N * sizeof(float);

    // 1. 分配主机内存
    float* h_A = (float*)malloc(size);
    float* h_B = (float*)malloc(size);
    float* h_C = (float*)malloc(size);

    // 初始化输入数据
    for (int i = 0; i < N; i++) {
        h_A[i] = (float)i;
        h_B[i] = (float)(i * 2);
    }

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

    // 3. 将输入数据从 Host 拷贝到 Device
    cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

    // 4. 配置并启动 Kernel
    int blockSize = 256;
    int gridSize = (N + blockSize - 1) / blockSize;
    vectorAdd<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);

    // 5. 将结果从 Device 拷回 Host
    cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

    // 6. 验证结果
    for (int i = 0; i < N; i++) {
        if (fabsf(h_C[i] - (h_A[i] + h_B[i])) > 1e-5f) {
            printf("Verification FAILED at index %d!\n", i);
            break;
        }
    }
    printf("Verification PASSED!\n");

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

    return 0;
}

3.2 编译与运行

1
2
3
$ nvcc vectorAdd.cu -o vectorAdd
$ ./vectorAdd
Verification PASSED!

4. CUDA 错误检查

CUDA API 调用可能失败（内存不足、参数错误等），生产代码中必须检查每个调用的返回值。

4.1 宏定义

1
2
3
4
5
6
7
8
9
#define CUDA_CHECK(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)

4.2 使用方式

1
2
3
4
5
6
7
CUDA_CHECK(cudaMalloc(&d_A, size));
CUDA_CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));

// Kernel 启动的错误检查需要两步
vectorAdd<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);
CUDA_CHECK(cudaGetLastError());          // 检查启动参数是否正确
CUDA_CHECK(cudaDeviceSynchronize());     // 检查执行过程是否有错误

4.3 常见错误码

错误码 含义 常见原因
cudaErrorMemoryAllocation 内存分配失败 显存不足
cudaErrorInvalidValue 参数无效 size=0 或空指针
cudaErrorInvalidConfiguration 启动配置无效 Block 超过 1024 线程
cudaErrorIllegalAddress 非法内存访问 越界/野指针
cudaErrorLaunchTimeout Kernel 执行超时 死循环或计算量过大

💡 提示：开发阶段建议在每个 CUDA API 后都加上错误检查。发布时可以用编译宏控制，在 Release 模式下去掉检查开销。

5. 性能测量：CPU vs GPU

5.1 CPU 基线实现

1
2
3
4
5
void vectorAddCPU(const float* A, const float* B, float* C, int N) {
    for (int i = 0; i < N; i++) {
        C[i] = A[i] + B[i];
    }
}

5.2 使用 CUDA Event 计时

CUDA Event 是 GPU 端精确计时的标准方法，精度达微秒级：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);

// 计时开始
cudaEventRecord(start);

// 要计时的操作
vectorAdd<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);

// 计时结束
cudaEventRecord(stop);
cudaEventSynchronize(stop);

// 计算耗时
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
printf("GPU Kernel time: %.4f ms\n", milliseconds);

// 清理
cudaEventDestroy(start);
cudaEventDestroy(stop);

5.3 包含数据传输的端到端计时

1
2
3
4
5
6
7
8
9
10
11
12
13
cudaEventRecord(start);

cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
vectorAdd<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);
cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

cudaEventRecord(stop);
cudaEventSynchronize(stop);

float totalTime = 0;
cudaEventElapsedTime(&totalTime, start, stop);
printf("GPU total (with memcpy): %.4f ms\n", totalTime);

5.4 性能结果对比（参考值）

在 RTX 4090 上，N = $2^{24}$（约 1600 万元素）的典型结果：

📊 指标 CPU (单核) GPU (仅 Kernel) GPU (含传输)
耗时 ~25 ms ~0.3 ms ~8 ms
有效带宽 ~5 GB/s ~400 GB/s ~15 GB/s

⚠️ 注意：向量加法是内存密集型操作（每个元素只做 1 次加法但需要读 2 个 float、写 1 个 float）。包含数据传输时，PCIe 带宽成为瓶颈。GPU 的优势在数据已经在显存上时最明显。

6. 完整可运行代码

将以上所有部分整合为一个完整的、可直接编译运行的程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <chrono>
#include <cuda_runtime.h>

// 错误检查宏
#define CUDA_CHECK(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)

// GPU Kernel
__global__ void vectorAdd(const float* A, const float* B, float* C, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < N) {
        C[idx] = A[idx] + B[idx];
    }
}

// CPU 基线
void vectorAddCPU(const float* A, const float* B, float* C, int N) {
    for (int i = 0; i < N; i++) {
        C[i] = A[i] + B[i];
    }
}

int main() {
    const int N = 1 << 24;  // 约 1600 万元素
    const size_t size = N * sizeof(float);
    printf("Vector size: %d elements (%.1f MB)\n", N, size / (1024.0 * 1024.0));

    // ========== 分配主机内存 ==========
    float* h_A = (float*)malloc(size);
    float* h_B = (float*)malloc(size);
    float* h_C_cpu = (float*)malloc(size);
    float* h_C_gpu = (float*)malloc(size);

    // 初始化数据
    for (int i = 0; i < N; i++) {
        h_A[i] = sinf(i * 0.001f);
        h_B[i] = cosf(i * 0.001f);
    }

    // ========== CPU 计算 ==========
    auto cpu_start = std::chrono::high_resolution_clock::now();
    vectorAddCPU(h_A, h_B, h_C_cpu, N);
    auto cpu_end = std::chrono::high_resolution_clock::now();
    float cpu_ms = std::chrono::duration<float, std::milli>(cpu_end - cpu_start).count();
    printf("CPU time: %.4f ms\n", cpu_ms);

    // ========== GPU 计算 ==========
    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc(&d_A, size));
    CUDA_CHECK(cudaMalloc(&d_B, size));
    CUDA_CHECK(cudaMalloc(&d_C, size));

    // 创建 CUDA Event 计时器
    cudaEvent_t start, stop, kernel_start, kernel_stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventCreate(&kernel_start);
    cudaEventCreate(&kernel_stop);

    // 端到端计时（含传输）
    cudaEventRecord(start);

    CUDA_CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice));

    int blockSize = 256;
    int gridSize = (N + blockSize - 1) / blockSize;

    // 仅 Kernel 计时
    cudaEventRecord(kernel_start);
    vectorAdd<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);
    cudaEventRecord(kernel_stop);

    CUDA_CHECK(cudaGetLastError());
    CUDA_CHECK(cudaMemcpy(h_C_gpu, d_C, size, cudaMemcpyDeviceToHost));

    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    float gpu_total_ms = 0, kernel_ms = 0;
    cudaEventElapsedTime(&gpu_total_ms, start, stop);
    cudaEventElapsedTime(&kernel_ms, kernel_start, kernel_stop);

    printf("GPU kernel time: %.4f ms\n", kernel_ms);
    printf("GPU total time (with memcpy): %.4f ms\n", gpu_total_ms);
    printf("Speedup (kernel only): %.1fx\n", cpu_ms / kernel_ms);
    printf("Speedup (end-to-end): %.1fx\n", cpu_ms / gpu_total_ms);

    // ========== 结果验证 ==========
    int errors = 0;
    for (int i = 0; i < N; i++) {
        if (fabsf(h_C_gpu[i] - h_C_cpu[i]) > 1e-5f) {
            if (errors < 5) {
                printf("Mismatch at %d: CPU=%.6f, GPU=%.6f\n",
                       i, h_C_cpu[i], h_C_gpu[i]);
            }
            errors++;
        }
    }
    if (errors == 0) {
        printf("Verification PASSED!\n");
    } else {
        printf("Verification FAILED: %d mismatches\n", errors);
    }

    // ========== 计算有效带宽 ==========
    // 向量加法读 2 个 float、写 1 个 float = 12 bytes/element
    float bandwidth_gb = (3.0f * size) / (kernel_ms * 1e-3f) / 1e9f;
    printf("Effective bandwidth: %.1f GB/s\n", bandwidth_gb);

    // ========== 清理 ==========
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    cudaEventDestroy(kernel_start);
    cudaEventDestroy(kernel_stop);
    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));
    free(h_A);
    free(h_B);
    free(h_C_cpu);
    free(h_C_gpu);

    return 0;
}
1
2
3
4
5
# 编译
nvcc -O3 -arch=sm_89 vectorAdd.cu -o vectorAdd

# 运行
./vectorAdd

7. 性能分析与讨论

7.1 向量加法的算术强度

算术强度（Arithmetic Intensity）= 计算量 / 内存访问量：

$$
\text{AI} = \frac{\text{FLOPs}}{\text{Bytes}} = \frac{1 \text{ FLOP}}{12 \text{ Bytes}} \approx 0.083 \text{ FLOP/Byte}
$$

这个值极低，说明向量加法是典型的内存密集型操作，性能受限于内存带宽而非计算能力。

7.2 带宽利用率分析

如果 Kernel 时间为 0.3 ms，处理 $3 \times 2^{24} \times 4$ Bytes = 201 MB 数据：

$$
\text{有效带宽} = \frac{201 \text{ MB}}{0.3 \text{ ms}} \approx 670 \text{ GB/s}
$$

对比 RTX 4090 的理论峰值带宽 1008 GB/s，利用率约 66%——对于如此简单的 Kernel 来说已经不错。

7.3 PCIe 传输的开销

📊 环节 典型耗时 占比
H2D 传输 (64MB × 2) ~5 ms ~60%
Kernel 执行 ~0.3 ms ~4%
D2H 传输 (64MB) ~3 ms ~36%

📌 关键点：对于简单的逐元素操作，数据传输往往是主要瓶颈。实际应用中应尽量让数据”留在 GPU 上”——多个 Kernel 串联处理同一批数据，避免反复搬运。

7.4 何时值得用 GPU

场景 CPU 更好 GPU 更好
数据量小（< 10K 元素） ✅ ❌
单次简单操作 + 数据不在 GPU ✅ ❌
大数据 + 数据已在 GPU ❌ ✅
大数据 + 多步串联计算 ❌ ✅
计算密集型（如矩阵乘法） ❌ ✅

8. 进阶练习

完成向量加法后，尝试以下练习来巩固所学：

8.1 练习 1：向量点积

实现两个向量的点积 $\text{result} = \sum_{i=0}^{N-1} A[i] \times B[i]$。

提示：这需要规约（Reduction）——每个线程计算部分乘积，然后在 Block 内用共享内存累加。

8.2 练习 2：SAXPY

实现 SAXPY 操作：$Y[i] = \alpha \times X[i] + Y[i]$

这是 BLAS 的基础操作，比向量加法多一个标量乘法。

8.3 练习 3：Grid-stride Loop 版本

将向量加法改为 Grid-stride loop 写法，启动固定数量的线程处理任意大小的数组：

1
2
3
4
5
6
7
8
9
10
11
12
__global__ void vectorAddStride(const float* A, const float* B,
                                float* C, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    for (int i = idx; i < N; i += stride) {
        C[i] = A[i] + B[i];
    }
}

// 启动时不再需要 gridSize = (N + 255) / 256
// 可以选择任意合理的 Grid 大小
vectorAddStride<<<128, 256>>>(d_A, d_B, d_C, N);

8.4 练习 4：异步传输与计算重叠

使用 CUDA Stream 让数据传输和计算并行执行：

1
2
3
4
5
6
7
8
9
10
11
cudaStream_t stream;
cudaStreamCreate(&stream);

// 异步传输和计算可以重叠
cudaMemcpyAsync(d_A, h_A, size, cudaMemcpyHostToDevice, stream);
cudaMemcpyAsync(d_B, h_B, size, cudaMemcpyHostToDevice, stream);
vectorAdd<<<gridSize, blockSize, 0, stream>>>(d_A, d_B, d_C, N);
cudaMemcpyAsync(h_C, d_C, size, cudaMemcpyDeviceToHost, stream);

cudaStreamSynchronize(stream);
cudaStreamDestroy(stream);

💡 提示：要实现真正的重叠，主机端内存必须使用 cudaMallocHost()（Pinned Memory）分配。

📝 总结

通过向量加法这个完整案例，我们掌握了 CUDA 编程的核心流程：

  1. 五步标准流程：cudaMalloc → cudaMemcpy(H2D) → Kernel Launch → cudaMemcpy(D2H) → cudaFree
  2. Kernel 编写三要素：全局索引计算 + 边界检查 + 并行计算
  3. 错误检查：用 CUDA_CHECK 宏包裹每个 API 调用
  4. 性能度量：CUDA Event 计时 + 有效带宽计算
  5. 性能认知：简单操作的瓶颈在内存传输，GPU 优势在数据驻留和计算密集场景

向量加法虽然简单，但它建立起了 CUDA 编程的完整框架。后续学习矩阵乘法、卷积、规约等高级 Kernel 时，都是在这个框架上增加复杂度。

🎯 自我检验清单

  • 能独立写出 CUDA 向量加法的完整代码（含内存管理和错误检查）
  • 能解释为什么 Kernel 中需要 if (idx < N) 边界检查
  • 能使用 CUDA Event 精确测量 Kernel 执行时间
  • 能计算给定 Kernel 的有效内存带宽并与理论峰值对比
  • 能判断一个操作是计算密集型还是内存密集型
  • 能解释为什么包含数据传输的端到端时间远大于纯 Kernel 时间
  • 能使用 Grid-stride loop 模式处理任意大小的数据
  • 能使用 CUDA_CHECK 宏编写健壮的 CUDA 代码

📚 参考资料