cuda编程入门（1）

简介

这是这个系列的第一篇入门文章，这个系列的博客不会讲解太多的东西，毕竟官方文档就是最好的教程，这个系列的博客主要是写一些 cuda 代码跑起来试试，记录下自己的学习收获。
官方文档就是最好最权威的学习资源
​​CUDA Toolkit Documentation v11.7.0​​

此外作为入门，个人觉得这篇文章是很不错的 ，强烈建议先看完下面这篇文章。
​​CUDA编程入门极简教程​​ 它主要是翻译自这篇文章
​​even-easier-introduction-cuda​​

GPU简介

简而言之，GPU加速就是利用多线程加速。相比CPU，GPU有更多的核心，能够同时跑更多的线程。同时线程之间的切换代价也非常小。

可以看出，GPU的core数远大于CPU，因此，我们把一个计算任务，分为很多小份，多线程同时跑在许多核心上，就可以实现计算加速。

CUDA编程模型简介

典型的CUDA程序的执行流程如下：

1. 分配host内存，并进行数据初始化；
2. 分配device内存，并从host将数据拷贝到device上；
3. 调用CUDA的核函数在device上完成指定的运算；
4. 将device上的运算结果拷贝到host上；
5. 释放device和host上分配的内存。

上面流程中最重要的一个过程是调用CUDA的核函数来执行并行计算，kernel是CUDA中一个重要的概念，kernel是在device上线程中并行执行的函数，核函数用__global__符号声明，在调用时需要用<<<grid, block>>>来指定kernel要执行的线程数量，在CUDA中，每一个线程都要执行核函数，并且每个线程会分配一个唯一的线程号thread ID，这个ID值可以通过核函数的内置变量threadIdx来获得。

要深刻理解kernel，必须要对kernel的线程层次结构有一个清晰的认识。首先GPU上很多并行化的轻量级线程。kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格（grid），同一个网格上的线程共享相同的全局内存空间，grid是线程结构的第一层次，而网格又可以分为很多线程块（block），一个线程块里面包含很多线程，这是第二个层次。线程两层组织结构如下图所示，这是一个gird和block均为2-dim的线程组织。grid和block都是定义为dim3类型的变量，dim3可以看成是包含三个无符号整数（x，y，z）成员的结构体变量，在定义时，缺省值初始化为1。因此grid和block可以灵活地定义为1-dim，2-dim以及3-dim结构，对于图中结构（主要水平方向为x轴），定义的grid和block如下所示，kernel在调用时也必须通过​​执行配置​​<<<grid, block>>>来指定kernel所使用的线程数及结构。

dim3 grid(3, 2);
dim3 block(5, 3);
kernel_fun<<< grid, block >>>(prams…);

更详细的内容请直接参考博客​​CUDA编程入门极简教程​​

向量加法实战

完整代码如下：

// 两个向量加法kernel，grid和block均为一维
__global__ void add(float* x, float * y, float* z, int n)
{
    // 获取该线程的全局索引
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    // 步长(线程总数)
    int stride = blockDim.x * gridDim.x;
    for (int i = index; i < n; i += stride)
    {
        z[i] = x[i] + y[i];
    }
}

#include <iostream>
int main()
{
    int N = 1 << 20; //(1048576)
    int nBytes = N * sizeof(float);
    // 申请host内存
    float *x, *y, *z;
    x = (float*)malloc(nBytes);
    y = (float*)malloc(nBytes);
    z = (float*)malloc(nBytes);

    // 初始化数据
    for (int i = 0; i < N; ++i)
    {
        x[i] = 10.0;
        y[i] = 20.0;
    }

    // 申请device内存
    float *d_x, *d_y, *d_z;
    cudaMalloc((void**)&d_x, nBytes);
    cudaMalloc((void**)&d_y, nBytes);
    cudaMalloc((void**)&d_z, nBytes);

    // 将host数据拷贝到device
    cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);
    cudaMemcpy((void*)d_y, (void*)y, nBytes, cudaMemcpyHostToDevice);
    // 定义kernel的执行配置
    dim3 blockSize(256);
    dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
    // 执行kernel
    add << < gridSize, blockSize >> >(d_x, d_y, d_z, N);

    // 将device得到的结果拷贝到host
    cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyDeviceToHost);

    // 检查执行结果
    float maxError = 0.0;
    for (int i = 0; i < N; i++)
        maxError = fmax(maxError, fabs(z[i] - 30.0));
    std::cout << "最大误差: " << maxError << std::endl;

    // 释放device内存
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);
    // 释放host内存
    free(x);
    free(y);
    free(z);

    return 0;
}

上述代码中 blockSize = 256； gridSize = 4906

运行分析：

既然是多线程加速，那线程是不是越多越好呢。很明显并不是，线程太多了，线程之间切换的开销也会变大。不如测试一下：

上面的测试设置的线程数为：

blockSize = 256； gridSize = 4906 耗时57.280us

改为 blockSize = 1024； gridSize = 4906 试试

可见耗时增加了，因为线程太多了，线程之间的切换也增加了耗时。

但最佳线程数是多少呢？

从官方文档中可知，每个 block 的线程数一般限制在 1024

There is a limit to the number of threads per block, since all threads of a block are expected to reside on the same processor core and must share the limited memory resources of that core. On current GPUs, a thread block may contain up to 1024 threads.

同时 block 在SM（string multiprocessors : 流式多处理器）之间调度

Indeed, each block of threads can be scheduled on any of the available multiprocessors within a GPU, in any order, concurrently or sequentially, so that a compiled CUDA program can execute on any number of multiprocessors as illustrated by Figure 3, and only the runtime system needs to know the physical multiprocessor count.

所以 GPU的核心数 为 SM * 1024。 当线程数 <=核心数时，线程之间一般就不存在切换，减小了切换开销，同时最大程度上利用了GPU的资源。

可以查询GPU硬件的一些信息，代码如下：

int dev = 0;
    cudaDeviceProp devProp;
    cudaGetDeviceProperties(&devProp, dev);
    std::cout << "使用GPU device " << dev << ": " << devProp.name << std::endl;
    std::cout << "SM的数量：" << devProp.multiProcessorCount << std::endl;
    std::cout << "每个线程块的共享内存大小：" << devProp.sharedMemPerBlock / 1024.0 << " KB" << std::endl;
    std::cout << "每个线程块的最大线程数：" << devProp.maxThreadsPerBlock << std::endl;
    std::cout << "每个SM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << std::endl;

重新设置线程的层次结构：

// 定义kernel的执行配置
dim3 blockSize(devProp.maxThreadsPerMultiProcessor);
dim3 gridSize(devProp.multiProcessorCount);

运入如下：

有趣的是：

如果线程的层次结构设置如下，性能会出现大幅下降：

// 定义kernel的执行配置
dim3 blockSize(devProp.maxThreadsPerMultiProcessor);
dim3 gridSize(devProp.multiProcessorCount + 1);

设置为devProp.multiProcessorCount的倍数，性能又会恢复：

// 定义kernel的执行配置
dim3 blockSize(devProp.maxThreadsPerMultiProcessor);
dim3 gridSize(devProp.multiProcessorCount * 2);