《OpenACC并行编程实战》—— 1.3 CUDA C

本节书摘来自华章出版社《OpenACC并行编程实战》一 书中的第1章，第1.3节，作者何沧平，更多章节内容可以访问云栖社区“华章计算机”公众号查看。

1.3　CUDA C

本节简要介绍CUDA C编程的相关概念，使读者能够看懂OpenACC编译过程中出现的CUDA内置变量，理解并行线程的组织方式。如果读者已有CUDA编程经验，请跳过。
CPU用得好好的，为什么要费心费力地改写程序去到GPU上运行呢？只有一个理由：跑得更快。小幅的性能提升吸引力不够，必须有大幅提升才值得采购新设备、学习新工具、设计新算法。从图1.19可以看出，在双精度浮点峰值和内存带宽这两个关键指标上，GPU的性能都达到同时期主力型号CPU的5～7倍。如果利用得当，可以预期获得5～7的性能提升。以前只在CPU上运行，计算方法的数学理论和程序代码实现已经迭代发展多年，花很大力气才能提速10%～20%，提速50%已经很厉害了。简单粗暴地更换硬件设备就能立刻提速几倍，全世界的科学家、工程师一拥而上，GPU加速的应用遍地开花。注意，评价GPU应用性能的时候，至少要和2颗中高端CPU相对，并且两种代码都优化到最好。任何超过硬件潜能的加速结果都是有问题的。

那么问题来了。GPU的芯片面积与CPU差不多，价格也接近，为什么性能这么强悍呢？图1.20是CPU和GPU芯片的组成示意图，左边是一个单核超标量CPU，4个算术逻辑单元（ALU）承担着全部计算任务，却只占用一小部分芯片面积。“控制”是指分支预测、乱序执行等功能，占用芯片面积大而且很费电。服务器CPU通常有三级缓存，占用的芯片面积最大，有的型号甚至高达70%。ALU、控制、缓存都在CPU内部，大量内存条插在主板上，与CPU通过排线相连。GPU中绝大部分芯片面积都是计算核心（4行紧挨着的小方块，每行12个），带阴影的水平小块是控制单元，控制单元上面的水平条是缓存。

通用CPU对追踪链表这样拥有复杂逻辑控制的程序运行得很好，但大规模的科学与工程计算程序的流程控制都比较简单，CPU的长处难以施展。为了解释GPU如何获得极高的性能，需要先了解一下CPU中的控制、缓存、多线程的作用。
ALU承担最终的计算工作，越多越好。“控制”的目标是预取到正确的指令和数据以保证流水线不中断，挖掘指令流里的并行度，让尽量多的部件都在忙碌工作，从而提高性能。缓存的作用是为了填补CPU频率与内存条频率的差距、减小CPU与内存条之间数据延时。目前中高端CPU的频率在2.0～3.2GHz，而内存条的频率还处于1600MHz、1866MHz、2133MHz，内存条供应、承接数据的速度赶不上CPU处理数据的速度。由于ALU到主板上内存条的路径较长，延时高，而如果需要的数据已经在缓存中，那么就能有效降低延时，提高数据处理速度。缓存没有命中怎么办？只能到内存条上取，延时高。为了进一步降低延时，英特尔CPU有超线程功能，开启后，一个CPU物理核心就变成了两个逻辑核心，两个逻辑核心分时间片轮流占用物理核心资源。当然了，按时间片切换是有代价的：换出时要保留正在运行的程序的现场，换入时再恢复现场以便接着上次继续运行。在缓存命中率比较低的情况下，超线程功能能够提高性能。
GPU天生是为并行计算设计的：处理图像的大量像素，像素之间相互独立，可以同时计算，而且没有复杂的流程跳转控制。正如图1.19所示，GPU的大部分芯片面积都是计算核心，缓存和控制单元很小，那么它是怎么解决分支预测、乱序执行、数据供应速度、存取数据延时这些问题的呢？
GPU的设计目标是大批量的简单计算，没有复杂的跳转，因此直接取消分支预测、乱序执行等高级功能。更进一步，多个计算核心（例如32个）共用一个控制单元再次削减控制单元占用的芯片面积。这样做的效果就是：发射一条指令，例如加法，32个计算核心步调一致地做加法，只是每个计算核心操作不同的数据。如果只让第1个计算核心做加法，那么在第1个计算核心做加法运算的时候，剩余的计算核心空闲等待。这种情形下资源浪费，性能低下，要尽量避免。让大量计算核心空转的应用程序不适合GPU，用CPU计算性能更好。
计算核心与显存之间的频率差异如何填补？特别简单，降低计算核心的频率。考虑到芯片功耗与频率的平方近似成正比，降低频率不但能解决数据供应速度问题，而且能降低GPU的功耗，一举两得。从表1.1可以看出GPU产品的频率在562～875MHz，远低于主力CPU的2.0GHz～3.2GHz。
最重要是延时，GPU的缓存那么小，怎么解决访问显存的巨大延时呢？答案是多线程，每个计算核心分摊10个以上的线程。执行每条指令之前都要从就绪队列中挑选出一组线程，每组线程每次只执行一条指令，执行完毕立即到后面排队。如果恰巧碰上了延时较多的访存操作，那么该线程进入等待队列，访存操作完成后再转入就绪队列。只要线程足够多，计算核心总是在忙碌，隐藏了访存延时。有人立刻会问，这么频繁地切换线程、保存现场、恢复现场也需消耗不少时间吧，会不会得不偿失呢？实际上GPU线程切换瞬间完成，这是因为每个线程都有一份独占资源（例如寄存器），不需要保存、恢复现场，线程切换只是计算核心使用权的转移。

1.3.1　线程组织方式

一块GPU上有几千个核心，每个核心都能运行10个以上线程，可见线程数量庞大，需要按照一定结构组织起来，方便使用和管理。所有的线程合在一起称为一个网格（grid），网格再剖分成线程块（block），线程块包含若干线程。图1.21中的线程按照二维形式组织，网格包含2×3个线程块，每个线程块又包含3×4个线程。实际上，线程还可以按照一维、三维形式组织。

既然线程能够以不同的形式组织起来，那么每个线程都要有一个唯一的编号。为此CUDA C引入了一个新的数据类型dim3。dim3相当于一个结构体，3个成员分别为：
        unsigned int x;
        unsigned int y;
        unsigned int z;
dim3类型变量的3个成员的默认值都是1。网格尺寸用内置变量gridDim表示，gridDim.x、gridDim.y、gridDim.z分别表示x、y、z方向上的线程块数量；网格中每个线程块的编号用内置变量blockIdx表示，blockIdx.x、blockIdx.y、blockIdx.z分别表示当前线程块在x、y、z方向上的编号，从0开始编号；线程块的尺寸用内置变量blockDim表示，blockDim.x、blockDim.y、blockDim.z分别表示当前线程块在x、y、z方向上拥有的线程数量；任意一个线程块内的线程编号用内置变量threadIdx来表示，threadIdx.x、threadIdx.y、threadIdx.z分别表示当前线程在x、y、z方向上的编号，从0开始编号。以图1.21中的网格、线程块（1，1）、线程块（1，2）为例，这些内置变量的值如表1.3：

1.3.2　运行过程

在GPU编程话语体系里，称CPU为主机，称GPU为设备。图1.22演示了CUDA C程序的执行过程：在带有设备的计算机上，与C语言程序一样，从主机开始执行，主机上执行串行代码，并为设备上的并行计算做准备，包括数据初始化、开辟设备内存、将数据复制到设备内存中。准备工作完成之后，在主机上以特殊形式调用一个在设备上执行的函数（称为内核，调用时比C函数多了一对三尖号），然后设备执行内核中的并行代码。内核代码执行完以后，控制权交还主机，主机从设备上取回内核的并行计算结果，程序继续向下执行。图1.22中只画出一个内核，实际上一个CUDA程序可以包含多个内核。

下面以实际例子演示CUDA C代码的编写方法和执行过程。两个长度为N的向量a和b对应元素相加，将结果存入向量c。从图1.23可以看出，N个加法操作之间没有依赖关系，可以并行计算。实现代码见例1.1。

例1.1中第10行定义3个主机向量a、b、c，第11行定义3个指针用于存放设备向量，第12～14行为3个设备向量分配设备内存空间。第15～19行的循环为主机向量a、b赋初值，第20～21行使用内置函数cudaMemcpy将主机向量a和b中的元素值复制到设备向量a_d和b_d之中，即从主机内存复制到设备内存。第22行定义了2个dim3变量block和grid。block用于指定每个线程块的形状：一维，x方向长度为32；grid用于指定线程网格的形状：一维，x方向的尺寸用block.x和N计算出来，以适应N不能被32整除的情形。至此，准备工作完毕。
第24行从主机调用内核add，三尖号＜＜＜＞＞＞里的参数称为执行配置，第1个参数指定线程网格的形状，第2个参数指定线程块的形状，紧跟着的圆括号里面是和C函数一样的实参。执行配置参数要求启动2个线程块共64个线程来执行内核add。内核add在设备上运行，它将设备向量a_d和b_d并行相加，结果存入设备向量c_d。内核add的定义在第4～7行，第4行上的修饰符__global__表示该函数需要在主机上调用且在设备上执行。第5行计算线程的全局编号，N为64，每个线程块有32个线程，因此网格中有2个线程块。在每个线程块中，线程的本地编号threadIdx.x分别是0，1，2，…，31，blockDim.x的值为32，所以执行内核的64个线程的tid分别为0，1，2，...，63，见图1.24。第6行也被64个线程同时执行，每个线程执行1次加法，共同完成两个向量的对应相加。

第25行将设备上的计算结果复制回主机内存，即把向量c_d的元素值复制到向量c中。第27～28行输出计算结果以便检验正确性，可以预见是64行1＋2＝3。第29～31行释放设备内存。
在已经部署CUDA C开发工具的Linux环境上编译、运行：
        $ nvcc -o addvec.exe addvec.cu
        $ ./addvec.exe
         1 + 2 = 3
         1 + 2 = 3
         1 + 2 = 3
    【共64行，后面省略】

1.3.3　内存层级

从1.1.1节的硬件架构图中已经看到，GPU中有多种内存：处于芯片外部的全局内存（Global Memory)，芯片内部的共享内存（Shared Meory）、寄存器（Register）、纹理内存、常量内存、L1缓存、L2缓存。每种内存都有不同的特性，有不同的使用技巧。对开发CUDA程序最重要的三种内存分别是寄存器、共享内存和全局内存。
如图1.25所示，每个线程都有自己专用的寄存器，从内核开始时，一旦拥有某个寄存器的使用权，就一直独占，直到内核结束才释放，从而线程之间无法通过寄存器交换数据。虽然有大量的寄存器，但也有大量的线程，平均下来每个线程只能分配到几十个至几百个寄存器，复杂程序仍然要控制线程消耗的寄存器数量。每个线程块都能分配一块共享内存，本块内的线程可以访问这块共享内存的任意位置，因此可以用共享内存来交换数据。一个线程块不能访问其他线程块的共享内存，因而线程块之间不能用共享内存交换数据。共享内存容量比寄存器要大，例如Tesla P100的每个流式多处理器拥有64KB共享内存，每个线程块最多可以拥有32KB。所有的线程块、线程网格都能访问全局内存，只要不显式地释放或者程序结束，全局内存中的数据会一直存在，因此可以用于线程块之间、线程网格之间的数据交换。全局内存更大，以GB为单位。
不同内存的访问延时差别很大，寄存延时最小，共享内存次之，全局内存最大。对Pascal之前的架构，全局内存与GPU芯片相互分离，通过板卡上的排线相连，访问延时达到几百个时钟周期。Pascal架构中，全局内存与GPU芯片距离很近，延时应该有大幅减小，

具体数值还需要等待官方公布。
不同构件下的内存层级多少都有些变化，要想使CUDA程序达到最好性能，必须做针对性优化。

1.3.4　性能优化技术

CUDA程序编写容易，调优不易。程序员能够掌控很多事情，包括但不限于分配全局内存：全局内存中的数据对齐、维数，为每个线程块分配的共享内存大小，将哪些数据以什么样的组织方式放入共享内存，哪些数据放入纹理内存，哪些数据放入常量内存，线程网格如何划分，线程块是一维、二维还是三维，线程块每个维度的大小是多少，线程与数据元素的对应关系，不同线程访问的数据是否有冲突，不同线程同时访问的数据是否会走相同的通道；单个内核是否能够用满资源，如何同时运行多个内核以提高设备利用率，有几个数据复制引擎，如何安排异步队列来重叠数据的来往传输，如何重叠数据传输与计算，如何填补PCIe带宽与全局内存带宽之间的差异，数据复制操作是否需要锚定主机内存；计算密度够不够大，计算核心要等待数据多久，一个Warp内的线程的流程分支有多少，多少个线程才能隐藏延时；GPU上的算术指令与CPU上对应指令的差异，双精度操作、单精度操作、半精度操作、三角函数等特殊操作的计算资源分配。
管事多，操心就多。每个问题都有相应的优化方法和一定的约束条件，具体技巧请参考英伟达官方文档《CUDA C BEST PRACTICES GUIDE》。需要注意，不同架构下的优化技术会有一些差别。
影响最大的优化技巧是主机与设备间的数据传输。从图1.4可以看出，设备与主机通过PCIe×16通道相连，在采用2016年发布的最新CPU的服务器上，PCIe 3.0×16的理论带宽为16GB/s，与表1.1中几百GB/s的显存（全局内存）带宽差别可达30倍，与Tesla P100的差别会更大。因此，应尽量减少主机与设备间的数据传输量与传输次数。