Garden | GPU 与 CUDA

随着近年来深度学习的爆发，原来被用于图形渲染的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

尽管现在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…

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 程序的执行流程如下：

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

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

Streaming Multiprocessor

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

将 SM 放大，单个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。

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 工具，你需要通过配置环境变量来使用他们：

编译器：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次。

实际上，线程（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 %}

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 计算的结果。