随着近年来深度学习的爆发，原来被用于图形渲染的 GPU 被大量用于并行加速深度学习的模型训练中，在这个过程中 CUDA 作为 NVIDIA 推出的基于 GPU 的一个通用并行计算平台和编程模型也得到了广泛的使用。或许你已经十分了解 现代 CPU 的体系架构，但是对于 GPU 还不甚清晰，GPU 的体系架构到底和 CPU 有何区别，CUDA 模型是什么，我们该如何使用 CUDA 实现并行计算，本文将为你扫盲祛魅，本文中使用到的所有代码可以在我的 GitHub 中找到。

GPU 体系架构

为什么我们需要 GPU

如前所述，GPU （Graphics Processing Unit）最开始只是用于游戏、视频中的图形渲染，而现在最热门的一个应用领域是在深度学习的加速计算上。为什么需要 GPU 来加速计算呢？我们知道，随着摩尔定律的发展，在过去五十年间 CPU 的性能获得了巨大的提升，不论是从芯片上晶体管数目，还是时钟频率，到后来的从单核处理器发展到后来的多核多处理器。

下图是过去五十年间各款 CPU 处理器上晶体管数目的变化，基本上满足每 18 个月提升一倍的规律，虽然现在看起来 50 十年后摩尔定律对 CPU 来说有停滞的迹象（这是另一个话题，此处不表）

在 CPU 算力快速提升的这五十年，人们需要的计算量也同时在迅猛发展着，从最开始的桌面互联网，到后来的移动互联网，以及 5 年前爆发的深度学习，无一不需要庞大的计算力。在这个过程中，仅仅依靠 CPU 的算力开始力有不逮，这个过程中像 GPU、FPGA、DSP 等异构计算单元开始得到广泛的应用。下面，我回归计算的本质，以 GPU 为例来分析为什么我们需要这些异构计算单元。

无论是 CPU 还是 GPU，我们可以把计算模型抽象为下面这张图，这也是典型的冯诺伊曼体系架构。

影响计算能力的 4 个主要因素如下：

  • Parallel Processing：Amount of data processed at one time
  • Clock Frequency：Processing speed on each data element
  • Memory Bandwidth：Amount of data transferred at one time
  • Memory Lantency：Time for each data element to be transferred

对于 CPU，依次分析这几个因素：

  • 为了提供并行处理能力，我们从单核单处理器发展到多核多处理器，每个时钟周期 CPU 也能够处理多条指令
  • 因为 CPU 时钟频率和功率的关系 $ Power = k _ ClockFrequency _ Voltage^2 $ ，在 CPU 过去的发展历史中，通过提高 CPU 时钟频率可以变得更快，与此同时为了保持 CPU 功耗的正常，也需要不断降低电压。但是当主频逐渐逼近到 4GHz 时，电压已经不能再降低了，因为这已经到达了晶体管高低电平反转的极限，关于这部分的更多内容可以参考 摩尔定律 。
  • 现在 CPU 用的是常规的 DDR 内存，明显存在着内存带宽限制
  • 从 CPU 到 DDR 内存的延时很高，2020 年的时候大概有 100ns，具体可以参考 Key Numbers Every Programmer Should Know。CPU 通过其他的方式隐藏了这个问题：
    • Large On-Chip Low-Latency Cache，大概 1ns
    • MultiThreading
    • Out-of-order execution

Credit to https://queue.acm.org/detail.cfm?id=2181798
Credit to https://queue.acm.org/detail.cfm?id=2181798

尽管现在 CPU 的能力还在发展，但是以上的问题极大的限制了其算力的提高，当前仅靠 CPU 已经不能够满足人们对庞大算力的需求了。因此我们需要其他的专用芯片来帮助 CPU 一起计算，这就是异构计算的来源。GPU 等专用计算单元虽然工作频率较低，但具有更多的内核数和并行计算能力，总体性能/芯片面积比和性能/功耗比都很高。随着人工智能时代的降临，GPU 从游戏走进了人们的视野。

无论是 CPU 还是 GPU，在进行计算时都需要用核心（Core）来做算术逻辑运算。核心中有 ALU（逻辑运算单元）和寄存器等电路。在进行计算时，一个核心只能顺序执行某项任务。CPU 作为通用计算芯片，不仅仅做算术逻辑计算，其很重要的一部分功能是做复杂的逻辑控制，一般而言 CPU 上的 Core 数目相对较少，数据中心的服务器一般也就 40 左右个 CPU 核心。但是 GPU 动辄有上千个核心，这些核心可以独立的进行算术逻辑计算，大大提高了并行计算处理能力。

GPU 时代的最大获益者是 NVIDIA，当然 AMD 他们家也有 GPU 产品，但是因为 AMD 并没有形成 CUDA 这样的软件生态导致深度学习中主要用的都是 NVIDIA 的 GPU，后面的分析都将基于 NVIDIA 的 GPU 产品。NVIDIA 不同时代产品的芯片设计不同，每代产品背后有一个架构代号，架构均以著名的物理学家为名，以向先贤致敬，对于消费者而言，英伟达主要有两条产品线：

  • 消费级产品 GeForce 系列：GeForce 2080 Ti…
  • 高性能计算产品 Telsa 系列：Telsa V100、Telsa P100、Telsa P40…

NVIDIA GPU产品体系
NVIDIA GPU产品体系

GPU 硬件模型

Host and Device

GPU 并不是一个独立运行的计算平台，而是需要与 CPU 的协同工作，可以看作是 CPU 的协处理器，因此当我们说 GPU 并行计算的时候，实质上是指的 CPU+GPU 的异构计算架构。由于 CPU 和 GPU 是分开的，在 NVIDIA 的设计理念里，CPU 和主存被称为 Host，GPU 和显存被称为 Device。Host 和 Device 概念会贯穿整个 NVIDIA GPU 编程。

基于 CPU + GPU 的异构计算平台可以优势互补，CPU 负责处理逻辑复杂的串行程序，GPU 重点处理数据密集型的并行计算程序，从而发挥最大功效。CUDA 程序中既包含 Host 程序，又包含 Device 程序，它们分别在 CPU 和 GPU 上运行。

同时， Host 与 Device 之间通过 PCIe 总线交互进行数据拷贝，典型的 CUDA 程序的执行流程如下：

  1. 初始化后，将数据从 Main Memory 拷贝到 GPU Memory
  2. CPU 调用 CUDA 的核函数
  3. GPU 的 CUDA Core 并行执行核函数
  4. 将 Device 上的运算结果拷贝到 Host 上

GPU 核心在做计算时，只能直接从显存中读写数据，程序员需要在代码中指明哪些数据需要从内存和显存之间相互拷贝。这些数据传输都是在总线上，因此总线的传输速度和带宽成了部分计算任务的瓶颈。当前最新的总线技术是 NVLink，IBM 的 Power CPU 和 NVIDIA 的高端显卡可以通过 NVLink 直接通信，Intel 的 CPU 目前不支持 NVLink，只能使用 PCIe 技术。同时，单台机器上的多张英伟达显卡也可以使用 NVLink 相互通信，适合多 GPU 卡并行计算的场景。

NVLink可以连接CPU和GPU
NVLink可以连接CPU和GPU

Streaming Multiprocessor

在 NVIDIA 的设计里，一张 GPU 卡有多个 Streaming Multiprocessor（SM），每个 SM 中有多个计算核心，SM 是运算和调度的基本单元。下图为当前计算力最强的显卡 Tesla V100，密密麻麻的绿色小格子就是 GPU 小核心，多个小核心一起组成了一个 SM。

Tesla V100 with 84 SM Units
Tesla V100 with 84 SM Units

将 SM 放大，单个 SM 的结构如图所示：

Tesla V100 Streaming Multiprocessor(SM)
Tesla V100 Streaming Multiprocessor(SM)

可以看到一个 SM 中包含了计算核心和存储部分，SM 的核心组件包括 CUDA 核心，共享内存，寄存器等，SM 可以并发地执行数百个线程，并发能力就取决于 SM 所拥有的资源数。

  • 针对不同计算的小核心（绿色小格子），包括优化深度学习的 TENSOR CORE，32 个 64 位浮点核心（FP64），64 个整型核心(INT)，64 个 32 位浮点核心(FP32)
  • 计算核心直接从寄存器（Register）中读写数据
  • 调度和分发器（Scheduler 和 Dispatch Unit）
  • L0 和 L1 级缓存

具体而言，SM 中的 FP32 进行 32 位浮点加乘运算，INT 进行整型加乘运算，SFU（Special Functional Unit）执行一些倒数和三角函数等运算。Tensor Core 是 NVIDIA 新的微架构中提出的一种混合精度的计算核心。我们知道，当前深度神经网络中使用到最频繁的矩阵运算是： $ D = A \times B + C $。Tensor Core 可以对 $ 4 \times 4 $ 的矩阵做上述运算。其中：

  • 涉及乘法的 A 和 B 使用 FP16 的 16 位浮点运算，精度较低
  • 涉及加法的 C 和 D 使用 FP16 或 FP32 精度

Tensor Core 是在 Volta 架构开始提出的，使用 Volta 架构的 V100 在深度学习上的性能远超 Pascal 架构的 P100。

Tensor Core是一种为优化深度学习计算核心
Tensor Core是一种为优化深度学习计算核心

CUDA 编程模型

前面提到，NVIDIA 相对于 AMD 的一个巨大优势是它的 CUDA 软件生态，下图是 NVIDIA GPU 编程的软件栈，从底层的 GPU 驱动和 CUDA 工具包，上面还提供了科学计算所必需的 cuBLAS 线性代数库，cuFFT 快速傅里叶变换库以及 cuDNN 深度神经网络加速库，当前常见的 TensorFlow 和 PyTorch 深度学习框架底层大多都基于 cuDNN 库。

Hello World

在进一步学习 CUDA 编程模型之前，我们首先配置好 CUDA 的运行环境，跑通 Hello World 从而对 CUDA 编程有一个直观的认识，这里使用的是腾讯云的 GPU 服务器，机器安装的是 CentOS 7 系统，CUDA 环境配置可以参考 CUDA Installation Guide Linux 。

根据上图的 NVIDIA GPU 软件栈，有了一个插上了 GPU 的服务器之后，我们首先查看机器上的 GPU，可以看到当前机器上装 GPU 是 Tesla P40：

1
2

$ lspci | grep -i nvidia
00:08.0 3D controller: NVIDIA Corporation GP102GL [Tesla P40] (rev a1)

接下来在 这里下载 CUDA Toolkit，这里选择的是 rpm local 的安装方式：

1
2
3
4
5

$ wget https://developer.download.nvidia.com/compute/cuda/11.1.1/local_installers/cuda-repo-rhel7-11-1-local-11.1.1_455.32.00-1.x86_64.rpm
$ sudo rpm -i cuda-repo-rhel7-11-1-local-11.1.1_455.32.00-1.x86_64.rpm
$ sudo yum clean all
$ sudo yum -y install nvidia-driver-latest-dkms cuda
$ sudo yum -y install cuda-drivers

执行上面的安装操作之后，我们可以看到在 /usr/lib64/ 看到 libcuda.so ：

1
2
3
4

$ ls /usr/lib64 -al | grep cuda
lrwxrwxrwx   1 root root        20 Nov 21 15:05 libcuda.so -> libcuda.so.455.32.00
lrwxrwxrwx   1 root root        20 Nov 21 15:05 libcuda.so.1 -> libcuda.so.455.32.00
-rwxr-xr-x   1 root root  21074296 Oct 15 06:58 libcuda.so.455.32.00

下面是一些我们会经常用到的 CUDA 工具，你需要通过配置环境变量来使用他们：

1
2
3
4

编译器：nvcc (C/C++)
调试器：nvcc-gdb
性能分析：nsight, nvprof
函数库：cublas, nvblas, cusolver, cufftw, cusparse, nvgraph

设置环境变量如下：

1
2
3
4
5
6
7

$ export PATH=/usr/local/cuda-11.1/bin${PATH:+:${PATH}}
$ nvcc --version
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2020 NVIDIA Corporation
Built on Mon_Oct_12_20:09:46_PDT_2020
Cuda compilation tools, release 11.1, V11.1.105
Build cuda_11.1.TC455_06.29190527_0

除此之外，对于 64 位系统，需要设置 LD_LIBRARY_PATH：

1

$ export LD_LIBRARY_PATH=/usr/local/cuda-11.1/lib64${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}

这个时候可以确认驱动的版本：

1
2
3

$ cat /proc/driver/nvidia/version
NVRM version: NVIDIA UNIX x86_64 Kernel Module  455.32.00  Wed Oct 14 22:46:18 UTC 2020
GCC version:  gcc version 4.8.5 20150623 (Red Hat 4.8.5-39) (GCC)

可以使用 nvidia-smi命令查看显卡情况，比如这台机器上几张显卡，CUDA 版本，显卡上运行的进程等。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

$ nvidia-smi
Sat Nov 21 17:09:13 2020
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 455.32.00    Driver Version: 455.32.00    CUDA Version: 11.1     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  Tesla P40           Off  | 00000000:00:08.0 Off |                    0 |
| N/A   27C    P0    49W / 250W |      0MiB / 22919MiB |      3%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+

+-----------------------------------------------------------------------------+
| Processes:                                                                  |
|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
|        ID   ID                                                   Usage      |
|=============================================================================|
|  No running processes found                                                 |
+-----------------------------------------------------------------------------+

CUDA 自己提供了一系列的代码示例，可以通过下面的方法安装：

1

$ cuda-install-samples-11.1.sh <dir>

在对应目录下，我们可以看到 CUDA 提供的源代码：

1
2
3

$ ls NVIDIA_CUDA-11.1_Samples
0_Simple     2_Graphics  4_Finance      6_Advanced       bin     EULA.txt  Makefile
1_Utilities  3_Imaging   5_Simulations  7_CUDALibraries  common  LICENSE

直接在这个目录下执行 make，可以在 bin目录下得到所有代码的二进制程序，选择其中的 deviceQuery 执行：

 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

$ ./deviceQuery
./deviceQuery Starting...

 CUDA Device Query (Runtime API) version (CUDART static linking)

Detected 1 CUDA Capable device(s)

Device 0: "Tesla P40"
  CUDA Driver Version / Runtime Version          11.1 / 11.1
  CUDA Capability Major/Minor version number:    6.1
  Total amount of global memory:                 22919 MBytes (24032378880 bytes)
  (30) Multiprocessors, (128) CUDA Cores/MP:     3840 CUDA Cores
  GPU Max Clock rate:                            1531 MHz (1.53 GHz)
  Memory Clock rate:                             3615 Mhz
  Memory Bus Width:                              384-bit
  L2 Cache Size:                                 3145728 bytes
  Maximum Texture Dimension Size (x,y,z)         1D=(131072), 2D=(131072, 65536), 3D=(16384, 16384, 16384)
  Maximum Layered 1D Texture Size, (num) layers  1D=(32768), 2048 layers
  Maximum Layered 2D Texture Size, (num) layers  2D=(32768, 32768), 2048 layers
  Total amount of constant memory:               65536 bytes
  Total amount of shared memory per block:       49152 bytes
  Total shared memory per multiprocessor:        98304 bytes
  Total number of registers available per block: 65536
  Warp size:                                     32
  Maximum number of threads per multiprocessor:  2048
  Maximum number of threads per block:           1024
  Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
  Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)
  Maximum memory pitch:                          2147483647 bytes
  Texture alignment:                             512 bytes
  Concurrent copy and kernel execution:          Yes with 2 copy engine(s)
  Run time limit on kernels:                     No
  Integrated GPU sharing Host Memory:            No
  Support host page-locked memory mapping:       Yes
  Alignment requirement for Surfaces:            Yes
  Device has ECC support:                        Enabled
  Device supports Unified Addressing (UVA):      Yes
  Device supports Managed Memory:                Yes
  Device supports Compute Preemption:            Yes
  Supports Cooperative Kernel Launch:            Yes
  Supports MultiDevice Co-op Kernel Launch:      Yes
  Device PCI Domain ID / Bus ID / location ID:   0 / 0 / 8
  Compute Mode:
     < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >

deviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 11.1, CUDA Runtime Version = 11.1, NumDevs = 1
Result = PASS

到现在，CUDA Toolkit 安装完毕，接下来通过编写一个简单的 hello world 来直观感受 CUDA 编程：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

#include <stdio.h>

__global__ void hello_from_gpu()
{
    printf( "\"Hello, world!\", says the GPU.\n" );
}

void hello_from_cpu()
{
    printf( "\"Hello, world!\", says the CPU.\n" );
}

// host code entrance
int main( int argc, char **argv )
{
    hello_from_cpu();
    hello_from_gpu <<< 2, 4>>>();
    cudaDeviceReset();
    return 0;
}

可以看到，CUDA 程序基本上和标准 C 语言程序一样，主要的区别在于 __global__ 限定词 和 <<<... >>> 符号。其中 __global__ 标记用来告诉编译器这段代码会运行在 Device （GPU）上，它会被运行在 Host 上的代码调用，也被称作是在 Device 上线程中并行执行的核函数（Kernel），是在 Device 上线程中并行执行的函数。

当一个核函数被调用时，需要通过 <<<grid, block>>> 符号 来设置核函数执行时的配置，在 CUDA 的术语中，这称作 kernel lauch，在后面我们将深入介绍这部分。

hello world 程序写完，我们以 hello.cu 这样的后缀名来保存，接下来使用 nvcc 来编译，整体上用法与 gcc 几乎一样：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

$ nvcc hello.cu -o hello
$./hello
"Hello, world!", says the CPU.
"Hello, world!", says the GPU.
"Hello, world!", says the GPU.
"Hello, world!", says the GPU.
"Hello, world!", says the GPU.
"Hello, world!", says the GPU.
"Hello, world!", says the GPU.
"Hello, world!", says the GPU.
"Hello, world!", says the GPU.

可以看到，来自 CPU 的 Hello World 执行了一次，来自 GPU 的 Hello World 执行了 8 次。

核函数与线程模型

上文提到，为了实现 GPU 并行加速计算，我们需要在 Host 上执行 kernel launch，让 核函数 在 Device 上的多个线程并发执行。具体的方式就是在调用核函数的时候通过 <<<grid, block>>> 来指定核函数要执行的线程数量 N，之后 GPU 上的 N 个 Core 会并行执行核函数，并且每个线程会分配一个唯一的线程号 threadID，这个 ID 值可以通过核函数的内置变量 threadIdx来获得。

CUDA 将核函数所定义的运算称为线程（Thread），多个线程组成一个块（Block），多个块组成网格（Grid）。这样一个 Grid 可以定义成千上万个线程，也就解决了并行执行上万次操作的问题。 <<<grid, block>>> 中括号中第一个数字表示整个 Grid 有多少个 Block，括号中第二个数字表示一个 Block 有多少个 Thread。前面 Hello World 用 2 个 Block，每个 Block 中有 4 个 Thread，所以总共执行了 8 次。

Grid of Thread Blocks
Grid of Thread Blocks

实际上，线程（Thread）是一个编程上的软件概念。从硬件来看，Thread 运行在一个 CUDA 核心上，多个 Thread 组成的 Block 运行在 Streaming Multiprocessor（SM），多个 Block 组成的 Grid 运行在一个 GPU 显卡上。当一个 kernel 被执行时，它的 gird 中的线程块被分配到 SM 上，一个线程块只能在一个 SM 上被调度。SM 一般可以调度多个线程块，这要看 SM 本身的能力。那么有可能一个 kernel 的各个线程块被分配多个 SM，所以 grid 只是逻辑层，而 SM 才是执行的物理层。

grid 和 block都是定义为 dim3类型的变量，dim3可以看成是包含三个无符号整数（x，y，z）成员的结构体变量，在定义时，缺省值初始化为 1。因此 grid 和 block 可以灵活地定义为 1-dim，2-dim 以及 3-dim 结构，对于上图中结构（主要水平方向为 x 轴），定义的 grid和 block 如下所示， kernel 在调用时也必须通过执行配置 <<<grid, block>>>来指定 kernel 所使用的线程数及结构。

1
2
3

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

所以，一个线程需要两个内置的坐标变量 （blockIdx，threadIdx）来唯一标识，它们都是 dim3类型变量，其中 blockIdx 指明线程所在 grid 中的位置，而 threaIdx 指明线程所在 block 中的位置，如图中的 Thread (1,1) 满足：

1
2
3
4

threadIdx.x = 1
threadIdx.y = 1
blockIdx.x = 1
blockIdx.y = 1

不同的执行配置会影响 GPU 程序的速度，一般需要多次调试才能找到较好的执行配置，在实际编程中，执行配置 <<<grid, block>>>应参考下面的方法：

  • Block 运行在 SM 上，不同硬件架构（Turing、Volta、Pascal…）的 CUDA 核心数不同，一般需要根据当前硬件来设置 Block 的大小 block（执行配置中第二个参数）。一个 Block 中的 Thread 数最好是 32、128、256 的倍数。注意，限于当前硬件的设计，Block 大小不能超过 1024。
  • Grid 的大小 grid（执行配置中第一个参数），即一个 Grid 中 Block 的个数可以由总次数 N除以 block，并向上取整。

例如，我们想并行启动 1000 个 Thread，可以将 blockDim 设置为 128，1000 ÷ 128 = 7.8，向上取整为 8。使用时，执行配置可以写成 gpuWork<<<8, 128>>>()，CUDA 共启动 8 * 128 = 1024个 Thread，实际计算时只使用前 1000 个 Thread，多余的 24 个 Thread 不进行计算。

{% note info%}

这几个变量比较容易混淆，再次明确一下：block是 Block 中 Thread 的个数，一个 Block 中的 threadIdx最大不超过 block；grid是 Grid 中 Block 的个数，一个 Grid 中的 blockIdx最大不超过 grid。

{% endnote %}

这几个变量比较容易混淆，再次明确一下：block是 Block 中 Thread 的个数，一个 Block 中的 threadIdx最大不超过 block；grid是 Grid 中 Block 的个数，一个 Grid 中的 blockIdx最大不超过 grid。

kernel 的这种线程组织结构天然适合 vector，matrix 等运算，我们将在后面实现向量加法和矩阵乘法。如我们将利用上图 2-dim 结构实现两个矩阵的加法，每个线程负责处理每个位置的两个元素相加，代码如下所示。线程块大小为(16, 16)，然后将 $ N*N $ 大小的矩阵均分为不同的线程块来执行加法运算。

SM 采用的是SIMT (Single-Instruction, Multiple-Thread，单指令多线程)架构，基本的执行单元是 线程束（wraps)，线程束包含 32 个线程，这些线程同时执行相同的指令，但是每个线程都包含自己的指令地址计数器和寄存器状态，也有自己独立的执行路径。

当线程块被划分到某个 SM 上时，它将进一步划分为多个线程束，因为这才是 SM 的基本执行单元，但是一个 SM 同时并发的线程束数是有限的。这是因为资源限制，SM 要为每个线程块分配共享内存，而也要为每个线程束中的线程分配独立的寄存器。所以 SM 的配置会影响其所支持的线程块和线程束并发数量。由于 SM 的基本执行单元是包含 32 个线程的线程束，所以 block 大小一般要设置为 32 的倍数。(16, 16)的二维 Block 是一个常用的配置，共 256 个线程。之前也曾提到过，每个 Block 的 Thread 个数最好是 128、256 或 512，这与 GPU 的硬件架构高度相关。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

// Kernel定义
__global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N])
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    if (i < N && j < N)
        C[i][j] = A[i][j] + B[i][j];
}
int main()
{
    ...
    // Kernel 线程配置
    dim3 threadsPerBlock(16, 16);
    dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
    // kernel调用
    MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
    ...
}

线程块中的线程数是有限制的，现代 GPUs 的线程块可支持的线程数可达 1024 个。有时候，我们要知道一个线程在 blcok 中的全局 ID，此时就必须还要知道 block 的组织结构，这是通过线程的内置变量 blockDim来获得。它获取线程块各个维度的大小。

  • 对于一个 2-dim 的 block $ (D_x, D_y) $ ，线程 $ (x, y) $ 的 ID 值为 $ (x + y * D_x) $
  • 对于一个 3-dim 的 block $ (D_x, D_y, D_z) $，线程 $(x, y, z)$ 的 ID 值为 $ (x + y _ D_z + z _ D_z * D_y) $

另外线程还有内置变量 gridDim，用于获得网格块各个维度的大小。

内存模型与管理

此外这里简单介绍一下 CUDA 的内存模型，如下图所示。可以看到，

  • 每个 Thread 有自己的私有本地内存（Local Memory）
  • 每个 Block 有包含共享内存（Shared Memory），可以被线程块中所有线程共享，其生命周期与线程块一致
  • 所有的 Thread 都可以访问全局内存（Global Memory）
  • 访问一些只读内存块：常量内存（Constant Memory）和纹理内存（Texture Memory）
  • L1 Cache，L2 Cache

{% gp group-layout %}

{% endgp %}

下面简单介绍一下 CUDA 编程中内存管理常用的 API。首先是在 Device 上分配内存的 cudaMalloc 、cudaFree 和 cudaMemcpy函数，分别对应 C 语言中的 malloc、free和 memcpy函数：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

// 在 Device 上申请一定字节大小的显存，其中 `devPtr` 是指向所分配内存的指针
cudaError_t cudaMalloc(void** devPtr, size_t size);

// 在 Device 上释放一定大小的现存， `devPtr` 是指向所释放内存的指针
cudaError_t cudaFree(void* devPtr);

// 负责 Host 和 Device 之间数据通信，src指向数据源，dst是目标区域，count是复制的字节数，kind控制复制的方向
// 这里的 kind 有四种类型：
// - cudaMemcpyHostToHost
// - cudaMemcpyHostToDevice
// - cudaMemcpyDeviceToHost
// - cudaMemcpyDeviceToDevice
cudaError_t cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind)

CUDA 编程实战

知道了 CUDA 编程基础，接下来我们以两个向量的加法为例，介绍如何利用 CUDA 编程来实现 GPU 加速计算。

CPU 向量加法：传统计算方法

我们首先来看利用 CPU 来计算向量加法该如何编程：

 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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#define N 10000000
#define MAX_ERR 1e-6

void vector_add(float *out, float *a, float *b, int n) {
    for(int i = 0; i < n; i++){
        out[i] = a[i] + b[i];
    }
}

int main(){
    float *a, *b, *out;

    // Allocate memory
    a   = (float*)malloc(sizeof(float) * N);
    b   = (float*)malloc(sizeof(float) * N);
    out = (float*)malloc(sizeof(float) * N);

    // Initialize array
    for(int i = 0; i < N; i++){
        a[i] = 1.0f;
        b[i] = 2.0f;
    }

    // Main function
    vector_add(out, a, b, N);

    // Verification
    for(int i = 0; i < N; i++){
        assert(fabs(out[i] - a[i] - b[i]) < MAX_ERR);
    }

    printf("out[0] = %f\n", out[0]);
    printf("PASSED\n");
}

GPU 向量加法：一个 Block 一个 Thread

我们将 CPU 的向量加法转换成 CUDA 程序，使用 GPU 来计算，下面这段代码演示了如何使用 CUDA 编程规范来编写程序。实际上仍然只是使用一个 core 来进行计算，不仅没有提高并行度，反而还增加了数据拷贝的成本，显然相比原来的计算是会更慢的，这里主要作为演示。

 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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <cuda.h>
#include <cuda_runtime.h>

#define N 10000000
#define MAX_ERR 1e-6

__global__ void vector_add(float *out, float *a, float *b, int n) {
    for(int i = 0; i < n; i ++){
        out[i] = a[i] + b[i];
    }
}

int main(){
    float *a, *b, *out;
    float *d_a, *d_b, *d_out;

    // Allocate host memory
    a   = (float*)malloc(sizeof(float) * N);
    b   = (float*)malloc(sizeof(float) * N);
    out = (float*)malloc(sizeof(float) * N);

    // Initialize host arrays
    for(int i = 0; i < N; i++){
        a[i] = 1.0f;
        b[i] = 2.0f;
    }

    // Allocate device memory
    cudaMalloc((void**)&d_a, sizeof(float) * N);
    cudaMalloc((void**)&d_b, sizeof(float) * N);
    cudaMalloc((void**)&d_out, sizeof(float) * N);

    // Transfer data from host to device memory
    cudaMemcpy(d_a, a, sizeof(float) * N, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, b, sizeof(float) * N, cudaMemcpyHostToDevice);

    // Executing kernel
    vector_add<<<1,1>>>(d_out, d_a, d_b, N);

    // Transfer data back to host memory
    cudaMemcpy(out, d_out, sizeof(float) * N, cudaMemcpyDeviceToHost);

    // Verification
    for(int i = 0; i < N; i++){
        assert(fabs(out[i] - a[i] - b[i]) < MAX_ERR);
    }

    printf("PASSED\n");

    // Deallocate device memory
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_out);

    // Deallocate host memory
    free(a);
    free(b);
    free(out);
}

GPU 向量加法：一个 Block 多个 Thread

为了提高并行度，我们设置一个 Block 多个 Thread 同时进行计算，如下图所示总共有 256 个 Thread，每个 Thread 负责处理 Vector 中的一部分。每一次迭代中，256 个 Thread 分别计算 Vector 的这 256 个数，然后在下一次迭代中每个 Thread 往后推进 256 个数，继续计算。

 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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <cuda.h>
#include <cuda_runtime.h>

#define N 10000000
#define MAX_ERR 1e-6

__global__ void vector_add(float *out, float *a, float *b, int n) {
    int index = threadIdx.x;
    int stride = blockDim.x;

    for(int i = index; i < n; i += stride){
        out[i] = a[i] + b[i];
    }
}

int main(){
    float *a, *b, *out;
    float *d_a, *d_b, *d_out;

    // Allocate host memory
    a   = (float*)malloc(sizeof(float) * N);
    b   = (float*)malloc(sizeof(float) * N);
    out = (float*)malloc(sizeof(float) * N);

    // Initialize host arrays
    for(int i = 0; i < N; i++){
        a[i] = 1.0f;
        b[i] = 2.0f;
    }

    // Allocate device memory
    cudaMalloc((void**)&d_a, sizeof(float) * N);
    cudaMalloc((void**)&d_b, sizeof(float) * N);
    cudaMalloc((void**)&d_out, sizeof(float) * N);

    // Transfer data from host to device memory
    cudaMemcpy(d_a, a, sizeof(float) * N, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, b, sizeof(float) * N, cudaMemcpyHostToDevice);

    // Executing kernel
    vector_add<<<1,256>>>(d_out, d_a, d_b, N);

    // Transfer data back to host memory
    cudaMemcpy(out, d_out, sizeof(float) * N, cudaMemcpyDeviceToHost);

    // Verification
    for(int i = 0; i < N; i++){
        assert(fabs(out[i] - a[i] - b[i]) < MAX_ERR);
    }

    printf("PASSED\n");

    // Deallocate device memory
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_out);

    // Deallocate host memory
    free(a);
    free(b);
    free(out);
}

相比 CPU 程序，这里的并行度显著提高，GPU 计算的时间也大大减小。

GPU 向量加法：多个 Block 多个 Thread

在上一个方案中，我们的 256 个 Thread 仍然需要计算多个数字，如果我们将并行度继续扩大，让每个 Thread 只需要计算 Vector 中的一个数，那么计算消耗时间将会更短。如下图所示，我们使用多个 Block 多个 Thread，其中每个 Block 还是 256 个 Thread，但是我们现在的 Grid 有多个 Block，Block 数字由 Vector 的长度除以 BlockSize 得到。

 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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <cuda.h>
#include <cuda_runtime.h>

#define N 10000000
#define MAX_ERR 1e-6

__global__ void vector_add(float *out, float *a, float *b, int n) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    // Handling arbitrary vector size
    if (tid < n){
        out[tid] = a[tid] + b[tid];
    }
}

int main(){
    float *a, *b, *out;
    float *d_a, *d_b, *d_out;

    // Allocate host memory
    a   = (float*)malloc(sizeof(float) * N);
    b   = (float*)malloc(sizeof(float) * N);
    out = (float*)malloc(sizeof(float) * N);

    // Initialize host arrays
    for(int i = 0; i < N; i++){
        a[i] = 1.0f;
        b[i] = 2.0f;
    }

    // Allocate device memory
    cudaMalloc((void**)&d_a, sizeof(float) * N);
    cudaMalloc((void**)&d_b, sizeof(float) * N);
    cudaMalloc((void**)&d_out, sizeof(float) * N);

    // Transfer data from host to device memory
    cudaMemcpy(d_a, a, sizeof(float) * N, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, b, sizeof(float) * N, cudaMemcpyHostToDevice);


    // Executing kernel
    int block_size = 256;
    int grid_size = ((N + block_size - 1) / block_size);
    vector_add<<<grid_size,block_size>>>(d_out, d_a, d_b, N);

    // Transfer data back to host memory
    cudaMemcpy(out, d_out, sizeof(float) * N, cudaMemcpyDeviceToHost);

    // Verification
    for(int i = 0; i < N; i++){
        assert(fabs(out[i] - a[i] - b[i]) < MAX_ERR);
    }

    printf("PASSED\n");

    // Deallocate device memory
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_out);

    // Deallocate host memory
    free(a);
    free(b);
    free(out);
}

GPU 向量加法：Unified Memory

在上面的实现中，我们需要单独在 Host 和 Device 上进行内存分配，并且要进行数据拷贝，这是很容易出错的。好在 CUDA 6.0 引入统一内存（Unified Memory）来避免这种麻烦，简单来说就是统一内存使用一个托管内存来共同管理 Host 和 Device 中的内存，并且自动在 Host 和 Device 中进行数据传输。CUDA 中使用 cudaMallocManaged 函数分配托管内存：

1

cudaError_t cudaMallocManaged(void **devPtr, size_t size, unsigned int flag=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
49
50
51
52

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <cuda.h>
#include <cuda_runtime.h>

#define N 10000000
#define MAX_ERR 1e-6

__global__ void vector_add(float *out, float *a, float *b, int n) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    // Handling arbitrary vector size
    if (tid < n){
        out[tid] = a[tid] + b[tid];
    }
}

int main(){
    // Allocate managed memory
    float *x, *y, *z;
    cudaMallocManaged((void**)&x, sizeof(float) * N);
    cudaMallocManaged((void**)&y, sizeof(float) * N);
    cudaMallocManaged((void**)&z, sizeof(float) * N);

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

    // Executing kernel
    int block_size = 256;
    int grid_size = ((N + block_size - 1) / block_size);
    vector_add<<<grid_size,block_size>>>(z, x, y, N);

    // 同步 Device 保证结果能正确访问
    cudaDeviceSynchronize();

    // Verification
    for(int i = 0; i < N; i++){
        assert(fabs(out[i] - a[i] - b[i]) < MAX_ERR);
    }

    printf("PASSED\n");

    // Deallocate managed memory
    cudaFree(x);
    cudaFree(y);
    cudaFree(z);
}

相比之前的代码，使用统一内存更简洁了，值得注意的是 kernel 执行是与 Host 异步的，由于托管内存自动进行数据传输，这里要用 cudaDeviceSynchronize() 函数保证 Device 和 Host 同步，这样后面才可以正确访问 kernel 计算的结果。

参考资料