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

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

本节简要介绍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线程切换瞬间完成，这是因为每个线程都有一份独占资源（例如寄存器），不需要保存、恢复现场，线程切换只是计算核心使用权的转移。

一块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：

在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

    【共64行，后面省略】

从1.1.1节的硬件架构图中已经看到，gpu中有多种内存：处于芯片外部的全局内存（global memory)，芯片内部的共享内存（shared meory）、寄存器（register）、纹理内存、常量内存、l1缓存、l2缓存。每种内存都有不同的特性，有不同的使用技巧。对开发cuda程序最重要的三种内存分别是寄存器、共享内存和全局内存。

如图1.25所示，每个线程都有自己专用的寄存器，从内核开始时，一旦拥有某个寄存器的使用权，就一直独占，直到内核结束才释放，从而线程之间无法通过寄存器交换数据。虽然有大量的寄存器，但也有大量的线程，平均下来每个线程只能分配到几十个至几百个寄存器，复杂程序仍然要控制线程消耗的寄存器数量。每个线程块都能分配一块共享内存，本块内的线程可以访问这块共享内存的任意位置，因此可以用共享内存来交换数据。一个线程块不能访问其他线程块的共享内存，因而线程块之间不能用共享内存交换数据。共享内存容量比寄存器要大，例如tesla p100的每个流式多处理器拥有64kb共享内存，每个线程块最多可以拥有32kb。所有的线程块、线程网格都能访问全局内存，只要不显式地释放或者程序结束，全局内存中的数据会一直存在，因此可以用于线程块之间、线程网格之间的数据交换。全局内存更大，以gb为单位。

不同内存的访问延时差别很大，寄存延时最小，共享内存次之，全局内存最大。对pascal之前的架构，全局内存与gpu芯片相互分离，通过板卡上的排线相连，访问延时达到几百个时钟周期。pascal架构中，全局内存与gpu芯片距离很近，延时应该有大幅减小，

具体数值还需要等待官方公布。

不同构件下的内存层级多少都有些变化，要想使cuda程序达到最好性能，必须做针对性优化。

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的差别会更大。因此，应尽量减少主机与设备间的数据传输量与传输次数。