CUDA编程入门——一个将CPU上函数转移到GPU上计算的Case_cpu程序到gpu程序的迁移

这篇文章将简单地介绍NVIDIA公司推出的并行计算平台和编程模型 CUDA。 CUDA C++ 是使用 CUDA
创建大规模并行应用程序的方法之一。它使得我们可以使用C++编程语言来开发在GPU
上运行的数千个并行线程加速的高性能算法。下面我们将以两个float数组（数组大小为2^20）相加的场景为例，说明一段原本运行在CPU上的C++程序如何迁移到GPU上执行。

前置工作

以下相关代码的学习需要一台带有支持 CUDA 的 GPU 的计算机，或者一个带有 GPU 的云实例（腾讯云、阿里云、华为云、AutoDL等），推荐操作系统为Linux。

Simple Case on CPU

如下是一段在CPU上运行两数组相加的完整代码。

#include<iostream>
#include <math.h>

// 执行两个float数组相加的函数
void add(int n , float*x, float*y){
    for(int i = 0; i < n; ++i){
        y[i] = x[i] + y[i];
    }
}

int main(void){
    int N = 1 << 20; // 1M个元素

    // 申请内存：在CPU Memory（Host）上申请数组占用的内存空间
    float *x = new float[N];
    float *y = new float[N];

    // 内存初始化：在Host上初始化两个数组 x, y
    for(int i = 0; i < N; ++i){
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    // 执行核函数：在CPU上运行两个包含有1M个元素的数组相加的函数add
    add(N, x, y);

    //检验运算的正确与否：经过运行之后，两者加和的结果应该存储在y处且元素值均为3.0f
    float maxError = 0.0f;
    for(int i = 0; i < N; ++i){
        maxError = fmax(maxError, fabs(y[i]-3.0f));
    }
    std::cout << "Max error: " << maxError << std::endl;

    // 释放内存：释放在Host上动态申请的内存空间
    delete[] x;
    delete[] y;
    return 0;
}
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

如何编译并运行该段代码？

1. 将该段代码拷贝到文件中并保存为add.cpp；
2. 编译代码：g++ add.cpp -o add，如果是在Mac机器，编译器需修改为：clang++
3. 运行代码并得到结果：

$ ./add 
Max error: 0
1
2

Running on GPU

现在，当我们想在 GPU 运行add函数时，我们需要对代码做以下几个修改。

1. 为add函数添加指定标识符__global__

__global__ // 表明这是一个能从CPU代码中调用并运行在GPU上的函数
void add(int n , float*x, float*y){
    for(int i = 0; i < n; ++i){
        y[i] = x[i] + y[i];
    }
}
1
2
3
4
5
6

通过为函数添加**__global__**标识符，将会告诉CUDA C++编译器，这是一个可以调用CPU处代码并在GPU上运行的函数，这个函数在GPU编程中被称为核函数(kernel)。同时为了将GPU上运行的代码和内存与CPU区分开来，前者通常称为device code，后者通常称为host code。

2. GPU上的内存分配

在GPU上运行核函数进行计算，就要求对于核函数的输入所在内存空间能被CUDA程序所访问，因此对于add函数的输入x,y，我们需要为其在GPU上申请一段可访问的内存空间。
而CUDA 中的unified memory提供了一个系统中所有 GPU 和 CPU 都能访问的单一内存空间。要在unified memory中分配空间，需要调用 cudaMallocManaged()，它将返回一个指针，而后可以通过host（CPU）代码或device（GPU）代码访问该指针。要释放数据，只需将指针传递给cudaFree()。

// 申请内存：在CPU Memory（Host）上申请数组占用的内存空间
float *x = new float[N];
float *y = new float[N];

// ======================================> 替换

// 申请device memory：从Unified Memory分配数组需要占用的内存空间
// Unified Memory对于CPU和GPU都是可访问的。
float *x, *y;
cudaMallocManaged(&x, N*sizeof(float));
cudaMallocManaged(&y, N*sizeof(float));
1
2
3
4
5
6
7
8
9
10
11

实际操作：将上面代码中对 new 的调用替换为对 cudaMallocManaged() 的调用，并将对 delete [] 的调用替换为对cudaFree()的调用。

3. CUDA核函数的调用

根据CUDA语法，在GPU上调用核函数，需要采用三角括号语法<<<>>>来指定。调用时我们需要将其加入到函数的参数列表之前。

// Run kernel on 1M elements on CPU
add<<<1, 1>>>(N, x, y);
1
2

启动 add() 内核，在 GPU 上调用它。CUDA 内核启动使用三角括号语法 <<< >>>指定。我只需将其添加到参数列表之前的 add 调用中即可。

4. 同步Host(CPU)与device(GPU)之间的操作

cudaDeviceSynchronize() 是 CUDA 编程中的一个函数，它用于同步主机（CPU）和设备（GPU）之间的操作。这个函数会等待设备上的所有流（stream）中的任务完成，然后返回。
在CUDA编程中，GPU可以并行执行多个操作，这些操作通常是异步的，同时CUDA内核启动时也不会阻塞调用CPU的线程。这意味着当在Host（CPU）上执行代码时，GPU可能仍在执行之前的某些任务。cudaDeviceSynchronize() 的作用是等待GPU上的所有任务完成，确保主机上的代码不会继续执行，直到GPU上的所有任务都完成。
因此，在如上代码中，为了保证最后结果检验代码的正确性，我们需要在 CPU 上进行最后的错误检查之前调用cudaDeviceSynchronize()以保证在GPU上运算结果全部完成之后再进行错误检查。

// 等待 GPU 运算完成后再访问主机
cudaDeviceSynchronize();
1
2

以下为Host code迁移到GPU上后的完整代码：

#include<iostream>
#include <math.h>

__global__ // 表明这是一个能从CPU代码中调用并运行在GPU上的函数
void add(int n , float*x, float*y){
    for(int i = 0; i < n; ++i){
        y[i] = x[i] + y[i];
    }
}

int main(void){
    int N = 1 << 20;

    // Allocate Unified Memory —— accessible from CPU or GPU
    float *x, *y;
    cudaMallocManaged(&x, N*sizeof(float));
    cudaMallocManaged(&y, N*sizeof(float));

    // initialize x, y
    for(int i = 0; i < N; ++i){
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    // Run kernel on 1M elements on CPU
    add<<<1, 1>>>(N, x, y);

    // Wait for GPU to finish before accessing on host
    cudaDeviceSynchronize();

    // Check for errors(all values should be 3.0f)
    float maxError = 0.0f;
    for(int i = 0; i < N; ++i){
        maxError = fmax(maxError, fabs(y[i]-3.0f));
    }
    std::cout << "Max error: " << maxError << std::endl;

    cudaFree(x);
    cudaFree(y);
    return 0;
}
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

如何编译并运行该段代码？

1. 将该段代码拷贝到文件中并保存为add.cu；【CUDA文件通常以.cu为后缀；】
2. 编译代码：nvcc add.cu -o add_cuda；【CUDA文件编译时通常采用nvcc编译器】
   • 可能出现的问题：如果执行该指令时，出现sh: 1: cicc: not found的报错，需要设置环境变量export PATH=/usr/local/cuda/bin:$PATH。细节请参考https://github.com/ShuangLI59/person_search/issues/56。
3. 运行代码并得到结果：

$ ./add_cuda
Max error: 0
1
2

以上是由CPU迁移到GPU的过程，但当前的代码仍然只是单线程，即同时只有一个线程会执行该核函数，对数组进行操作。倘若更改为多线程，那么由于多个并行线程都会读写相同的位置，因此存在Race Condition。
【Race Condition: 指在多线程或并发环境中，程序的执行结果依赖于线程执行的相对时间或执行顺序的不确定性，从而导致程序在不同的执行情境下产生不同的结果。这种情况通常发生在多个线程或进程尝试同时访问和修改共享资源的情况下。（这里对Race Condition不作进一步解释，有兴趣者可自行检索）】

多线程并行

如何多线程并行运行核函数呢？关键在于 CUDA 的 <<<1, 1>>>语法，此处称之为核函数的执行配置（execution configuration），通过该配置我们可以设定在GPU上运行核函数时要启动多少个并行线程。这里有两个参数，此处跳过第一个，先关注并修改第二个参数：线程块中的线程数（the number of threads in a thread block）。对于CUDA GPU 在设定thread block的threads number时，需要设定为32的倍数，这里我们选择设定为256 个线程。

// Run kernel on 1M elements on CPU
add<<<1, 1>>>(N, x, y);

// ======================> Modify
add<<<1, 256>>>(N, x, y);
1
2
3
4
5

如果仅仅只做这样的改动，那么每个线程均会对x,y所在的整个内存空间进行一次计算，而不是将不同内存地址的数据的计算分散到并行线程中，由不同线程执行同一个核函数但接受不同内存地址的输入。为了正确执行，我们需要修改核函数，使得不同的线程和不同的内存空间的输入对应起来。
CUDA C++ 提供的关键字threadIdx使得我们可以在核函数中获取当前运行该核函数的线程的索引Id，有了索引Id后，我们就可以让不同id对应的线程和不同内存地址的数据对应起来并执行计算。具体来说，threadIdx.x 即为当前线程的索引，而 blockDim.x 则为一个threadBlock中所设定的线程数，即为256。下面我们通过修改核函数中的循环，将256个线程和它们对应的内存数据分别对应起来。

__global__ // 表明这是一个能从CPU代码中调用并运行在GPU上的函数
void add(int n , float*x, float*y){
    int index = threadIdx.x; // 当前线程索引
    int stride = blockDim.x; // 一个block中线程的数量
    for(int i = index; i < n; i+=stride){
        y[i] = x[i] + y[i];
    }
}
1
2
3
4
5
6
7
8

根据上述代码，线程id即index将在[0, 256)范围内，那么对于0号线程，所要执行的计算应为：

y[0] = x[0] + y[0]
y[256] = x[256] + y[256]
...
y[0+n*256] = x[0+n*256] + y[0+n*256]
...
0 =< n < 2^12
1
2
3
4
5
6

如此一来，每个线程都将执行2^12次加法计算。事实上，将 index 设置为 0、stride 设置为 1 后，它在语义上与第一个版本add.cu完全相同。以下是完整代码：

#include<iostream>
#include <math.h>
#include <chrono>

__global__ // 表明这是一个能从CPU代码中调用并运行在GPU上的函数
void add(int n , float*x, float*y){
    int index = threadIdx.x; // 当前线程索引
    int stride = blockDim.x; // 一个block中线程的数量
    for(int i = index; i < n; i+=stride){
        y[i] = x[i] + y[i];
    }
}

int main(void){
    int N = 1 << 20;

    // Allocate Unified Memory —— accessible from CPU or GPU
    float *x, *y;
    cudaMallocManaged(&x, N*sizeof(float));
    cudaMallocManaged(&y, N*sizeof(float));

    // initialize x, y
    for(int i = 0; i < N; ++i){
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    add<<<1, 256>>>(N, x, y);
    cudaDeviceSynchronize();

    // Check for errors(all values should be 3.0f)
    float maxError = 0.0f;
    for(int i = 0; i < N; ++i){
        maxError = fmax(maxError, fabs(y[i]-3.0f));
    }
    std::cout << "Max error: " << maxError << std::endl;

    cudaFree(x);
    cudaFree(y);
    return 0;
}
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

如何编译并运行该段代码？

1. 将该段代码拷贝到文件中并保存为add_block.cu；
2. 编译代码：nvcc add_block.cu -o add_block；
3. 运行代码并得到结果：

$ ./add_block
Max error: 0
1
2

多Block

CUDA GPU 有许多并行处理器，它们被分组组成Streaming Multiprocessors(SMs, 流式多处理器)。每个 SM 可以运行多个并发线程块。例如，基于Pascal GPU 架构的 Tesla P100 GPU 有 56 个 SM，每个SM有64个CUDA Core，每个 SM 最多可支持 2048 个活动线程并发；基于Volta GPU架构的Tesla V100 GPU有80个SM，每个SM有64个CUDA Core，每个SM最多可支持2048个活动线程并发。为了充分利用所有这些线程，此处我们应该用多个线程块启动核函数。
由上文可推测知，Execution Configuration的第二个参数指定了每个线程块中线程的数量，而第一个参数指定了线程块的数量。CUDA中多个线程组成一个线程块thread block，多个线程块共同组成所谓的网格grid。
对于数组相加这个场景，由于有 N 个元素要处理，而每个区块有 256 个线程，我们希望启动N个线程，每个线程负责一次加法处理，因此对于要设定的线程块的数量，只需用 N 除以线程块大小（每个线程块中包含的线程数）即可获得（注意四舍五入，以防 N 不是 blockSize 的倍数）。

// numBlocks: the number of thread blocks
// blockSize: the number of threads in a thread block
int blockSize = 256;
int numBlocks = (N + blockSize - 1) / blockSize; // 4096
add<<<numBlocks, blockSize>>>(N, x, y);
1
2
3
4
5

同样地，在更改了启动线程数后，我们仍然需要核函数以保证线程id和需要负责计算内存数据之间对应起来。CUDA 提供了 gridDim.x（网格中的线程块数）和 blockIdx.x（包含网格中当前线程块的索引）。下图展示了在 CUDA 中使用 blockDim.x、gridDim.x 和 threadIdx.x 对数组（一维）进行索引的方法。其原理是，每个线程通过计算其区块开头的偏移量（区块索引乘以区块大小：blockIdx.x * blockDim.x）并加上线程在区块中的索引（threadIdx.x）来获得索引。代码 blockIdx.x * blockDim.x + threadIdx.x 是惯用的 CUDA 代码。

__global__ // 表明这是一个能从CPU代码中调用并运行在GPU上的函数
void
add(int n, float *x, float *y)
{
    int index = blockIdx.x * blockDim.x + +threadIdx.x; // 当前线程索引
    int stride = blockDim.x * gridDim.x;                // 一个grid即所有的block中线程的总数
    for (int i = index; i < n; i += stride)
        {
            y[i] = x[i] + y[i];
        }
}
1
2
3
4
5
6
7
8
9
10
11

更新后的内核还会将 stride 设置为网格中的线程总数（blockDim.x * gridDim.x）。CUDA 内核中的这种循环通常称为grid-stride循环。
完整代码如下：

#include <iostream>
#include <math.h>
#include <chrono>

__global__ // 表明这是一个能从CPU代码中调用并运行在GPU上的函数
void
add(int n, float *x, float *y)
{
    int index = blockIdx.x * blockDim.x + +threadIdx.x; // 当前线程索引
    int stride = blockDim.x * gridDim.x;                // 偏移量为总的线程数
    for (int i = index; i < n; i += stride)
        {
            y[i] = x[i] + y[i];
        }
}

int main(void)
{
    int N = 1 << 20;

    // Allocate Unified Memory —— accessible from CPU or GPU
    float *x, *y;
    cudaMallocManaged(&x, N * sizeof(float));
    cudaMallocManaged(&y, N * sizeof(float));

    // initialize x, y
    for (int i = 0; i < N; ++i)
        {
            x[i] = 1.0f;
            y[i] = 2.0f;
        }

    add<<<numBlocks, blockSize>>>(N, x, y);
    // Wait for GPU to finish before accessing on host
    cudaDeviceSynchronize();

    // Check for errors(all values should be 3.0f)
    float maxError = 0.0f;
    for (int i = 0; i < N; ++i)
        {
            maxError = fmax(maxError, fabs(y[i] - 3.0f));
        }
    std::cout << "Max error: " << maxError << std::endl;

    cudaFree(x);
    cudaFree(y);
    return 0;
}
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

如何编译并运行该段代码？

1. 将该段代码拷贝到文件中并保存为add_grid.cu；
2. 编译代码：nvcc add_grid.cu -o add_grid；
3. 运行代码并得到结果：

$ ./add_grid
Max error: 0
1
2

速度对比

这一部分我们比较一下add.cu， add_block.cu， add_grid.cu中核函数的运行时间，为了使得对比更加明显，我们运行1024次然后比较总的时间，此处需要用循环包住核函数以及同步操作：

for(int i = 0; i < 1024; ++i){
    // add.cu
    add<<<1, 1>>>(N, x, y);
    // // add_block
    // add<<<1, 256>>>(N, x, y);
    // // add_grid
    // add<<<1, 256>>>(N, x, y);
    // Wait for GPU to finish before accessing on host
    cudaDeviceSynchronize();
}
1
2
3
4
5
6
7
8
9
10

CUDA提供了nvprof工具帮助我们快速获得核函数运行所耗费的时间，只需在命令行中输入 nvprof ./add_cuda即可。单线程、单block多线程以及多block多线程在Tesla V100上分别运行的时间统计如下：

参考链接

An Even Easier Introduction to CUDA | NVIDIA Technical Blog