1 GPGPU简介

在上两篇中，我们使用多种技术提高了Blur算子的速度，效果最好的技术是利用多核和SIMD指令集，通过并行计算来提升速度。这给我们一个很大的启发，如果我们的CPU能扩展到数百数千甚至数万的核，不就能进一步加速了吗？事实上，传统超算就是这么干的，比如我国的神威-太湖之光。但是传统超算的做法是一种分布式计算的思路，即用数十个机柜，每个机柜内放置数百个计算单元（看成一台服务器），每个计算单元有几个CPU，以此组成一个巨大的CPU集群。集群间通过高速网络，配合上分布式的编程来使用整个集群的数以万记的CPU执行计算任务。

之所以用分布式计算的思路提高整个系统的并行计算的能力，主要原因是CPU的发展已面临物理定律的限制，人们无法在有限的成本下创造出核心数更多的CPU。此时人们发现，原本用来做图形渲染的专用硬件，GPU（Graphic Processing Unit），非常适合用于并行计算。

我个人认为原因有三点：首先，GPU在渲染画面时需要同时渲染数以百万记的顶点或三角形，因此GPU硬件设计时就利用多种技术优化并行计算/渲染的效率；第二，GPU作为专用的硬件，只负责并行计算/并行渲染这一类任务，不承担复杂的逻辑控制和分时切换等任务，可更专注于计算本身，通过专用指令集和流水线等技术，优化计算效率和吞吐；最后，GPU作为一类外插设备，在尺寸、功率、散热、兼容性等方面的限制远远小于CPU，让GPU可以有更大的缓存、显存和带宽，可以放置更多的计算核心，而这进一步体现到成本上的优势（成本对于技术的普及非常重要）。

CPU和GPU对比示意图

上图是CPU和GPU的一个对比示意图。由于CPU除了计算任务外，还需要负责大量的逻辑控制，分时调度，兼容多种指令集等等“包袱”，真正用来计算的部分其实很少。而GPU则不然，GPU可以在更大的自由度下，放置更多的计算核心。

其实从最本质的角度来讲，GPU之所以适合并行计算场景，是因为GPU使用的是SIMT（Single Instruction Multiple Threads）模型。SIMT可以看成是SIMD（Single Instruction Multiple Data）模型的一种增强。

Nvidia公司在2000年推出第一款GPU，但是直到2006年才把GPU用于通用计算任务，其转折点就是CUDA（Compute Unified Device Architecture）的出现。CUDA是一整套软件技术，把Nvidia GPU的计算能力封装成一套编程规范和编程语言接口，同时提供相应的编译器、调试器等工具，这使得程序员能方便的使用GPU并行计算能力，执行通用的计算任务，因此又被称为GPGPU（General-purpose computing on graphics processing units）。在CUDA出现6年之后，Alex Krizhevsky和其老师Geoffrey Hinton在CUDA的帮助下成功地实现了AlexNet深度神经网络，一举打开了此轮人工智能浪潮的大门。AlexNet论文中提到：模型的深度对于提高算法性能至关重要，但是其计算成本很高，GPU的使用让深度神经网络的训练具有可行性。

甚至可以武断地说，如果没有GPU和CUDA，尤其是后者，以深度神经网络为主要代表的人工智能算法很可能是另外一种命运。当然，Nvidia并不是唯一一家生产GPU和提供上层编程框架的公司，除了Nvidia GPU和CUDA之外，Intel、AMD和最近几年爆火的如寒武纪、海思等公司，也推出了不同的硬件和软件平台，但是无论完善度、丰富性尤其是生态上，Nvidia仍然是毫无疑问的王者。因此这篇文章介绍的还是Nvidia的GPU，下文中提到的GPU，没有特殊说明都代指Nvidia GPU。

Nvidia的软硬件生态

上图是Nvidia提供的一整套软硬件平台。Nvidia针对嵌入式端，桌面级消费卡、专业工作站和数据中心提供了不同的硬件产品，但是在软件层，通过CUDA抽象成了统一的编程接口并提供C/C++/Python/Java等多种编程语言的支持。CUDA的这一层抽象非常重要，因为有了这层抽象，向下开发者的代码可在不同的硬件平台上快速迁移；向上，在CUDA基础上封装了用于科学计算的cuBLAS（兼容BLAS接口），用于深度学习的cuDNN等中间件和代码库。这些中间件对于Tensorflow，PyTorch这一类的深度学习框架非常重要，可以更容易地使用CUDA和底层硬件进行机器学习的计算任务。

2 GPU硬件架构

对于一般开发者而言，大多数情况下接触到的都是CUDA层。但是由于GPU的特殊性，为了能正确高效使用CUDA，非常有必要学习一下GPU的硬件架构。由于Nvidia的GPU发展非常迅速，平均1-2年就会推出一款新GPU或核心架构，因此这里简单介绍下，不涉及硬件细节（主要是我也不太懂硬件）。

Nvidia A100 GPU的硬件架构

我们以Nvidia 2020年推出的Ampere架构和A100 GPU为例（A10、A30、RTX3090等型号也是Ampere架构）为例。上图是整个A100的硬件架构，从上到下，其主要组成部分：

PCIe接口层，目前大部分GPU都是通过PCIe以外部设备的方式集成到服务器上；

GigaThread Engine with MIG是GPU用于调度资源的引擎，Ampere架构还支持了MIG（Multiple Instance GPU）控制即允许对GPU进行硬件隔离切分；

中间占据面积最大的绿色部分是GPU的核心计算单元，共有6912个CUDA Cores，这是重点部分，我们在后面详细介绍；

中间蓝色部分是L2缓存，在GPU层面提供大于40MB的共享缓存；

最下层的NVLink是用于多个GPU间快速通信的模块（在分布式训练中非常重要）；

两侧的HBM2（High Bandwidth Memory）即我们通常说的显存。A100提供了高达1.56TB带宽，40GB/80GB的显存；

其中最重要的绿色部分是计算核心，在A100中这些计算单元又进一步被“分组”成了GPC（Graphic Processing Cluster）- TPC（Texture Processing Cluster） - SM（Streaming Multiprocessor）几个概念。最需要了解的是SM的架构，把其中一个SM放大：

Nvidia A100的SM硬件架构

SM可以看成是GPU计算调度的一个基本单位，只有对其架构和硬件规格有一定认识，才能更高效地利用GPU的硬件能力。SM中自上而下分别有：

L1 Shared Memory，在SM层面提供192KB的共享缓存

接下来是四个PB（Process Block），每个PB内部分别有

L0缓存

Warp Scheduler用于计算任务调度

Dispatch Unit用于计算任务发射

Register File寄存器

16个专用于INT32计算的核心

16个专用于FP32计算的核心

8个专用于FP64计算的核心

4个Tensor Core，这是受Google TPU的影响，专门设计用于执行A*B+C这类矩阵计算的硬件单元，甚至还支持稀疏矩阵计算，Tensor Core对深度学习的加速非常明显

8个用于LOAD/STORE数据的硬件单元

4个SFU（Special Function Unit）用于执行超越函数、插值等特殊计算

因此一个SM内就有64个用于计算的核心，Nvidia称之为CUDA Core（但是我没搞明白64是怎么算出来的）。而一整块A100 GPU则密密麻麻地放置了108个SM，总计6912个CUDA Core，可以简单理解成这是一个拥有6912个核心的巨型CPU！

至此，我们应该能搞懂为什么GPU比CPU更适合执行并行计算了，就是“专业人干专业事，大力就会出奇迹“嘛。

3 CUDA基本概念

为了方便地使用GPU的硬件能力，Nvidia在软件层做了抽象，形成软件上的一系列逻辑概念，这些逻辑概念跟GPU的硬件有一定的对应关系。

GPU软件和硬件的对应关系

如上图，CUDA编程中的最小单元称之为Thread，可以简单认为一个软件Thread会在一个硬件CUDA Core中执行，而Thread中执行的内容或函数称之为Kernel。多个相同的Thread组成一个Thread Block，软件Thread Block会被调度到一个硬件SM上执行，同一个Thread Block内的多个Thread执行相同的kernel并共享SM内的硬件资源。而多个Thread Block又可以进一步组成一个Grid，一个软件Grid可以看成一次GPU的计算任务，被提交到一整个GPU硬件中执行。这几个概念非常重要，简单总结下：

kernel：Thread执行的内容/代码/函数

Thread：执行kernel的最小单元，被SM调度到CUDA Core中执行（其实还有一个Warp的概念，为了简单，这里先略过）

Thread Block：多个Thread组合，GPU任务调度的最小单元（这个描述不太准确，应该是Warp，为了简单暂时先不细究），被调度到SM中执行。一个SM可以同时执行多个Thread Block，但是一个Thread Block只能被调度到一个SM上。

Grid：多个Thread Block的组合，被调度到整个GPU中执行

同时，Thread、Thread Block和Grid由于所处层次不同，他们可以访问的存储资源也不同。如Thread只能访问自身的寄存器，Thread Block可以访问SM中的L1缓存，而Grid则可以访问L2缓存和更大的HBM显存。在第一篇文章中我们就介绍过不同层次的存储其访问速度往往是数量级的差别，GPU也不例外，在后续的文章中我们会看到，针对CUDA的优化很大一部分就是如何正确高效的使用GPU中的多级存储来提高GPU的方寸比，从而进一步提高GPU的计算效率。

Nvidia GPU的内存模型

到底如何使用CUDA？由于GPU是作为一类外挂设备，通过PCIe之类的接口插到服务器上，因此在CUDA编程中，称GPU为Device，而服务器上的CPU和内存统称为Host。在编写CUDA代码时，最主要的工作就是编写kernel函数，然后利用CUDA提供的接口，把kernel函数从Host发射到Device中执行，Device则从 Grid → Thread Block → Warp → Thread 一层层的调度下最终完成kernel的计算。

在计算开始前，还需要用CUDA接口把计算要用到的数据从Host拷贝到Device中。Device完成计算后，则需要利用CUDA接口把计算结果从Device拷贝回Host中，因此总结下来，CUDA编程分为三步：

从Host拷贝数据到Device

把需要Device执行的kernel函数发射给Device

从Device拷贝计算结果到Host

如果有点抽象，我们看个例子。

4 CUDA HelloWorld

CUDA HelloWorld代码很简单，求两个数组A和B的和。代码分两部分，首先实现上面介绍的三步流程。

#include <cuda.h> 
// Host code
int main() {
  int N = 1024;
  size_t size = N * sizeof(float);
 
  // Host端申请内存
  float* h_A = (float*)malloc(size);
  float* h_B = (float*)malloc(size);
  float* h_C = (float*)malloc(size);
 
  // Device端申请显存，分别用来存放A、B和结果C
  float* d_A, d_B, d_C;
  cudaMalloc(&d_A, size);
  cudaMalloc(&d_B, size);
  cudaMalloc(&d_C, size);
 
  // 把需要计算的结果从Host拷贝到Device
  cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
  cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
 
  // 执行kernel，需要指定两个参数
  int threadsPerBlock = 256;
  int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
  VecAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);
 
  // 把计算结果拷贝回Host
  cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);
 
  // 释放Device中的显存
  cudaFree(d_A);
  cudaFree(d_B);
  cudaFree(d_C);
 
  // 释放Host的内存
  ...
}

代码比较直白，其中cuda开头的函数都是CUDA库中提供的，比如cudaMalloc会在Device中申请一段显存，cudaMemcpy则可以在Host和Device之间进行数据拷贝等，不再赘述。

核心是kernel函数部分，先看看kernel函数的实现：

__global__ void VecAdd(float* A, float* B, float* C, int N) {
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  if (i < N) C[i] = A[i] + B[i];

Kernel函数只有两行并不复杂，但是需要好好解释。

这个kernel函数把数组A和数组B逐元素相加，结果写到数组C中。即C[i] = A[i] + B[i]。如果在CPU上编程，一个循环遍历就搞定了。但是当用CUDA编程时，我们希望利用GPU中数千个计算核心同时计算，每个核心只负责A和B的一个元素。但是，如何分配任务？

如果GPU中的每个核心都有唯一的ID，比如0号核心执行C[0] = A[0] + B[0]，127号核心执行C[127] = A[127] + B[127]就好了。但是如果我们要操作的数据大于GPU的核心数怎么办？因此CUDA做了一层逻辑抽象。当一个计算任务在GPU中执行时，这个任务相关的Grid、Thread Block和Thread都有一系列的身份标识，即几个全局变量：

blockDim：表明kernel所在的Thread Block的尺寸，包含x, y, z三个维度

blockIdx：表明kernel所在的Thread Block的index，包含x, y, z三个维度

threadIdx：表明kernel在在的Thead在Thread Block内的idx，包括x, y, z三个维度

只要引入了cuda.h的头文件，在代码中可以直接使用上面几个全局变量，这几个全局变量的值会在执行时由CUDA自动维护更新。因此，kernel函数 i 的计算方法就是用blockDim.x乘以blockIdx.x再加上threadIdx.x，就可以算出当前一个唯一的、逻辑ID被当前的thread使用。那为什么只用到了x维度，没有用到y和z维度？

答案在于Host端代码中，在启动kernel时，代码中指定了两个模板类参数：blocksPerGrid和threadsPerBlock。其含义通过名字就能看懂。这两个参数被我们指定成两个int型数字。

int threadsPerBlock = 256;
int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
 
VecAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);

其实这两个参数还可以是dim2或dim3类型。如果是dim2类型，那么逻辑ID可以通过逻辑x和逻辑y唯一确定，其计算方法如下。

int y = blockIdx.y * blockDim.y + threadIdx.y;
int x = blockIdx.x * blockDim.x + threadIdx.x;
int i = blockDim.y * y + x;

原理可以参考上图。这与二维数组索引转换成一维数组索引的过程是类似的，相信有基本编程经验的同学很快就能搞懂。当然这里我们只是举了两个最简单的例子，在实际编程中这块会更复杂一些。网上关于这部分的内容比较多，大家可以自行学习。