CUDA 编程学习手册（一）

原文：annas-archive.org/md5/f6da79e769f988319eb178273ecbf55b

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

不要在第一次胜利后就休息。因为如果你在第二次失败，会有更多的人说你的第一次胜利只是运气。

A. P. J. Abdul Kalam

传统上，计算需求与中央处理单元（CPU）相关联，CPU 已经从单核发展到现在的多核。每一代新的 CPU 都提供了更多的性能，但科学和高性能计算社区每年都要求更多的性能，导致应用程序需求与硬件/软件堆栈提供的计算之间存在差距。与此同时，传统上用于视频图形的新架构也进入了科学领域。图形处理单元（GPU）——基本上是用于加速计算机图形的并行计算处理器——在 2007 年 CUDA 推出时进入了 HPC 领域。CUDA 成为了使用 GPU 进行通用计算的事实标准；即非图形应用程序。

自 CUDA 诞生以来，已经发布了许多版本，现在 CUDA 的版本为 10.x。每个版本都提供支持新硬件架构的新功能。本书旨在帮助您学习 GPU 并行编程，并指导您在现代应用中的应用。在它的帮助下，您将能够发现现代 GPU 架构的 CUDA 编程方法。本书不仅将指导您了解 GPU 功能、工具和 API，还将帮助您了解如何使用示例并行编程算法来分析性能。本书将确保您获得丰富的优化经验和对 CUDA 编程平台的洞察，包括各种库、开放加速器（OpenACC）和其他语言。随着您的进步，您将发现如何在一个盒子或多个盒子中利用多个 GPU 生成额外的计算能力。最后，您将探索 CUDA 如何加速深度学习算法，包括卷积神经网络（CNN）和循环神经网络（RNN）。

本书旨在成为任何新手或初学者开发者的入门点。但到最后，您将能够为不同领域编写优化的 CUDA 代码，包括人工智能。

如果您符合以下情况，这本书将是一个有用的资源：

您是 HPC 或并行计算的新手
您有代码，并希望通过将并行计算应用于 GPU 来提高其性能
您是深度学习专家，想利用 GPU 加速深度学习算法，如 CNN 和 RNN
您想学习优化代码和分析 GPU 应用性能的技巧和窍门，并发现优化策略
您想了解最新的 GPU 功能，以及高效的、分布式的多 GPU 编程。

如果您觉得自己属于以上任何一类，请加入我们一起踏上这段旅程。

这本书适合谁

这本初学者级别的书适用于希望深入研究并行计算、成为高性能计算社区的一部分并构建现代应用程序的程序员。假定具有基本的 C 和 C++编程经验。对于深度学习爱好者，本书涵盖了 Python InterOps、DL 库以及性能估计的实际示例。

为了充分利用这本书

本书适用于完全初学者和刚开始学习并行计算的人。除了计算机体系结构的基础知识外，不需要任何特定的知识，假定具有 C/C++编程经验。对于深度学习爱好者，在[第十章]（d0e9e8ff-bc17-4031-bb0e-1cfd310aff6f.xhtml）中，还提供了基于 Python 的示例代码，因此该章节需要一些 Python 知识。

本书的代码主要是在 Linux 环境中开发和测试的。因此，熟悉 Linux 环境是有帮助的。任何最新的 Linux 版本，如 CentOS 或 Ubuntu，都可以。代码可以使用 makefile 或命令行进行编译。本书主要使用免费软件堆栈，因此无需购买任何软件许可证。本书中将始终使用的两个关键软件是 CUDA 工具包和 PGI 社区版。

由于本书主要涵盖了利用 CUDA 10.x 的最新 GPU 功能，为了充分利用所有培训材料，最新的 GPU 架构（Pascal 及更高版本）将是有益的。虽然并非所有章节都需要最新的 GPU，但拥有最新的 GPU 将有助于您重现本书中实现的结果。每一章都有一个关于首选或必备 GPU 架构的部分，位于技术要求部分。

下载示例代码文件

您可以从您在www.packt.com的帐户中下载本书的示例代码文件。如果您在其他地方购买了本书，您可以访问www.packtpub.com/support并注册，以便将文件直接发送到您的邮箱。

您可以按照以下步骤下载代码文件：

在www.packt.com上登录或注册。
选择支持选项卡。
点击“代码下载”。
在搜索框中输入书名，然后按照屏幕上的说明操作。

下载文件后，请确保您使用最新版本的解压缩或提取文件夹：

WinRAR/7-Zip for Windows
Zipeg/iZip/UnRarX for Mac
7-Zip/PeaZip for Linux

该书的代码包也托管在 GitHub 上，网址为github.com/PacktPublishing/Learn-CUDA-Programming。如果代码有更新，将在现有的 GitHub 存储库上进行更新。

我们还提供了来自我们丰富书籍和视频目录的其他代码包，可在**github.com/PacktPublishing/**上找到。快来看看吧！

下载彩色图像

我们还提供了一个 PDF 文件，其中包含本书中使用的屏幕截图/图表的彩色图像。您可以在这里下载：static.packt-cdn.com/downloads/9781788996242_ColorImages.pdf。

使用的约定

本书中使用了许多文本约定。

CodeInText：表示文本中的代码词、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 句柄。以下是一个例子：“请注意，cudaMemcpy有一个异步的替代方案。”

代码块设置如下：

#include<stdio.h>
#include<stdlib.h>

__global__ void print_from_gpu(void) {
    printf("Hello World! from thread [%d,%d] \
        From device\n", threadIdx.x,blockIdx.x);
}

当我们希望引起您对代码块的特定部分的注意时，相关行或项目将以粗体显示：

int main(void) {
    printf("Hello World from host!\n");
    print_from_gpu<<<1,1>>>();
    cudaDeviceSynchronize();
    return 0;
}

任何命令行输入或输出都以如下形式编写：

$ nvcc -o hello_world hello_world.cu

粗体：表示新术语、重要单词或屏幕上看到的单词。例如，菜单中的单词或对话框中的单词会以这种形式出现在文本中。以下是一个例子：“对于 Windows 用户，在 VS 项目属性对话框中，您可以在 CUDA C/C++ | Device | Code Generation 中指定您的 GPU 的计算能力。”

警告或重要说明如下。

提示和技巧如下。

CUDA 编程简介

自 2007 年首次发布以来，统一计算设备架构（CUDA）已经成长为使用图形计算单元（GPU）进行通用计算的事实标准，即非图形应用程序。那么，CUDA 到底是什么？有人可能会问以下问题：

这是一种编程语言吗？
这是一个编译器吗？
这是一种新的计算范式吗？

在本章中，我们将揭开有关 GPU 和 CUDA 的一些神秘之谜。本章通过提供对高性能计算（HPC）历史的简化视图，并用摩尔定律和丹纳德缩放等法则加以证实，为异构计算奠定基础，这些法则一直在推动半导体行业，因此也推动了处理器架构本身。您还将了解 CUDA 编程模型，并了解 CPU 和 GPU 架构之间的根本区别。通过本章的学习，您将能够使用 C 语言中的 CUDA 编程构造编写和理解Hello World!程序。

尽管本章主要使用 C 语言来演示 CUDA 构造，但我们将在其他章节中涵盖其他编程语言，如 Python、Fortran 和 OpenACC。

本章将涵盖以下主题：

高性能计算的历史
来自 CUDA 的 Hello World
使用 CUDA 进行矢量加法
CUDA 的错误报告
CUDA 中的数据类型支持

高性能计算的历史

高性能计算一直在不断突破极限，以实现科学发现。处理器架构和设计的根本转变有助于跨越 FLOP 障碍，从百万浮点运算（MFLOPs）开始，现在能够在一秒内进行 PetaFLOP 计算。

每秒浮点运算（FLOPs）是衡量任何计算处理器理论峰值的基本单位。MegaFLOP 代表 FLOPS 的 10 的 6 次方。PetaFLOP 代表 FLOPS 的 10 的 15 次方。

指令级并行性（ILP）是一个概念，其中独立于代码的指令可以同时执行。为了使指令并行执行，它们需要彼此独立。所有现代 CPU 架构（甚至 GPU 架构）都提供了五到 15 个以上的阶段，以实现更快的时钟频率：

Instr 1: add = inp1 + inp2

Instr 2: mult = inp1 * inp2

Instr 3: final_result = mult / add

用于计算mult和add变量的操作不相互依赖，因此可以在计算final_result时同时计算，而final_result依赖于Instr 1和Instr 2操作的结果。因此，在计算add和mult之前无法计算它。

当我们从技术变革的角度看高性能计算的历史，这些变革导致了新处理器设计的根本转变，以及对科学界的影响，有三个主要的变革可以被称为时代：

时代 1：超级计算机的历史可以追溯到 CRAY-1，它基本上是一个提供峰值 160 MegaFLOP/MFLOP 计算能力的单一矢量 CPU 架构。
时代 2：通过从单核设计转向 CRAY-2 的多核设计，跨越了 MegaFLOP 障碍，CRAY-2 是一个 4 核矢量 CPU，提供了 2 GigaFLOPs 的峰值性能。
时代 3：跨越 GigaFLOP 计算性能是一个根本性的转变，需要计算节点相互协作，并通过网络进行通信，以提供更高的性能。 Cray T3D 是第一批提供 1 TeraFLOP 计算性能的机器之一。网络是 3D Torus，提供 300 MB/s 的带宽。这是标准微处理器周围丰富shell的第一个重要实现。

此后，将近 20 年没有根本性的创新。技术创新主要集中在三个架构创新上：

从 8 位到 16 位再到 32 位，现在是 64 位指令集
增加 ILP
增加核心数量

这得到了时钟频率的增加，目前为 4 GHz。由于驱动半导体行业的基本定律，这是可能的。

摩尔定律：这个定律观察到密集集成电路中的晶体管数量每两年翻一番。

摩尔的预测在几十年来一直准确无误。摩尔定律是对历史趋势的观察和预测。

Dennard 缩放：这是一个使摩尔定律保持活力的缩放定律。Dennard 观察到晶体管尺寸和功率密度之间的关系，并用以下公式总结了这一观察：

P = QfCV² + V I[leakage]

在这个方程中，Q是晶体管数量，f是操作频率，C是电容，V是操作电压，*I[leakage]*是泄漏电流。

Dennard 缩放和摩尔定律彼此相关，因为可以推断出，减小晶体管的尺寸可以在成本效益方面导致芯片上的晶体管数量越来越多。

根据 Dennard 缩放规则，对于给定尺寸的芯片，多个处理器世代的总芯片功率保持不变。晶体管数量翻倍，而尺寸不断缩小（1/S速率），并且每两年以 40%的速度增加频率。当特征尺寸达到 65 纳米以下时，这种情况停止了，因为泄漏电流呈指数增长，这些规则不再能够持续。为了减少泄漏电流的影响，新的创新被强制执行在开关过程中。然而，这些突破仍然不足以恢复电压的缩放。电压在许多处理器设计中保持在 1V 恒定。不再可能保持功率包络恒定。这也被称为 Powerwall。

Dennard 缩放从 1977 年一直持续到 1997 年，然后开始衰退。因此，从 2007 年到 2017 年，处理器从 45 纳米变为 16 纳米，但导致每芯片能耗增加了三倍。

同时，管线阶段从五个阶段发展到了最新架构的 15+个阶段。为了保持指令管线的充分，使用了先进的技术，比如推测。推测单元涉及预测程序的行为，比如预测分支和内存地址。如果预测准确，它可以继续；否则，它会撤销已经完成的工作并重新开始。深层管线阶段和传统软件的编写方式导致了未使用的晶体管和浪费的时钟周期，这意味着应用性能没有改善。

然后出现了 GPU，最初主要用于图形处理。研究人员马克·哈里斯首次利用 GPU 进行了非图形任务，并创造了新术语使用 GPU 进行通用计算（GPGPU）。GPU 在某些数据并行类任务方面被证明是有效的。毫不奇怪，许多 HPC 应用程序中的大部分计算密集型任务在性质上都是数据并行的。它们主要是矩阵乘法，这在基本线性代数规范（BLAS）中是常规且广泛使用的。

用户在适应和使用 GPU 时唯一的问题是他们必须了解图形管线以利用 GPU。提供给 GPU 上任何计算工作的唯一接口围绕着着色器执行。需要提供一个更通用的接口，让在 HPC 社区工作的开发人员熟悉。这个问题在 2007 年引入 CUDA 时得到解决。

虽然 GPU 架构也受到相同的定律约束（摩尔定律和 Dennard 缩放），但处理器的设计采用了不同的方法，为不同的用途专门分配晶体管，并实现了比传统的同质架构更高的性能。

以下图表显示了计算机体系结构从顺序处理到分布式内存的演变及其对编程模型的影响：

随着 GPU 被添加到现有服务器上，应用程序在两种处理器（CPU 和 GPU）上运行，引入了异构的概念。这是我们将在下一节介绍的内容。

异构计算

围绕 GPU 的一个常见误解是它是 CPU 的替代品。GPU 用于加速代码中并行的部分。加速器是一个常用术语，用于描述 GPU，因为它们通过更快地运行代码的并行部分来加速应用程序，而 CPU 运行另一部分代码，即延迟绑定的部分。因此，高效的 CPU 与高吞吐量的 GPU 相结合，可以提高应用程序的性能。

以下图表代表了在多种处理器类型上运行的应用程序：

这个概念可以很好地用阿姆达尔定律来定义。阿姆达尔定律用于定义当应用程序的一部分被并行化时可以实现的最大加速。为了演示这一点，前面的图表显示了代码的两个部分。一个部分是延迟绑定的，而另一个是吞吐量绑定的。我们将在下一节中介绍这两个术语，区分 CPU 和 GPU 体系结构。

关键点是，CPU 对于某些延迟绑定的代码部分很好，而 GPU 擅长并行运行代码的单指令多数据（SIMD）部分。如果在优化后只有其中一个，即 CPU 代码或 GPU 代码，运行速度更快，这不一定会导致整体应用程序的速度提升。需要的是，当两个处理器都得到最佳利用时，性能方面才能获得最大的好处。这种从处理器上卸载某些类型的操作到 GPU 的方法被称为异构计算。

以下图表描述了所有应用程序具有的两种部分，即延迟绑定和吞吐量绑定：

在这里，使用阿姆达尔定律演示了改进两个部分的重要性。

编程范式

计算机体系结构的分类是使用弗林分类法进行的，描述了四类体系结构。弗林的分类之一 SIMD 用于描述 GPU 体系结构。然而，两者之间存在微妙的差异。SIMD 用于描述同一指令并行应用于多个数据点的体系结构。这种描述适用于具有矢量化能力的处理器。相比之下，在单指令多线程（SIMT）中，不是单个线程发出指令，而是多个线程向不同的数据发出相同的指令。与 SIMD 相比，GPU 体系结构更适合 SIMT 类别。

让我们看一个例子，将两个数组相加并将数据存储在第三个数组中。这个操作的数据集包括数组A、B和C。用于加法的相同操作被用于数组的每个元素：

Cx = Ax + Bx

显然，每个任务都是独立的，但所有线程都在应用相同的操作。

以下截图显示了矢量加法，展示了这种范例的一个例子：

低延迟与高吞吐量

正如我们在前一节中提到的，CPU 架构被优化用于低延迟访问，而 GPU 架构被优化用于数据并行吞吐量计算。如下截图所示，与 GPU 相比，CPU 架构具有大量缓存并且具有许多类型。我们越高，即从 L3 到 L1，缓存的数量就越少，但延迟就越低。CPU 架构旨在实现对缓存数据集的低延迟访问。大量晶体管用于实现推测执行和乱序执行。由于 CPU 以非常高的时钟速度运行，因此有必要通过频繁地将使用的数据存储在缓存中并预测下一条要执行的指令来隐藏获取数据的延迟。可以最佳地利用 CPU 缓存的应用程序可以探索这种时间局部性。此外，可以利用填充指令管线的应用程序，例如代码中没有if和else语句的应用程序，通过隐藏获取指令的延迟来受益。因此，CPU 架构是一种减少延迟的架构。

以下截图显示了 CPU 和 GPU 架构如何为不同的内存和计算单元分配芯片芯片区域。GPU 使用大量晶体管进行计算 ALUs，而 CPU 使用它来减少延迟。

另一方面，GPU 架构被称为“减少延迟”或“高吞吐量架构”。GPU 架构通过来自其他线程的计算来隐藏延迟。当一个线程在等待数据可用进行计算时，其他线程可以开始执行，因此不会浪费任何时钟周期。如果您熟悉 CUDA，那么您可能已经了解到 warp 的概念。我们将在接下来的章节中介绍 warp 的概念。（在 CUDA 中，执行单元是 warp 而不是线程。因此，上下文切换发生在 warp 而不是线程之间）。

有些人可能已经在想为什么我们不能在 CPU 中创建这些线程并做同样的事情来隐藏延迟。原因是 GPU 有大量的寄存器，并且所有线程上下文切换信息已经存在于其中。这是最快的内存。然而，在 CPU 中，寄存器集是有限的，因此线程相关的信息通常存储在较低的内存层次结构中，比如缓存。例如，Volta 包含 20MB 的寄存器存储。因此，与 GPU 相比，CPU 中线程之间的上下文切换时间要长得多。

现在，让我们来看看在 GPU 编程方面的不同方法。

GPU 的编程方法

现在，让我们回到我们最初的问题，即 CUDA 是什么？CUDA 是由 NVIDIA 开发的并行计算平台和编程模型架构，它将 GPU 上的通用计算作为一流能力进行暴露。与任何其他处理器一样，GPU 架构可以使用各种方法进行编码。提供快速加速的最简单方法是利用现有库。另外，开发人员可以选择使用 OpenACC 指令以获得快速加速结果和可移植性。另一种选择是选择通过使用 C、C++、Fortran、Python 等语言构造来深入研究 CUDA，以获得最高的性能和灵活性。我们将在接下来的章节中详细介绍所有这些方法。

以下截图表示了我们可以进行 GPU 编程的各种方式。

在本节中，我们为您提供了处理器和高性能计算随时间演变的视角。我们为您提供了异构编程模型对于从应用程序中获得最佳性能的关键性概述，以及 GPU 编程的方法。在下一节中，我们将开始在 GPU 上编写一个 Hello World 程序。

技术要求

本章需要一台安装了现代 NVIDIA GPU（Pascal 架构及以上）的 Linux/Windows PC，以及所有必要的 GPU 驱动程序和安装了 CUDA Toolkit（10.0 及以上版本）。如果您不确定您的 GPU 架构，请访问 NVIDIA 的 GPU 网站（developer.nvidia.com/cuda-gpus）并确认您的 GPU 架构。本章的代码也可以在 GitHub 上找到：github.com/PacktPublishing/Learn-CUDA-Programming。

本章的代码示例是使用 CUDA Toolkit 的 10.1 版本开发和测试的，但建议尽可能使用最新的 CUDA 版本。

来自 CUDA 的 Hello World

CUDA 是一个包括 CPU 和 GPU 在内的异构编程模型。CUDA C/C++编程接口由 C 语言扩展组成，以便您可以将源代码的部分目标定为在设备（GPU）上并行执行。它基于行业标准的 C/C++，并提供了一系列 C 函数库，可以在主机（CPU）上执行，以便它可以与设备进行交互。

在 CUDA 中，有两个相互配合的处理器。主机通常被称为 CPU，而设备通常被称为 GPU。主机负责调用设备函数。正如我们已经提到的，运行在 GPU 上的代码的一部分被称为设备代码，而在 CPU 上运行的串行代码被称为主机代码。

让我们从在 C 中编写我们的第一个 CUDA 代码开始。我们的意图是采取一个系统化的逐步方法，从一些顺序代码开始，通过添加一些额外的关键字将其转换为 CUDA 感知代码。正如我们之前提到的，没有必要学习一门新语言，我们只需要在现有语言中添加一些关键字，以便在 CPU 和 GPU 的异构环境中运行它。

让我们来看看我们的第一段代码。这段代码的作用只是从主机和设备上打印 Hello World!

#include<stdio.h>
#include<stdlib.h>

__global__ void print_from_gpu(void) {
    printf("Hello World! from thread [%d,%d] \
        From device\n", threadIdx.x,blockIdx.x);
}

int main(void) {
    printf("Hello World from host!\n");
    print_from_gpu<<<1,1>>>();
    cudaDeviceSynchronize();
    return 0;
}

让我们尝试编译和运行前面的代码片段：

编译代码：将前面的代码放入一个名为hello_world.cu的文件中，并使用NVIDIA C Compiler（nvcc）进行编译。请注意，文件的扩展名是.cu，这告诉编译器这个文件里面有 GPU 代码：

$ nvcc -o hello_world hello_world.cu

执行 GPU 代码：在执行 GPU 代码后，我们应该收到以下输出：

到目前为止，您可能已经注意到 CUDA C 代码的使用方式并没有太大不同，只需要学习一些额外的构造来告诉编译器哪个函数是 GPU 代码，以及如何调用 GPU 函数。这并不像我们需要完全学习一门新语言。

在前面的代码中，我们添加了一些构造和关键字，如下：

__global__：在函数之前添加此关键字，告诉编译器这是一个将在设备上而不是主机上运行的函数。但请注意，它是由主机调用的。这里另一个重要的事情是设备函数的返回类型始终是"void"。算法的数据并行部分在设备上作为内核执行。
<<<,>>>: 这个关键字告诉编译器这是对设备函数的调用，而不是对主机函数的调用。此外，1,1参数基本上决定了在内核中启动的线程数。我们将稍后介绍尖括号内的参数。目前，1,1参数基本上意味着我们只启动一个线程的内核，也就是说，除了打印之外，我们在代码中没有做任何重要的事情。
threadIdx.x, blockIdx.x: 这是给所有线程的唯一 ID。我们将在下一节更详细地介绍这个主题。
cudaDeviceSynchronize(): CUDA 中的所有内核调用都是异步的。在调用内核后，主机变得空闲，并在之后开始执行下一条指令。这应该不足为奇，因为这是一个异构环境，因此主机和设备都可以并行运行，以利用可用的处理器类型。如果主机需要等待设备完成，CUDA 编程提供了 API 使主机代码等待设备函数完成。其中一个 API 是cudaDeviceSynchronize，它会等待所有先前对设备的调用完成。

尝试删除cudaDeviceSynchronize()调用，看看设备输出是否可见。或者，尝试在打印主机代码之前放置这个调用。

线程层次结构

现在，让我们开始玩弄两个参数，即threadIdx.x和blockIdx.x。

实验 1：首先，将参数从<<<1,1>>>更改为<<<2,1>>>并查看输出。运行多个线程-单个块的 Hello World 代码的输出应该如下：

正如我们所看到的，现在我们不是一个线程，而是两个线程打印值。请注意，它们的唯一 ID 是不同的。

实验 2：现在，不要更改第一个参数，而是更改第二个参数，即将<<<1,1>>>更改为<<<1,2>>>，并观察运行多个单线程块的 Hello World 代码的输出，如下所示：

如您所见，被启动到内核中的线程总数是两个，就像以前一样——唯一的区别是它们的 ID 不同。那么，这些线程和块的概念是什么？为了解决这个问题，让我们更深入地了解 GPU 架构。

GPU 架构

CUDA 变得如此受欢迎的一个关键原因是因为硬件和软件被设计和紧密绑定，以获得应用程序的最佳性能。因此，有必要展示软件 CUDA 编程概念与硬件设计本身之间的关系。

以下截图显示了 CUDA 的两个方面：

我们可以看到，CUDA 软件已经映射到了 GPU 硬件。

根据前面的截图，以下表解释了 CUDA 编程模型的软件和硬件映射：

软件	执行于/作为	硬件
CUDA 线程	CUDA 核心/SIMD 代码
CUDA 块	流多处理器
网格/内核	GPU 设备

让我们详细看一下前表的组件：

CUDA 线程：CUDA 线程在 CUDA 核心上执行。CUDA 线程不同于 CPU 线程。CUDA 线程非常轻量级，并提供快速的上下文切换。快速上下文切换的原因是由于 GPU 中有大量的寄存器和硬件调度器。线程上下文存在于寄存器中，而不是像 CPU 中那样在较低的内存层次结构中，比如缓存中。因此，当一个线程处于空闲/等待状态时，另一个准备好的线程几乎可以立即开始执行。每个 CUDA 线程必须执行相同的内核，并且独立地处理不同的数据（SIMT）。
CUDA 块：CUDA 线程被组合成一个称为 CUDA 块的逻辑实体。CUDA 块在单个流多处理器（SM）上执行。一个块在一个 SM 上运行，也就是说，一个块内的所有线程只能在一个 SM 的核心上执行，不会在其他 SM 的核心上执行。每个 GPU 可能有一个或多个 SM，因此为了有效地利用整个 GPU，用户需要将并行计算划分为块和线程。
GRID/核心：CUDA 块被组合成一个称为 CUDA GRID 的逻辑实体。然后在设备上执行 CUDA GRID。

乍一看，这可能听起来有些复杂。在接下来的部分，我们将以向量加法的示例来解释这个问题。希望事情会变得更清晰。

使用 CUDA 进行向量加法

我们要解决的问题是向量加法。正如我们所知，向量加法是一种数据并行操作。我们的数据集包括三个数组：A、B和C。每个元素执行相同的操作：

Cx = Ax + Bx

每个加法是相互独立的，但所有 CUDA 线程都执行相同的操作。要开始，根据以下步骤配置您的环境：

准备你的 GPU 应用程序。这段代码将放在01_cuda_introduction/01_vector_addition中。
使用nvcc编译器编译您的应用程序，命令如下：

$nvcc -o vector_addition vector_addition.cu

上述代码是顺序代码。我们将按照以下步骤逐步将此代码转换为可以在 GPU 上运行的代码：

#include<stdio.h>
#include<stdlib.h>

#define N 512

void host_add(int *a, int *b, int *c) {
    for(int idx=0;idx<N;idx++)
        c[idx] = a[idx] + b[idx];
}

//basically just fills the array with index.
void fill_array(int *data) {
    for(int idx=0;idx<N;idx++)
        data[idx] = idx;
}

void print_output(int *a, int *b, int*c) {
    for(int idx=0;idx<N;idx++)
        printf("\n %d + %d = %d", a[idx] , b[idx], c[idx]);
}

int main(void) {
    int *a, *b, *c;
    int size = N * sizeof(int);
   // Alloc space for host copies of a, b, c and setup input values
    a = (int *)malloc(size); fill_array(a);
    b = (int *)malloc(size); fill_array(b);
    c = (int *)malloc(size);
    host_add(a,b,c);
    print_output(a,b,c);
    free(a); free(b); free(c);
    return 0;
}

在转换顺序代码之前，让我们看一下 CUDA 代码和顺序代码之间所采取的基本变化或步骤：

顺序代码	CUDA 代码
步骤 1	在 CPU 上分配内存，即`malloc new`。	步骤 1	在 CPU 上分配内存，即`malloc new`。
步骤 2	填充/初始化 CPU 数据。	步骤 2	在 GPU 上分配内存，即`cudaMalloc`。
步骤 3	调用处理数据的 CPU 函数。在这种情况下，实际算法是向量加法。	步骤 3	填充/初始化 CPU 数据。
步骤 4	处理数据，这里是打印出来的。	步骤 4	用`cudaMemcpy`将数据从主机传输到设备。
步骤 5	用`<<<,>>>`括号调用 GPU 函数。
步骤 6	用`cudaDeviceSynchronize`同步设备和主机。
步骤 7	用`cudaMemcpy`将数据从设备传输到主机。
步骤 8	处理数据，这里是打印出来的。

本书不是 CUDA API 指南的替代品，也不涵盖所有 CUDA API。如需广泛使用 API，请参考 CUDA API 指南。

正如我们所看到的，CUDA 处理流程有一些额外的步骤需要添加到顺序代码中。具体如下：

**在 GPU 上分配内存：**CPU 内存和 GPU 内存是物理上分开的内存。malloc在 CPU 的 RAM 上分配内存。GPU 核心/设备函数只能访问已分配/指向设备内存的内存。要在 GPU 上分配内存，我们需要使用cudaMalloc API。与malloc命令不同，cudaMalloc不会返回指向已分配内存的指针；相反，它以指针引用作为参数，并更新相同的已分配内存。
**将数据从主机内存传输到设备内存：**主机数据然后被复制到使用前一步中使用的cudaMalloc命令分配的设备内存。用于在主机和设备之间复制数据的 API 是cudaMemcpy。与其他memcopy命令一样，此 API 需要目标指针、源指针和大小。它额外需要一个参数，即复制的方向，也就是说，我们是从主机到设备复制，还是从设备到主机复制。在 CUDA 的最新版本中，这是可选的，因为驱动程序能够理解指针是指向主机内存还是设备内存。请注意，cudaMemcpy有一个异步的替代方案。这将在其他章节中更详细地介绍。
**调用和执行 CUDA 函数：**如同在 Hello World CUDA 程序中所示，我们通过<<<,>>>括号调用一个核函数，它提供了块和线程大小的参数。在完成所有步骤后，我们将更详细地介绍这一点。
**同步：**正如我们在 Hello World 程序中提到的，核函数调用是异步的。为了确保主机确保核执行已完成，主机调用cudaDeviceSynchronize函数。这确保了之前启动的所有设备调用都已完成。
**将数据从主机内存传输到设备内存：**使用相同的cudaMemcpy API 将数据从设备复制回主机，用于后处理或验证任务，如打印。与第一步相比，唯一的变化是我们颠倒了复制的方向，也就是说，目标指针指向主机，而源指针指向在内存中分配的设备。
**释放分配的 GPU 内存：**最后，使用cudaFree API 释放分配的 GPU 内存。

更改顺序向量加法代码的main函数以反映这些新步骤。main函数将如下所示：

int main(void) {
    int *a, *b, *c;
    int *d_a, *d_b, *d_c; // device copies of a, b, c
    int size = N * sizeof(int);

    // Alloc space for host copies of a, b, c and setup input values
    a = (int *)malloc(size); fill_array(a);
    b = (int *)malloc(size); fill_array(b);
    c = (int *)malloc(size);

    // Alloc space for device copies of vector (a, b, c)
    cudaMalloc((void **)&d_a, N* * sizeof(int));
    cudaMalloc((void **)&d_b, N* *sizeof(int));
    cudaMalloc((void **)&d_c, N* * sizeof(int));

    // Copy from host to device
    cudaMemcpy(d_a, a, N * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, b, N* sizeof(int), cudaMemcpyHostToDevice);

    device_add<<<1,1>>>(d_a,d_b,d_c);

    // Copy result back to host
    cudaMemcpy(c, d_c, N * sizeof(int), cudaMemcpyDeviceToHost);

    print_output(a,b,c);
    free(a); free(b); free(c);

    //free gpu memory
    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);

    return 0;
}

现在，让我们看一下如何编写核心代码并管理线程和块大小。为此，我们将进行多个实验。

实验 1 - 创建多个块

在这一部分，我们将利用 CUDA 块在 GPU 上并行运行向量加法代码。将会暴露与我们如何索引 CUDA 块相关的附加关键字。更改对device_add函数的调用如下：

//changing from device_add<<<1,1>>> to
device_add<<<N,1>>>

这将使device_add函数并行执行N次，而不是一次。device_add函数的每个并行调用称为一个块。现在，让我们添加一个__global__设备函数，如下所示：

__global__ void device_add(int *a, int *b, int *c) {
 c[blockIdx.x] = a[blockIdx.x] + b[blockIdx.x];
}

通过使用blockIdx.x来索引数组，每个块处理数组的不同元素。在设备上，每个块可以并行执行。让我们看一下以下的屏幕截图：

前面的屏幕截图代表了向量加法 GPU 代码，其中每个块显示了多个单线程块的索引。

实验 2 - 创建多个线程

在这一部分，我们将利用 CUDA 线程在 GPU 上并行运行向量加法代码。将会暴露与我们如何索引 CUDA 线程相关的附加关键字。

一个块可以分成多个线程。更改对device_add函数的调用如下：

//changing from device_add<<<1,1>>> to
device_add<<<1,N>>>

这将并行执行device_add函数N次，而不是一次。device_add函数的每个并行调用被称为一个线程。更改设备例程以反映内核，如下所示：

__global__ void device_add(int *a, int *b, int *c) {
     c[threadIdx.x] = a[threadIdx.x] + b[threadIdx.x];
}

一个显著的区别是，我们使用threadIdx.x而不是blockIdx.x，如下面的截图所示：

前面的截图代表了向量加法 GPU 代码，其中每个块显示了单个块-多个线程的索引。

实验 3 - 结合块和线程

到目前为止，我们已经通过在实验 1 - 创建多个块部分使用多个块和一个线程，以及在实验 2 - 创建多个线程部分使用一个块和多个线程来查看并行向量加法。在这个实验中，我们将使用多个块以及包含多个线程的单独块。在如何找到索引方面，这变得更具挑战性，因为我们需要结合threadIdx和blockIdx来生成一个唯一的 ID。

让我们看一下两种不同组合的场景，开发人员可以从中选择：

场景 1： 假设向量元素的总数是 32。每个块包含八个线程，总共有四个块。
场景 2： 假设向量元素的总数是 32。每个块包含四个线程，总共有八个块。

在这两种情况下，并行执行的数量都是 32，所有 32 个元素都会并行填充。开发人员根据问题的大小和每个硬件的限制选择块内的线程和块的数量。我们将在另一章节中详细介绍基于架构的正确尺寸选择的细节。

以下截图显示了不同块和线程配置的向量加法 GPU 索引代码：

现在，让我们看看如何改变内核代码以结合线程和块来计算全局索引：

__global__ void device_add(int *a, int *b, int *c) {
     int index = threadIdx.x + blockIdx.x * blockDim.x;
     c[index] = a[index] + b[index];
}

在从main()函数中调用内核时，开发人员选择了块和线程的配置，如前面提到的两种情况所示的代码：

场景 1： 以下是用于计算每个块八个线程的向量加法 GPU 网格和块大小的代码：

threads_per_block = 8;
no_of_blocks = N/threads_per_block;
device_add<<<no_of_blocks,threads_per_block>>>(d_a,d_b,d_c);

场景 2： 以下是用于计算每个块四个线程的向量加法 GPU 网格和块大小的代码：

threads_per_block = 4;
no_of_blocks = N/threads_per_block;
device_add<<<no_of_blocks,threads_per_block>>>(d_a,d_b,d_c);

通过线程和块的组合，可以计算出线程的唯一 ID。如前面的代码所示，所有线程都被赋予另一个变量。这被称为blockDim。这个变量包含了块的维度，也就是每个块的线程数。让我们看一下下面的截图：

在这里，我们可以看到场景 1 的向量加法 GPU 索引计算。

为什么要费心处理线程和块？

可能不明显为什么我们需要这种额外的线程和块的层次结构。它增加了开发人员需要找到正确的块和网格大小的复杂性。全局索引也变得具有挑战性。这是因为 CUDA 编程模型设置了限制。

与并行块不同，线程有有效的通信和同步机制。现实世界的应用程序需要线程之间进行通信，并且可能希望在继续之前等待某些数据进行交换。这种操作需要线程进行通信，CUDA 编程模型允许同一块内的线程进行通信。属于不同块的线程在内核执行期间无法进行通信/同步。这种限制允许调度程序独立地在 SM 上调度块。其结果是，如果发布了具有更多 SM 的新硬件，并且代码具有足够的并行性，则代码可以线性扩展。换句话说，这允许硬件根据 GPU 的能力并行运行块的数量。

线程之间使用一种称为共享内存的特殊内存进行通信。我们将在第二章中广泛介绍共享内存，即CUDA 内存管理，在那里我们将介绍 GPU 中的其他内存层次结构及其最佳使用方法。以下屏幕截图演示了在不同 GPU 上扩展块，这些 GPU 包含不同数量的 SM：

现在，让我们更多地了解在多个维度中启动内核。

在多个维度中启动内核

到目前为止，我们一直在一维中启动线程和块。这意味着我们只使用一个维度的索引；例如，我们一直在使用threadIdx.x，其中x表示我们只使用一个x维度的线程索引。同样，我们一直在使用blockIdx.x，其中x表示我们只使用一个x维度的块索引。我们可以在一、二或三个维度中启动线程和块。在二维中启动线程和块的一个例子是当我们在图像上进行并行操作，例如使用滤波器模糊图像。开发人员可以选择在二维中启动线程和块，这是一个更自然的选择，因为图像在本质上是二维的。

重要的是要了解每个 GPU 架构也对线程和块的维度施加了限制。例如，NVIDIA Pascal 卡允许在x和y维度中每个线程块最多有 1,024 个线程，而在z维度中，您只能启动 64 个线程。同样，在 Pascal 架构中，网格中的最大块数限制为y和z维度中的 65,535 个，x维度中为2³¹ -1。如果开发人员使用不受支持的维度启动内核，应用程序会抛出运行时错误。

到目前为止，我们一直假设我们编写的代码是没有错误的。但在现实世界中，每个程序员都会写有错误的代码，必须捕捉这些错误。在下一节中，我们将看看 CUDA 中的错误报告是如何工作的。

CUDA 中的错误报告

在 CUDA 中，主机代码管理错误。大多数 CUDA 函数调用cudaError_t，它基本上是一个枚举类型。cudaSuccess（值 0）表示0错误。用户还可以使用cudaGetErrorString()函数，该函数返回描述错误条件的字符串。

 cudaError_t e;
 e = cudaMemcpy(...);
 if(e)
     printf("Error: %sn", cudaGetErrorString(err));

内核启动没有返回值。我们可以在这里使用cudaGetLastError()这样的函数，它返回最后一个 CUDA 函数（包括内核启动）的错误代码。在多个错误的情况下，只报告最后一个错误：

MyKernel<<< ... >>> (...);
cudaDeviceSynchronize();
e = cudaGetLastError();

在生产代码中，建议在逻辑检查点处使用错误检查代码，因为即使 GPU 内核崩溃，CPU 代码也会继续正常执行，导致结果不正确。

在下一节中，我们将向您介绍 CUDA 编程模型中支持的数据类型。

CUDA 中的数据类型支持

与任何处理器架构一样，GPU 也有不同类型的内存，每种用于不同的目的。我们将在第二章中更详细地介绍它们，CUDA 内存管理。然而，重要的是要理解支持的不同数据类型及其对性能和精度的影响。CUDA 编程支持开发人员在其各自语言中熟悉的所有标准数据类型。除了不同大小的标准数据类型（char为 1 字节，float为 4 字节，double为 8 字节等）之外，它还支持矢量类型，如float2和float4。

建议数据类型自然对齐，因为对于大小为 1、2、4、8 或 16 字节的数据类型的对齐数据访问，可以确保 GPU 调用单个内存指令。如果它们没有对齐，编译器将生成多个交错的指令，导致内存和指令总线的低效利用。因此，建议在 GPU 内存中使用自然对齐的类型。对于char、short、int、long、long long、float和double等内置类型，如float2和float4，对齐要求会自动满足。

此外，CUDA 编程还支持复杂的数据结构，如结构和类（在 C 和 C++的上下文中）。对于复杂的数据结构，开发人员可以利用对齐说明符来强制编译器满足对齐要求，如下面的代码所示：

struct __align__(16) {
    float r;
    float g;
    float b;
};

每个 GPU 都有一组有限的核心，因此 FLOPS 是不同的。例如，具有 Volta 架构的 Tesla V100 卡具有 2560 个 FP64 核心（双精度），而具有两倍数量的 32 位单精度核心。很明显，根据算法的精度要求使用正确的数据类型是至关重要的。现在正在开发混合精度算法，以利用不同类型的核心，其中算法的某些部分以更高精度运行，而某些部分以较低精度运行。我们将在即将到来的章节中更多地涵盖这个主题。目前，重要的是要理解 GPU 内存层次结构是不同的，因此使用正确的数据类型很重要。

虽然这是对 GPU 支持的数据类型的一般介绍，但有关所有支持的数据类型的更多详细信息可以在docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#built-in-vector-types找到。

总结

在本章中，我们通过历史和高性能计算为您提供了异构计算的视角。我们详细介绍了两个处理器，即 CPU 和 GPU 的不同之处。我们还在 GPU 上编写了一个 Hello World 和矢量加法 CUDA 程序。最后，我们看了如何检测 CUDA 中的错误，因为对 CUDA API 的大多数调用都是异步的。

在下一章中，我们将看一下不同类型的 GPU 内存以及如何最优地利用它们。

CUDA 内存管理

正如我们在第一章中所描述的，CUDA 编程简介，CPU 和 GPU 架构在根本上是不同的，它们的内存层次结构也是不同的。它们不仅在大小和类型上有所不同，而且在目的和设计上也有所不同。到目前为止，我们已经学习了每个线程如何通过索引（blockIdx和threadIdx）访问自己的数据。我们还使用了诸如cudaMalloc之类的 API 在设备上分配内存。GPU 中有许多内存路径，每个路径的性能特征都不同。启动 CUDA 核心可以帮助我们实现最大性能，但只有在以最佳方式使用正确类型的内存层次结构时才能实现。将数据集映射到正确的内存类型是开发人员的责任。

根据经验，如果我们绘制一个图表，概述 GPU 上的顶级应用性能约束，它将看起来像以下图表：

上述饼图粗略地分解了大多数基于 CUDA 的应用程序中出现的性能问题。很明显，大多数情况下，应用程序的性能将受到与内存相关的约束的限制。根据应用程序和采取的内存路径，内存相关的约束进一步划分。

让我们以不同的方式来看待这种方法，并了解有效使用正确类型的内存的重要性。最新的 NVIDIA GPU 采用 Volta 架构，提供了 7,000 GFLOP 的峰值性能，其设备内存带宽为 900 GB/s。您将首先注意到的是 FLOP 与内存带宽的比率，大约为 7:1。这是假设所有线程都访问 4 字节（浮点数）数据执行操作。执行此操作所需的总带宽是47,000 = 28,000* GB/s，即达到峰值性能所需的带宽。900 GB/s 将执行限制为 225 GFLOP。这将执行速率限制为峰值的 3.2%（225 GFLOP 是设备的 7,000 GFLOP 峰值的 3.2%）。正如您现在所知，GPU 是一种隐藏延迟的架构，有许多可用于执行的线程，这意味着它在理论上可以容忍长的内存访问延迟。然而，对内存的过多调用可能会导致一些 SMs 空闲，导致一些线程停顿或等待。CUDA 架构提供了其他几种方法，我们可以使用这些方法来访问内存，以解决内存瓶颈问题。

从 CPU 内存到被 SM 用于处理的数据路径在下图中展示。在这里，我们可以看到数据元素在到达 SM 核心进行计算之前的旅程。每个内存带宽的数量级都不同，访问它们的延迟也不同：

在上图中，我们可以看到从 CPU 到达寄存器的数据路径，最终计算是由 ALU/核心完成的。

下图显示了最新 GPU 架构中存在的不同类型的内存层次结构。每种内存可能具有不同的大小、延迟、吞吐量和应用程序开发人员的可见性：

上图显示了最新 GPU 架构中存在的不同类型的内存及其在硬件中的位置。

在本章中，您将学习如何最佳地利用不同类型的 GPU 内存。我们还将研究 GPU 统一内存的最新特性，这使得程序员的生活变得更简单。本章将详细介绍以下内存主题：

全局内存/设备内存
共享内存
只读数据/缓存
固定内存
统一内存

但在我们查看内存层次结构之前，我们将遵循优化周期，如下所示：

步骤 1：分析
步骤 2：并行化
步骤 3：优化

对应用程序的分析要求我们不仅要了解我们应用程序的特性，还要了解它在 GPU 上的有效运行方式。为此，我们将首先向您介绍 Visual Profiler，然后再进入 GPU 内存。由于我们在这里使用了一些最新的 CUDA 功能，请在继续本章之前阅读以下部分。

技术要求

本章需要一台带有现代 NVIDIA GPU（Pascal 架构或更高版本）的 Linux PC，以及安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（10.0 或更高版本）。如果您不确定您的 GPU 架构，请访问 NVIDIA GPU 网站developer.nvidia.com/cuda-gpus进行确认。本章的代码也可以在 GitHub 上找到github.com/PacktPublishing/Learn-CUDA-Programming。

本章的示例代码示例是使用 CUDA Toolkit 的 10.1 版本开发和测试的。但是，建议使用最新的 CUDA 版本或更高版本。

在下一节中，我们将向您介绍 Visual Profiler，它将帮助我们分析我们的应用程序。我们还将看一下它在 GPU 上的运行情况。

NVIDIA Visual Profiler

为了了解不同内存层次结构的有效利用，重要的是在运行时分析应用程序的特性。分析器是非常方便的工具，可以测量和显示不同的指标，帮助我们分析内存、SM、核心和其他资源的使用方式。 NVIDIA 决定提供一个 API，供分析器工具的开发人员用于连接到 CUDA 应用程序，随着时间的推移，一些分析工具已经发展出来，如 TAU 性能系统、Vampir Trace 和 HPC Toolkit。所有这些工具都利用CUDA 分析器工具接口（CUPTI）为 CUDA 应用程序提供分析信息。

NVIDIA 本身开发并维护作为 CUDA Toolkit 的一部分提供的分析工具。本章使用这两个分析工具（NVPROF 和 NVVP）来演示不同内存类型的有效使用，并不是分析工具的指南。

我们将使用 NVPROF 或 NVVP 来演示 CUDA 应用程序的特性。NVPROF 是一个命令行工具，而nvvp具有可视化界面。nvvp有两种格式，一种是独立版本，另一种是集成在 Nsisght Eclipse 中的版本。

我们将广泛使用的 NVVP 分析器窗口如下所示：

这是在 macOS 上拍摄的 NVVP 9.0 版本窗口快照。

窗口中有四个视图可用：时间轴、指南、分析结果和摘要。时间轴视图显示了随时间发生的 CPU 和 GPU 活动。Visual Profiler 显示了 CUDA 编程模型的内存层次结构的摘要视图。分析视图显示了分析结果。Visual Profiler 提供了两种分析模式：

**引导分析：**顾名思义，它通过逐步方法指导开发人员了解关键性能限制器。我们建议初学者在成为了解不同指标的专家之前先使用此模式，然后再转到无引导模式。
**无引导分析：**开发人员必须手动查看此模式下的结果，以了解性能限制器。

CUDA Toolkit 提供了两个 GPU 应用程序性能分析工具，NVIDIA Profiler（NVPROF）和NVIDIA Visual Profiler（NVVP）。为了获得性能限制器信息，我们需要进行两种类型的分析：时间线分析和度量分析。此代码可在 02_memory_overview/04_sgemm 中访问。分析命令可以执行如下：

$ nvcc -o sgemm sgemm.cu
$ nvprof -o sgemm.nvvp ./sgemm
$ nvprof --analysis-metrics -o sgemm-analysis.nvvp ./sgemm

让我们打开 Visual Profiler。如果你使用的是 Linux 或 OSX，你可以在终端中执行 nvvp。或者，你可以从安装了 CUDA Toolkit 的二进制文件中找到 nvvp 可执行文件。如果你使用的是 Windows，你可以使用 Windows 搜索框执行该工具，命令为 nvvp。

要打开两个分析数据，我们将使用“文件”|“导入...”菜单，如下所示：

然后，我们将继续点击底部的“下一步”按钮：

。

我们的 CUDA 应用程序使用一个进程。因此，让我们继续点击底部的“下一步”按钮：

现在，让我们将收集的分析数据放入 Visual Profiler 中。以下截图显示了一个示例。通过右侧的“浏览...”按钮，将时间线数据放入第二个文本框。然后，以相同的方式将度量分析数据放入下一个文本框中：

有关性能分析工具的详细使用，请参阅 CUDA Profiling 指南，该指南作为 CUDA Toolkit 的一部分提供（相应的网页链接为 docs.nvidia.com/cuda/profiler-users-guide/index.html）。

在基于 Windows 的系统中，在安装了 CUDA Toolkit 后，你可以从“开始”菜单中启动 Visual Profiler。在具有 X11 转发的 Linux 系统中，你可以通过运行 nvvp 命令来启动 Visual Profiler，nvvp 代表 NVIDIA Visual Profiler：

$ ./nvvp

既然我们现在对将要使用的分析工具有了一个公平的理解，让我们进入第一个也是绝对最关键的 GPU 内存——全局内存/设备内存。

全局内存/设备内存

本节将详细介绍如何使用全局内存，也称为设备内存。在本节中，我们还将讨论如何高效地将数据从全局内存加载/存储到缓存中。由于全局内存是一个暂存区，所有数据都从 CPU 内存中复制到这里，因此必须充分利用这种内存。全局内存或设备内存对于内核中的所有线程都是可见的。这种内存也对 CPU 可见。

程序员使用 cudaMalloc 和 cudaFree 显式地管理分配和释放。数据使用 cudaMalloc 分配，并声明为 __device__。全局内存是从 CPU 使用 cudaMemcpy API 传输的所有内存的默认暂存区。

全局内存上的矢量加法

我们在第一章中使用的矢量加法示例演示了全局内存的使用。让我们再次查看代码片段，并尝试理解全局内存的使用方式：

__global__ void device_add(int *a, int *b, int *c) {
     int index = threadIdx.x + blockIdx.x * blockDim.x;
     c[index] = a[index] + b[index];
}
int main (void) {
...
    // Alloc space for device copies of a, b, c
    cudaMalloc((void **)&d_a, size);
    cudaMalloc((void **)&d_b, size);
    cudaMalloc((void **)&d_c, size);
...

   // Free space allocated for device copies
   cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
...

}

cudaMalloc 在设备内存上分配数据。内核中的参数指针（a、b 和 c）指向这个设备内存。我们使用 cudaFree API 释放这个内存。正如你所看到的，块中的所有线程都可以在内核中访问这个内存。

此代码可在 02_memory_overview/01_vector_addition 中访问。要编译此代码，可以使用以下命令：

$ nvcc -o vec_addition ./vector_addition_gpu_thread_block.cu

这是一个使用全局内存的简单示例。在下一节中，我们将看看如何最优地访问数据。

合并与未合并的全局内存访问

为了有效使用全局内存，了解 CUDA 编程模型中 warp 的概念是非常重要的，这是我们到目前为止忽略的。warp 是 SM 中的线程调度/执行单位。一旦一个块被分配给一个 SM，它被划分为一个 32 个线程的单位，称为warp。这是 CUDA 编程中的基本执行单位。

为了演示 warp 的概念，让我们看一个例子。如果两个块被分配给一个 SM，每个块有 128 个线程，那么块内的 warp 数量是128/32 = 4个 warp，SM 上的总 warp 数量是4 * 2 = 8个 warp。以下图表显示了 CUDA 块如何在 GPU SM 上被划分和调度：

块和 warp 在 SM 和其核心上的调度更多地是与体系结构相关的，对于 Kepler、Pascal 和最新的 Volta 等不同的架构，情况会有所不同。目前，我们可以忽略调度的完整性。在所有可用的 warp 中，具有下一条指令所需操作数的 warp 变得可以执行。根据运行 CUDA 程序的 GPU 的调度策略，选择要执行的 warp。当被选择时，warp 中的所有线程执行相同的指令。CUDA 遵循单指令，多线程（SIMT）模型，也就是说，warp 中的所有线程在同一时间实例中获取和执行相同的指令。为了最大程度地利用全局内存的访问，访问应该合并。合并和未合并之间的区别如下：

合并的全局内存访问： 顺序内存访问是相邻的。
未合并的全局内存访问： 顺序内存访问不是相邻的。

以下图表更详细地展示了这种访问模式的示例。图表的左侧显示了合并访问，其中 warp 中的线程访问相邻数据，因此导致了一个 32 位宽的操作和 1 次缓存未命中。图表的右侧显示了一种情况，即 warp 内的线程访问是随机的，可能导致调用 32 次单个宽度的操作，因此可能有 32 次缓存未命中，这是最坏的情况：

为了进一步理解这个概念，我们需要了解数据如何通过缓存行从全局内存到达。

情况 1： warp 请求 32 个对齐的、4 个连续的字节

地址落在 1 个缓存行内和一个 32 位宽的操作内。总线利用率为 100%，也就是说，我们利用从全局内存中获取的所有数据到缓存中，并没有浪费任何带宽。如下图所示：

上图显示了合并的访问，导致了总线的最佳利用。

情况 2： warp 请求 32 个分散的 4 字节单词

虽然 warp 需要 128 字节，但在未命中时执行了 32 次单个宽度的获取，导致32 * 128字节在总线上移动。如下图所示，总线利用率实际上低于 1%：

上图显示了未合并的访问，导致了总线带宽的浪费。

正如我们在前面的图表中看到的，warp 内的线程如何从全局内存中访问数据非常重要。为了最大程度地利用全局内存，改善合并是非常重要的。有多种可以使用的策略。其中一种策略是改变数据布局以改善局部性。让我们看一个例子。将滤波器应用于图像或将掩模应用于图像的计算机视觉算法需要将图像存储到数据结构中。当开发人员声明图像类型时，有两种选择。

以下代码片段使用Coefficients_SOA数据结构以数组格式存储数据。Coefficients_SOA结构存储与图像相关的数据，如 RGB、色调和饱和度值：

//Data structure representing an image stored in Structure of Array Format
struct Coefficients_SOA {
 int r;
 int b;
 int g;
 int hue;
 int saturation;
 int maxVal;
 int minVal;
 int finalVal;
};

以下图表显示了关于Coefficients_SOA存储数据的数据布局，以及在内核中由不同线程访问数据的情况：

通过这样做，我们可以看到 AOS 数据结构的使用导致了不连续的全局内存访问。

同样的图像可以以数组结构格式存储，如下面的代码片段所示：

//Data structure representing an image stored in Array of Structure Format
struct Coefficients_AOS {
 int* r;
 int* b;
 int* g;
 int* hue;
 int* saturation;
 int* maxVal;
 int* minVal;
 int* finalVal;
};

以下图表显示了关于Coefficients_AOS存储数据的数据布局，以及在内核中由不同线程访问数据的情况：

通过这样做，我们可以看到使用 SOA 数据结构导致了不连续的全局内存访问。

虽然 CPU 上的顺序代码更喜欢 AOS 以提高缓存效率，但在单指令多线程（SIMT）模型（如 CUDA）中，SOA 更受欢迎，以提高执行和内存效率。

让我们尝试通过使用分析器来分析这一方面。根据以下步骤配置你的环境：

准备好你的 GPU 应用程序。例如，我们将使用两段代码来演示全局内存的有效使用。aos_soa.cu文件包含了使用 AOS 数据结构的朴素实现，而aos_soa_solved.cu则使用了 SOA 数据结构，可以有效地利用全局内存。这段代码可以在02_memory_overview/02_aos_soa中找到。
使用nvcc编译器编译你的应用程序，然后使用nvprof编译器对其进行分析。以下命令是对此的一个nvcc命令的示例。然后我们使用nvprof命令对应用程序进行分析。还传递了--analysis-metrics标志，以便我们可以获得内核的指标。
生成的分析文件，即aos_soa.prof和aos_soa_solved.prof，然后加载到 NVIDIA Visual Profiler 中。用户需要从“文件|打开”菜单中加载分析输出。此外，不要忘记在文件名选项中选择“所有文件”。

$ nvcc -o aos_soa ./aos_soa.cu
$ nvcc -o aos_soa_solved ./aos_soa_solved.cu
$ nvprof --analysis-metrics --export-profile aos_soa.prof ./aos_soa
$ nvprof --analysis-metrics --export-profile aos_soa_solved.prof ./aos_soa_solved

以下是分析输出的屏幕截图。这是一个使用 AOS 数据结构的朴素实现：

前面的图表显示了在引导分析模式下分析器的输出。

你将看到的第一件事是分析器明确指出应用程序受到内存限制。正如你所看到的，分析器不仅显示指标，还分析了这些指标的含义。在这个例子中，由于我们使用 AOS，分析器明确指出访问模式不高效。但编译器是如何得出这个结论的呢？让我们看一下下面的屏幕截图，它提供了更多的细节：

正如我们所看到的，它清楚地说明了访问数据的理想事务数为四，而实际运行时进行了 32 次事务/访问。

当我们将数据结构从 AOS 更改为 SOA 时，瓶颈得到了解决。当你运行aos_soa_solved可执行文件时，你会发现内核时间减少了，这对我们的计时来说是一个改进。在 V100 16 GB 卡上，时间从 104 微秒减少到 47 微秒，这是一个2.2x的加速因子。分析输出aos_soa_solved.prof将显示内核仍然受到内存限制，这是非常明显的，因为我们读写的内存数据比进行计算时要多。

内存吞吐量分析

对于应用程序开发人员来说，了解应用程序的内存吞吐量非常重要。这可以通过两种方式来定义：

从应用程序的角度来看： 计算应用程序请求的字节数
从硬件的角度来看： 计算硬件传输的字节数

这两个数字完全不同。其中许多原因包括未协调的访问导致未利用所有事务字节，共享内存银行冲突等。我们应该从内存角度使用两个方面来分析应用程序：

地址模式：在实际代码中确定访问模式是非常困难的，因此使用诸如性能分析器之类的工具变得非常重要。性能分析器显示的指标，如全局内存效率和每次访问的 L1/L2 事务，需要仔细观察。
**飞行中的并发访问数量：**由于 GPU 是一种隐藏延迟的架构，饱和内存带宽变得非常重要。但是确定并发访问的数量通常是不够的。此外，从硬件的角度来看，吞吐量与理论值相比要不同得多。

下图演示了每个 SM 中飞行的约 6KB 数据可以达到 Volta 架构峰值带宽的 90%。在以前的一代架构上进行相同的实验会得到不同的图表。一般来说，建议了解特定架构的 GPU 内存特性，以便从硬件中获得最佳性能：

本节为我们提供了全局内存的示例用法以及如何以最佳方式利用它。有时，全局内存的协调数据访问很困难（例如，在 CFD 领域，对于非结构化网格，相邻单元格的数据可能不会相邻存储在内存中）。为了解决这样的问题或减少对性能的影响，我们需要利用另一种形式的内存，称为共享内存。

共享内存

共享内存一直在 CUDA 内存层次结构中扮演着重要角色，被称为用户管理的缓存。这为用户提供了一种机制，可以以协调的方式从全局内存中读取/写入数据并将其存储在内存中，这类似于缓存但可以由用户控制。在本节中，我们不仅将介绍利用共享内存的步骤，还将讨论如何有效地从共享内存中加载/存储数据以及它在银行中的内部排列。共享内存只对同一块中的线程可见。块中的所有线程看到共享变量的相同版本。

共享内存具有类似于 CPU 缓存的好处；然而，CPU 缓存无法明确管理，而共享内存可以。共享内存的延迟比全局内存低一个数量级，带宽比全局内存高一个数量级。但共享内存的关键用途来自于块内线程可以共享内存访问。CUDA 程序员可以使用共享变量来保存在内核执行阶段中多次重复使用的数据。此外，由于同一块内的线程可以共享结果，这有助于避免冗余计算。直到 9.0 版本，CUDA Toolkit 没有提供可靠的通信机制来在不同块的线程之间进行通信。我们将在后续章节中更详细地介绍 CUDA 9.0 通信机制。目前，我们将假设在 CUDA 中只能通过使用共享内存来实现线程之间的通信。

共享内存上的矩阵转置

用于演示共享内存的最原始的例子之一是矩阵转置。矩阵转置是一个内存绑定的操作。以下代码片段使用matrix_transpose_naive内核，展示了矩阵转置内核的示例实现：

__global__ void matrix_transpose_naive(int *input, int *output) {
     int indexX = threadIdx.x + blockIdx.x * blockDim.x;
     int indexY = threadIdx.y + blockIdx.y * blockDim.y;
     int index = indexY * N + indexX;
     int transposedIndex = indexX * N + indexY;
     output[index] = input[transposedIndex];
}

上述代码展示了使用全局内存的矩阵转置的朴素实现。如果以朴素的方式实现，这将导致在读取矩阵或写入矩阵时出现未协调的访问。在 V100 PCIe 16 GB 卡上，内核的执行时间约为 60 微秒。

根据以下步骤配置您的环境：

准备您的 GPU 应用程序。此代码可以在02_memory_overview/02_matrix_transpose中找到。
使用nvcc编译器编译您的应用程序，然后使用nvprof编译器对其进行分析。以下命令是对此的nvcc命令的一个示例。然后，我们使用nvprof命令对应用程序进行分析。还传递了--analysis-metrics标志以获取内核的指标。
生成的配置文件，即matrix_transpose.prof，然后加载到 NVIDIA Visual Profiler 中。用户需要从“文件|打开”菜单中加载分析输出。还要记得选择“所有文件”作为文件名选项的一部分：

$ nvcc -o matrix_transpose ./matrix_transpose.cu
$ nvcc -o conflict_solved ./conflict_solved.cu
$ nvprof --analysis-metrics --export-profile matrix_transpose.prof ./matrix_transpose
$ nvprof --analysis-metrics --export-profile conflict_solved.prof ./conflict_solved

以下截图显示了性能分析的输出。输出清楚地表明对全局内存的访问是未协调的，这是需要解决的关键指标，以便我们可以提高性能：

解决这个问题的一种方法是利用高带宽和低延迟的内存，比如共享内存。这里的诀窍是以协调的方式从全局内存读取和写入。在这里，对共享内存的读取或写入可以是未协调的模式。使用共享内存会带来更好的性能，时间缩短到 21 微秒，这是 3 倍的加速时间：

__global__ void matrix_transpose_shared(int *input, int *output) {

    __shared__ int sharedMemory [BLOCK_SIZE] [BLOCK_SIZE];

    //global index
     int indexX = threadIdx.x + blockIdx.x * blockDim.x;
     int indexY = threadIdx.y + blockIdx.y * blockDim.y;

    //transposed global memory index
     int tindexX = threadIdx.x + blockIdx.y * blockDim.x;
     int tindexY = threadIdx.y + blockIdx.x * blockDim.y;

    //local index
     int localIndexX = threadIdx.x;
     int localIndexY = threadIdx.y;
     int index = indexY * N + indexX;
     int transposedIndex = tindexY * N + tindexX;

    //transposed the matrix in shared memory. 
    // Global memory is read in coalesced fashion
     sharedMemory[localIndexX][localIndexY] = input[index];
     __syncthreads();

    //output written in global memory in coalesed fashion.
     output[transposedIndex] = sharedMemory[localIndexY][localIndexX];
}

上述代码片段显示了使用共享内存的矩阵转置的实现。全局内存读取/写入是协调的，而转置发生在共享内存中。

银行冲突及其对共享内存的影响

与使用全局内存相比的良好加速并不一定意味着我们有效地使用了共享内存。如果我们转换从引导分析到未引导分析的分析器输出，即matrix_transpose.prof，我们将看到共享内存访问模式显示出对齐问题，如下截图所示：

我们可以看到分析器显示了共享内存的非最佳使用，这是银行冲突的一个迹象。

为了有效地理解这个对齐问题，重要的是要理解bank的概念。共享内存被组织成 bank 以实现更高的带宽。每个 bank 可以在一个周期内服务一个地址。内存可以为它有的 bank 提供多个同时访问。Volta GPU 有 32 个 bank，每个 bank 宽度为 4 字节。当一个数组存储在共享内存中时，相邻的 4 字节单词会进入连续的 bank，如下图所示：

上述图中的逻辑视图显示了数据在共享内存中的存储方式。

warp 内的线程对 bank 的多个同时访问会导致 bank 冲突。换句话说，当 warp 内的两个或多个线程访问同一个 bank 中的不同 4 字节单词时，就会发生 bank 冲突。从逻辑上讲，这是当两个或多个线程访问同一个 bank 中的不同行时。以下图示例展示了不同n-way bank 冲突的例子。最坏的情况是 32-way 冲突 | 31 次重播 - 每次重播都会增加一些延迟：

上述情景显示了来自同一个 warp 的线程访问驻留在不同 bank 中的相邻 4 字节元素，导致没有 bank 冲突。看一下下图：

这是另一个没有银行冲突的场景，同一 warp 的线程访问随机的 4 字节元素，这些元素位于不同的银行中，因此没有银行冲突。由于共享内存中的 2 路银行冲突，顺序访问如下图所示：

前面的图表显示了一个场景，其中来自同一 warp 的线程 T0 和 T1 访问同一银行中的 4 字节元素，因此导致了 2 路银行冲突。

在前面的矩阵转置示例中，我们利用了共享内存来获得更好的性能。然而，我们可以看到 32 路银行冲突。为了解决这个问题，可以使用一种称为填充的简单技术。所有这些都是在共享内存中填充一个虚拟的，即一个额外的列，这样线程就可以访问不同的银行，从而获得更好的性能：

__global__ void matrix_transpose_shared(int *input, int *output) {

     __shared__ int sharedMemory [BLOCK_SIZE] [BLOCK_SIZE + 1];

    //global index
     int indexX = threadIdx.x + blockIdx.x * blockDim.x;
     int indexY = threadIdx.y + blockIdx.y * blockDim.y;

    //transposed index
     int tindexX = threadIdx.x + blockIdx.y * blockDim.x;
     int tindexY = threadIdx.y + blockIdx.x * blockDim.y;
     int localIndexX = threadIdx.x;
     int localIndexY = threadIdx.y;
     int index = indexY * N + indexX;
     int transposedIndex = tindexY * N + tindexX;

    //reading from global memory in coalesed manner 
    // and performing tanspose in shared memory
     sharedMemory[localIndexX][localIndexY] = input[index];

    __syncthreads();

    //writing into global memory in coalesed fashion 
    // via transposed data in shared memory
     output[transposedIndex] = sharedMemory[localIndexY][localIndexX];
}

前面的代码片段中，我们使用了matrix_transpose_shared内核，展示了填充的概念，这样可以消除银行冲突，从而更好地利用共享内存带宽。像往常一样，运行代码并借助可视化分析器验证这种行为。通过这些更改，您应该看到内核的时间减少到 13 微秒，这进一步提高了 60%的速度。

在本节中，我们看到了如何最大限度地利用共享内存，它提供了读写访问作为一个临时存储器。但有时，数据只是只读输入，不需要写访问。在这种情况下，GPU 提供了一种称为纹理内存的最佳内存。我们将在下一章中详细介绍这一点，以及它为开发人员提供的其他优势。我们将在下一节中介绍只读数据。

只读数据/缓存

根据内存名称，只读缓存适合存储只读数据，并且在内核执行过程中不会发生更改。该缓存针对此目的进行了优化，并且根据 GPU 架构，释放并减少了其他缓存的负载，从而提高了性能。在本节中，我们将详细介绍如何利用只读缓存，以及如何使用图像处理代码示例进行图像调整。

GPU 中的所有线程都可以看到只读数据。对于 GPU 来说，这些数据被标记为只读，这意味着对这些数据的任何更改都会导致内核中的未指定行为。另一方面，CPU 对这些数据具有读写访问权限。

传统上，这个缓存也被称为纹理缓存。虽然用户可以显式调用纹理 API 来利用只读缓存，但是在最新的 GPU 架构中，开发人员可以在不显式使用 CUDA 纹理 API 的情况下利用这个缓存。使用最新的 CUDA 版本和像 Volta 这样的 GPU，标记为const __restrict__的内核指针参数被视为只读数据，通过只读缓存数据路径传输。开发人员还可以通过__ldg内在函数强制加载这个缓存。

只读数据在算法要求整个 warp 读取相同地址/数据时理想地使用，这主要导致每个时钟周期对所有请求数据的线程进行广播。纹理缓存针对 2D 和 3D 局部性进行了优化。随着线程成为同一 warp 的一部分，从具有 2D 和 3D 局部性的纹理地址读取数据往往会获得更好的性能。纹理在要求随机内存访问的应用程序中已被证明是有用的，特别是在 Volta 架构之前的显卡中。

纹理支持双线性和三线性插值，这对于图像处理算法如缩放图像特别有用。

下图显示了一个 warp 内的线程访问空间位置在 2D 空间中的元素的示例。纹理适用于这类工作负载：

现在，让我们看一个关于缩放的小型实际算法，以演示纹理内存的使用。

计算机视觉-使用纹理内存进行图像缩放

我们将使用图像缩放作为示例来演示纹理内存的使用。图像缩放的示例如下截图所示：

图像缩放需要在 2 维中插值图像像素。纹理提供了这两个功能（插值和对 2D 局部性的高效访问），如果直接通过全局内存访问，将导致内存访问不连续。

根据以下步骤配置您的环境：

准备您的 GPU 应用程序。此代码可以在02_memory_overview/03_image_scaling中找到。
使用nvcc编译器编译您的应用程序，使用以下命令：

$nvcc -c scrImagePgmPpmPackage.cpp 
$nvcc -c image_scaling.cu
$nvcc -o image_scaling image_scaling.o scrImagePgmPpmPackage.o

scrImagePgmPpmPackage.cpp文件包含了读取和写入.pgm扩展名图像的源代码。纹理代码位于image_scaling.cu中。

用户可以使用 IrfanView（www.irfanview.com/main_download_engl.htm）等查看器来查看pgm文件，这些查看器是免费使用的。

主要有四个步骤是必需的，以便我们可以使用纹理内存：

声明纹理内存。
将纹理内存绑定到纹理引用。
在 CUDA 内核中使用纹理引用读取纹理内存。
从纹理引用中解绑纹理内存。

以下代码片段显示了我们可以使用的四个步骤来使用纹理内存。从 Kepler GPU 架构和 CUDA 5.0 开始，引入了一项名为无绑定纹理的新功能。这暴露了纹理对象，它基本上是一个可以传递给 CUDA 内核的 C++对象。它们被称为无绑定，因为它们不需要手动绑定/解绑，这是早期 GPU 和 CUDA 版本的情况。纹理对象使用cudaTextureObject_t类 API 声明。现在让我们通过这些步骤：

首先，声明纹理内存：

texture<unsigned char, 2, cudaReadModeElementType> tex;

创建一个通道描述，我们在链接到纹理时将使用：

cudaArray* cu_array;
cudaChannelFormatKind kind = cudaChannelFormatKindUnsigned;
cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(8, 0, 0, 0, kind);

然后，指定纹理对象参数：

struct cudaTextureDesc texDesc;
memset(&texDesc, 0, sizeof(texDesc)); 
//set the memory to zero
texDesc.addressMode[0] = cudaAddressModeClamp; 
// setting the x dimension addressmode to Clamp
texDesc.addressMode[1] = cudaAddressModeClamp; 
//Setting y dimension addressmode to Clamp
texDesc.filterMode = cudaFilterModePoint; 
// Filter mode set to Point
texDesc.readMode = cudaReadModeElementType; 
// Reading element type and not interpolated
texDesc.normalizedCoords = 0;

接下来，在 CUDA 内核中从纹理引用中读取纹理内存：

imageScaledData[index] = tex2D<unsigned char>(texObj,(float)(tidX*scale_factor),(float)(tidY*scale_factor));

最后，销毁纹理对象：

cudaDestroyTextureObject(texObj);

纹理内存的重要方面，它们像配置一样由开发人员设置，如下所示：

**纹理维度：**这定义了纹理是作为 1D、2D 还是 3D 数组寻址。纹理中的元素也称为纹素。深度、宽度和高度也被设置以定义每个维度。请注意，每个 GPU 架构都定义了可接受的每个维度的最大尺寸。
**纹理类型：**这定义了基本整数或浮点纹素的大小。
**纹理读取模式：**纹理的读取模式定义了元素的读取方式。它们可以以NormalizedFloat或ModeElement格式读取。标准化浮点模式期望在[0.0 1.0]和[-1.0 1.0]范围内的索引，对于无符号整数和有符号整数类型。
**纹理寻址模式：**纹理的一个独特特性是它如何寻址超出范围的访问。这听起来可能很不寻常，但实际上在许多图像算法中非常常见。例如，如果您正在通过平均相邻像素来应用插值，那么边界像素的行为应该是什么？纹理为开发人员提供了这个选项，以便他们可以选择将超出范围视为夹紧、包裹或镜像。在调整大小的示例中，我们已将其设置为夹紧模式，这基本上意味着超出范围的访问被夹紧到边界。
**纹理过滤模式：**设置模式定义在获取纹理时如何计算返回值。支持两种类型的过滤模式：cudaFilterModePoint和cudaFilterModeLinear。当设置为线性模式时，可以进行插值（1D 的简单线性，2D 的双线性和 3D 的三线性）。仅当返回类型为浮点类型时，线性模式才有效。另一方面，ModePoint不执行插值，而是返回最近坐标的纹素。

在本节中引入纹理内存的关键目的是为您提供其用法的示例，并向您展示纹理内存的有用之处。它提供了不同配置参数的良好概述。有关更多信息，请参阅 CUDA API 指南（docs.nvidia.com/cuda/cuda-runtime-api/index.html）。

在本节中，我们通过示例描述了使用纹理内存的目的。在下一节中，我们将研究 GPU 内存中最快（最低延迟）的可用内存（寄存器）。与 CPU 相比，GPU 中富裕地存在这种内存。

GPU 中的寄存器

CPU 和 GPU 架构之间的一个基本区别是 GPU 中寄存器的丰富性相对于 CPU。这有助于线程将大部分数据保存在寄存器中，从而减少上下文切换的延迟。因此，使这种内存达到最佳状态也很重要。

寄存器的范围是单个线程。在 GRID 中为所有启动的线程创建变量的私有副本。每个线程都可以访问其变量的私有副本，而其他线程的私有变量则无法访问。例如，如果使用 1,000 个线程启动内核，那么作为线程范围的变量将获得其自己的变量副本。

作为内核的一部分声明的局部变量存储在寄存器中。中间值也存储在寄存器中。每个 SM 都有一组固定的寄存器。在编译期间，编译器（nvcc）尝试找到每个线程的最佳寄存器数量。如果寄存器数量不足，通常发生在 CUDA 内核较大且具有许多局部变量和中间计算时，数据将被推送到本地内存，该内存可以位于 L1/L2 缓存中，甚至更低的内存层次结构中，例如全局内存。这也被称为寄存器溢出。每个线程的寄存器数量在 SM 上可以激活多少个块和线程方面起着重要作用。这个概念在下一章中有详细介绍，该章节专门讨论了占用率。一般来说，建议不要声明大量不必要的局部变量。如果寄存器限制了可以在 SM 上调度的线程数量，那么开发人员应该考虑通过将内核拆分为两个或更多个（如果可能）来重新构建代码。

作为vecAdd内核的一部分声明的变量存储在寄存器内存中。传递给内核的参数，即A、B和C，指向全局内存，但变量本身存储在基于 GPU 架构的共享内存或寄存器中。以下图显示了 CUDA 内存层次结构和不同变量类型的默认位置：

到目前为止，我们已经看到了关键内存层次结构（全局、纹理、共享和寄存器）的目的和最佳用法。在下一节中，我们将看一些可以提高应用程序性能并增加开发人员在编写 CUDA 程序时的生产力的 GPU 内存的优化和特性。

固定内存

现在是时候回想一下数据的路径，即从 CPU 内存到 GPU 寄存器，最终由 GPU 核心用于计算。尽管 GPU 具有更高的计算性能和更高的内存带宽，但由于 CPU 内存和 GPU 内存之间的传输，应用程序获得的加速效果可能会变得规范化。数据传输是通过总线/链接/协议进行的，例如 PCIe（对于英特尔和 AMD 等 CPU 架构）或 NVLink（对于 OpenPower Foundation 的power等 CPU 架构）。

为了克服这些瓶颈，建议采用以下技巧/指南：

首先，建议在可能的情况下尽量减少主机和设备之间传输的数据量。这甚至可能意味着将顺序代码的一部分作为 GPU 上的内核运行，与在主机 CPU 上顺序运行相比，几乎没有或没有加速。
其次，通过利用固定内存，重要的是在主机和设备之间实现更高的带宽。
建议将小的传输批量成一个大的传输。这有助于减少调用数据传输 CUDA API 所涉及的延迟，根据系统配置，这可能从几微秒到几毫秒不等。
最后，应用程序可以利用异步数据传输来重叠内核执行和数据传输。

我们将在本节中更详细地介绍固定内存传输。异步传输将在第四章中更详细地介绍，内核执行模型和优化策略，在那里我们将使用一个称为 CUDA 流的概念。

带宽测试-固定与分页

默认情况下，称为malloc()的内存分配 API 分配的是可分页的内存类型。这意味着，如果需要，作为页面映射的内存可以被其他应用程序或操作系统本身交换出去。因此，大多数设备，包括 GPU 和其他设备（如 InfiniBand 等），也位于 PCIe 总线上，都希望在传输之前将内存固定。默认情况下，GPU 将不访问可分页内存。因此，当调用内存传输时，CUDA 驱动程序会分配临时固定内存，将数据从默认可分页内存复制到此临时固定内存，然后通过设备内存控制器（DMA）将其传输到设备。

这个额外的步骤不仅会增加延迟，还有可能会将请求的页面传输到已经被交换并需要重新传输到 GPU 内存的 GPU 内存。

为了了解使用固定内存的影响，让我们尝试编译和运行一段示例代码。这已作为 CUDA 示例的一部分提供。根据以下步骤配置您的环境：

准备您的 GPU 应用程序。此代码位于<CUDA_SAMPLES_DIR>/1_Utilities/bandwidthTest中。
使用make命令编译您的应用程序。
以分页和固定两种模式运行可执行文件，如下所示：

$make
$./bandwidthTest --mode=shmoo --csv --memory=pageable > pageable.csv
$./bandwidthTest --mode=shmoo --csv --memory=pinned >  pinned.csv

注意，CUDA_SAMPLES_DIR是 CUDA 安装所在目录的路径。

正如我们所看到的，与之前的代码相比，关键的变化是我们迄今为止编写的是数据分配 API。以下代码片段显示了使用cudaMallocHost API 而不是malloc来分配内存：

cudaError_t status = cudaMallocHost((void**)&h_aPinned, bytes);
if (status != cudaSuccess)
 printf("Error allocating pinned host memory\n");

cudaMallocHost API 使内存成为固定内存而不是分页内存。虽然分配 API 已更改，但我们仍然可以使用相同的数据传输 API，即cudaMemcpy()。现在，重要的问题是，*固定内存是什么，为什么它提供更好的带宽？*我们将在下一节中介绍这个问题。

性能的影响可以从带宽测试的输出中看出。我们已经将结果绘制成图表，以便您可以轻松理解影响。x轴显示了以 KB 为单位传输的数据，而y轴显示了以 MB/sec 为单位的实现带宽。

第一个图是主机到设备的传输，而第二个图是设备到主机的传输。您将看到的第一件事是可实现的最大带宽约为 12 GB/sec。PCIe Gen3 的理论带宽为 16 GB/sec，但实际可实现的范围在 12 GB/sec 左右。可实现的带宽高度取决于系统（主板、CPU、PCIe 拓扑等）：

如您所见，对于固定内存，在较小的传输大小时带宽始终更高，而在可分页内存中，随着数据大小的增加，带宽变得相等，因为驱动程序和 DMA 引擎开始通过应用诸如重叠的概念来优化传输。尽管建议使用固定内存，但过度使用也有缺点。为应用程序分配整个系统内存作为固定内存可能会降低整体系统性能。这是因为它会占用其他应用程序和操作系统任务可用的页面。应该固定的正确大小非常依赖于应用程序和系统，并且没有通用的公式可用。我们能做的最好的事情是在可用系统上测试应用程序并选择最佳的性能参数。

此外，重要的是要了解新的互连技术，如 NVLink，为受这些数据传输限制的应用程序提供了更高的带宽和更低的延迟。目前，CPU 和 GPU 之间的 NVLink 仅与 Power CPU 一起提供。

在本节中，我们将看看如何提高 CPU 和 GPU 之间的数据传输速度。现在我们将继续利用 CUDA 的一个新特性，称为统一内存，这有助于提高编写 CUDA 程序的开发人员的生产力。

统一内存

随着每一次新的 CUDA 和 GPU 架构发布，都会添加新的功能。这些新功能提供了更高的性能和更便捷的编程，或者允许开发人员实现新的算法，否则无法使用 CUDA 在 GPU 上进行移植。从 CUDA 6.0 开始发布的一个重要功能是统一内存，从 Kepler GPU 架构开始实现。在本章中，我们将统一内存称为 UM。

以更简单的话来说，UM 为用户提供了一个单一内存空间的视图，所有 GPU 和 CPU 都可以访问该空间。下图对此进行了说明：

在本节中，我们将介绍如何使用 UM，优化它，并突出利用它的关键优势。与全局内存访问一样，如果以不连续的方式进行，会导致性能不佳，如果未正确使用 UM 功能，也会导致应用程序整体性能下降。我们将采取逐步的方法，从一个简单的程序开始，并在此基础上构建，以便我们可以理解 UM 及其对性能的影响。

让我们尝试编译和运行一些示例代码。根据以下步骤配置您的环境：

准备您的 GPU 应用程序。此代码可以在02_memory_overview/unified_memory中找到。
使用以下nvcc命令编译您的应用程序：

$nvcc -o unified_simple.out unified_memory.cu
$nvcc -o unified_initialized.out unified_memory_initialized.cu
$nvcc -o unified_prefetch.out unified_memory_prefetch.cu
$nvcc -o unified_64align.out unified_memory_64align.cu

请注意，本节中显示的结果是针对 Tesla P100 卡的。当在其他架构（如 Kepler）上运行相同的代码时，预计会产生不同的结果。本节的重点是最新的架构，如 Pascal 和 Volta。

了解统一内存页面分配和传输

让我们从 UM 的朴素实现开始。代码的第一部分unified_memory.cu演示了这个概念的基本用法。代码中的关键更改是使用cudaMallocManaged()API 来分配内存，而不是使用malloc，如下面的代码片段所示：

float *x, *y;
int size = N * sizeof(float);
...
cudaMallocManaged(&x, size);
cudaMallocManaged(&y, size);
...

 for (int ix = 0; ix < N; ix++) {
    x[ix] = rand()%10;
    y[ix] = rand()%20;
  }
...

 add<<<numBlocks, blockSize>>>(x, y, N);

如果我们仔细查看源代码，我们会发现x和y变量只被分配一次并指向统一内存。同一个指针被发送到 GPU 的add<<<>>>()内核，并且在 CPU 中使用for循环进行初始化。这对程序员来说非常简单，因为他们不需要跟踪指针是指向 CPU 内存还是 GPU 内存。但这是否意味着我们能获得良好的性能或传输速度呢？不一定，所以让我们尝试通过对这段代码进行性能分析来深入了解，如下面的屏幕截图所示：

我们使用以下命令来获取性能分析输出：

$ nvprof ./unified_simple.out

正如预期的那样，大部分时间都花在了add<<<>>>内核上。让我们尝试理论计算带宽。我们将使用以下公式来计算带宽：

带宽 = 字节/秒 = (3 * 4,194,304 字节 * 1e-9 字节/GB) / 2.6205e-3 秒 = 5 GB/s

如您所见，P100 提供了 720 GB/s 的理论带宽，而我们只能实现 5 GB/s，这实在是太差了。您可能想知道为什么我们只计算内存带宽。这是因为应用程序受内存限制，因为它完成了三次内存操作和仅一次加法。因此，只集中在这个方面是有意义的。

从 Pascal 卡开始，cudaMallocManaged()不再分配物理内存，而是基于首次触摸的基础上分配内存。如果 GPU 首次触摸变量，页面将被分配并映射到 GPU 页表；否则，如果 CPU 首次触摸变量，它将被分配并映射到 CPU。在我们的代码中，x和y变量在 CPU 中用于初始化。因此，页面被分配给 CPU。在add<<<>>>内核中，当访问这些变量时，会发生页面错误，并且页面迁移的时间被添加到内核时间中。这是内核时间高的根本原因。现在，让我们深入了解页面迁移的步骤。

页面迁移中完成的操作顺序如下：

首先，我们需要在 GPU 和 CPU 上分配新页面（首次触摸）。如果页面不存在并且映射到另一个页面，会发生设备页表页错误。当在 GPU 中访问当前映射到 CPU 内存的page 2中的x时，会发生页面错误。请看下图：

接下来，CPU 上的旧页面被取消映射，如下图所示：

接下来，数据从 CPU 复制到 GPU，如下图所示：

最后，新页面在 GPU 上映射，旧页面在 CPU 上释放，如下图所示：

GPU 中的转换后备缓冲器（TLB）与 CPU 中的类似，执行从物理地址到虚拟地址的地址转换。当发生页面错误时，相应 SM 的 TLB 被锁定。这基本上意味着新指令将被暂停，直到执行前面的步骤并最终解锁 TLB。这是为了保持一致性并在 SM 内维护内存视图的一致状态。驱动程序负责删除这些重复项，更新映射并传输页面数据。正如我们之前提到的，所有这些时间都被添加到总体内核时间中。

所以，我们现在知道问题所在。但解决方案是什么呢？为了解决这个问题，我们将采用两种方法：

首先，我们将在 GPU 上创建一个初始化内核，以便在add<<<>>>内核运行期间没有页面错误。然后，我们将通过利用每页的 warp 概念来优化页面错误。
我们将预取数据。

我们将在接下来的部分中介绍这些方法。

使用每页 warp 优化统一内存

让我们从第一种方法开始，即初始化内核。如果你看一下unified_memory_initialized.cu文件中的源代码，我们在那里添加了一个名为init<<<>>>的新内核，如下所示：

__global__ void init(int n, float *x, float *y) {
 int index = threadIdx.x + blockIdx.x * blockDim.x;
 int stride = blockDim.x * gridDim.x;
 for (int i = index; i < n; i += stride) {
   x[i] = 1.0f;
   y[i] = 2.0f;
  }
}

通过在 GPU 本身添加一个初始化数组的内核，页面在init<<<>>>内核中首次被触摸时被分配和映射到 GPU 内存。让我们来看看这段代码的性能分析结果输出，其中显示了初始化内核的性能分析输出：

我们使用以下命令获取了性能分析输出

nvprof ./unified_initialized.out

正如你所看到的，add<<<>>>内核的时间减少到了 18 微秒。这有效地给了我们以下内核带宽：

带宽 = 字节/秒 = (3 * 4,194,304 字节 * 1e-9 字节/GB) / 18.84e-6 秒 = 670 GB/s

这个带宽是你在非统一内存场景中所期望的。正如我们从前面截图中的天真实现中所看到的，性能分析输出中没有主机到设备的行。然而，你可能已经注意到，即使add<<<>>>内核的时间已经减少，init<<<>>>内核也没有成为占用最长时间的热点。这是因为我们在init<<<>>>内核中首次触摸内存。此外，你可能想知道这些 GPU 错误组是什么。正如我们之前讨论的，个别页面错误可能会根据启发式规则和访问模式进行分组，以提高带宽。为了进一步深入了解这一点，让我们使用--print-gpu-trace重新对代码进行分析，以便我们可以看到个别页面错误。正如你从以下截图中所看到的，GPU 跟踪显示了错误的整体跟踪和发生错误的虚拟地址：

我们使用以下命令获取了性能分析输出：

$ nvprof --print-gpu-trace ./unified_initialized.out

第二行显示了相同页面的 11 个页面错误。正如我们之前讨论的，驱动程序的作用是过滤这些重复的错误并只传输每个页面一次。在复杂的访问模式中，通常驱动程序没有足够的信息来确定哪些数据可以迁移到 GPU。为了改善这种情况，我们将进一步实现每页 warp 的概念，这基本上意味着每个 warp 将访问位于相同页面中的元素。这需要开发人员额外的努力。让我们重新实现init<<<>>>内核。你可以在之前编译的unified_memory_64align.cu文件中看到这个实现。以下是内核的快照：

#define STRIDE_64K 65536
__global__ void init(int n, float *x, float *y) {
  int lane_id = threadIdx.x & 31;
  size_t warp_id = (threadIdx.x + blockIdx.x * blockDim.x) >> 5;
  size_t warps_per_grid = (blockDim.x * gridDim.x) >> 5;
  size_t warp_total = ((sizeof(float)*n) + STRIDE_64K-1) / STRIDE_64K;
  for(; warp_id < warp_total; warp_id += warps_per_grid) {
    #pragma unroll
    for(int rep = 0; rep < STRIDE_64K/sizeof(float)/32; rep++) {
      size_t ind = warp_id * STRIDE_64K/sizeof(float) + rep * 32 + lane_id;
      if (ind < n) {
        x[ind] = 1.0f;
        y[ind] = 2.0f;
      }
    }
  }
}

该内核显示索引是基于warp_id。 GPU 中的 warp 大小为 32，负责填充具有 64KB 范围的索引中的x和y变量，也就是说，warp 1 负责前 64KB 的部分，而 warp 2 负责接下来 64KB 的元素。warp 中的每个线程循环（最内层的for循环）以填充相同 64KB 内的索引。让我们来看看这段代码的性能分析结果。正如我们从以下截图中的性能分析输出中所看到的，init<<<>>>内核的时间已经减少，GPU 错误组也大大减少：

我们可以通过使用--print-gpu-trace运行分析器来重新确认这一点：

$ nvprof --print-gpu-trace ./unified_64align.out

以下截图清楚地显示了 GPU 每页的页面错误已经减少：

统一内存的优化使用数据预取

现在，让我们看一个更简单的方法，称为数据预取。CUDA 的一个关键特点是它为开发人员提供了不同的方法，从最简单的方法到需要忍者编程技能的方法。数据预取基本上是对驱动程序的提示，以在使用设备之前预取我们认为将在设备中使用的数据。CUDA 为此目的提供了一个名为cudaMemPrefetchAsync()的预取 API。要查看其实现，请查看我们之前编译的unified_memory_prefetch.cu文件。以下代码片段显示了此代码的快照：

// Allocate Unified Memory -- accessible from CPU or GPU
 cudaMallocManaged(&x, N*sizeof(float));  cudaMallocManaged(&y, N*sizeof(float));
// initialize x and y arrays on the host
 for (int i = 0; i < N; i++) {  x[i] = 1.0f;  y[i] = 2.0f;  } 
//prefetch the memory to GPU
cudaGetDevice(&device);
cudaMemPrefetchAsync(x, N*sizeof(float), device, NULL);
cudaMemPrefetchAsync(y, N*sizeof(float), device, NULL); 
...
 add<<<numBlocks, blockSize>>>(N, x, y);
//prefetch the memory to CPU
 cudaMemPrefetchAsync(y, N*sizeof(float), cudaCpuDeviceId, NULL);
 // Wait for GPU to finish before accessing on host
 cudaDeviceSynchronize();
...
for (int i = 0; i < N; i++)
 maxError = fmax(maxError, fabs(y[i]-3.0f));

代码非常简单，而且解释自己。概念相当简单：在已知将在特定设备上使用哪些内存的情况下，可以预取内存。让我们来看一下在以下截图中显示的分析结果。

正如我们所看到的，add<<<>>>内核提供了我们期望的带宽：

统一内存是一个不断发展的功能，随着每个 CUDA 版本和 GPU 架构的发布而发生变化。预计您通过访问最新的 CUDA 编程指南（docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#um-unified-memory-programming-hd）来保持自己的信息。

到目前为止，我们已经看到了 UM 概念的用处，它不仅提供了编程的便利（不需要使用 CUDA API 显式管理内存），而且在移植原本不可能移植到 GPU 上的应用程序或者原本移植困难的应用程序时更加强大和有用。使用 UM 的一个关键优势是超额订阅。与 CPU 内存相比，GPU 内存非常有限。最新的 GPU（Volta 卡 V100）每个 GPU 提供 32GB 的最大内存。借助 UM，多个 GPU 内存片段以及 CPU 内存可以被视为一个大内存。例如，拥有 16 个 Volta GPU 的 NVIDIA DGX2 机器，其内存大小为 323GB，可以被视为具有最大 512GB 大小的 GPU 内存集合。这对于诸如计算流体动力学（CFD）和分析等应用程序来说是巨大的优势。以前，很难将问题大小适应 GPU 内存，现在却是可能的。手动移动片段容易出错，并且需要调整内存大小。

此外，高速互连的出现，如 NVLink 和 NVSwitch，允许 GPU 之间进行高带宽和低延迟的快速传输。您实际上可以通过统一内存获得高性能！

数据预取，结合指定数据实际所在位置的提示，对于需要同时访问相同数据的多个处理器是有帮助的。在这种情况下使用的 API 名称是cudaMemAdvice()。因此，通过全面了解您的应用程序，您可以通过利用这些提示来优化访问。如果您希望覆盖某些驱动程序启发式方法，这些提示也是有用的。目前 API 正在采用的一些建议如下：

cudaMemAdviseSetReadMostly：顾名思义，这意味着数据大部分是只读的。驱动程序会创建数据的只读副本，从而减少页面错误。重要的是要注意，数据仍然可以被写入。在这种情况下，页面副本将变得无效，除了写入内存的设备：

// Sets the data readonly for the GPU
cudaMemAdvise(data, N, ..SetReadMostly, processorId); 
mykernel<<<..., s>>>(data, N);

cudaMemAdviseSetPreferredLocation：此建议将数据的首选位置设置为设备所属的内存。设置首选位置不会立即导致数据迁移到该位置。就像在以下代码中，mykernel<<<>>>将会出现页面错误并在 CPU 上生成数据的直接映射。驱动程序试图抵制将数据迁离设置的首选位置，使用cudaMemAdvise：

cudaMemAdvise(input, N, ..PreferredLocation, processorId); 
mykernel<<<..., s>>>(input, N);

cudaMemAdviseSetAccessedBy：这个建议意味着数据将被设备访问。设备将在 CPU 内存中创建输入的直接映射，不会产生页面错误：

cudaMemAdvise(input, N, ..SetAccessedBy, processorId); 
mykernel<<<..., s>>>(input, N);

在接下来的部分，我们将以整体的视角来看 GPU 中不同的内存是如何随着新的架构而发展的。

GPU 内存的演变

GPU 架构随着时间的推移发生了变化，内存架构也发生了相当大的变化。如果我们看一下过去四代的情况，会发现一些共同的模式，其中一些如下：

内存容量总体上已经提高了几个级别。
内存带宽和容量随着新一代架构的出现而增加。

以下表格显示了过去四代的属性：

内存类型	属性	Volta V100	Pascal P100	Maxwell M60	Kepler K80
寄存器	每个 SM 的大小	256 KB	256 KB	256 KB	256 KB
L1	大小	32...128 KiB	24 KiB	24 KiB	16...48 KiB
行大小	32	32 B	32 B	128 B
L2	大小	6144 KiB	4,096 KiB	2,048 KiB	1,536 Kib
行大小	64 B	32B	32B	32B
共享内存	每个 SMX 的大小	高达 96 KiB	64 KiB	64 KiB	48 KiB
每个 GPU 的大小	高达 7,689 KiB	3,584 KiB	1,536 KiB	624 KiB
理论带宽	13,800 GiB/s	9,519 GiB/s	2,410 GiB/s	2,912 GiB/s
全局内存	内存总线	HBM2	HBM2	GDDR5	GDDR5
大小	32,152 MiB	16,276 MiB	8,155 MiB	12,237 MiB
理论带宽	900 GiB/s	732 GiB/s	160 GiB/s	240 GiB/s

总的来说，前面的观察结果已经帮助 CUDA 应用在新的架构下运行得更快。但与此同时，CUDA 编程模型和内存架构也进行了一些根本性的改变，以便为 CUDA 程序员简化工作。我们观察到的一个这样的改变是纹理内存，之前在 CUDA 5.0 之前，开发人员必须手动绑定和解绑纹理，并且必须在全局声明。但在 CUDA 5.0 中，这是不必要的。它还取消了应用程序中开发人员可以拥有的纹理引用数量的限制。

我们还研究了 Volta 架构以及为简化开发人员编程而进行的一些根本性改变。Volta 的总容量是每个 SM 128 KB，比其上一代显卡 Pascal P100 多了七倍，这为开发人员提供了更大的缓存。此外，由于 Volta 架构中 L1 缓存的延迟要小得多，这使得它对频繁重用的数据具有高带宽和低延迟的访问。这样做的关键原因是让 L1 缓存操作获得共享内存性能的好处。共享内存的关键问题是需要开发人员显式控制。在使用 Volta 等新架构时，这种需求就不那么必要了。但这并不意味着共享内存变得多余。一些极客程序员仍然希望充分利用共享内存的性能，但许多其他应用程序不再需要这种专业知识。Pascal 和 Volta L1 缓存和共享内存之间的区别如下图所示：

前面的图表显示了与 Pascal 相比共享内存和 L1 缓存的统一。重要的是要理解，CUDA 编程模型从诞生以来几乎保持不变。尽管每个架构的内存容量、带宽或延迟都在变化，但相同的 CUDA 代码将在所有架构上运行。不过，随着这些架构变化，性能的影响肯定会发生变化。例如，在 Volta 之前利用共享内存的应用程序，与使用全局内存相比可能会看到性能提升，但在 Volta 中可能不会看到这样的加速，因为 L1 和共享内存的统一。

为什么 GPU 有缓存？

在这个演变过程中，还很重要的一点是要理解 CPU 和 GPU 缓存是非常不同的，而且有不同的用途。作为 CUDA 架构的一部分，我们通常在每个 SM 上启动数百到数千个线程。数万个线程共享 L2 缓存。因此，L1 和 L2 对每个线程来说都很小。例如，在每个 SM 上有 2,048 个线程，共有 80 个 SM，每个线程只能获得 64 字节的 L1 缓存和 38 字节的 L2 缓存。GPU 缓存中存储着许多线程访问的公共数据。这有时被称为空间局部性。一个典型的例子是当线程的访问是不对齐和不规则的时候。GPU 缓存可以帮助减少寄存器溢出和局部内存的影响，因为 CPU 缓存主要用于时间局部性。

总结

我们在本章开始时介绍了不同类型的 GPU 内存。我们详细讨论了全局、纹理和共享内存，以及寄存器。我们还看了 GPU 内存演变提供了哪些新功能，例如统一内存，这有助于提高程序员的生产力。我们看到了这些功能在最新的 GPU 架构（如 Pascal 和 Volta）中是如何实现的。

在下一章中，我们将深入讨论 CUDA 线程编程的细节，以及如何最优地启动不同的线程配置，以发挥 GPU 硬件的最佳性能。我们还将介绍新的 CUDA Toolkit 功能，例如用于灵活线程编程的协作组和 GPU 上的多精度编程。

CUDA 线程编程

CUDA 具有分层线程架构，因此我们可以控制 CUDA 线程的分组。了解它们在 GPU 上并行工作的方式有助于您编写并行编程代码并实现更好的性能。在本章中，我们将介绍 CUDA 线程操作及其与 GPU 资源的关系。作为实际经验，我们将研究并行减少算法，并看看如何通过使用优化策略来优化 CUDA 代码。

在本章中，您将学习 CUDA 线程在 GPU 中的操作：并行和并发线程执行，warp 执行，内存带宽问题，控制开销，SIMD 操作等等。

本章将涵盖以下主题：

层次化的 CUDA 线程操作
了解 CUDA 占用率
跨多个 CUDA 线程共享数据
识别应用程序的性能限制
最小化 CUDA warp 分歧效应
增加内存利用率和网格跨距循环
用于灵活线程处理的协作组
warp 同步编程
低/混合精度操作

技术要求

本章建议使用比 Pascal 架构更晚的 NVIDIA GPU 卡。换句话说，您的 GPU 的计算能力应等于或大于 60。如果您不确定您的 GPU 架构，请访问 NVIDIA 的 GPU 网站developer.nvidia.com/cuda-gpus，并确认您的 GPU 的计算能力。

在我们写这本书的时候，示例代码是使用 10.1 版本开发和测试的。一般来说，如果适用的话，建议使用最新的 CUDA 版本。

在本章中，我们将通过对代码进行性能分析来进行 CUDA 编程。如果您的 GPU 架构是图灵架构，建议安装 Nsight Compute 来对代码进行性能分析。它是免费的，您可以从developer.nvidia.com/nsight-compute下载。在我们写这本书的时候，这是性能分析工具的过渡时刻。您可以在第五章的使用 Nsight Compute 对内核进行性能分析部分了解其基本用法，CUDA 应用性能分析和调试。

CUDA 线程、块和 GPU

CUDA 编程中的基本工作单元是 CUDA 线程。基本的 CUDA 线程执行模型是单指令多线程（SIMT）。换句话说，内核函数的主体是单个 CUDA 线程的工作描述。但是，CUDA 架构执行具有相同操作的多个 CUDA 线程。

在概念上，多个 CUDA 线程以组的形式并行工作。CUDA 线程块是多个 CUDA 线程的集合。多个线程块同时运行。我们称线程块的组为网格。以下图表显示了它们之间的关系：

这些分层的 CUDA 线程操作与分层的 CUDA 架构相匹配。当我们启动 CUDA 内核时，每个流多处理器上会执行一个或多个 CUDA 线程块。此外，根据资源的可用性，一个流多处理器可以运行多个线程块。线程块中的线程数量和网格中的块数量也会有所不同。

流多处理器以任意和并发的方式执行线程块，执行尽可能多的 GPU 资源。因此，可并行执行的线程块数量取决于块需要的 GPU 资源量以及 GPU 资源的可用量。我们将在接下来的部分中介绍这一点。流多处理器的数量取决于 GPU 规格。例如，Tesla V100 为 80，RTX 2080（Ti）为 48。

CUDA 流多处理器以 32 个线程的组形式控制 CUDA 线程。一个组被称为warp。这样，一个或多个 warp 配置一个 CUDA 线程块。以下图显示了它们的关系：

小绿色框是 CUDA 线程，它们被 warp 分组。warp 是 GPU 架构的基本控制单元。因此，它的大小对 CUDA 编程具有隐式或显式的影响。例如，最佳线程块大小是在可以充分利用块的 warp 调度和操作的多个 warp 大小中确定的。我们称之为占用率，这将在下一节中详细介绍。此外，warp 中的 CUDA 线程并行工作，并具有同步操作。我们将在本章的Warp 级别基元编程部分讨论这一点。

利用 CUDA 块和 warp

现在，我们将研究 CUDA 线程调度及其使用 CUDA 的printf进行隐式同步。并行 CUDA 线程的执行和块操作是并发的。另一方面，从设备打印输出是一个顺序任务。因此，我们可以轻松地看到它们的执行顺序，因为对于并发任务来说输出是任意的，而对于并行任务来说是一致的。

我们将开始编写打印全局线程索引、线程块索引、warp 索引和 lane 索引的内核代码。为此，代码可以编写如下：

__global__ void index_print_kernel() {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int warp_idx = threadIdx.x / warpSize;
    int lane_idx = threadIdx.x & (warpSize - 1);

    if ((lane_idx & (warpSize/2 - 1)) == 0)
        //thread, block, warp, lane
        printf(" %5d\t%5d\t %2d\t%2d\n", idx, blockIdx.x, 
               warp_idx, lane_idx);
}

这段代码将帮助我们理解 warp 和 CUDA 线程调度的并发性。让我们让我们的代码从 shell 获取参数，以便轻松测试各种网格和线程块配置。

然后，我们将编写调用内核函数的主机代码：

int main() {
    int gridDim = 4, blockDim = 128;
    puts("thread, block, warp, lane");
    index_print_kernel<<< gridDim, blockDim >>>();
    cudaDeviceSynchronize();
}

最后，让我们编译代码，执行它，并查看结果：

nvcc -m64 -o cuda_thread_block cuda_thread_block.cu

以下结果是输出结果的一个示例。实际输出可能会有所不同：

$ ./cuda_thread_block.cu 4 128
thread, block, warp, lane
 64     0     2     0
 80     0     2    16
 96     0     3     0
 112     0     3    16
 0     0     0     0
 16     0     0    16
 ...
 352     2     3     0
 368     2     3    16
 288     2     1     0
 304     2     1    16

从结果中，您将看到 CUDA 线程以 warp 大小启动，并且顺序是不确定的。另一方面，lane 输出是有序的。从给定的结果中，我们可以确认以下事实：

**无序块执行：**第二列显示线程块的索引。结果表明，它不保证按照块索引的顺序执行。
**无序 warp 索引与线程块：**第三列显示块中 warp 的索引。warp 的顺序在块之间变化。因此，我们可以推断 warp 执行顺序没有保证。
**在 warp 中执行的分组线程：**第四列显示 warp 中的 lane。为了减少输出数量，应用程序限制只打印两个索引。从每个 warp 内的有序输出中，我们可以类比printf函数的输出顺序是固定的，因此没有倒置。

总之，CUDA 线程被分组为 32 个线程，它们的输出和 warp 的执行没有顺序。因此，程序员必须牢记这一点，以便进行 CUDA 内核开发。

理解 CUDA 占用率

CUDA 占用率是活动 CUDA warps 与每个流多处理器可以同时执行的最大 warps 的比率。一般来说，更高的占用率会导致更有效的 GPU 利用率，因为有更多的 warp 可用来隐藏停滞 warp 的延迟。然而，它也可能由于 CUDA 线程之间资源争用的增加而降低性能。因此，开发人员理解这种权衡是至关重要的。

找到最佳的 CUDA 占用率的目的是使 GPU 应用程序能够有效地使用 GPU 资源发出 warp 指令。GPU 在流多处理器上使用多个 warp 调度器调度多个 warp。当多个 warp 有效地调度时，GPU 可以隐藏 GPU 指令或内存延迟之间的延迟。然后，CUDA 核心可以执行连续从多个 warp 发出的指令，而未调度的 warp 必须等待，直到它们可以发出下一条指令。

开发人员可以使用两种方法确定 CUDA 占用率：

由 CUDA 占用率计算器确定的理论占用率：这个计算器是 CUDA 工具包提供的一个 Excel 表。我们可以从内核资源使用和 GPU 流多处理器理论上确定每个内核的占用率。
由 GPU 确定的实现占用率：实现占用率反映了在流多处理器上并发执行的 warp 的真实数量和最大可用 warp。这种占用率可以通过 NVIDIA 分析器进行度量分析来测量。

理论占用率可以被视为最大的上限占用率，因为占用率数字不考虑指令依赖性或内存带宽限制。

现在，让我们看看这个占用率和 CUDA C/C++之间的关系。

设置 NVCC 报告 GPU 资源使用

首先，我们将使用简单矩阵乘法（SGEMM）内核代码，如下所示：

__global__ void sgemm_gpu_kernel(const float *A, const float *B, 
        float *C, int N, int M, int K, alpha, float beta) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    float sum = 0.f;
    for (int i = 0; i < K; ++i) {
        sum += A[row * K + i] * B[i * K + col];
    }
    C[row * M + col] = alpha * sum + beta * C[row * M + col];
}

然后，我们将使用以下内核代码调用内核函数：

void sgemm_gpu(const float *A, const float *B, float *C,
            int N, int M, int K, float alpha, float beta) {
    dim3 dimBlock(BLOCK_DIM, BLOCK_DIM);
    dim3 dimGrid(M / dimBlock.x, N / dimBlock.y);
    sgemm_gpu_kernel<<< dimGrid, dimBlock >>>(A, B, C, N, M, K, alpha, beta);
}

您可能希望提供适当的 GPU 内存及其大小信息。我们将使用 2048 作为N，M和K。内存大小是该数字的平方。我们将把BLOCK_DIM设置为16。

现在，让我们看看如何使nvcc编译器报告内核函数的 GPU 资源使用情况。

Linux 设置

在 Linux 环境中，我们应该提供两个编译器选项，如下所示：

--resource-usage（--res-usage）：为 GPU 资源使用设置详细选项
-gencode：指定要编译和生成操作码的目标架构如下：
Turing：compute_75,sm_75
Volta：compute_70,sm_70
Pascal：compute_60,sm_60，compute_61,sm_61

如果您不确定您正在使用哪种架构，您可以从 CUDA GPU 网站上找到（developer.nvidia.com/cuda-gpus）。例如，nvcc编译命令可以有以下编译选项：

$ nvcc -m 64 --resource-usage \
 -gencode arch=compute_70,code=sm_70 \
 -I/usr/local/cuda/samples/common/inc \
 -o sgemm ./sgemm.cu

我们还可以编译代码以针对多个 GPU 架构，如下所示：

$ nvcc -m64 --resource-usage \
      -gencode arch=compute_70,code=sm_70 \
      -gencode arch=compute_75,code=sm_75 \
      -I/usr/local/cuda/samples/common/inc \
      -o sgemm ./sgemm.cu

如果您想使您的代码与新的 GPU 架构（Turing）兼容，您需要提供以下附加选项：

$ nvcc -m64 --resource-usage \
      -gencode arch=compute_70,code=sm_70 \
      -gencode arch=compute_75,code=sm_75 \
      -gencode arch=compute_75,code=compute_75 \
      -I/usr/local/cuda/samples/common/inc \
      -o sgemm ./sgemm.cu

如果您想了解更多关于这些选项的信息，您可以在这个文档中找到相关信息：docs.nvidia.com/cuda/turing-compatibility-guide/index.html#building-turing-compatible-apps-using-cuda-10-0。

现在，让我们编译源代码。我们可以从 NVCC 的输出中找到一个资源使用报告。以下结果是使用前面的命令生成的：

NVCC 为每个计算能力报告 CUDA 内核资源使用信息。在前面的输出截图中，我们可以看到每个线程的寄存器数量和常量内存使用情况。

Windows 设置

当我们开发 Windows 应用程序时，我们可以在 Visual Studio 项目属性对话框中设置这些设置。以下是该对话框的截图：

要打开此对话框，我们应该打开 debug_vs 属性页，然后在左侧面板上转到 CUDA C/C++ | 设备选项卡。然后，我们应该设置以下选项如下：

Verbose PTXAS Output: No | Yes
代码生成：更新选项以指定您的目标架构如下：
图灵：compute_75,sm_75
伏尔塔：compute_70,sm_70
帕斯卡：compute_60,sm_60;compute_61,sm_61

我们可以使用分号(;)指定多个目标架构。

现在，让我们构建源代码，我们将在 Visual Studio 的输出面板上看到 NVCC 的报告。然后，你会看到类似以下的输出：

这与 Linux 中 NVCC 的输出相同。

现在，让我们使用资源使用报告来分析内核的占用情况。

使用占用率计算器分析最佳占用率

实际上，我们可以使用 CUDA 占用率计算器，它是 CUDA 工具包提供的。使用这个，我们可以通过提供一些内核信息来获得理论上的占用率。计算器是一个 Excel 文件，你可以在以下位置找到它，根据你使用的操作系统：

Windows: C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\<cuda-version>\tools
Linux: /usr/local/cuda/tools
macOS: /Developer/NVIDIA/<cuda-version>/tools

以下是计算器的屏幕截图：

CUDA 占用率计算器

这个计算器有两部分：内核信息输入和占用信息输出。作为输入，它需要两种信息，如下所示：

GPU 的计算能力（绿色）
线程块资源信息（黄色）:
每个 CUDA 线程块的线程
每个 CUDA 线程的寄存器
每个块的共享内存

计算器在这里显示了 GPU 的占用信息：

GPU 占用数据（蓝色）
GPU 的 GPU 计算能力的物理限制（灰色）
每个块分配的资源（黄色）
每个流多处理器的最大线程块（黄色、橙色和红色）
根据三个关键的占用资源（线程、寄存器和每个块的共享内存），绘制占用限制图
图表上的红色三角形，显示当前的占用数据

现在，让我们把获得的信息放入计算器中。我们可以编辑 Excel 表格中的绿色和橙色区域：

输入你获得的内核资源信息，看看表格如何变化。

根据计算能力和输入数据，占用情况会发生变化，如下面的屏幕截图所示：

根据计算能力和输入数据的变化

蓝色区域显示了内核函数实现的占用率。在这个屏幕截图中，它显示了 100%的占用率。表格的右侧显示了 GPU 资源的占用率利用图：CUDA 线程、共享内存和寄存器。

一般来说，由于许多原因，内核代码不能达到 100%的理论占用率。然而，设置峰值占用率是有效利用 GPU 资源的开始。

占用率调整 - 限制寄存器使用

当内核算法复杂或处理数据类型为双精度时，CUDA 寄存器使用可能会增加。在这种情况下，由于活动 warp 大小有限，占用率会下降。在这种情况下，我们可以通过限制寄存器使用来增加理论上的占用率，并查看性能是否提高。

调整 GPU 资源使用的一种方法是在内核函数中使用__launch_bound__限定符。这告诉 NVCC 保证每个流多处理器的最大块大小的最小线程块。然后，NVCC 找到实现给定条件的最佳寄存器大小。如果你在编译时知道使你的算法有效运行的大小，你可以使用这个。标识符可以如下使用：

int maxThreadPerBlock = 256;
int minBlocksPerMultiprocessor = 2;
__global__ void
__launch_bound__ (maxThreadPerBlock, minBlocksPerMultiprocessor) foo_kernel() {
    ...
}

然后，编译器检查上限资源并减少每个块的限制资源使用。如果其资源使用没有超过上限，编译器会调整寄存器使用，如果 CUDA 可以调度额外的多处理器线程块，如果没有给出第二个参数。或者，编译器会增加寄存器使用以隐藏单线程指令延迟。

此外，我们可以简单地在应用程序级别限制占用寄存器的数量。--maxrregcount标志到NVCC将指定数量，编译器将重新排列寄存器使用。以下编译命令显示了如何在 Linux 终端中使用该标志：

$ nvcc -m64 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 --resource-usage --maxrregcount 24 -o sgemm ./sgemm.cu

但是，请记住，以这种方式限制寄存器使用可能会引入由寄存器限制引起的线程性能下降。即使编译器无法将其设置在限制之下，也可以将寄存器分割为本地内存，并且本地变量放置在全局内存中。

从分析器获取实现的占用

现在，我们可以使用 Visual Profiler 从分析的度量数据中获取实现的占用。单击目标内核时间轴条。然后，我们可以在属性面板中看到理论和实现的占用。我们还可以从内核延迟菜单中获取更多详细信息。以下屏幕截图显示了我们使用的示例代码的实现性能：

显示实现和理论占用的性能

通过这种占用调整，我们可以设计 CUDA 块大小，充分利用流多处理器中的 warp 调度。然而，这并没有解决我们在上一节中发现的 54.75%的内存限制问题。这意味着多处理器可能会停顿，无法掩盖由于受阻内存请求而产生的内存访问延迟。我们将在本章讨论如何优化这一点，并且在第七章《CUDA 中的并行编程模式》中，我们将讨论矩阵乘法优化。

理解并行归约

归约是一种简单但有用的算法，可以获得许多参数的公共参数。这个任务可以按顺序或并行完成。当涉及到并行处理到并行架构时，并行归约是获得直方图、均值或任何其他统计值的最快方式。

以下图表显示了顺序归约和并行归约之间的差异：

通过并行进行归约任务，可以将并行归约算法的总步骤减少到对数级别。现在，让我们开始在 GPU 上实现这个并行归约算法。首先，我们将使用全局内存实现一个简单的设计。然后，我们将使用共享内存实现另一个归约版本。通过比较这两种实现，我们将讨论是什么带来了性能差异。

使用全局内存的天真并行归约

归约的第一种基本方法是使用并行的 CUDA 线程，并使用全局内存共享归约输出。对于每次迭代，CUDA 内核通过将其大小减少两倍来从全局内存获取累积值。归约的工作如下图所示，显示了使用全局内存数据共享的天真并行归约：

这种方法在 CUDA 中很慢，因为它浪费了全局内存的带宽，并且没有利用任何更快的片上内存。为了获得更好的性能，建议使用共享内存来节省全局内存带宽并减少内存获取延迟。我们将讨论这种方法如何浪费带宽。

现在，让我们实现这个归约。首先，我们将编写归约内核函数，如下所示：

__global__ void naive_reduction_kernel
     (float *data_out, float *data_in, int stride, int size) {
     int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
     if (idx_x + stride < size)
         data_out[idx_x] += data_in[idx_x + stride];
}

我们将在迭代过程中减半减小步长大小，直到stride大小为 1 时调用内核函数。

void naive_reduction(float *d_out, float *d_in, int n_threads, int size) {
    int n_blocks = (size + n_threads - 1) / n_threads;
    for (int stride = 1; stride < size; stride *= 2)
        naive_reduction_kernel<<<n_blocks, n_threads>>>(d_out, d_in, stride, size);
}

在这个实现中，内核代码使用跨距寻址获取设备内存并输出一个减少结果。主机代码触发每个步骤的减少内核，并且参数大小减半。我们不能有内部内核循环，因为 CUDA 不能保证线程块和流多处理器之间的同步操作。

使用共享内存减少内核

在这种减少中，每个 CUDA 线程块减少输入值，并且 CUDA 线程使用共享内存共享数据。为了进行适当的数据更新，它们使用块级内在同步函数__syncthreads()。然后，下一个迭代操作上一个减少结果。其设计如下图所示，显示了使用共享内存的并行减少：

黄点框表示 CUDA 线程块的操作范围。在这个设计中，每个 CUDA 线程块输出一个减少结果。

块级减少允许每个 CUDA 线程块进行减少，并输出单个减少输出。由于它不需要我们将中间结果保存在全局内存中，CUDA 内核可以将过渡值存储在共享内存中。这种设计有助于节省全局内存带宽并减少内存延迟。

与全局减少一样，我们将实现这个操作。首先，我们将编写内核函数，如下所示：

__global__ void reduction_kernel(float* d_out, float* d_in, 
                                 unsigned int size) {
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;

    extern __shared__ float s_data[];
    s_data[threadIdx.x] = (idx_x < size) ? d_in[idx_x] : 0.f;

    __syncthreads();

    // do reduction
    for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
        // thread synchronous reduction
        if ( (idx_x % (stride * 2)) == 0 )
            s_data[threadIdx.x] += s_data[threadIdx.x + stride];

        __syncthreads();
    }

    if (threadIdx.x == 0)
        d_out[blockIdx.x] = s_data[0];
}

然后，我们将调用内核函数，如下所示：

void reduction(float *d_out, float *d_in, int n_threads, int size)
{
    cudaMemcpy(d_out, d_in, size * sizeof(float), cudaMemcpyDeviceToDevice);
    while(size > 1) {
        int n_blocks = (size + n_threads - 1) / n_threads;
        reduction_kernel
            <<< n_blocks, n_threads, n_threads * sizeof(float), 0 >>>
            (d_out, d_out, size);
        size = n_blocks;
    }
}

在这段代码中，我们提供了n_threads * sizeof (float)字节，因为每个 CUDA 线程将共享每个字节的单个变量。

编写性能测量代码

为了测量每个版本的性能，我们将使用 CUDA 示例timer辅助函数：

// Initialize timer
StopWatchInterface *timer;
sdkCreateTimer(&timer);
sdkStartTimer(&timer);

... Execution code ...

// Getting elapsed time
cudaDeviceSynchronize(); // Blocks the host until GPU finishes the work
sdkStopTimer(&timer);

// Getting execution time in micro-secondes
float execution_time_ms = sdkGetTimerValue(&timer)

// Termination of timer
sdkDeleteTimer(&timer);

这个函数集有助于在微秒级别测量执行时间。此外，建议在性能测量之前调用内核函数，以消除设备初始化开销。有关更详细的实现，请访问global_reduction.cu和reduction.cu文件中的实现代码。这些代码集在本章中用于评估优化效果以及分析器。

两种减少-全局和共享内存的性能比较

现在，我们可以比较两个并行减少操作的执行时间。性能可能会因 GPU 和实现环境而异。分别运行以下命令进行全局减少和使用共享内存进行减少：

# Reduction with global memory
$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction_global ./reduction_global.cpp reduction_global_kernel.cu

# Reduction using shared memory
$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction_shared ./reduction_shared.cpp reduction_shared_kernel.cu

使用我的 Tesla V100 PCIe 卡，两种减少的估计性能如下。元素数量为2²⁴个：

操作	估计时间（毫秒）	加速
原始方法（使用全局内存进行减少）	4.609	1.0x
使用共享内存的减少	0.624	7.4x

从这个结果中，我们可以看到在减少中使用共享内存共享数据如何快速返回输出。第一个实现版本在global_reduction.cu中，第二个版本在shared_reduction.cu中，所以您可以自行比较实现。

通过将减少与共享内存结合，我们可以显著提高性能。然而，我们无法确定这是否是我们可以获得的最大性能，并且不知道我们的应用程序有什么瓶颈。为了分析这一点，我们将在下一节中涵盖性能限制器。

识别应用程序的性能限制器

之前，我们看到了如何通过保存全局内存来使 CUDA 内核的性能受益。一般来说，使用片上缓存比使用片外内存更好。但是，我们无法确定这种简单类比是否还有很多优化空间。

性能限制因素显示了限制应用程序性能的因素，它最显著地限制了应用程序的性能。根据其分析信息，它分析了计算和内存带宽之间的性能限制因素。根据这些资源的利用率，应用程序可以被分类为四种类型：计算受限，带宽受限，延迟受限和计算和延迟受限。以下图表显示了这些类别与计算和内存利用率的关系：

在确定了限制因素之后，我们可以使用下一个优化策略。如果任一资源的利用率很高，我们可以专注于优化该资源。如果两者都未充分利用，我们可以从系统的 I/O 方面应用延迟优化。如果两者都很高，我们可以调查是否存在内存操作停顿问题和与计算相关的问题。

现在让我们看看如何获得利用率信息。

找到性能限制因素并进行优化

现在，让我们将此分析应用于两个减少实现。我们将对它们进行比较，并讨论共享内存如何有助于性能限制因素分析以改善性能。首先，让我们使用以下命令对基于全局内存的减少应用程序进行度量分析：

$ nvprof -o reduction_global.nvvp ./reduction_global 
$ nvprof --analysis-metrics -o reduction_global_metric.nvvp ./reduction_global

然后，我们将从 NVIDIA 分析器获得以下图表，显示了基于全局内存的第一个减少性能的限制因素：

在这张图表上，我们需要查看性能执行比来查看是否通过检查内核延迟分析来平衡。因为，如前图表所示，计算和内存之间的利用率差距很大，这可能意味着由于内存瓶颈，计算中会有很多延迟。以下图表显示了基于采样的分析结果，我们可以确定 CUDA 核心由于内存依赖而饥饿：

如您所见，由于内存等待，内核执行被延迟。现在，让我们基于共享内存对减少进行分析。我们可以使用以下命令来做到这一点：

$ nvprof -o reduction_shared.nvvp ./reduction_shared 
$ nvprof --analysis-metrics -o reduction_shared_metric.nvvp ./reduction_shared

然后，我们将获得以下图表，显示了基于共享内存的第二个减少性能的限制因素：

我们可以确定它是计算受限的，内存不会使 CUDA 核心饥饿。

现在让我们回顾我们的核心操作，以优化计算操作。以下代码显示了内核函数中的并行减少部分：

for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
     if ( (idx_x % (stride * 2)) == 0 )
         s_data[threadIdx.x] += s_data[threadIdx.x + stride];
     __syncthreads();
 }

作为算术操作，模运算是一种重型操作。由于stride变量是2的指数倍数，因此可以用位操作替换，如下所示：

for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
     if ( (idx_x & (stride * 2 - 1)) == 0 )
         s_data[threadIdx.x] += s_data[threadIdx.x + stride];
     __syncthreads();
 }

运行以下命令以查看优化后的输出：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction_shared ./reduction_shared.cpp reduction_shared_kernel.cu

然后，新的估计时间为0.399 毫秒，我们可以实现更优化的性能，如下表所示：

操作	估计时间（毫秒）	加速比
原始方法（使用全局内存进行减少）	4.609	1.0 倍
使用共享内存进行减少	0.624	7.4 倍
将条件操作从`%`更改为`&`	0.399	11.55 倍

以下图表显示了更新后的性能限制因素：

我们可以确定其操作是计算和延迟受限。因此，我们可以确定我们可以通过优化计算效率来增加内存利用率。

最小化 CUDA warp 分歧效应

在单指令，多线程（SIMT）执行模型中，线程被分组成 32 个线程的集合，每个集合称为warp。如果一个 warp 遇到条件语句或分支，其线程可以分歧并串行执行每个条件。这称为分支分歧，它会显著影响性能。

CUDA warp 分歧是指在 warp 中 CUDA 线程的分歧操作。如果条件分支具有if-else结构，并且 warp 具有此 warp 分歧，所有 CUDA 线程对于分支代码块都有活动和非活动操作部分。

下图显示了 CUDA warp 中的 warp 分歧效应。不处于空闲状态的 CUDA 线程会降低 GPU 线程的有效使用：

随着分支部分的增加，GPU 调度吞吐量变得低效。因此，我们需要避免或最小化这种 warp 分歧效应。您可以选择几种选项：

通过处理不同的 warp 来避免分歧效应
通过合并分支部分来减少 warp 中的分支
缩短分支部分；只有关键部分进行分支
重新排列数据（即转置，合并等）
使用协作组中的tiled_partition来对组进行分区

确定分歧作为性能瓶颈

从先前的减少优化中，您可能会发现由于计算分析中的分歧分支而导致内核效率低下的警告，如下所示：

73.4％的分歧意味着我们有一个低效的操作路径。我们可以确定减少寻址是问题所在，如下所示：

__global__ void reduction_kernel(float* d_out, float* d_in, unsigned int size) {
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;

    extern __shared__ float s_data[];
    s_data[threadIdx.x] = (idx_x < size) ? d_in[idx_x] : 0.f;

    __syncthreads();

    // do reduction
    for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
        // thread synchronous reduction
        if ( (idx_x % (stride * 2 - 1)) == 0 )
            s_data[threadIdx.x] += s_data[threadIdx.x + stride];

        __syncthreads();
    }

    if (threadIdx.x == 0)
        d_out[blockIdx.x] = s_data[0];
}

在减少寻址方面，我们可以选择以下 CUDA 线程索引策略之一：

交错寻址
顺序寻址

让我们回顾一下它们，并通过实施这些策略来比较它们的性能。由于我们只会修改减少内核，因此我们可以重用主机代码进行下两个实现。

交错寻址

在这种策略中，连续的 CUDA 线程使用交错寻址策略获取输入数据。与之前的版本相比，CUDA 线程通过增加步幅值来访问输入数据。以下图表显示了 CUDA 线程如何与减少项交错：

可以实现以下交错寻址：

__global__ void
 interleaved_reduction_kernel(float* g_out, float* g_in, unsigned int size) {
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;

    extern __shared__ float s_data[];
    s_data[threadIdx.x] = (idx_x < size) ? g_in[idx_x] : 0.f;
    __syncthreads();

    // do reduction
    // interleaved addressing
    for (unsigned int stride = 1; stride < blockDim.x; stride *= 2) {
        int index = 2 * stride * threadIdx.x;
        if (index < blockDim.x)
            s_data[index] += s_data[index + stride];
        __syncthreads();
    }
    if (threadIdx.x == 0)
        g_out[blockIdx.x] = s_data[0];
}

运行以下命令来编译上述代码：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction ./reduction.cpp ./reduction_kernel_interleaving.cu

在 Tesla V100 上，测得的内核执行时间为 0.446 毫秒。这比之前的版本慢，因为在这种方法中每个线程块都没有完全利用。通过对其指标进行分析，我们可以得到更多细节。

现在我们将尝试另一种寻址方法，该方法旨在使每个线程块计算更多数据。

顺序寻址

与之前的版本相比，这具有高度合并的索引和寻址。这种设计更有效，因为当步幅大小大于 warp 大小时就没有分歧。以下图表显示了合并的线程操作：

现在，让我们编写一个内核函数，以在减少项上使用顺序寻址。

__global__ void
 sequantial_reduction_kernel(float *g_out, float *g_in, 
                             unsigned int size)
{
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;

    extern __shared__ float s_data[];

    s_data[threadIdx.x] = (idx_x < size) ? g_in[idx_x] : 0.f;

    __syncthreads();

    // do reduction
    // sequential addressing
    for (unsigned int stride = blockDim.x / 2; stride > 0; 
         stride >>= 1)
    {
        if (threadIdx.x < stride)
            s_data[threadIdx.x] += s_data[threadIdx.x + stride];

        __syncthreads();
    }

    if (threadIdx.x == 0)
        g_out[blockIdx.x] = s_data[0];
}

运行以下命令来编译上述代码：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction ./reduction.cpp ./reduction_kernel_sequential.cu

在 Tesla V100 GPU 上，测得的执行时间为 0.378 毫秒，略快于之前的策略（0.399 毫秒）。

由于避免 warp 分歧，我们可以在原始计算上获得 12.2 倍的性能提升。以下图表显示了更新后的性能限制器分析：

与之前的性能限制器相比，我们可以看到减少了控制流操作并增加了内存利用率。

性能建模和平衡限制器

根据性能限制器分析，我们当前的减少性能受到计算延迟的限制，这是由于内存带宽，尽管限制器分析显示每个资源的充分利用。让我们讨论为什么这是一个问题，以及如何通过遵循 Roofline 性能模型来解决这个问题。

Roofline 模型

Roofline 模型是一种直观的视觉性能分析模型，用于为并行处理单元上的给定计算内核提供估计性能。根据这个模型，并行编程中的开发人员可以确定算法应该受到什么限制，并确定哪些应该进行优化。

以下图表显示了 Roofline 模型的一个示例：

倾斜部分表示内存受限，平坦部分表示算术受限。每个并行算法和实现都有自己的 Roofline 模型，因为它们具有不同的计算能力和内存带宽。有了这个模型，算法可以根据它们的操作密度（flops/bytes）进行放置。如果一个实现不符合这个模型的预期性能，我们可以确定这个版本受到延迟的限制。

考虑到我们并行减少的复杂性，它必须是内存受限的。换句话说，它的操作密度低，因此我们的策略应尽可能最大化内存带宽。

因此，我们需要确认我们的减少内核函数如何使用性能分析器中的内存带宽。以下图表显示了全局内存的带宽使用情况：

如图所示，我们没有充分利用内存带宽。Tesla V100 GPU 的总带宽为 343.376 GB/s，利用了大约三分之一的带宽，因为这款 GPU 具有 900 GB/s 带宽的 HBM2 内存。因此，下一步是通过让每个 CUDA 线程处理更多数据来增加带宽使用率。这将解决延迟限制的情况，并使我们的应用程序受限于内存带宽。

现在，让我们讨论如何增加内存带宽。

通过网格跨步循环最大化内存带宽

我们可以通过一个简单的想法实现这一点。减少问题允许我们使用 CUDA 线程累积输入数据并开始减少操作。以前，我们的减少实现是从输入数据大小开始的。但现在，我们将迭代到一组 CUDA 线程的输入数据，并且该大小将是我们内核函数的网格大小。这种迭代风格称为网格跨步循环。这种技术有许多好处，可以控制多个 CUDA 核心，并在本文中介绍：devblogs.nvidia.com/cuda-pro-tip-write-flexible-kernels-grid-stride-loops。

以下代码显示了更新后的减少内核函数：

__global__ void reduction_kernel(float *g_out, float *g_in, 
                                 unsigned int size) {
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
    extern __shared__ float s_data[];

    // cumulates input with grid-stride loop 
 // and save to the shared memory
 float input = 0.f;
 for (int i = idx_x; i < size; i += blockDim.x * gridDim.x)
 input += g_in[i];
 s_data[threadIdx.x] = input;
 __syncthreads();

    // do reduction
    for (unsigned int stride = blockDim.x / 2; stride > 0; 
         stride >>= 1) {
        if (threadIdx.x < stride)
            s_data[threadIdx.x] += s_data[threadIdx.x + stride];
        __syncthreads();
    }
    if (threadIdx.x == 0)
        g_out[blockIdx.x] = s_data[0];
}

您会发现这个内核函数首先专注于累积输入数据，然后减少加载的数据。

现在，我们需要确定网格大小。为了使我们的 GPU 代码在各种 GPU 目标上运行，我们必须在运行时确定它们的大小。此外，我们需要利用 GPU 中的所有多处理器。CUDA C 提供了相关函数。我们可以使用cudaOccpancyMaxActiveBlocksPerMultiprocessor()函数获得占用率感知的每个多处理器的最大活动块数。此外，我们可以使用cudaDeviceGetAttribte()函数获得目标 GPU 上的多处理器数量。以下代码显示了如何使用这些函数并调用内核函数：

int reduction(float *g_outPtr, float *g_inPtr, int size, int n_threads)
{
    int num_sms;
    int num_blocks_per_sm;
    cudaDeviceGetAttribute(&num_sms, 
                           cudaDevAttrMultiProcessorCount, 0);
    cudaOccupancyMaxActiveBlocksPerMultiprocessor(&num_blocks_per_sm, 
           reduction_kernel, n_threads, n_threads*sizeof(float));
    int n_blocks = min(num_blocks_per_sm * num_sms, (size 
                       + n_threads - 1) / n_threads);

    reduction_kernel<<<n_blocks, n_threads, n_threads * 
                       sizeof(float), 0>>>(g_outPtr, g_inPtr, size);
    reduction_kernel<<<1, n_threads, n_threads * sizeof(float), 
                       0>>>(g_outPtr, g_inPtr, n_blocks);
    return 1;
}

这个函数还有一个额外的修改。为了节省占用率计算开销，它再次启动reduction_kernel()函数，这次只使用一个块。运行以下命令：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction ./reduction.cpp ./reduction_kernel.cu

更新后的减少性能为0.278 ms，在 Tesla V100 上比以前的方法快了大约 100 ms。

现在，让我们回顾一下我们如何利用内存带宽。以下图表显示了在 Visual Profiler 中的内存利用分析，并显示了我们如何将内存带宽增加了两倍：

尽管它显示出了增加的带宽，但我们仍然有进一步增加的空间。让我们来看看如何实现更多的带宽。

平衡 I/O 吞吐量

从分析器得到的结果来看，局部变量 input 有大量的加载/存储请求。这样大量的 I/O 会影响线程块的调度，因为存在操作依赖。当前数据累积中最糟糕的是它对设备内存有依赖。因此，我们将使用额外的寄存器来发出更多的加载指令以减轻依赖。以下代码显示了我们如何做到这一点：

#define NUM_LOAD 4
__global__ void
 reduction_kernel(float *g_out, float *g_in, unsigned int size)
{
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;

    extern __shared__ float s_data[];

    // cumulates input with grid-stride loop 
    // and save to the shared memory
    float input[NUM_LOAD] = {0.f};
    for (int i = idx_x; i < size; i += blockDim.x * 
         gridDim.x * NUM_LOAD)
    {
        for (int step = 0; step < NUM_LOAD; step++)
            input[step] += (i + step * blockDim.x * gridDim.x < size) ? 
                g_in[i + step * blockDim.x * gridDim.x] : 0.f;
    }
    for (int i = 1; i < NUM_LOAD; i++)
        input[0] += input[i];
    s_data[threadIdx.x] = input[0];

    __syncthreads();

    // do reduction
    for (unsigned int stride = blockDim.x / 2; stride > 0; 
         stride >>= 1)
    {
        if (threadIdx.x < stride)
            s_data[threadIdx.x] += s_data[threadIdx.x + stride];

        __syncthreads();
    }

    if (threadIdx.x == 0) {
        g_out[blockIdx.x] = s_data[0];
    }
}

这段代码使用了三个额外的寄存器来收集全局内存数据。NUM_LOAD的值可能会因 GPU 的不同而有所不同，因为它受 GPU 的内存带宽和 GPU 中 CUDA 核心数量的影响：

运行以下命令时，使用 Tesla V100 卡的性能达到了0.264毫秒：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction ./reduction.cpp ./reduction_kernel_opt.cu

warp 级原语编程

CUDA 9.0 引入了新的 warp 同步编程。这一重大变化旨在避免 CUDA 编程依赖隐式 warp 同步操作，并明确处理同步目标。这有助于防止 warp 级同步操作中的疏忽竞争条件和死锁。

从历史上看，CUDA 只提供了一个显式同步 API，即__syncthreads()用于线程块中的 CUDA 线程，并依赖于 warp 的隐式同步。下图显示了 CUDA 线程块操作的两个级别的同步：

然而，最新的 GPU 架构（Volta 和 Turing）具有增强的线程控制模型，其中每个线程可以执行不同的指令，同时保持其 SIMT 编程模型。下图显示了它是如何改变的：

直到 Pascal 架构（左图），线程是在 warp 级别进行调度的，并且它们在 warp 内部隐式同步。因此，CUDA 线程在 warp 中隐式同步。然而，这可能会导致意外的死锁。

Volta 架构对此进行了改进，并引入了独立线程调度。这种控制模型使每个 CUDA 线程都有自己的程序计数器，并允许 warp 中的一组参与线程。在这个模型中，我们必须使用显式的同步 API 来指定每个 CUDA 线程的操作。

因此，CUDA 9 引入了显式的 warp 级原语函数：

	warp 级原语函数
识别活动线程	`__activemask()`
屏蔽活动线程	`__all_sync()`, `__any_sync()`, `__uni_sync()`, __ballot_sync()``__match_any_sync(), `__match_all_sync()`
同步数据交换	`__shfl_sync()`, `__shfl_up_sync()`, `__shfl_down_sync()`, `__shfl_xor_sync()`
线程同步	`__syncwarp()`

有三类 warp 级原语函数，分别是 warp 识别、warp 操作和同步。所有这些函数都隐式地指定了同步目标，以避免意外的竞争条件。

使用 warp 原语进行并行归约

让我们看看这如何有益于我们的并行归约实现。这个示例将使用 Cooperative Groups 中的shfl_down()函数和 warp 原语函数中的shfl_down_sync()。下图显示了shfl_down_sync()如何与 shift down 操作一起工作：

在这个集体操作中，warp 中的 CUDA 线程可以将指定的寄存器值移动到同一个 warp 中的另一个线程，并与其同步。具体来说，集体操作有两个步骤（第三个是可选的）：

识别、屏蔽或投票源 CUDA 线程在一个 warp 中将进行操作。
让 CUDA 线程移动数据。
warp 中的所有 CUDA 线程都在同步（可选）。

对于并行归约问题，我们可以使用__shfl_down_sync()进行 warp 级别的归约。现在，我们可以通过以下图来增强我们的线程块级别的归约：

每个 warp 的归约结果都存储在共享内存中，以与其他 warp 共享。然后，通过再次进行 warp-wise 收集，可以获得最终的块级归约。

我们使用__shfl_down_sync()，因为我们只需要一个线程进行 warp 级别的归约。如果您需要让所有 CUDA 线程都进行 warp 级别的归约，可以使用__shfl_xor_sync()。

第一个块级别的归约的数量是网格的维度，输出存储在全局内存中。通过再次调用，我们可以使用 warp 级别的同步函数构建一个并行归约核。

现在，让我们使用 warp 级别的原始函数来实现 warp 级别的归约。首先，我们将编写一个使用 warp-shifting 函数进行 warp 级别归约的函数。以下代码显示了如何实现这一点：

__inline__ __device__ float warp_reduce_sum(float val) {
    for (int offset = warpSize / 2; offset > 0; offset >>= 1) {
        unsigned int mask = __activemask();
        val += __shfl_down_sync(mask, val, offset);
    }
    return val;
}

对于 warp-shifting，我们需要让 CUDA 调度程序识别活动线程，并让 warp-shifting 函数进行归约。

第二步是使用先前的 warp 级别归约编写一个块级别的归约函数。我们将在共享内存中收集先前的结果，并从结果中进行第二次归约。以下代码显示了如何实现这一点：

__inline__ __device__ float block_reduce_sum(float val) {
    // Shared mem for 32 partial sums
    static __shared__ float shared[32]; 
    int lane = threadIdx.x % warpSize;
    int wid = threadIdx.x / warpSize;

    val = warp_reduce_sum(val); // Warp-level partial reduction
    if (lane == 0)
        shared[wid] = val; // Write reduced value to shared memory
    __syncthreads(); // Wait for all partial reductions

    //read from shared memory only if that warp existed
    if (wid == 0) {
        val = (threadIdx.x < blockDim.x / warpSize) ? shared[lane] : 0;
        val = warp_reduce_sum(val); //Final reduce within first warp
    }
    return val;
}

现在，我们将实现归约核函数，累积输入数据，并从我们实现的块级归约中进行归约。由于我们只关注优化 warp 级别的优化，因此整体设计与之前的版本相同。以下代码显示了核函数：

__global__ void
reduction_kernel(float *g_out, float *g_in, unsigned int size) {
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
    // cumulates input with grid-stride loop and save to share memory
    float sum[NUM_LOAD] = { 0.f };
    for (int i = idx_x; i < size; i += blockDim.x * gridDim.x * NUM_LOAD) {
        for (int step = 0; step < NUM_LOAD; step++)
            sum[step] += (i + step * blockDim.x * gridDim.x < size) ? g_in[i + step * blockDim.x * gridDim.x] : 0.f;
    }
    for (int i = 1; i < NUM_LOAD; i++)
        sum[0] += sum[i];
    // warp synchronous reduction
    sum[0] = block_reduce_sum(sum[0]);

    if (threadIdx.x == 0)
        g_out[blockIdx.x] = sum[0];
}

然后，让我们使用以下命令编译代码：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction ./reduction.cpp ./reduction_wp_kernel.cu

以下屏幕截图显示了执行时间的减少：

在主机代码修改上，没有从 warp 原语切换到协作组。因此，我们可以对两种归约实现使用相同的主机代码。

我们已经介绍了 CUDA 中的 warp 同步编程。它的应用不仅限于归约，还可以用于其他并行算法：扫描、双调排序和转置。如果您需要了解更多信息，可以查看以下文章：

灵活处理线程的协作组

CUDA 9.0 引入了一个名为协作组的新 CUDA 编程特性。这通过指定组操作来引入了一种新的 CUDA 编程设计模式。使用这个特性，程序员可以编写显式控制 CUDA 线程的 CUDA 代码。

首先，让我们看看协作组是什么，以及它的编程优势。

CUDA 线程块中的协作组

协作组提供了显式的 CUDA 线程分组对象，帮助程序员更清晰、更方便地编写集体操作。例如，我们需要获取一个掩码来控制 warp 中活动的 CUDA 线程以进行 warp-shifting 操作。另一方面，协作组对象将可用的线程绑定为一个瓦片，并将它们作为一个对象进行控制。这为 CUDA C 编程带来了 C++语言的好处。

协作组的基本类型是thread_group。这使得 C++类样式的类型thread_group能够提供其配置信息，使用is_valid()、size()和thread_rank()函数。此外，这提供了可以应用于组中所有 CUDA 线程的集体函数。这些函数如下：

	thread_group 集体函数
识别活动线程	`tiled_partition()`, `coalesced_threads()`
屏蔽活动线程	`any()`, `all()`, ballot()``match_any(), `match_all()`
同步数据交换	`shfl()`, `shfl_up()`, `shfl_down()`, `shfl_xor()`
线程同步	`sync()`

这些函数列表类似于 warp 级别的原始函数。因此，warp 级别的原始操作可以用协作组替换。thread_group可以被较小的thread_group、thread_block_tile或coalesced_group分割。

协作组还提供了线程块编程的灵活性。使用以下代码行，我们可以处理一个线程块：

thread_block block = this_thread_block();

thread_block提供了 CUDA 内置关键字包装函数，我们使用它来获取块索引和线程索引：

dim3 group_index();  // 3-dimensional block index within the grid
dim3 thread_index(); // 3-dimensional thread index within the block

我们可以使用this_thread_block()来获取一个线程块对象，如下所示：

thread_block block = this_thread_block();

现在，让我们看看协作组的好处与传统的 CUDA 内置变量相比有什么好处。

协作组的好处

使用协作组提供了更多的 C++可编程性，而不是使用传统的 CUDA 内置变量。使用thread_block组，您可以将您的内核代码从使用内置变量切换到协作组的索引。但是，协作组的真正力量不仅仅如此。让我们在以下部分介绍其优势。

模块化

使用协作组，程序员可以将集体操作的内核代码模块化，对应于屏障目标。这有助于避免假设所有线程都在同时运行而导致的疏忽，从而引发死锁和竞争条件。以下是 CUDA 线程同步的死锁和正常操作的示例：

对于左侧的示例，内核代码意图同步 CUDA 线程块中的一部分线程。通过指定屏障目标，此代码最小化了同步开销。然而，它引入了死锁情况，因为__syncthreads()调用了一个屏障，等待所有 CUDA 线程到达屏障。然而，__synchthroead()无法满足其他线程的要求并等待。右侧的示例显示了良好的操作，因为它没有任何死锁点，因为线程块中的所有线程都可以满足__syncthreads()。

在协作组 API 中，CUDA 程序员指定线程组进行同步。协作组使得显式同步目标成为可能，因此程序员可以让 CUDA 线程显式同步。这个项目也可以被视为一个实例，因此我们可以将实例传递给设备函数。

以下代码显示了协作组如何提供显式同步对象并将它们作为实例处理：

__device__ bar(thread_group block, float *x) {
    ...
    block.sync();
    ...
}

__global__ foo() {
    bar(this_thread_block(), float *x);
}

正如在前面的示例代码中所示，内核代码可以指定同步组并将它们作为thread_group参数传递。这有助于我们在子例程中指定同步目标。因此，程序员可以通过使用协作组来防止意外死锁。此外，我们可以将不同类型的组设置为thread_group类型并重用同步代码。

显式分组线程的操作和避免竞争条件

协作组通过在 warp 中平铺线程来支持 warp 级协作操作。如果 tile 大小与 warp 大小匹配，CUDA 可以省略 warps 的隐式同步，确保正确的内存操作以避免竞争条件。通过消除隐式同步，可以增强 GPU 的性能。从历史上看，有经验的 CUDA 程序员使用分离的 warps 进行 warp 级同步。这意味着 warp 中的协作操作不必与其他 warp 操作同步。这释放了 GPU 的性能。但是，这是有风险的，因为它引入了协作操作之间的竞争条件。

动态活动线程选择

CUDA Cooperative Groups 的另一个好处是程序员可以选择 warp 中的活动线程，以避免分支分歧效应。由于 CUDA 是 SIMT 架构，一个指令单元发出一组线程，并且如果它们遇到分支，就无法禁用分歧。但是，从 CUDA 9.0 开始，程序员可以使用coalesced_threads()选择在分支块中活动的线程。这通过禁用不参与分支的线程返回聚合的线程。然后，SM 的指令单元发出下一个活动线程组中的活动线程。

应用于并行减少

我们将更新以前的减少内核代码以使用协作组。从以前的内核代码中，您可以轻松应用协作组的thread_block，如下所示：

__global__ void
 reduction_kernel(float* g_out, float* g_in, unsigned int size)
{
    unsigned int idx_x = blockIdx.x * blockDim.x + threadIdx.x;

    thread_block block = this_thread_block();

    extern __shared__ float s_data[];

我们不必更新数据输入累积部分，因此让我们为每个线程块更新减少部分。以下代码显示了一个块大小的减少的示例：

    // do reduction
    for (unsigned int stride = block.group_dim().x / 2; stride > 0; 
         stride >>= 1) {
        if (block.thread_index().x < stride) {
            s_data[block.thread_index().x] += 
                s_data[block.thread_index().x + stride];
            block.sync(); // threads synchronization in a branch
        }
    }
}

使用以下命令的估计操作性能为 0.264 毫秒：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction_cg -rdc=true ./reduction.cpp ./reduction_cg_kernel.cu

前面的命令显示了与以前版本相同的性能。

协作组以避免死锁

协作组可以支持独立的 CUDA 线程调度。因此，我们可以使用一个组单独控制 CUDA 线程，并显式地对它们进行同步。目标组可以是预定义的 tile，但也可以根据条件分支确定，如下所示：

// do reduction
for (unsigned int stride = block.group_dim().x / 2; stride > 0; 
     stride >>= 1) {
    // scheduled threads reduce for every iteration
    // and will be smaller than a warp size (32) eventually.
    if (block.thread_index().x < stride) { 
        s_data[block.thread_index().x] += s_data[
                       block.thread_index().x + stride];

        // __syncthreads(); // (3) Error. Deadlock.
        // block.sync();    // (4) Okay. Benefit of Cooperative Group
    }
    // __syncthreads();     // (1) Okay
    block.sync();           // (2) Okay
}

这段代码有四种线程块同步选项。选项(1)和(2)是具有不同 API 的等效操作。另一方面，选项(3)和(4)则不是。选项(3)引入了 CUDA 线程的死锁，主机无法返回 CUDA 内核，因为活动的 CUDA 线程无法与未激活的 CUDA 线程同步。另一方面，选项(4)由于协作组的自动活动线程识别而起作用。这有助于我们避免意外错误并轻松开发复杂的算法。

NVIDIA 提供了有关协作组的详细描述，可以在以下文档中找到：

您还可以从cooperative_groups.h本身了解其架构和完整的 API 列表。

CUDA 内核中的循环展开

CUDA 也可以像其他编程语言一样受益于循环展开。通过这种技术，CUDA 线程可以减少或消除循环控制开销，例如循环结束测试每次迭代，分支惩罚等。

如果 CUDA 可以识别循环的迭代次数，它会自动展开小循环。程序员还可以使用#pragma unroll指令向编译器提供提示，或者将循环代码重写为一组独立的语句。应用循环展开很简单，因此您可以轻松应用到当前的工作代码中。

让我们将这应用到我们的并行减少实现中。就像 C/C++中的普通循环展开指令一样，我们可以在for循环的顶部放置#pragma循环展开指令。NVCC 编译器可以展开循环，因为编译器可以自行获得group.size()的确切大小：

template <typename group_t>
__inline__ __device__ float
 warp_reduce_sum(group_t group, float val)
{
    #pragma unroll
    for (int offset = group.size() / 2; offset > 0; offset >>= 1)
        val += group.shfl_down(val, offset);
    return val;
}

使用以下命令，估计的操作性能为 0.263 毫秒：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction_cg -rdc=true ./reduction.cpp ./reduction_cg_kernel.cu

如果您更喜欢使用 warp 原始函数，可以像下面这样编写warp_reduce_sum。循环代码可以通过用warpSize替换group.size()来重用，但在这种情况下稍微更快：

#define FULL_MASK 0xFFFFFFFF
__inline__ __device__ float
warp_reduce_sum(float val) {
#pragma unroll 5
    for (int offset = 1; offset < 6; offset++)
        val += __shfl_down_sync(FULL_MASK, val, warpSize >> offset);
    return val;
}

运行以下命令来编译上述代码：

nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction_wp -rdc=true ./reduction.cpp ./reduction_wp_kernel.cu

其结果是 0.263 毫秒，与先前的结果相同。

使用循环展开存在一个陷阱。展开的代码执行可能导致寄存器使用增加而降低占用率。此外，由于代码执行大小增加，可能会出现更高的指令缓存未命中惩罚。

原子操作

在 CUDA 编程中，程序员可以使用原子 API 从多个 CUDA 线程更新共享资源。这些原子 API 保证消除对共享资源的竞争条件，因此我们可以期望并行执行产生一致的输出。这个操作对于获取统计参数（如直方图、均值、总和等）特别有用。我们还可以简化代码实现。例如，可以使用以下代码中的atomicAdd()函数编写减少操作：

__global__ void
 atomic_reduction_kernel(float *data_out, float *data_in, int size)
 {
     int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
     atomicAdd(&data_out[0], data_in[idx_x]);
 }

正如您所看到的，原子函数简化了所需的操作。但是，由于原子操作将所有请求串行化到共享资源，因此其性能较慢。运行以下命令查看执行时间：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o mixed_precision_single ./mixed_precision.cu

这个显示的内核函数在我的 Tesla V100 上花费了 39 毫秒，比原始版本（4.609 毫秒）慢得多。因此，建议的原子操作使用是只在必要时限制请求。例如，对于并行减少问题，我们可以在某个级别并行减少项目，并使用原子操作输出最终结果。

下图显示了另一种可能的方法。这将块内减少替换为atomicAdd：

在上图中，我们可以看到有两个减少点：一个是warp，一个是线程块，并且，块内减少结果通过单个全局内存变量原子地累积。因此，我们可以消除第二次减少迭代。以下截图显示了第二次减少迭代的内核优化优先级（左侧）和性能限制分析（右侧）：

内核优化优先级与性能限制分析（第二次迭代）

换句话说，第二次迭代的性能受其较小的网格大小的延迟限制。因此，通过删除这一点，我们将能够减少执行时间。

现在让我们实现该设计并看看性能如何改变。我们只需要更新减少内核函数的最后部分：

__global__ void
 reduction_kernel(float* g_out, float* g_in, unsigned int size)
{
    unsigned int idx_x = blockIdx.x * (2 * blockDim.x) + threadIdx.x;

    thread_block block = this_thread_block();

    // cumulates input with grid-stride loop and save to share memory
    float sum[NUM_LOAD] = { 0.f };
    for (int i = idx_x; i < size; i += blockDim.x 
         * gridDim.x * NUM_LOAD)
    {
        for (int step = 0; step < NUM_LOAD; step++)
            sum[step] += (i + step * blockDim.x * gridDim.x < size) ? 
                         g_in[i + step * blockDim.x * gridDim.x] : 0.f;
    }
    for (int i = 1; i < NUM_LOAD; i++)
        sum[0] += sum[i];
    // warp synchronous reduction
    sum[0] = block_reduce_sum(block, sum[0]);

    sum[0] = block_reduce_sum(sum[0]);

    // Performing Atomic Add per block
    if (block.thread_rank() == 0) {
        atomicAdd(&g_out[0], sum);
    }
}

然后，我们将删除第二次迭代函数调用。因此，如果原子操作的延迟短于那个，我们可以消除内核调用延迟并实现更好的性能。运行以下命令：

$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o reduction_atomic_block ./reduction.cpp ./reduction_blk_atmc_kernel.cu

幸运的是，在 Tesla V100 上估计的执行时间为 0.259 毫秒，因此我们可以获得稍微增强的结果。

如果您想了解 CUDA C 中原子操作的更多信息，请查看此链接的编程指南：docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomic-functions。

低/混合精度操作

混合精度是一种探索低精度并获得高精度结果的技术。这种技术使用低精度计算核心操作，并使用高精度操作生成输出。与高精度计算相比，低精度操作计算具有减少内存带宽和更高的计算吞吐量的优势。如果低精度足以从具有高精度的应用程序中获得目标精度，这种技术可以通过这种权衡来提高性能。NVIDIA 开发者博客介绍了这种可编程性：devblogs.nvidia.com/mixed-precision-programming-cuda-8。

在这些情况下，CUDA 将其支持扩展到低于 32 位数据类型的低精度工具，例如 8/16 位整数（INT8/INT16）和 16 位浮点数（FP16）。对于这些低精度数据类型，GPU 可以使用一些特定的 API 进行单指令，多数据（SIMD）操作。在本节中，我们将研究这两种用于混合精度目的的低精度操作的指令。

要从中受益，您需要确认您的 GPU 是否支持低混合精度操作和支持的数据类型。特定 GPU 支持低精度计算是可能的，精度取决于 GPU 芯片组。具体来说，GP102（Tesla P40 和 Titan X），GP104（Tesla P4）和 GP106 支持 INT8；GP100（Tesla P100）和 GV100（Tesla V100）支持 FP16（半精度）操作。Tesla GV100 兼容 INT8 操作，没有性能下降。

CUDA 具有一些特殊的内置函数，可以为低精度数据类型启用 SIMD 操作。

半精度操作

CUDA 为半精度浮点数据类型（FP16）提供了内置函数，并且开发人员可以选择 CUDA 是否为每条指令计算一个或两个值。CUDA 还提供了单精度和半精度之间的类型转换函数。由于 FP16 的精度限制，您必须使用转换内置函数来处理单精度值。

现在，让我们实现和测试 GPU 的 FP16 操作。GPU 可以支持高于计算能力 5.3 的这种类型的本机计算。但是有些 GPU 不支持这一点，因此请仔细检查您的 GPU 是否支持这种半精度操作。

CUDA C 中的半精度数据类型是half，但您也可以使用__half类型。对于 API，CUDA 提供了与此数据类型相关的内置函数，例如__hfma()、__hmul()和__hadd()。这些内置函数还提供了使用__hfma2()、__hmul2()和__hadd2()一次处理两个数据的本机操作。使用这些函数，我们可以编写混合精度操作的核心代码：

__global__ void hfma_kernel(half *d_x, half *d_y, float *d_z, int size)
 {
     int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
     int stride = gridDim.x * blockDim.x;

     half2 *dual_x = reinterpret_cast<half2*>(d_x);
     half2 *dual_y = reinterpret_cast<half2*>(d_y);
     float2 *dual_z = reinterpret_cast<float2*>(d_z);

     extern __shared__ float2 s_data[];

 #if __CUDA_ARCH__ >= 530
     for (int i = idx_x; i < size; i+=stride) {
         s_data[threadIdx.x] = __half22float2(__hmul2(dual_y[i], 
                                                      dual_x[i]));
         __syncthreads();
         dual_z[i] = s_data[threadIdx.x];
     }
     #else
     for (int i = idx_x; i < size; i+=stride) {
         s_data[threadIdx.x] = __half22float2(dual_x[i]) * 
                               __half22float2(dual_y[i]);
         __syncthreads();
         dual_z[i] = s_data[threadIdx.x];
     }
     #endif
 }

对于那些不支持本机半精度操作的 GPU，我们的代码在编译时检查 CUDA 的计算能力，并确定应采取哪种操作。

以下代码调用了半精度网格大小的核函数，因为每个 CUDA 线程将操作两个数据：

int n_threads = 256;
int num_sms;
int num_blocks_per_sm;
cudaDeviceGetAttribute(&num_sms, cudaDevAttrMultiProcessorCount, 0);
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&num_blocks_per_sm,   
    hfma_kernel, n_threads, n_threads*sizeof(float2));
int n_blocks = min(num_blocks_per_sm * num_sms, 
                   (size/2 + n_threads - 1) / n_threads);
hfma_kernel<<< n_blocks, n_threads, n_threads * sizeof(float2) 
               >>>(X.d_ptr_, Y.d_ptr_, Z.d_ptr_, size/2);

其他初始化代码和基准代码在示例配方代码中实现，因此请查看它。

我们已经涵盖了 FP16 精度操作中的 FMA 操作。CUDA C 提供了各种半精度操作（docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH__INTRINSIC__HALF.html）。请查看其他操作。

8 位整数和 16 位数据的点积运算和累加（DP4A 和 DP2A）

对于 8 位/16 位整数，CUDA 提供了矢量化的点积操作。这些是 DP4A（四元素点积累加）和 DP2A（两元素点积累加）。使用这些函数，CUDA 开发人员可以进行更快的操作。CUDA 8.0 开发博客通过直观的图示介绍了这些函数 (devblogs.nvidia.com/mixed-precision-programming-cuda-8/)。以下显示了 GPU 的点积和累加操作的工作原理：

使用这个，你可以编写只有 8 位或 8 位/16 位混合操作的 32 位整数累加。其他操作，如求和、加法和比较，也可以使用 SIMD 内在函数。

如前所述，有特定的 GPU 可以支持 INT8/INT16 操作，具有特殊功能（dp4a和dp2a）。支持的 GPU 的计算能力必须高于 6.1。

现在，让我们实现一个使用dp4aAPI 的内核函数，如下所示：

__global__ void dp4a_kernel(char *d_x, char *d_y, int *d_z, int size)
 {
     int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
     int stride = gridDim.x * blockDim.x;

 #if __CUDA_ARCH__ >= 610
     char4 *quad_x = (char4 *)d_x;
     char4 *quad_y = (char4 *)d_y;

     for (int i = idx_x; i < size; i+=stride)
         d_z[i] = __dp4a(quad_y[i], quad_x[i], 0);
 #else
     for (int i = idx_x; i < size; i+=4*stride) {
         int sum = 0;
         for (int j = 0; j < 4; j++)
             sum += d_y[4 * i + j] * d_x[4 * i + j];
         d_z[i] = sum + 0;
     }
 #endif
 }

在这个函数中，__dp4a获取两个字符数组，合并四个项目，并输出其点积输出。自帕斯卡以来，这个 API 就得到了支持，具有 CUDA 计算能力（版本 6.1）。但是旧的 GPU 架构，低于版本 6.1，需要使用原始操作。

以下代码显示了我们将如何调用实现的内核函数。由于每个 CUDA 线程将操作四个项目，其网格大小减小了四倍：

int n_threads = 256;
int num_sms;
int num_blocks_per_sm;
cudaDeviceGetAttribute(&num_sms, cudaDevAttrMultiProcessorCount, 0);
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&num_blocks_per_sm, 
    dp4a_kernel, n_threads, n_threads*sizeof(int));
int n_blocks = min(num_blocks_per_sm * num_sms, (size/4 + n_threads 
                                                  - 1) / n_threads);
dp4a_kernel<<< n_blocks, n_threads, n_threads * sizeof(int) >>>  
    (X.d_ptr_, Y.d_ptr_, Z.d_ptr_, size/4);

其他初始化代码和基准代码都在示例代码中实现，就像前面的示例代码一样。

我们已经介绍了 INT8 的点操作，但 CUDA C 还提供了其他 INT8 类型的 SIMD 内在函数（docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH__INTRINSIC__SIMD.html）。请查阅此文档以了解其他操作。

性能测量

示例代码有三个混合精度操作的版本：单精度、半精度和 INT8。随着精度的降低，我们可以为每个 CUDA 线程添加更多的操作。

运行以下命令进行单精度、半精度和 INT8 操作：

# Single-precision
$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o mixed_precision_single ./mixed_precision.cu

# Half-precision
$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o mixed_precision_half ./mixed_precision_half.cu

# INT8 
$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc -o mixed_precision_int ./mixed_precision_int.cu

以下表格显示了每种精度操作的估计性能：

精度	测得的性能
FP32	59.441 GFlops
FP16	86.037 GFlops
INT8	196.225 Gops

由于我们的实现没有经过优化，所以测得的性能比 Tesla V100 的理论性能要低得多。当你对它们进行分析时，它们会报告它们的内存绑定性很高。换句话说，我们需要优化它们，使它们在算术上受限，以接近理论性能。

总结

在本章中，我们介绍了如何配置 CUDA 并行操作并对其进行优化。为了做到这一点，我们必须了解 CUDA 的分层体系结构线程块和流多处理器之间的关系。通过一些性能模型——占用率、性能限制分析和 Roofline 模型——我们可以优化更多性能。然后，我们介绍了一些新的 CUDA 线程可编程性，合作组，并学习了如何简化并行编程。我们优化了并行减少问题，并在元素中实现了 0.259 毫秒，这是与相同 GPU 相比速度提高了 17.8。最后，我们了解了 CUDA 的半精度（FP16）和 INT8 精度的 SIMD 操作。

本章我们的经验集中在 GPU 的并行处理级别编程上。然而，CUDA 编程包括系统级编程。基本上，GPU 是额外的计算资源，独立于主机工作。这增加了额外的计算能力，但另一方面也可能引入延迟。CUDA 提供了可以利用这一点并隐藏延迟并实现 GPU 的全面性能的 API 函数。我们将在下一章中介绍这一点。

CUDA-编程学习手册（一）

CUDA 编程学习手册（一）

前言

这本书适合谁

为了充分利用这本书

下载示例代码文件

下载彩色图像

使用的约定

CUDA 编程简介

高性能计算的历史

异构计算

编程范式

低延迟与高吞吐量

GPU 的编程方法

技术要求

来自 CUDA 的 Hello World

线程层次结构

GPU 架构

使用 CUDA 进行向量加法

实验 1 - 创建多个块

实验 2 - 创建多个线程

实验 3 - 结合块和线程

为什么要费心处理线程和块？

在多个维度中启动内核

CUDA 中的错误报告

CUDA 中的数据类型支持

总结

CUDA 内存管理

技术要求

NVIDIA Visual Profiler

全局内存/设备内存

全局内存上的矢量加法

合并与未合并的全局内存访问

内存吞吐量分析

共享内存

共享内存上的矩阵转置

银行冲突及其对共享内存的影响

只读数据/缓存

计算机视觉-使用纹理内存进行图像缩放

GPU 中的寄存器

固定内存

带宽测试-固定与分页

统一内存

了解统一内存页面分配和传输

使用每页 warp 优化统一内存

统一内存的优化使用数据预取

GPU 内存的演变

为什么 GPU 有缓存？

总结

CUDA 线程编程

技术要求

CUDA 线程、块和 GPU

利用 CUDA 块和 warp

理解 CUDA 占用率

设置 NVCC 报告 GPU 资源使用

Linux 设置

Windows 设置

使用占用率计算器分析最佳占用率

占用率调整 - 限制寄存器使用

从分析器获取实现的占用

理解并行归约

使用全局内存的天真并行归约

使用共享内存减少内核

编写性能测量代码

两种减少-全局和共享内存的性能比较

识别应用程序的性能限制器

找到性能限制因素并进行优化

最小化 CUDA warp 分歧效应

确定分歧作为性能瓶颈

交错寻址

顺序寻址

性能建模和平衡限制器

Roofline 模型

通过网格跨步循环最大化内存带宽

平衡 I/O 吞吐量

warp 级原语编程

使用 warp 原语进行并行归约

灵活处理线程的协作组

CUDA 线程块中的协作组

协作组的好处