当今科学计算越来越复杂，各种模型越来越大，使用SVM做分类的时代早已一去不复返，通过GPU加速应用的重要性不言而喻。我在MRI领域工作的这几年，经常碰高度复杂的MRI应用，Matlab计算时间在几小时到几天不等，即使高度优化的C++的程序也需要10至30分钟左右。这些应用的CPU利用率已然接近100%，因此通过GPU进一步优化程序成为最可行的手段。

GPU市场依然是NVIDIA一家独大，AMD想要与之抗衡还需要时间。尤其在生态方面，NVIDIA基于CUDA的护城河太深，以至于科学计算大部分都使用N家的GPU。如果你已经拥有N家的GPU,恐怕使用CUDA编写核函数加速应用仍然是最好的选择。

在开始之前，介绍下我的基本配置和驱动情况：

  • CPU: AMD Ryzen5 1600
  • RAM: 8G x 4
  • GPU: NVIDIA GTX 1080 x 2 (Driver: 550.67 CUDA: 12.4)
  • System: Manjaro Linux

基本编程模型

CUDA为GPU的硬件结构做了和CPU类似的抽象。以我的CPU配置为例，从硬件抽象来说，Ryzen5 1600有8个核心，每个核心有2个线程，因此可以同时运行16个线程。从软件抽象来说，我们则可以创建上百个逻辑线程，每一时刻最多有16个线程同时交给CPU运行，线程间的调度由操作系统负责，保证所有线程都有执行的可能性。

CUDA的GPU编程模型和CPU类似，thread是模型中的最小运行单位。从软件抽象来说，CUDA划分了grid、block和thread三个不同的层级，block包含一组固定数目的thread，而grid包含一组固定数目的block。比如以64个threads组成一个block,以10个blocks组成一个grid，那么这个grid实际上有10*64=640个thread。每个thread都会基于其所在的block有一个唯一的索引号（从0开始索引），而每个block都会基于所在的grid同样有一个唯一的索引号（从0开始索引）。因此每个thread在grid中的全局索引号是可以唯一确定下来的。block和grid总共可以有xyz三个维度（可以把它们想象成立方体）。

这种grid、block、thread的设计主要是为了耦合科学计算，因为科学计算的本质是去操作向量、矩阵等高维度的数组。比如我们希望相加两个32×3232 \times 32的矩阵A,B\mathbf{A},\mathbf{B}，那么可以设计为1个2维block，每个block在xy方向各含有32个threads，每个thread只负责一次简单标量加法Ci,j=Ai,j+Bi,j\mathbf{C}_{i,j} = \mathbf{A}_{i,j}+\mathbf{B}_{i,j}。理想情况下，我们可以让GPU同时执行这32×3232 \times 32个线程，在一次运算周期内完成矩阵加法的操作（相比让CPU每个线程执行标量加法的模型，我们总共需要循环32∗32/16=6432 * 32 / 16=64次运算周期来覆盖所有的元素，所以说多就是快，多就是好）。

GPU在硬件上的抽象则没有CPU那么好懂，有几个NVIDIA的基本概念需要澄清。streaming multiprocessor(SM) 类比于CPU的核心，是互相独立的组件；cuda core是物理上运行1个thread的计算单元，类比于CPU的线程；一个SM可以包含很多的cuda cores（类比于CPU核心有两个线程）。与CPU不同的地方在于，GPU对SM里的cuda cores做了进一步的划分，每32个cuda cores视作一个warp，它是启动thread计算的最小单元。需要16个threads来做计算，GPU会启动一个warp（32个threads）。需要48个threads？GPU会启动两个warps（64个threads）。多出来的threads会执行同样的计算，只不过运行的结果会被丢弃掉。

实际运行时，一个或多个blocks由GPU负责调度交给一个SM去运行，硬件会将block中的threads按32划分成warps并执行。一个SM能同时执行多少个blocks取决于程序的设计方式和对资源的占用多少，这也是CUDA优化的重点之一，即尽量让SM的利用率变高。比如GTX 1080总共有20个SMs，每个SM包含128个cuda cores，因此每个SM可以最多同时运行4个warps。如果我的程序需要40个blocks，每个block包含48个threads，并且每个thread占用的资源都很少，理论上来说，GTX 1080可以同时在一个SM上启动2个block的计算，所有的SMs都能有效利用，只不过每个SM里有32个threads是在做无用的计算。相比之下，RTX 4080有76个SMs，因此可以同时处理更多的blocks。

CPU和GPU有独立的内存空间，因此数据必须先从CPU移动到GPU上，处理结束后再从GPU移动到CPU上。类比于熟知的DRAM、L3 cache、L2 cache和L1 cache，GPU分为global memory、local memory、L2 cache、shared memory、L1 cache和registers。最主要的区别当然是存取速度有着数量级上的提升。global memory是最慢的，但是所有thread都能访问，我们使用cudaMalloc分配内存也都在global memory上，host（CPU）也能够访问global memory。L2 cache充当global memory和SMs的桥梁，所有的SMs都能访问。shared memory只有同一个block中的thread能够访问，它的速度大概是global memory的100倍以上。L1 cache和shared memory也类似。registers则是每个thread独有的，速度最快。local memory则有点特别，它是在register放不下时候的缓存区，本质上就是global memory，每个thread只能访问自己的部分。

NVCC编译

CUDA的核函数使用NVIDIA自定义的C++扩展语法编写，不属于C++标准的一部分，因此常见的编译器是不支持CUDA代码的编译，必须使用NVIDIA的nvcc编译器。nvcc编译器做的事情类似于对常规编译器的扩展：

  • 先将CUDA核函数与常规C++代码分离开
  • 将核函数编译成汇编（PTX）或者二进制（cubin）
  • 用CUDA runtime函数替换原始代码中的核函数调用部分
  • 调用常规C++编译器处理剩下的编译工作

最基本的nvcc命令如下所示，将vecadd.cu编译成可执行二进制文件vecadd：

nvcc vecadd.cu
        -gencode arch=compute_50,code=sm_50
        -gencode arch=compute_61,code=sm_61
        -gencode arch=compute_70,code=\"compute_70,sm_70\"
        -o vecadd

其中arch用来控制基于何种虚拟架构生成汇编代码，而code控制基于汇编如何生成和优化目标架构的二进制文件。所有的虚拟架构都以compute_xy表示，而真实架构都以sm_xy表示。我们主要重点关注code的部分即可。上述命令的意思是生成目标架构为sm_50的二进制，然后生成目标架构为sm_61的二进制，最后生成目标架构为sm_70的二进制和虚拟架构为compute_70的PTX代码。这些二进制和代码都被封装在同一个文件中，在实际执行时由CUDA运行时负责选择最接近的架构。

搞这么多复杂的架构主要还是为了兼容性的问题。NVIDIA有以下规定：

  • code=sm_xy生成的二进制只能在符合compute capability的x.y设备上运行，跨版本的二进制文件是不允许的
  • 基于arch=compute_ab生成的PTX可以编译更高版本的二进制文件code=sm_xy
  • 基于arch=compute_xy生成的PTX不可以编译低版本的二进制文件code=sm_ab
  • 如果code=compute_xy，则可以依赖JIT编译在compute capability的x.z设备上运行，其中z>=y
  • 如果code=compute_ab，则可以依赖JIT编译在compute capability的x.y设备上运行，其中x.y的版本要大于a.b版本
  • 基于arch=compute_ab生成的PTX不能编译高版本code=compute_xy，因为此时已经没办法附加更高版本的PTX代码进入最终文件，同理高版本arch=compute_xy也没有办法编译低版本code=compute_ab

GTX 1080的Compute Capability只到6.1，只支持最基本硬件特性。因此，上述编译的可执行文件可以在我的GPU上运行（因为存在和我的GPU架构匹配的二进制）。根据规则，下面命令编译的文件也可以通过JIT在我的GPU上运行，尽管二进制的目标架构是sm_70：

nvcc vecadd.cu -gencode arch=compute_61,code=\"compute_61,sm_70\" -o vecadd

在实施最简单的demo前，首先让我们看下GTX 1080的基本参数：

  • SMs: 20
  • Warps(CUDA Cores)/SM: 4(128)
  • Global Memory(GB): 8
  • L2 Cache(KB): 2048
  • Shared Memory(KB)/SM: 96
  • L1 Cache(KB)/SM: 48

矢量加法Demo

一个最简单的矢量加法demo如下：

#include <iostream>
#include <cuda_runtime.h>

#define N 10000

__global__ void addvec(int *a, int *b, int *c)
{
    int iThreadIdx = threadIdx.x + blockIdx.x * blockDim.x;
    int iStride = blockDim.x * gridDim.x;
    for (int i = iThreadIdx; i < N; i += iStride)
    {
        c[i] = a[i] + b[i];
    }
}

void initArray(int *p, int value)
{
    for (int i = 0; i < N; ++i)
    {
        p[i] = value;
    }
}

int main()
{
    int *a, *b, *c;
    a = (int *)malloc(N * sizeof(int));
    b = (int *)malloc(N * sizeof(int));
    c = (int *)malloc(N * sizeof(int));

    initArray(a, 1);
    initArray(b, 2);
    initArray(c, 0);

    int *da, *db, *dc;
    cudaMalloc(&da, N * sizeof(int));
    cudaMalloc(&db, N * sizeof(int));
    cudaMalloc(&dc, N * sizeof(int));

    cudaMemcpy(da, a, N * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(db, b, N * sizeof(int), cudaMemcpyHostToDevice);

    int iThreads = 64;
    int iBlocks = (N - 1) / iThreads + 1;
    std::cout << "Threads: " << iThreads << " Blocks: " << iBlocks << std::endl;

    addvec<<<iBlocks, iThreads>>>(da, db, dc);

    cudaDeviceSynchronize();

    cudaMemcpy(c, dc, N * sizeof(int), cudaMemcpyDeviceToHost);

    for (int i = 0; i < 10; ++i)
    {
        std::cout << c[i] << " ";
    }
    std::cout << std::endl;

    cudaFree(da);
    cudaFree(db);
    cudaFree(dc);
    free(a);
    free(b);
    free(c);

    return 0;
}

带有__global__的返回值为空的函数addvec就是在GPU上每个thread都运行的核函数，这里的threadIdx、blockIdx、blockDim和gridDim都是CUDA自己定义的内部变量，是uint的三元组，分别提供前面编程模型中提到的thread、block的索引以及block中thread的数量、block的数量。核函数就是利用这些信息让每个thread做相同的标量加法运算，但作用于不同的元素位置。此外，因为我们的元素数目要大于实际可用的cuda cores，所以要利用循环分批次处理数据。

main函数中分别在host（CPU）和device（GPU）上申请内存，并将CPU的内容复制到GPU内存中。addvec<<<iBlocks, iThreads>>>(da, db, dc)调用核函数并立即返回（异步执行），通过cudaDeviceSynchronize来等待所有thread计算完毕。最后将计算完成的内容copy回host的内存上。考虑到GTX 1080的架构，每个block我选择64个threads（2个warps），这样理想情况下每个SM可以塞进两个blocks的计算。

用Nsight分析程序性能

NVIDIA提供Nsight Systems工具来帮助分析程序性能，执行以下命令：

nsys profile --stats=true ./vecadd

会对我们的vecadd可执行文件做详细的分析，输出各种运行时间、调用次数的报告。需要注意的是，nsys执行过程相当漫长，terminal上会显示“The target application terminated. One or more process it created re-parented. Waiting for termination of re-parented processes.”的提示，我一开始还以为是程序出错，其实是正常现象，只要耐心等待就好。执行完毕后，会在当前文件夹下生成.nsys-rep和.sqlite文件，以便后续追踪分析。

Time (%)  Total Time (ns)  Instances  Avg (ns)  Med (ns)  Min (ns)  Max (ns)  StdDev (ns)             Name            
--------  ---------------  ---------  --------  --------  --------  --------  -----------  ---------------------------
   100.0             2881          1    2881.0    2881.0      2881      2881          0.0  addvec(int *, int *, int *)

在上例中，我们的addvec核函数显示总运行时间是2.8us；如果稍作修改，比如iThreads=48，运行时间则为：

Time (%)  Total Time (ns)  Instances  Avg (ns)  Med (ns)  Min (ns)  Max (ns)  StdDev (ns)             Name            
--------  ---------------  ---------  --------  --------  --------  --------  -----------  ---------------------------
   100.0             3007          1    3007.0    3007.0      3007      3007          0.0  addvec(int *, int *, int *)

可以看出时间稍长了一些，符合我们的预期，因为有些warp里的thread计算被浪费掉了。