第1章 并行编程概览

对绝大多数人而言，编程语言只是一个工具，讲究简单高效。科学家的主要精力应该用在科研创新活动上，编程工作仅仅是用来验证创新的理论，编程水平再高也不可能获得诺贝尔奖。对学生和企业程序员而言，技术无穷尽，永远学不完，不用即忘，应该认清技术发展方向，学习有前途的技术，不浪费青春年华。

OpenACC语言专为超级计算机设计，因此读者需要了解超级计算机的技术演进方向，特别是主流加速器的体系架构、编程模型，看清OpenACC的应用场景，有的放矢。普通读者虽然不会用到大型机群，但小型机群甚至单台服务器、普通显卡的计算模式都是相同的。

最近几年的著名超级计算机（见附录A）均采用加速器作为主要计算部件，可预见未来几年的上层应用仍将围绕加速器展开。

1.1 加速器产品

超级计算机的加速器历史上有多种，本节只介绍当前流行的两种：英伟达GPU和英特尔融核处理器。加速器的物理形态是PCIe板卡，样子大致如图1.1所示，图1.2是拆掉外壳后的样子，正中央的是GPU芯片，芯片周围的小黑块是显存颗粒，金黄色的边缘处是与PCIe连接的金手指，通过PCIe插槽与CPU相连。图1.3中的机架式服务器左下部装有4块GPU卡，图1.4是服务器的主板俯视图，箭头处就是4个PCIe插槽。

图1.1 英伟达Tesla GPU

图1.2 英伟达GPU内部构造

图1.3 GPU服务器

图1.4 GPU服务器的主板

从逻辑关系上来看，目前市场在售的英特尔服务器CPU中已经集成了内存控制器和PCIe通道。2颗CPU通过1个或2个QPI接口连接，CPU通过内置的内存控制器连接DDR3或DDR4内存条，GPU加速器通过PCIe总线与CPU相连（图1.5）。2016年6月主流的英特尔至强E5-2600 v3和E5-2600 v4 CPU每颗拥有40个PCIe Lan，而每块GPU的接口需要16个PCIe Lan，因此每颗CPU上最多全速挂载2块GPU。有些服务器的PCIe ×16插槽上只安排×8，甚至×4的信息速率，可以挂载更多的GPU。

图1.5 加速器在服务器上的逻辑位置

英特尔融核处理器在服务器上的位置与GPU一样，不再赘述。

1.1.1 英伟达GPU

通用计算GPU是英伟达公司发明的，每隔2～3年就会升级硬件架构。通用计算GPU的第一代架构代号是G80，每8个核心封装在一起称为一个流式多处理器（Streaming Multiprocessor, SM）,2个流多处理器共用一块L1缓存，所有的核心共用L2缓存，任意核心都能访问芯片外部的GPU显存（图1.6中的FB）。从图1.6中可以看出，G80架构GPU产品最多可以有128个核心。

图1.6 通用计算GPU的G80架构

接下来的架构代号是GT200，然后是Fermi。Fermi架构（图1.7）最多包含14个流式多处理器，每个流式多处理器包含32个核心。6个DRAM接口，1个PCIe接口（图1.16中的Host Interface）。

图1.7 Fermi架构图

流式多处理器内部组件也相当多（图1.8）:1个指令缓存（Instruction Cache）,2组Warp调度器（Warp Scheduler）、分发单元（Dispatch Unit），一堆寄存器；最重要的是32个核心，每个核心拥有一个单精度浮点单元和一个整数单元；16个Load/Store单元；4个特殊功能单元（Special Function Unit），负责计算双精度浮点数、三角函数、超越函数、倒数、平方根等。共享内存+L1缓存共64KB，可以灵活配置。

图1.8 Fermi架构的流式多处理器

接下来是Kepler架构（图1.9），这一代架构最大可以包含15个流式多处理器（Streamong Multiprocessor extreme, SMX)和6个64位内存控制器。不同的产品型号将使用不同的配置，可以包含13或14个SMX，多种参数都有升级和更改。Kepler架构的主要设计目标是提高用电效率，台积电的28nm制造工艺在降低功耗方面起着重要的作用。Kepler架构还提高了双精度计算能力。

图1.9 Kepler架构

流式多处理器SMX（图1.10）包含192个单精度核心、64个双精度单元、32个特殊功能单元（SFU）和32个加载/存储单元（LD/ST）、3个Warp调度器、6个分发器、一大堆寄存器。每个核心由1个浮点计算单元和1个整数算术逻辑单元组成，支持融加（FMA）运算。每个SMX拥有64KB的片上存储器，可配置为48KB的共享存储器和16KB的L1缓存，或配置为16KB的共享内存和48KB的L1缓存。

图1.10 kepler架构中的流式多处理SMX

这里简要介绍市面的主力GPU产品型号，见表1.1。

表1.1 2014～2016年市场主力GPU产品规格

2016年4月，英伟达在GPU技术会议（GPU Technology Conference）上发布了新一代Pascal架构和旗舰产品Tesla P100，引入一些新特性。

❑ 极致性能：为高性能计算、深度学习等计算领域设计。双精度浮点峰值5.3 Tflops，单精度浮点峰值10.6 Tflops，专为深度学习设计的半精度浮点峰值达到惊人的21.2 Tflops。按双精度峰值对比，Tesla P100是同期主力高端CPU英特尔至强E5-2680 v4的10倍。按照深度学习应用性能对比，Tesla P100是E5-2680 v4的20倍。

❑ NVLink：为应用扩展性全新设计的高速、高带宽互连协议。

❑ HBM2：快速、大容量、高效片上堆叠式内存。

❑ 统一内存：用统一的代码来管理主机CPU内存和GPU内存，方便开发代码。

Pascal最强劲的是GP100硬件架构（图1.11），包含60个Pascal流式多处理器和8个512位内存控制器（共4096位）。每个流式多处理器拥有64个CUDA核心和4个纹理单元。GP100共有3840个单精度核和240个纹理单元。每个内存控制器附带512KB的L2缓存，每个HBM2堆叠内存由一对内存控制器管理。L2缓存共计4096KB。Tesla P100共配置有56个流式多处理器。

图1.11 GP100架构全景图

GP100的第6代流式多处理器架构提升了CUDA核心的利用率和能源效率，可以运行在更高的频率上。每个流多处理器包含64个单精度（FP32）CUDA核心，分成两部分，每部分32个；含有1个指令缓冲器，1个Warp调度器，2个分发单元。与前面的Kepler、Maxwell流式多处理器相比，每个核心分到的寄存器数量增多，从而可以运行更多的线程。由于增加了流式多处理器，共享内存总量也相应增加，合并带宽也翻番。每个流式多处理器中更多的共享内存、寄存器、Warp，使代码可以更高效地执行。Pascal架构中的FP32 CUDA核心新增一项功能：既能处理16位精度又能处理32位精度的指令和数据，FP16操作的吞吐率是FP32的2倍。

如图1.12所示，每个流式多处理器拥有32个双精度浮点单元，是单精度核心的一半。共享内存独占64KB，不再与L1缓存相互调配。纹理和L1缓存共用一块缓存，可以灵活配置。

图1.12 Pascal GP100流式多处理器

Tesla P100的高带宽内存（High Bandwidth Memory 2, HBM2）跟前面几代相比有显著变化：从GPU芯片封装的外部移到了内部（图1.13），安装在同一块基板上，距离更近，传输更快。由于在竖直方向进行堆叠，占用的芯片面积更小，16GB的容量原来需要十几个GDDR5颗粒，现在只需要4个HBM2颗粒。

HBM2内存带来的另一个好处是原生支持错误校验码（Error Correcting Code, ECC）功能。有些应用会用到很多块GPU，计算时间也很长，对数据差错十分敏感，一旦中间步骤有微小误差，后续计算就会将误差迅速放大，污染最终结果。ECC技术能够检测并纠正1位差错。

图1.13 Tesla P100正视图

先前使用的GDDR5显存内部不提供ECC功能，仅能侦测GDDR5总线上的错误，内存控制器中的错误或DRAM自身的错误都不能侦测。Kepler架构GPU的ECC功能实现方法是划出6.25%的容量专门用来存放校验数据，从而内存带宽将损失12%～15%。

由于HBM2自带ECC功能，Tesla P100就不再损失内存容量和内存带宽。Tesla P100的寄存器、共享内存、L1缓存、L2缓存和HBM2内存都能侦测2位错误并纠正1位错误。

1.1.2 英特尔至强融核处理器

2012年，英特尔推出至强Phi融核处理器，代号Knights Corner，融核处理器又称为MIC（Many Integrated Core）。物理形态上跟GPU一样是PCIe板卡（图1.14），分为被动散热和主动散热（带风扇），具体分为6个型号（表1.2），市面上常用的型号是5110P和7120。

图1.14 Intel Xeon Phi协处理器

表1.2 英特尔融核处理器技术规格

此代架构MIC拥有最多63个核心（图1.15），每个核心有512KB的L2缓存，所有核心使用双向环互连。每个核心有一个512位向量处理单元，支持4个硬件线程，线程并发隐藏延时。

图1.15 Knights Corner架构

MIC上编程模式有4种（图1.16）：主机模式（Multicore Only）、卸载模式（Multicore Hosted with Many-core Offload）、对称模式（Symmetric）、原生模式（Many-Core Only）。主机模式中只使用主机CPU，不在MIC上运行代码；卸载模式跟GPU的模式类似，将部分代码卸载到MIC上运行；对称模式将CPU核心和MIC核心一视同仁，但是两种核心的运算能力不一样，极易导致负载不均；原生模式只使用MIC运行程序，CPU空闲。实际上，常用的还是卸载模式和原生模式。

图1.16 MIC的运行模式

此代MIC的理论性能是同时期主流CPU的3.2～3.45倍，实测性能2.2～2.9倍。因此，使用MIC编写运行程序不要期待过高的加速比。

第二代MIC的代号为Knights Landing，开发版于2016年4月开始发货，已经不是PCIe板卡，而是可以独立运行的CPU。

Knights Landing由36片组成（图1.17），每片包含2个CPU核心，每个核心包含2个VPU（Vector Processing Units，向量处理器单元）。物理形态上，不再是协处理器，而是真正的中央处理器，与通常的CPU一样使用。每个CPU核心可以运行4个线程，芯片的乱序执行也有实质性提升，每个线程的性能提高约3倍，整个处理器的理论峰值约为3Tflops。芯片上有16GB MCDRAM（Multi-Channel DRAM）内存，6个内存通道支持高达384GB DDR4外部内存。与Knights Corner中的环状拓扑结构不同，新款处理器中使用了网状的拓扑结构。尽管图1.17显示为6×7，但是由于封装方面的一些问题，布局上可能会调整为4×9。

图1.17 英特尔至强融核处理器Knights Landing

与2016年2月发布的英特尔至强系列CPU的主力型号E5-2680 v4（主频2.4GHz,14核心，理论峰值537.6Gflops）相比，Knights Landing的理论峰值提高至约5.6倍。

1.2 并行编程语言

在并行计算发展史上出现过多种并行编程语言，至今仍在使用的只剩几种，它们各有特色。

（1）Pthreads

20世纪70年代，贝尔实验室发明了UNIX，并于20世纪80年代向美国各大高校分发V7版的源码以做研究。加利福尼亚大学伯克利分校在V7的基础上开发了BSD UNIX。后来很多商业厂家意识到UNIX的价值也纷纷以贝尔实验室的System V或BSD为基础来开发自己的UNIX，较著名的有Sun OS、AIX、VMS。随着操作系统的增多，应用程序的适配性工作越来越繁重。为了提高UNIX环境下应用程序的可迁移性，电气和电子工程师协会（Institute of Electrical and Electronics Engineers, IEEE）设计了POSIX标准。然而，POSIX并不局限于UNIX，许多其他的操作系统也支持POSIX标准。POSIX.1已经被国际标准化组织所接受，POSIX已发展成为一个非常庞大的标准族，一直处在发展之中。

POSIX线程（POSIX threads, Pthreads），是线程的POSIX标准。该标准定义了创建和操纵线程的一整套接口。在类UNIX操作系统（UNIX、Linux、Mac OS X等）中，都使用Pthreads作为操作系统的线程。Pthreads用来开发与操作系统紧密相关的应用程序，管理粒度很细，例如线程的创建与销毁、线程锁、线程属性、线程优先级、线程间通信等琐碎操作均需要程序员安排。对科学与工程类计算程序来说，程序员的精力应集中在业务模型和代码算法上，不应浪费在底层代码细节上。

虽然Pthreads不适合编写高性能计算程序，但它多线程并发的设计理念启发了其他并行语言。

（2）OpenMP

OpenMP是由一些大型IT厂商和一些学术机构组成的非盈利组织，官网是www.openmp.org。永久成员包括AMD、CAPS-Entreprise、Convey Computer、Cray、Fujitsu、HP、IBM、Intel、NEC、NVIDIA、Oracle Corporation、Red Hat、ST Microelectronics、Texas Instruments；正式成员是对OpenMP标准感兴趣，但不生产销售相关产品的组织，例如ANL、ASC/LLNL、BSC、cOMPunity、EPCC、LANL、NASA、ORNL、RWTH Aachen University、SNL-Sandia National Lab、Texas Advanced Computing Center、University of Houston。

计算热点都是在循环上，OpenMP的并行化思路是将循环的迭代步分摊到多个线程上，每个线程只承担一部分计算任务，循环运行的墙上时间（从开始到结束的自然流逝时间）自然也就减少了。分割方法是在循环上面添加一些预处理标记（图1.18），编译器识别到这些标记以后，将关联的循环翻译成并行代码，然后再与剩余的串行代码合并起来编译、链接成可执行程序。

图1.18 OpenMP的并行化模式

使用OpenMP编译并行程序时，程序员需要先保证串行代码正确，找出热点循环，然后在循环上添加OpenMP预处理标记。打开编译器的openmp选项可以编译成并行版本，否则编译器将忽略预处理器标记，仍然编译成串行版本。既不破坏原有代码，开发速度又快，省时省力。

（3）CUDA

CUDA（Compute Unified Device Architecture）是英伟达公司设计的GPU并行编程语言，一经推出就引发了GPU通用计算研究热潮。CUDA是闭源的，只能运行在英伟达的产品上。CUDA C/C++是对C/C++语言的扩展，添加了一些数据类型、库函数，并定义一种新的函数调用形式。CUDA起初支持C和少量C++特性，后来逐渐提高对C++的支持度。从CUDA 3.0开始与PGI合作支持Fortran。CUDA语言还可以细分为CUDA C/C++、CUDA Fortran，本书成稿时的CUDA C/C++最新版本是7.5,8.0版本即将发布。CUDA Fortran没有版本号，只是随着PGI编译器的升级而增加新特性。1.3节会介绍CUDA C/C++的编程模型和一些示例代码，CUDA Fortran的详细情况可以参考官网https://developer.nvidia.com/cuda-fortran，网络上有英伟达工程师撰写的图书《Best Practices for Efficient CUDA Fortran Programming》和中文翻译版《CUDA-Fortran高效编程实践》，此处不展开介绍。

（4）OpenCL

OpenCL（Open Computing Language，开放计算语言）是一个面向异构系统并行编程的免费标准，支持多种多样的设备，包括但不限于CPU、GPU、数字信号处理器（DSP）。OpenCL的优势是一套代码多处运行，只要为新的设备重新编译代码就可以运行，移植方便。

OpenCL由苹果公司首先提出，随后Khronos Group成立相关工作组，以苹果草案为基础，联合业界各大企业共同完成了标准制定工作，工作组的成员来自各行各业，且都是各自领域的领导者，成员名单请参见官网www.khronos.org/opencl/。

（5）OpenACC

这里不多说，后面会全面、详细讲解。

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相对，并且两种代码都优化到最好。任何超过硬件潜能的加速结果都是有问题的。

图1.19 同时期主力CPU与GPU的性能对比

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

图1.20 CPU（左）和GPU（右）的芯片面积占用情况

通用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个线程。实际上，线程还可以按照一维、三维形式组织。

图1.21 线程网格与线程块

既然线程能够以不同的形式组织起来，那么每个线程都要有一个唯一的编号。为此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 内置变量的取值

1.3.2 运行过程

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

图1.22 CUDA程序运行过程

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

图1.23 两个向量的对应元素相加

【例1.1】addvec.cu：向量并行相加。

      1  #include<stdio.h>
      2  #define N 64
      3
      4  __global__ void add( int ＊a_d, int ＊b_d, int ＊c_d ) {
      5      int tid = blockIdx.x ＊ blockDim.x + threadIdx.x;
      6      if (tid < N) c_d[tid] = a_d[tid] + b_d[tid];
      7  }
      8  int main()
      9  {
      10    int a[N], b[N], c[N];
      11    int ＊a_d, ＊b_d, ＊c_d;
      12    cudaMalloc((void＊＊)&a_d, N ＊ sizeof(int));
      13    cudaMalloc((void＊＊)&b_d, N ＊ sizeof(int));
      14    cudaMalloc((void＊＊)&c_d, N ＊ sizeof(int));
      15    for(int i=0; i<N; i++)
      16    {
      17      a[i] = 1;
      18      b[i] = 2;
      19    }
      20    cudaMemcpy(a_d, a, N＊sizeof(int), cudaMemcpyHostToDevice);
      21    cudaMemcpy(b_d, b, N＊sizeof(int), cudaMemcpyHostToDevice);
      22    dim3 block(32,1,1), grid;
      23    grid.x = (N+block.x-1)/block.x;
      24    add<<<grid, block>>>(a_d, b_d, c_d);
      25    cudaMemcpy(c, c_d, N＊sizeof(int), cudaMemcpyDeviceToHost);
      26
      27    for(int i=0; i<N; i++)
      28      printf("%2d +%2d =%2d\n", a[i], b[i], c[i]);
      29    cudaFree( a_d );
      30    cudaFree( b_d );
      31    cudaFree( c_d );
      32  return 0;
      33  }

例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次加法，共同完成两个向量的对应相加。

图1.24 线程在线程块内的编号和全局编号

第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芯片距离很近，延时应该有大幅减小，具体数值还需要等待官方公布。

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