背景

    近年来，随着人工智能、高性能数据分析和金融分析等计算密集型领域的兴起，传统通用计算已经无法满足对计算性能的需求，异构计算越来越引起学术界和产业界的重视。

    异构计算是指采用不同类型的指令集和体系架构的计算单元组成系统的计算方式。相比传统CPU，异构计算可以实现更高的效率和更低的延迟。目前的异构计算引擎主要有图形处理器（GPU，Graphics Processing Unit）、现场可编程门阵列（FPGA，Field Programming Gate Array）、专用集成电路（ASIC）等。

    当前的通用CPU设计得已经很复杂，配有几十个核心，运行频率高达几GHz，每个核心有自己的独立缓存。通常CPU已具备一级、二级、三级缓存。而GPU是目前科研领域比较常用的硬件计算工具。GPU的计算核心数通常是CPU的上百倍，运行频率尽管比CPU的低，但是核心数量多，整体性能好。所以，GPU比较适合计算密集型应用，比如视频处理、人工智能等，现在传统的科学计算、工程计算等也开始越来越适合在GPU上运行。相比来说，CPU的缺点就是太通用了，数据读写、计算、逻辑等各种功能都得照顾，反而影响了计算性能。

并行计算

    在过去的几十年间，人们对并行计算产生了越来越多的兴趣。并行计算的主要目标是提高运算速度。

    从纯粹的计算视角来看，并行计算可以被定义为计算的一种形式，在这种形式下，计算机可以同时进行许多运算，计算原则是一个大的问题往往可以被划分为很多可以同时解决的小问题。并行计算其实设计到两个不同的技术领域：

• 计算机架构（硬件）

• 并行程序设计（软件）

    计算机架构关注的是在结构级别上支持并行性，而并行编程设计关注的是充分使用计算机架构的计算能力来并发地解决问题。为了在软件中实现并行执行，硬件必须提供一个支持并行执行多进程或多线程的平台。
    大多数现代处理器都应用了哈佛体系结构（Harvard architecture），如下图所示，它主要由 3 个部分组成：

• 内存（指令内存，数据内存）

• 中央处理单元（控制单元和算数逻辑单元）

• 输入、输出接口

    从程序员的角度来看，一个程序应包含两个基本的组成部分：指令和数据。当一个计算问题被划分成许多小的计算单元后，每个计算单元都是一个任务。在一个任务中，单独的指令负责处理输入和调用一个函数并产生输出。当一个指令处理前一个指令产生的数据时，就有了数据相关性的概念。因此，你可以区分任何两个任务之间的依赖关系，如果一个任务处理的是另一个任务的输出，那么它们就是相关的，否则就是独立的。

    在并行算法的实现中，分析数据的相关性是最基本的内容，因为相关性是限制并行性的一个主要因素，而且在现代编程环境下，为了提高应用程序的运行速度，理解这些是很有必要的。在大多数情况下，具有依赖关系的任务之间的独立的关系链为并行化提供了很好的机会。

并行性

    如今，并行性的应用非常广泛，在编程领域，并行编程设计正在成为主流。多层次的并行性设计是架构设计的驱动力。在应用程序中有两种基本的并行类型。

• 任务并行

• 数据并行

    通常，数据是在一维空间中存储的。即便是多维逻辑数据，仍然要被映射到一维物理地址空间中。如何在线程中分配数据不仅与数据的物理储存方式密切相关，并且与每个线程的执行次序也有很大关系。组织线程的方式对程序的性能有很大的影响。

    SISD指的是传统计算机：一种串行架构。在这种计算机上只有一个核心。在任何时间点上只有一个指令流在处理一个数据流。

    SIMD是一种并行架构类型。在这种计算机上有多个核心。在任何时间点上所有的核心只有一个指令流处理不同的数据流。向量机是一种典型的SIMD类型的计算机，现在大多数计算机都采用了SIMD架构。SIMD最大的优势或许就是，在CPU上编写代码时，程序员可以继续按串行逻辑思考但对并行数据操作实现并行加速，而其他细节则由编译器来负责。

    MISD类架构比较少见，在这种架构中，每个核心通过使用多个指令流处理同一个数据流。

    MIMD是一种并行架构，在这种架构中，多个核心使用多个指令流来异步处理多个数据流，从而实现空间上的并行性。许多MIMD架构还包括SIMD执行的子组件。

    为了实现以下目的，在架构层次上已经取得了许多进展。

• 降低延迟

• 提高带宽

• 提高吞吐量

    延迟是一个操作从开始到完成所需要的时间，常用微秒来表示。带宽是单位时间内可处理的数据量，通常表示为MB/s或GB/s。吞吐量是单位时间内成功处理的运算数量，通常表示为gflops（即每秒十亿次的浮点运算数量），特别是在重点使用浮点计算的科学计算领域经常用到。延迟用来衡量完成一次操作的时间，而吞吐量用来衡量在给定的单位时间内处理的操作量。

    计算机架构也能根据内存组织方式进行进一步划分，一般可以分成下面两种类型。

• 分布式内存的多节点系统

• 共享内存的多处理器系统

    在多节点系统中，大型计算引擎是由许多网络连接的处理器构成的。每个处理器有自己的本地内存，而且处理器之间可以通过网络进行通信。下图所示为一个典型的分布式内存的多节点系统，这种系统常被称作集群。

    多处理器架构的大小通常是从双处理器到几十个或几百个处理器之间。这些处理器要么是与同一个物理内存相关联（如下图所示），要么共用一个低延迟的链路（如PCI-Express或PCIe）。尽管共享内存意味着共享地址空间，但并不意味着它就是一个独立的物理内存。这样的多处理器不仅包括由多个核心组成的单片机系统，即所谓的多核系统，而且还包括由多个芯片组成的计算机系统，其中每一个芯片都可能是多核的。目前，多核架构已经永久地取代了单核架构。

    “众核”（many-core）通常是指有很多核心（几十或几百个）的多核架构。近年来，计算机架构正在从多核转向众核。

    GPU代表了一种众核架构，几乎包括了前文描述的所有并行结构：多线程、MIMD（多指令多数据）、SIMD（单指令多数据），以及指令级并行。NVIDIA公司称这种架构为SIMT（单指令多线程）。

    GPU和CPU的来源并不相同。历史上，GPU是图形加速器。直到最近，GPU才演化成一个强大的、多用途的、完全可编程的，以及任务和数据并行的处理器，它非常适合解决大规模的并行计算问题。

GPU核心和CPU核心

尽管可以使用多核和众核来区分CPU和GPU的架构，但这两种核心是完全不同的。CPU核心比较重，用来处理非常复杂的控制逻辑，以优化串行程序执行。GPU核心较轻，用于优化具有简单控制逻辑的数据并行任务，注重并行程序的吞吐量。

GPGPU

    通用图形处理器（GPGPU，General Purpose Graphics Processing Unit）最早由NVIDIA公司的Mark J. Harris于2002年提出。基于图形渲染管线的流水线特征，GPU本质上是一个可同时处理多个计算任务的硬件加速器。由于GPU中包含了大量的计算资源，Mark J. Harris自2002年就开始尝试在GPU上做通用并行计算方面的研究。在此阶段，由于架构及编程平台的限制，研究人员采用将目标计算算法转换为图形运算算法的方式，使用GPU来实现通用并行计算需求。

    NVIDIA公司提出Tesla统一渲染架构以及CUDA（Compute Unified Device Architecture，计算统一设备架构）编程模型后，NVIDIA公司的GPU开始了对通用并行计算的全面支持。在CUDA提出近两年之后，开放计算语言标准OpenCL 1.0发布，这标志着利用GPU进行通用并行计算已基本成熟。目前市场上应用甚广的GPU芯片除了完成高质量的图形渲染之外，通用并行计算也已经成为一个主流应用。GPGPU在各个方面得到了不同GPU厂家为GPU通用计算提供的编程模型与平台，如CUDA和OpenCL，这些编程模型在C/C++基础之上做了面向大规模通用并行计算的语法扩展，为程序员提供了更好的、面向GPU的编程接口。

    GPGPU 通常由成百上千个架构相对简易的基本运算单元组成。在这些基本运算单元中，一般不提供复杂的诸如分支预测、寄存器重命名、乱序执行等处理器设计技术来提高单个处理单元性能，而是采用极简的流水线进行设计。每个基本运算单元可同时执行一至多个线程，并由GPGPU中相应的调度器控制。GPGPU作为一个通用的众核处理器，凭借着丰富的高性能计算资源以及高带宽的数据传输能力在通用计算领域占据了重要的席位。虽然各个GPGPU厂商的芯片架构各不相同，但几乎都是采用众核处理器阵列架构，在一个GPU芯片中包含成百上千个处理核心，以获得更高的计算性能和更大的数据带宽。

    GPU 中执行的线程对应的程序通常成为内核（kernel），这与操作系统中的内核是完全不同的两个概念。除此之外，GPU中执行的线程与CPU或者操作系统中定义的线程也有所区别，GPU中的线程相对而言更为简单，所包含的内容也更为简洁。在GPU众核架构中，多个处理核心通常被组织成一个线程组调度执行单位，线程以组的方式被调度在执行单元中执行，如NVIDIA的流多处理器、AMD的SIMD执行单元。同一个线程组中的线程执行相同的程序指令，并以同步的方式执行，每个线程处理不同的数据，实现数据级并行处理。不同GPU架构对线程组的命名各不一样，如NVIDIA将线程组称为warp，AMD将线程组称为wavefront。线程组中包含的线程数量各不相同，从4个到128个不等。除此以为，线程组的组织执行模式也各不相同，常见的执行模式有SIMT（Single Instruction Multiple Threads，单指令多线程）执行模式和SIMD（Single Instruction Multiple Data，单指令流多数据流）执行模式两种。

    在一个GPU程序中，避免不了对数据的加载和存储，同时也避免不了条件分支跳转指令。这两类指令通常会引起程序以不可预测的情况执行。对于前者，在第一级高速缓存命中缺失的情况下，指令的执行周期将不可预测。为了避免执行单元因为数据加载或者存储原因而造成运算资源的浪费，GPU的每个执行单元通常设置线程组缓冲区，以支持同时执行多个线程组。线程组之间的调度由线程组硬件调度器承担，与软件调度器不同的是，硬件调度过程一般为零负载调度。在执行单元中，即将执行的线程组首先被调度到缓冲区中，以队列的方式组织，当线程组被调度执行时，调度器从线程组队列中选择一个准备好的线程组启动执行。采用这种线程调度执行方式，可有效解决指令之间由于长延时操作所引起的停顿问题，更高效的应用执行单元中的计算资源。对于后者，在线程级并行执行过程中，条件分支指令的执行特点决定了程序执行的实际效率。无论是SIMD执行模式或是SIMT执行模式，当一组线程均执行相同的代码路径时可获得最佳性能。若一组线程中的每个线程各自执行不同的代码路径，为了确保所有线程执行的正确性，线程组中的多线程指令发送单元将串行地发送所有的指令代码，代码的执行效率将受到严重的影响。GPU架构采用各种控制方法来提高条件分支指令的执行效率。

    背景知识大部分内容引自 [陈国良, 吴俊敏. 并行计算机体系结构（第2版）[M]. 北京: 高等教育出版社, 2021.]。

GPU 异构计算

    CPU-GPU 协同是实现高性能计算的必要条件，称为CPU-GPU异构计算（HC，Heterogeneous Computing）。它通过将应用程序的计算密集型部分卸载到GPU来提供更高的性能，而其余代码仍然在CPU上运行，能智能地结合CPU和GPU的最佳特性以实现高计算增益，旨在将每个应用程序的需求与CPU/GPU架构的优势相匹配，并避免两个处理单元的空闲时间。需要新的优化技术来充分发挥HC的潜力并朝着百亿级性能的目标迈进。

    了解CPU和GPU之间差异的一种简单方法是比较它们处理任务的方式。CPU由几个针对顺序串行处理优化的内核组成，而GPU具有大规模并行架构，由数千个更小、更高效的内核组成，旨在同时处理多个任务。

    在GPU上解决计算问题原则上类似于使用多个CPU解决问题。手头的任务必须拆分为小任务，其中每个任务由单个GPU内核执行。GPU内核之间的通信由GPU芯片上的内部寄存器和内存处理。CUDA或OpenCL等特殊编程语言不是使用消息传递进行编程，而是提供主机CPU之间的数据交换和同步GPU内核的机制。

    一个现代超级计算系统实际上可能由大量节点组成，每个节点包含2到32颗常规CPU以及1到16个GPU。通常还会有一个高速网络和一个数据存储系统。该系统的软件可以使用传统编程语言（如C/C++、Fortran等）的组合编写，结合用于CPU并行化的消息传递系统以及用于GPU的CUDA或OpenCL。所有这些组件都必须进行调整和优化，以实现整个系统的最佳性能。

异构架构

    一个典型的异构计算节点包括两个多核CPU插槽和两个或更多个的众核GPU。GPU不是一个独立运行的平台而是CPU的协处理器。因此，GPU必须通过PCIe总线与基于CPU的主机相连来进行操作，如下图所示。这就是为什么CPU所在的位置被称作主机端而GPU所在的位置被称作设备端。

• 主机代码

• 设备代码

    以下是描述GPU容量的两个重要特征。

• CUDA核心数量

• 内存大小

    相应的，有两种不同的指标来评估GPU的性能。

• 峰值计算性能

• 内存带宽

    峰值计算性能是用来评估计算容量的一个指标，通常定义为每秒能处理的单精度或双精度浮点运算的数量。峰值性能通常用GFlops（每秒十亿次浮点运算）或TFlops（每秒万亿次浮点运算）来表示。内存带宽是从内存中读取或写入数据的比率。内存带宽通常用GB/s表示。下表所示为 Fermi 架构和 Kepler 架构的一些性能指标。

异构计算范例

    GPU 计算并不是要取代 CPU 计算。对于特定的程序来说，每种计算方法都有它自己的优点。CPU 计算适合处理控制密集型任务，GPU 计算适合处理包含数据并行的计算密集型任务。GPU 与 CPU 结合后，能有效提高大规模计算问题的处理速度与性能。CPU 针对动态工作负载进行了优化，这些动态工作负载是由短序列的计算操作和不可预测的控制流程标记的；而 GPU 在其他领域内的目的是：处理由计算任务主导的且带有简单控制流的工作负载。如下图所示，可以从两个方面来区分 CPU 和 GPU 应用的范围。

• 并行级

• 数据规模

    如果一个问题有较小的数据规模、复杂的控制逻辑和/或很少的并行性，那么最好选择 CPU 处理该问题，因为它有处理复杂逻辑和指令级并行性的能力。相反，如果该问题包含较大规模的待处理数据并表现出大量的数据并行性，那么使用 GPU 是最好的选择。因为 GPU 中有大量可编程的核心，可以支持大规模多线程运算，而且相比 CPU 有较大的峰值带宽。

因为 CPU 和 GPU 的功能互补性导致了 CPU＋GPU的 异构并行计算架构的发展，这两种处理器的类型能使应用程序获得最佳的运行效果。因此，为获得最佳性能，你可以同时使用CPU和GPU来执行你的应用程序，在 CPU 上执行串行部分或任务并行部分，在 GPU 上执行数据密集型并行部分，如下图所示。

    这种代码的编写方式能保证 GPU 与 CPU 相辅相成，从而使 CPU ＋ GPU 系统的计算能力得以充分利用。为了支持使用 CPU＋GPU 异构系统架构来执行应用程序，NVIDIA设计了一个被称为 CUDA 的编程模型。

CPU 线程与 GPU 线程

    CPU 上的线程通常是重量级实体。操作系统必须在 CPU 执行通道上和从 CPU 执行通道上交换线程以提供多线程功能。上下文切换缓慢且开销大。

    GPU 上的线程非常轻量级。在典型的系统中，数千个线程排队等待工作。如果 GPU 必须等待一组线程，它就会开始执行另一组线程的工作。

    CPU 内核旨在一次最小化一个或两个线程的延迟，而 GPU 内核旨在处理大量并发的轻量级线程以最大化吞吐量。

   现在，具有四个四核 CPU 可以同时运行 16 个线程，如果 CPU 支持超线程，则可以同时运行 32 个线程。

    现代 NVIDIA GPU 可以支持每个多处理器同时最多 1,536 个活动线程。在具有 16 个多处理器的 GPU 上，可以并发支持超过 24,000 个同时活跃的线程。

CUDA：一种异构计算平台

    CUDA是一种通用的并行计算平台和编程模型，它利用NVIDIA GPU中的并行计算引擎能更有效地解决复杂的计算问题。通过使用CUDA，你可以像在CPU上那样，通过GPU来进行计算。

    CUDA平台可以通过CUDA加速库、编译器指令、应用编程接口以及行业标准程序语言的扩展（包括C、C++、Fortran、Python，如图下图所示）来使用。

    CUDA C 是标准ANSI C 语言的一个扩展，它带有的少数语言扩展功能使异构编程成为可能，同时也能通过API来管理设备、内存和其他任务。CUDA 还是一个可扩展的编程模型，它使程序能对有不同数量核的 GPU 明显地扩展其并行性，同时对熟悉C编程语言的程序员来说也比较容易上手。

    CUDA 提供了两层 API 来管理 GPU 设备和组织线程，如下图所示：

• CUDA 驱动 API

• CUDA 运行时 API

    驱动 API 是一种低级API，它相对来说较难编程，但是它对于在 GPU 设备使用上提供了更多的控制。运行时 API 是一个高级 API，它在驱动 API 的上层实现。每个运行时 API 函数都被分解为更多传给驱动 API 的基本运算。

    运行时 API 和驱动程序 API 之间没有明显的性能差异。在设备端，内核是如何使用内存以及程序员是如何在设备上组织线程的，对性能有显著的影响。

    这两种 API 是相互排斥的。我们必须使用其中一种，混合使用两种函数调用是不可能的。

    一个CUDA程序包含了以下两个部分的混合。

• 在 CPU 上运行的主机代码

• 在 GPU 上运行的设备代码

更加确切的说，CUDA 的操作概括来说包含6 个步骤:

• CPU 在GPU 上分配内存: cudaMalloc

• CPU 把数据发送到GPU: cudaMemcpy

• CPU 在GPU 上启动kernel, 它是自己写的一段程序, 在每个线程上运行

• CPU 等待GPU 端完成之前CUDA 的任务: cudaDeviceSynchronize

• CPU 把数据从GPU 取回: cudaMemcpy

• CPU 释放GPU 上的内存: cudaFree

** * Vector addition: C = A + B. * * This sample is a very basic sample that implements element by element * vector addition. It is the same as the sample illustrating Chapter 2 * of the programming guide with some additions like error checking. */
#include <stdio.h>
// For the CUDA runtime routines (prefixed with "cuda_")#include <cuda_runtime.h>
#include <helper_cuda.h>/** * CUDA Kernel Device code * * Computes the vector addition of A and B into C. The 3 vectors have the same * number of elements numElements. */__global__ void vectorAdd(const float *A, const float *B, float *C,int numElements) {int i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < numElements) {C[i] = A[i] + B[i] + 0.0f;}}
/** * Host main routine */int main(void) {// Error code to check return values for CUDA callscudaError_t err = cudaSuccess;
// Print the vector length to be used, and compute its sizeint numElements = 50000;size_t size = numElements * sizeof(float);printf("[Vector addition of %d elements]\n", numElements);
// Allocate the host input vector Afloat *h_A = (float *)malloc(size);
// Allocate the host input vector Bfloat *h_B = (float *)malloc(size);
// Allocate the host output vector Cfloat *h_C = (float *)malloc(size);
// Verify that allocations succeededif (h_A == NULL || h_B == NULL || h_C == NULL) {fprintf(stderr, "Failed to allocate host vectors!\n");exit(EXIT_FAILURE);}
// Initialize the host input vectorsfor (int i = 0; i < numElements; ++i) {h_A[i] = rand() / (float)RAND_MAX;h_B[i] = rand() / (float)RAND_MAX;}
// Allocate the device input vector Afloat *d_A = NULL;err = cudaMalloc((void **)&d_A, size);
if (err != cudaSuccess) {fprintf(stderr, "Failed to allocate device vector A (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
// Allocate the device input vector Bfloat *d_B = NULL;err = cudaMalloc((void **)&d_B, size);
if (err != cudaSuccess) {fprintf(stderr, "Failed to allocate device vector B (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
// Allocate the device output vector Cfloat *d_C = NULL;err = cudaMalloc((void **)&d_C, size);
if (err != cudaSuccess) {fprintf(stderr, "Failed to allocate device vector C (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
// Copy the host input vectors A and B in host memory to the device input// vectors in// device memoryprintf("Copy input data from the host memory to the CUDA device\n");err = cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
if (err != cudaSuccess) {fprintf(stderr,"Failed to copy vector A from host to device (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
err = cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
if (err != cudaSuccess) {fprintf(stderr,"Failed to copy vector B from host to device (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
// Launch the Vector Add CUDA Kernelint threadsPerBlock = 256;int blocksPerGrid = (numElements + threadsPerBlock - 1) / threadsPerBlock;printf("CUDA kernel launch with %d blocks of %d threads\n", blocksPerGrid, threadsPerBlock);vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, numElements);err = cudaGetLastError();
if (err != cudaSuccess) {fprintf(stderr, "Failed to launch vectorAdd kernel (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
// Copy the device result vector in device memory to the host result vector// in host memory.printf("Copy output data from the CUDA device to the host memory\n");err = cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);
if (err != cudaSuccess) {fprintf(stderr,"Failed to copy vector C from device to host (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
// Verify that the result vector is correctfor (int i = 0; i < numElements; ++i) {if (fabs(h_A[i] + h_B[i] - h_C[i]) > 1e-5) {fprintf(stderr, "Result verification failed at element %d!\n", i);exit(EXIT_FAILURE);}}
printf("Test PASSED\n");
// Free device global memoryerr = cudaFree(d_A);
if (err != cudaSuccess) {fprintf(stderr, "Failed to free device vector A (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
err = cudaFree(d_B);
if (err != cudaSuccess) {fprintf(stderr, "Failed to free device vector B (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
err = cudaFree(d_C);
if (err != cudaSuccess) {fprintf(stderr, "Failed to free device vector C (error code %s)!\n",cudaGetErrorString(err));exit(EXIT_FAILURE);}
// Free host memoryfree(h_A);free(h_B);free(h_C);
printf("Done\n");return 0;}

    在后续的文章中，我们会介绍 CUDA 编程示例。

    NVIDIA 的 CUDA nvcc 编译器在编译过程中将设备代码从主机代码中分离出来。如下图所示，主机代码是标准的 C 代码，使用 C 编译器进行编译。设备代码，也就是核函数，是用扩展的带有标记数据并行函数关键字的 CUDA C 语言编写的。设备代码通过 nvcc 进行编译。在链接阶段，在内核程序调用和显示 GPU 设备操作中添加 CUDA 运行时库。