CUDA-编程学习手册（三）

布客飞龙
104 阅读23分钟

CUDA 编程学习手册（三）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

可扩展的多 GPU 编程

到目前为止，我们一直致力于在单个 GPU 上获得最佳性能。密集节点与多个 GPU 已成为即将到来的超级计算机的迫切需求，特别是自从 ExaFLOP（每秒千亿次操作）系统成为现实以来。 GPU 架构具有高能效，因此近年来，具有 GPU 的系统在 Green500 榜单（www.top500.org/green500）中占据了大多数前十名。在 2018 年 11 月的 Green500 榜单中，前十名中有七个基于 NVIDIA GPU。

NVIDIA 的 DGX 系统现在在一个服务器中有 16 个 V100 32GB。借助统一内存和诸如 NVLink 和 NvSwitch 之类的互连技术，开发人员可以将所有 GPU 视为一个具有 512GB 内存的大型 GPU（16 个 GPU *每个 32GB）。在本章中，我们将深入讨论编写 CUDA 代码的细节，并利用 CUDA-aware 库在多 GPU 环境中实现节点内和节点间的可伸缩性。

在本章中，我们将涵盖以下主题：

  • 使用高斯消元法解线性方程

  • GPUDirect 点对点

  • MPI 简介

  • GPUDirect RDMA

  • CUDA 流

  • 额外的技巧

技术要求

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

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

由于本章需要展示多 GPU 的交互，我们需要至少两个相同类型和架构的 GPU。还要注意，一些功能，如 GPUDirect RDMA 和 NVLink，仅支持 NVIDIA 的 Tesla 卡。如果您没有像 Tesla P100 或 Tesla V100 这样的 Tesla 卡，不要灰心。您可以安全地忽略其中一些功能。与我们在这里展示的情况相比，性能数字将会有所变化，但相同的代码将仍然有效。

在下一节中，我们将看一个示例，使用流行的高斯算法解决一系列线性方程，以演示如何编写多 GPU。

使用高斯消元法解线性方程

为了演示在节点内和节点间使用多个 GPU，我们将从一些顺序代码开始，然后将其转换为节点内和节点间的多个 GPU。我们将解决一个包含M个方程和N个未知数的线性方程组。该方程可以表示如下：

A × x = b

在这里，A是一个具有M行和N列的矩阵，x是一个列向量（也称为解向量），具有N行，b也是一个具有M行的列向量。找到解向量涉及在给定A和b时计算向量x。解线性方程组的标准方法之一是高斯消元法。在高斯消元法中，首先通过执行初等行变换将矩阵A减少为上三角矩阵或下三角矩阵。然后，通过使用回代步骤解决得到的三角形方程组。

以下伪代码解释了解线性方程所涉及的步骤：

1\. For iteration 1 to N (N: number of unknowns) 
    1.1 Find a row with non-zero pivot
    1.2 Extract the pivot row
    1.3 Reduce other rows using pivot row
2 Computing the solution vector through back substitution

让我们看一个示例，以便理解算法。假设方程组如下：

首先，我们将尝试设置基线系统，如下所示：

  1. 准备您的 GPU 应用程序。此代码可以在本书的 GitHub 存储库中的06_multigpu/gaussian文件夹中找到。

  2. 使用nvcc编译器编译您的应用程序，如下所示：

$ nvcc -o gaussian_sequential.out gaussian_sequential.cu
$ nvcc -o gaussian_single_gpu.out gaussian_single_gpu.cu
$ $ time ./gaussian_sequential.out
$ time ./gaussian_single_gpu.out

前面的步骤编译并运行了本章中存在的两个版本的代码：

  • 顺序运行的 CPU 代码

  • 在单个 GPU 上运行的 CUDA 代码

现在，让我们看看高斯消元的单 GPU 实现中的热点。

高斯消元的单 GPU 热点分析

让我们尝试理解和分析顺序和单 GPU 代码以建立基线。在此基础上，我们将增强并添加对多 GPU 运行的支持。

顺序 CPU 代码：以下代码显示了顺序实现的提取代码：

for( int n = 0; n < N; n++ ){
// M: number of equations, N: number of unknowns
    for( int pr = 0; pr < M; pr++ ){
        // finding the pivot row 
        //if pr satisfies condition for pivot i.e. is non zero 
        break; 
    }
    for( int r = 0; r < M; r++ ){
        // reduce all other eligible rows using the pivot row
        double ratio = AB[r*N+n]/AB[pr*N+n]
        for( int nn = n; nn < N + 1; nn++ ){
            AB[r * N + nn] -= (ratio*AB[pr * N + nn]);
        }
    }
}

从视觉上看，发生的操作如下：

在这里，高斯消元中的行数等于方程的数量，列数等于未知数的数量。在前面的图表中显示的pr行是主元行，将用于使用主元素减少其他行。

我们可以做出的第一个观察是，我们正在对增广矩阵进行操作，将A矩阵与b向量合并。因此，未知数的大小为N+1，因为增广矩阵的最后一列是b向量。创建增广矩阵有助于我们只处理一个数据结构，即矩阵。您可以使用以下命令对此代码进行分析。分析结果将显示guassian_elimination_cpu()函数完成所需的时间最长：

$ nvprof --cpu-profiling on ./guassian_sequential.out

CUDA 单 GPU 代码：通过前几章的学习，我们期望您已经熟悉了如何编写最佳的 GPU 代码，因此我们不会详细介绍单个 GPU 实现。以下摘录显示，在单个 GPU 实现中，三个步骤被称为三个用于找到N未知数的核心：

  • findPivotRowAndMultipliers<<<...>>>：该核心查找主元行和乘数，应用于行消除。

  • extractPivotRow<<<>>>：该核心提取主元行，然后用于执行行消除。

  • rowElimination<<<>>>：这是最终的核心调用，在 GPU 上并行进行行消除。

以下代码片段显示了数据在复制到 GPU 后迭代调用的三个核心：

<Copy input augmented matrix AB to GPU>
...
for( int n = 0; n < N; n++ ){
// M: number of equations, N: number of unknowns
    findPivotRowAndMultipliers<<<...>>>(); 
    extractPivotRow<<<...>>>(); 
    rowElimination<<<...>>>(); 

}

本章的重点是如何增强此单个 GPU 实现以支持多个 GPU。但是，为了填补 GPU 实现中的缺失部分，我们需要对单个 GPU 实现进行一些优化更改：

  • 高斯消元算法的性能受内存访问模式的影响很大。基本上，它取决于 AB 矩阵的存储方式：

  • 找到主元行更喜欢列主格式，因为如果矩阵以列主格式存储，则提供了合并访问。

  • 另一方面，提取主元行更喜欢行主格式。

  • 无论我们如何存储AB矩阵，内存访问中都无法避免一个合并和一个跨步/非合并的访问。

  • 列主格式对于行消除核心也是有益的，因此对于我们的高斯消元核心，我们决定存储 AB 矩阵的转置而不是 AB。AB 矩阵在代码开始时通过transposeMatrixAB()函数转置一次。

在下一节中，我们将启用多 GPU P2P 访问并将工作分配给多个 GPU。

GPU 直接点对点

GPUDirect 技术是为了允许 GPU 在节点内部和跨不同节点之间进行高带宽、低延迟的通信而创建的。该技术旨在消除一个 GPU 需要与另一个 GPU 通信时的 CPU 开销。GPUDirect 可以分为以下几个主要类别：

  • GPU 之间的点对点（P2P）传输：允许 CUDA 程序在同一系统中的两个 GPU 之间使用高速直接内存传输（DMA）来复制数据。它还允许对同一系统中其他 GPU 的内存进行优化访问。

  • 网络和存储之间的加速通信：这项技术有助于从第三方设备（如 InfiniBand 网络适配器或存储）直接访问 CUDA 内存。它消除了不必要的内存复制和 CPU 开销，从而减少了传输和访问的延迟。此功能从 CUDA 3.1 开始支持。

  • 视频的 GPUDirect：这项技术优化了基于帧的视频设备的流水线。它允许与 OpenGL、DirectX 或 CUDA 进行低延迟通信，并且从 CUDA 4.2 开始支持。

  • 远程直接内存访问（RDMA）：此功能允许集群中的 GPU 之间进行直接通信。此功能从 CUDA 5.0 及更高版本开始支持。

在本节中，我们将把我们的顺序代码转换为使用 GPUDirect 的 P2P 功能，以便在同一系统中的多个 GPU 上运行。

GPUDirect P2P 功能允许以下操作：

  • GPUDirect 传输：cudaMemcpy()启动了从 GPU 1 的内存到 GPU 2 的内存的 DMA 复制。

  • 直接访问：GPU 1 可以读取或写入 GPU 2 的内存（加载/存储）。

以下图表展示了这些功能：

要理解 P2P 的优势，有必要了解 PCIe 总线规范。这是为了通过 InfiniBand 等互连优化与其他节点进行通信而创建的。当我们想要从单个 GPU 优化地发送和接收数据时，情况就不同了。以下是一个样本 PCIe 拓扑，其中八个 GPU 连接到各种 CPU 和 NIC/InfiniBand 卡：

在前面的图表中，GPU0 和 GPU1 之间允许 P2P 传输，因为它们都位于同一个 PCIe 交换机中。然而，GPU0 和 GPU4 不能执行 P2P 传输，因为两个I/O Hub（IOHs）之间不支持 PCIe P2P 通信。IOH 不支持来自 PCI Express 的非连续字节进行远程对等 MMIO 事务。连接两个 CPU 的 QPI 链路的性质确保了如果 GPU 位于不同的 PCIe 域上，则不可能在 GPU 内存之间进行直接 P2P 复制。因此，从 GPU0 的内存到 GPU4 的内存的复制需要通过 PCIe 链路复制到连接到 CPU0 的内存，然后通过 QPI 链路传输到 CPU1，并再次通过 PCIe 传输到 GPU4。正如你所想象的那样，这个过程增加了大量的开销，无论是延迟还是带宽方面。

以下图表显示了另一个系统，其中 GPU 通过支持 P2P 传输的 NVLink 互连相互连接：

前面的图表显示了一个样本 NVLink 拓扑，形成了一个八立方网格，其中每个 GPU 与另一个 GPU 最多相连 1 跳。

更重要的问题是，*我们如何找出这个拓扑结构以及哪些 GPU 支持 P2P 传输？*幸运的是，有工具可以做到这一点。nvidia-smi就是其中之一，它作为 NVIDIA 驱动程序安装的一部分被安装。以下屏幕截图显示了在前面图表中显示的 NVIDIA DGX 服务器上运行nvidia-smi的输出：

前面的屏幕截图代表了在具有 8 个 GPU 的 DGX 系统上运行nvidia-smi topo -m命令的结果。如您所见，通过 SMP 互连（QPI/UPI）连接到另一个 GPU 的任何 GPU 都无法执行 P2P 传输。例如，GPU0将无法与GPU5、GPU6和GPU7进行 P2P 传输。另一种方法是通过 CUDA API 来找出这种传输，我们将在下一节中使用它来转换我们的代码。

现在我们已经了解了系统拓扑，我们可以开始将我们的应用程序转换为单个节点/服务器上的多个 GPU。

单节点-多 GPU 高斯消元

准备您的多 GPU 应用程序。此代码可以在本书的 GitHub 存储库中的06_multigpu/gaussian中找到。使用nvcc编译器编译您的应用程序，如下所示：

$ nvcc -o gaussian_multi_gpu_p2p.out gaussian_multi_gpu_p2p.cu
$ time ./gaussian_multi_gpu_p2p.out

从单 GPU 实现转换为多 GPU 实现，我们在上一小节中定义的三个内核将被原样使用。但是，线性系统被分成与 GPU 数量相等的部分。这些部分分配给每个 GPU 一个部分。每个 GPU 负责对分配给该 GPU 的部分执行操作。矩阵是按列分割的。这意味着每个 GPU 从所有行中获得相等数量的连续列。用于找到主元的内核在包含主元素的列上启动。主元元素的行索引被广播到其他 GPU。提取的主元行和行消除内核在所有 GPU 上启动，每个 GPU 都在矩阵的自己的部分上工作。以下图显示了行在多个 GPU 之间的分割以及主元行需要广播到其他进程的情况：

上述图表示了在多个 GPU 上的工作分配。目前，主元行属于GPU1，负责将主元行广播到其他 GPU。

让我们试着理解这些代码更改，以及用于启用 P2P 功能的 CUDA API：

  1. 在支持的 GPU 之间启用 P2P 访问。以下代码显示了这个步骤的第一步：启用 GPU 之间的 P2P 访问：
for( int i = 0; i < nGPUs; i++ ){   
    // setup P2P 
    cudaSetDevice(i);   
    for( int j = 0; j < nGPUs; j++ ) {      
        if (i == j) continue;      
        cudaDeviceCanAccessPeer(&canAccessPeer, i, j);
        if (canAccessPeer)      
            cudaDeviceEnablePeerAccess(j, 0);    
    } 
}

在上述代码中使用的关键 API 如下：

    • cudaDeviceCanAccessPeer(): 检查当前 GPU 是否可以对传递的 GPU ID 进行 P2P 访问

  • cudaDeviceEnablePeerAccess(): 如果cudaDeviceCanAccessPeer()返回True，则启用 P2P 访问

  1. 拆分并将内容传输到各自的 GPU：
for( int g = 0; g < nGPUs; g++ ){       
    cudaSetDevice(g);       
    //Copy  part ‘g’ of ABT to GPU ‘g’; 
}

在上述代码中使用的关键 API 是cudaSetDevice()。这将当前上下文设置为作为参数传递的 GPU ID。

  1. 找到主元行并通过 P2P 进行广播：
for( int n = 0; n < N; n++ ){        
    gp = GPU that holds n;        
    cudaSetDevice(gp);        
    findPivotRowAndMultipliers<<<...>>>();
    for( int g = 0; g < nGPUs; g++ ){ 
        if (g == gp) continue;
        cudaMemcpyPeer(pivotDatag, g, pivotDatagp, gp, numBytes);
     }  ...

用于将传输广播到 GPU 的 API 是cudaMemcpyPeer()。

  1. 提取主元行并执行行消除：
for( int n = 0; n < N; n++ ){
    ...
    for( int g = 0; g < nGPUs; g++ ){  
        cudaSetDevice(g); 
        extractPivotRow<<<...>>>(); 
        rowElimination<<<...>>>();   
    }  
}

如您所见，我们仍在重用相同的内核。唯一的区别是我们使用cudaSetDevice() API 告诉 CUDA 运行时内核应该在哪个 GPU 上启动。请注意，cudaSetDevice()是一个昂贵的调用，特别是在旧一代的 GPU 上。因此，建议您通过在 CPU 上并行调用nGPUs的 for 循环，利用OpenMP/OpenACC或 CPU 上的任何其他线程机制来调用。

  1. 从各自的 CPU 中复制数据回来：
for( int g = 0; g < nGPUs; g++ ){ 
    cudaSetDevice(g);  
    Copy  part ‘g’ of reduced ABT from GPU ‘g’ to Host; 
}

这五个步骤完成了将单个 GPU 实现转换为单个节点上的多个 GPU 的练习。

作为 CUDA 安装的一部分提供的 CUDA 示例包括一些测试 P2P 带宽性能的示例代码。它可以在samples/1_Utilities/p2pBandwidthLatencyTest文件夹中找到。建议您在系统上运行此应用程序，以便了解系统的 P2P 带宽和延迟。

现在我们已经在单个节点上实现了多 GPU，我们将改变方向并在多个 GPU 上运行此代码。但在将我们的代码转换为多个 GPU 之前，我们将提供一个关于 MPI 编程的简短介绍，这主要用于节点间通信。

MPI 的简要介绍

消息传递接口（MPI）标准是一种消息传递库标准，已成为在 HPC 平台上编写消息传递程序的行业标准。基本上，MPI 用于在多个 MPI 进程之间进行消息传递。相互通信的 MPI 进程可以驻留在同一节点上，也可以跨多个节点。

以下是一个 Hello World MPI 程序的示例：

#include <mpi.h> 
int main(int argc, char *argv[]) {     
    int rank,size;     
    /* Initialize the MPI library */     
    MPI_Init(&argc,&argv);     
    /* Determine the calling process rank and total number of ranks */
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);     
    MPI_Comm_size(MPI_COMM_WORLD,&size);     
    /* Compute based on process rank */     
    /* Call MPI routines like MPI_Send, MPI_Recv, ... */     
    ...     
    /* Shutdown MPI library */     
    MPI_Finalize();     
    return 0; 
}

正如您所看到的，MPI 程序涉及的一般步骤如下：

  1. 我们包括头文件mpi.h，其中包括所有 MPI API 调用的声明。

  2. 我们通过调用MPI_Init并将可执行参数传递给它来初始化 MPI 环境。在这个语句之后，多个 MPI 等级被创建并开始并行执行。

  3. 所有 MPI 进程并行工作，并使用诸如MPI_Send()、MPI_Recv()等消息传递 API 进行通信。

  4. 最后，我们通过调用MPI_Finalize()终止 MPI 环境。

我们可以使用不同的 MPI 实现库（如 OpenMPI、MVPICH、Intel MPI 等）来编译此代码：

$ mpicc -o helloWorldMPI helloWorldMPI.c
$ mpirun -n 4 --hostfile hostsList ./helloWorldMPI

我们使用mpicc编译器来编译我们的代码。mpicc基本上是一个包装脚本，它在内部扩展编译指令，以包括相关库和头文件的路径。此外，运行 MPI 可执行文件需要将其作为参数传递给mpirun。mpirun是一个包装器，它帮助在应用程序应该执行的多个节点上设置环境。-n 4参数表示我们要运行四个进程，并且这些进程将在主机名存储在文件主机列表中的节点上运行。

在本章中，我们的目标是将 GPU 内核与 MPI 集成，使其在多个 MPI 进程中运行。但我们不会涵盖 MPI 编程的细节。那些不熟悉 MPI 编程的人应该先查看computing.llnl.gov/tutorials/mpi/，了解分布式并行编程，然后再进入下一节。

GPUDirect RDMA

在集群环境中，我们希望在多个节点上利用 GPU。我们将允许我们的并行求解器将 CUDA 代码与 MPI 集成，以利用多节点、多 GPU 系统上的多级并行性。使用 CUDA-aware MPI 来利用 GPUDirect RDMA 进行优化的节点间通信。

GPUDirect RDMA 允许在集群中的 GPU 之间进行直接通信。它首先由 CUDA 5.0 与 Kepler GPU 卡支持。在下图中，我们可以看到 GPUDirect RDMA，即Server 1中的GPU 2直接与Server 2中的GPU 1通信：

GPUDirect RDMA 工作的唯一理论要求是网络卡和GPU共享相同的根复杂性。 GPU 和网络适配器之间的路径决定了是否支持 RDMA。让我们重新访问我们在上一节中运行的 DGX 系统上nvidia-smi topo -m命令的输出：

如果我们看一下GPU4行，它显示GPU4到mlx5_2连接类型为PIX（通过 PCIe 交换机遍历）。我们还可以看到GPU4到mlx_5_0连接类型为SYS（通过QPI遍历）。这意味着GPU4可以通过 Mellanox InfiniBand 适配器mlx_5_2执行 RDMA 传输，但如果需要从mlx_5_0进行传输，则无法进行 RDMA 协议，因为QPI不允许。

CUDA-aware MPI

所有最新版本的 MPI 库都支持 GPUDirect 功能。支持 NVIDIA GPUDirect 和统一虚拟寻址（UVA）的 MPI 库使以下功能可用：

  • MPI 可以将 API 传输直接复制到/从 GPU 内存（RDMA）。

  • MPI 库还可以区分设备内存和主机内存，无需用户提示，因此对 MPI 程序员透明。

  • 程序员的生产率提高了，因为少量应用代码需要更改以在多个 MPI 秩之间传输数据。

正如我们之前提到的，CPU 内存和 GPU 内存是不同的。没有 CUDA-aware MPI，开发人员只能将指向 CPU/主机内存的指针传递给 MPI 调用。以下代码是使用非 CUDA-aware MPI 调用的示例：

 //MPI rank 0:Passing s_buf residing in GPU memory 
 // requires it to be transferred to CPU memory
cudaMemcpy(s_buf_h,s_buf_d,size,cudaMemcpyDeviceToHost);
MPI_Send(s_buf_h,size,MPI_CHAR,1,100,MPI_COMM_WORLD);

//MPI rank 1: r_buf received buffer needs to be 
// transferred to GPU memory before being used in GPU
MPI_Recv(r_buf_h,size,MPI_CHAR,0,100,MPI_COMM_WORLD, &status);
cudaMemcpy(r_buf_d,r_buf_h,size,cudaMemcpyHostToDevice);

有了 CUDA-aware MPI 库，这是不必要的；GPU 缓冲区可以直接传递给 MPI，如下所示：

//MPI rank 0
MPI_Send(s_buf_d,size,MPI_CHAR,1,100,MPI_COMM_WORLD);

//MPI rank n-1
MPI_Recv(r_buf_d,size,MPI_CHAR,0,100,MPI_COMM_WORLD, &status);

例如，对于 Open MPI，CUDA-aware 支持存在于 Open MPI 1.7 系列及更高版本中。要启用此功能，需要在编译时配置 Open MPI 库以支持 CUDA，如下所示：

$ ./configure --with-cuda

拥有 CUDA-aware MPI 并不意味着总是使用 GPUDirect RDMA。如果数据传输发生在网络卡和 GPU 之间共享相同的根复杂，则使用 GPUDirect 功能。尽管如此，即使未启用 RDMA 支持，拥有 CUDA-aware MPI 也可以通过利用诸如消息传输之类的功能使应用程序更有效，如下图所示可以进行流水线处理：

！

上图显示了具有 GPUDirect 的 CUDA-aware MPI 与不具有 GPUDirect 的 CUDA-aware MPI。两个调用都来自 CUDA-aware MPI，但左侧是 GPUDirect 传输，右侧是没有 GPUDirect 传输。

非 GPUDirect 传输有以下阶段：

  • 节点 1：从 GPU1 传输到主机内存

  • 节点 1：从主机内存传输到网络适配器暂存区

  • 网络：通过网络传输

  • 节点 2：从网络暂存区传输到主机内存

  • 节点 2：从主机内存传输到 GPU 内存

如果支持 GPUDirect RDMA，则从 GPU 传输直接通过网络进行，涉及主机内存的额外副本都被删除。

现在我们已经掌握了这个概念，让我们开始将代码转换为使用 CUDA-aware MPI 编程启用多 GPU 支持。

多节点-多 GPU 高斯消元

准备您的 GPU 应用程序。此代码可以在本书的 GitHub 存储库中的06_multigpu/gaussian中找到。使用nvcc编译器编译和运行应用程序，如下所示：

$ mpicc-o gaussian_multi_gpu_rdma.out gaussian_multi_gpu_rdma.cu
$ mpirun -np 8 ./gaussian_multi_gpu_rdma.out

我们使用mpicc而不是nvcc来编译 MPI 程序。我们使用mpirun命令运行可执行文件，而不是直接运行已编译的可执行文件。本节中您将看到的结果是在同一系统上具有 8 个 V100 的 DGX 系统上运行的输出。我们利用 8 个最大 MPI 进程，将每个 GPU 映射为 1 个 MPI 进程。要了解如何将多个 MPI 进程映射到同一 GPU，请阅读本章后面的MPS子节。在本练习中，我们使用了已编译为支持 CUDA 的 Open MPI 1.10，如前一节所述。

多 GPU 实现涉及的步骤如下：

  1. MPI 进程的秩 0 生成线性系统（矩阵 A，B）的数据。

  2. 转置增广矩阵（AB^T）由根节点在 MPI 进程之间使用MPI_Scatterv()按行分割。

  3. 每个 MPI 进程并行计算其部分输入：

  • 三个内核的处理发生在 GPU 上。

  • 在findPivot操作后，通过MPI_Send()/Recv()实现了枢轴的共识。

  1. 减少的转置增广矩阵（ABT）使用MPI_Gatherv()在根节点上收集。

  2. 根节点执行回代以计算解 X。

展示前面代码的提取样本高斯代码如下：

void gaussianEliminationOnGPU() {
    cudaSetDevice(nodeLocalRank); //Set CUDA Device based on local rank
    //Copy  chuck of AB Transpose from Host to GPU; 
   for( int n = 0; n < N; n++ ){ 
       prank = MPI rank that holds n; 
       if (myRank == prank) 
           findPivotRowAndMultipliers<<<...>>>(); 
       bCastPivotInfo(); // from prank to other ranks 
       extractPivotRow<<<...>>>(); 
       rowElimination<<<...>>>(); 
   //Copy  myPartOfReducedTransposeAB from GPU to Host;
}

现在，让我们添加多 GPU 支持：

  1. 设置每个 MPI 等级的 CUDA 设备：在 Open MPI 中，您可以通过使用MPI_COMM_TYPE_SHARED作为MPI_Comm_split_type的参数来获得 MPI 进程的本地等级，如下面的代码所示：
MPI_Comm loc_comm;
MPI_Comm_split_type(MPI_COMM_WORLD, MPI_COMM_TYPE_SHARED, rank, MPI_INFO_NULL, &loc_comm);
int local_rank = -1;
MPI_Comm_rank(loc_comm,&local_rank);
MPI_Comm_free(&loc_comm);

现在我们有了本地等级，每个 MPI 进程都使用它来通过cudaSetDevice()设置当前 GPU，如下图所示：

  1. 使用MPI_Scatter将输入拆分并分发到不同的 MPI 进程：
void distributeInputs() {
    MPI_Scatterv(transposeAB, ..., myPartOfTransposeAB, recvCount, MPI_UNSIGNED, 0, MPI_COMM_WORLD); 
}
  1. 在 GPU 上执行高斯消元：
void gaussianEliminationOnGPU() { 
    cudaSetDevice(nodeLocalRank);
     for( int n = 0; n < N; n++ ){ 
        prank = MPI rank that holds n; 
        if (myRank == prank) 
            findPivotRowAndMultipliers<<<...>>>();
        MPI_Bcast(...); // from prank to other ranks 
        extractPivotRow<<<...>>>(); 
        rowElimination<<<...>>>(); 
}

在执行任何操作之前，基于本地等级设置当前 GPU。然后，由负责该行的进程提取枢轴行，然后将枢轴行广播到所有其他 MPI 等级，我们用于消除。

通过使用异步 MPI 调用而不是使用广播 API（如MPI_Bcast），可以提高传输时间的整体性能。实际上，不建议使用广播 API；它应该被替换为可以实现相同功能的MPI_Isend和MPI_Irecv，这些是异步版本。请注意，使调用异步会增加其他方面（如调试）的复杂性。因此，用户需要编写额外的代码来发送和接收数据。

本章提供了在向现有 MPI 程序添加 GPU 支持时的最佳编码实践，并不应被视为 MPI 编程的最佳编程实践的专家指南。

CUDA 流

流以 FIFO 方式工作，其中操作的顺序按照它们发出的顺序执行。从主机代码发出的请求被放入先进先出队列中。队列由驱动程序异步读取和处理，并且设备驱动程序确保队列中的命令按顺序处理。例如，内存复制在内核启动之前结束，依此类推。

使用多个流的一般想法是，在不同流中触发的 CUDA 操作可能会并发运行。这可能导致多个内核重叠或内核执行中的内存复制重叠。

为了理解 CUDA 流，我们将看两个应用程序。第一个应用程序是一个简单的矢量加法代码，添加了流，以便它可以重叠数据传输和内核执行。第二个应用程序是一个图像合并应用程序，也将在第九章中使用，使用 OpenACC 进行 GPU 编程。

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

  1. 准备您的 GPU 应用程序。例如，我们将合并两个图像。此代码可以在本书的 GitHub 存储库的06_multi-gpu/streams文件夹中找到。

  2. 使用nvcc编译器编译您的应用程序如下：

$ nvcc --default-stream per-thread -o vector_addition -Xcompiler -fopenmp -lgomp vector_addition.cu
$ nvcc --default-stream per-thread -o merging_muli_gpu -Xcompiler -fopenmp -lgomp scrImagePgmPpmPackage.cu image_merging.cu
$ ./vector addition
$ ./merging_muli_gpu

上述命令将创建两个名为vector_addition和merging_multi_gpu的二进制文件。正如您可能已经注意到的，我们在我们的代码中使用了额外的参数。让我们更详细地了解它们：

  • --default-stream per-thread：此标志告诉编译器解析代码中提供的 OpenACC 指令。

  • -Xcompiler -fopenmp -lgomp：此标志告诉nvcc将这些附加标志传递给 CPU 编译器，以编译代码的 CPU 部分。在这种情况下，我们要求编译器向我们的应用程序添加与 OpenMP 相关的库。

我们将把这一部分分为两部分。应用程序 1 和应用程序 2 分别演示了在单个和多个 GPU 中使用流。

应用程序 1-使用多个流来重叠数据传输和内核执行

我们需要遵循的步骤来重叠数据传输和内核执行，或者同时启动多个内核如下：

  1. 声明要固定的主机内存，如下面的代码片段所示：
cudaMallocHost(&hostInput1, inputLength*sizeof(float));
cudaMallocHost(&hostInput2, inputLength*sizeof(float));
cudaMallocHost(&hostOutput, inputLength*sizeof(float));

在这里，我们使用cudaMallocHost() API 来分配固定内存的向量。

  1. 创建一个Stream对象，如下面的代码片段所示：
for (i = 0; i < 4; i++) {
 cudaStreamCreateWithFlags(&stream[i],cudaStreamNonBlocking);

在这里，我们使用cudaStreamCreateWithFlags() API，传递cudaStreamNonBlocking作为标志，使此流非阻塞。

  1. 调用 CUDA 内核和内存复制时使用stream标志，如下面的代码片段所示：
for (i = 0; i < inputLength; i += Seglen * 4) {
    for (k = 0; k < 4; k++) {
        cudaMemcpyAsync(... , cudaMemcpyHostToDevice, stream[k]);
        cudaMemcpyAsync(... , cudaMemcpyHostToDevice, stream[k]);
        vecAdd<<<Gridlen, 256, 0, stream[k]>>>(...);
    }
}

如我们所见，我们不是通过一次复制整个数组来执行矢量加法，而是将数组分成段，并异步复制这些段。内核执行也是在各自的流中异步进行的。

当我们通过 Visual Profiler 运行这段代码时，我们可以看到以下特点：

前面的分析器截图显示，蓝色条（基本上是vector_addition内核）重叠了内存复制。由于我们在代码中创建了四个流，分析器中也有四个流。

每个 GPU 都有两个内存复制引擎。一个负责主机到设备的传输，另一个负责设备到主机的传输。因此，发生在相反方向的两个内存复制可以重叠。此外，内存复制可以与计算内核重叠。这可以导致n路并发，如下图所示：

每个 GPU 架构都有一定的约束和规则，根据这些规则，我们将在执行时看到这些重叠。一般来说，以下是一些指导方针：

  • CUDA 操作必须在不同的非 0 流中。

  • 使用cudaMemcpyAsync时，主机应该使用cudaMallocHost()或cudaHostAlloc()进行固定。

  • 必须有足够的资源可用。

  • 不同方向的cudaMemcpyAsyncs

  • 设备资源（SMEM、寄存器、块等）以启动多个并发内核

应用程序 2 - 使用多个流在多个设备上运行内核

为了在多个设备上运行内核并重叠内存传输，我们之前遵循的步骤保持不变，除了一个额外的步骤：设置 CUDA 设备以创建流。让我们看看以下步骤：

  1. 创建与系统中 CUDA 设备数量相等的流，如下面的代码片段所示：
cudaGetDeviceCount(&noDevices);
cudaStream_t *streams;
streams = (cudaStream_t*) malloc(sizeof(cudaStream_t) * noDevices);

我们使用cudaGetDeviceCount() API 来获取 CUDA 设备的数量。

  1. 在各自的设备中创建流，如下面的代码片段所示：
#pragma omp parallel num_threads(noDevices)
{
     int block = omp_get_thread_num();
    cudaSetDevice(block);
    cudaStreamCreate(&streams[block]);

我们启动与 CUDA 设备数量相等的 OpenMP 线程，以便每个 CPU 线程可以为其各自的设备创建自己的 CUDA 流。每个 CPU 线程执行cudaSetDevice()来根据其 ID 设置当前 GPU，然后为该设备创建流。

  1. 在该流中启动内核和内存复制，如下所示：
cudaMemcpyAsync(... cudaMemcpyHostToDevice,streams[block]);
cudaMemcpyAsync(..., cudaMemcpyHostToDevice, streams[block]);
merging_kernel<<<gridDim,blockDim,0,streams[block]>>>(...);
cudaMemcpyAsync(...,streams[block]);

在分析器中运行代码后的输出可以在下面的截图中看到，这代表了 Visual Profiler 的时间轴视图。这显示了一个 GPU 的内存复制与另一个 GPU 的内核执行重叠：

如您所见，我们在拥有四个 V100 的多 GPU 系统上运行了这段代码。不同 GPU 中的内存复制和内核重叠。在这段代码中，我们演示了利用 OpenMP 在不同设备上并行调用 CUDA 内核。这也可以通过利用 MPI 来启动利用不同 GPU 的多个进程来实现。

在下一节中，我们将看一些额外的主题，这些主题可以提高多 GPU 应用程序的性能，并帮助开发人员分析和调试他们的代码。

额外的技巧

在本节中，我们将涵盖一些额外的主题，这些主题将帮助我们了解多 GPU 系统的额外特性。

使用 InfiniBand 网络卡对现有系统进行基准测试

有不同的基准可用于测试 RDMA 功能。InfiniBand 适配器的一个这样的基准可以在www.openfabrics.org/找到。您可以通过执行以下代码来测试您的带宽：

$ git clone git://git.openfabrics.org/~grockah/perftest.git
$ cd perftest 
$ ./autogen.sh 
$ export CUDA_H_PATH=<<Path to cuda.h>> 
$ ./configure –prefix=$HOME/test 
$ make all install

然后，您可以运行以下命令来测试带宽：

For example host to GPU memory (H-G) BW test:
server$ ~/test/bin/ib_write_bw -n 1000 -O -a --use_cuda
client $ ~/test/bin/ib_write_bw -n 1000 -O -a server.name.org

//GPU to GPU memory (G-G) BW test:
server$ ~/test/bin/ib_write_bw -n 1000 -O -a --use_cuda
client $ ~/test/bin/ib_write_bw -n 1000 -O -a --use_cuda server.name.org

NVIDIA 集体通信库（NCCL）

NCCL 提供了常用于深度学习等领域的通信原语的实现。NCCL 1.0 从同一节点内多个 GPU 之间的通信原语实现开始，并发展到支持多个节点上的多个 GPU。NCCL 库的一些关键特性包括以下内容：

  • 支持来自多个线程和多个进程的调用

  • 支持多个环和树拓扑，以更好地利用节点内和节点间的总线

  • 支持 InfiniBand 节点间通信

  • 源代码包可以从 GitHub（github.com/nvidia/nccl）免费下载

NCCL 可以扩展到 24,000 个 GPU，延迟低于 300 微秒。请注意，尽管 NCCL 已被证明是深度学习框架中非常有用和方便的库，但在用于 HPC 应用时存在局限，因为它不支持点对点通信。NCCL 支持集体操作，这在深度学习应用中被使用，例如以下内容：

  • AllReduce

  • AllGather

  • ReduceScatter

  • Reduce

  • Broadcast

所有 NCCL 调用都作为 CUDA 内核运行，以更快地访问 GPU 内存。它使用较少的线程，实现为一个块。这最终只在一个 GPU SM 上运行，因此不会影响其他 GPU 的利用率。让我们看一下以下代码：

ncclGroupStart(); 
for (int i=0; i<ngpus; i++) 
{ 
    ncclAllGather(…, comms[i], streams[i]); 
} 
ncclGroupEnd();

正如我们所看到的，NCCL 调用简单，易于调用。

使用 NCCL 加速集体通信

**NVIDIA 集体通信库（NCCL）**提供了为多个 NVIDIA GPU 优化的性能集体通信原语。在本节中，我们将看到这个库是如何工作的，以及我们如何从中受益。

并不难找到使用多个 GPU 来训练网络的深度学习模型。由于两个 GPU 并行计算神经网络，我们很容易想象这种技术将随着 GPU 数量的增加而提高训练性能。不幸的是，世界并不那么简单。梯度应该在多个 GPU 之间共享，并且一个 GPU 中的权重更新过程应该等待其他 GPU 的梯度来更新其权重。这是使用多个 GPU 进行深度学习训练的一般过程，并在以下图表中显示：

集体通信有许多类型：全局归约、广播、归约、全局收集、归约散射等。在深度学习中，每个 GPU 在传输自己的数据的同时收集另一个 GPU 的数据。因此，我们可以确定深度学习在通信中需要所有类型的归约样式通信。

在 HPC 社区中，包括全局归约在内的集体通信是一个常见的话题。节点内和节点间处理器之间的通信是一个具有挑战性但至关重要的问题，因为它直接关系到可扩展性。正如我们在第六章中提到的，可扩展的多 GPU 编程，在多 GPU 编程部分，需要仔细考虑与每个 GPU 的通信。开发人员应该设计和实现 GPU 中的集体通信，即使 MPI 已经支持这样的通信模式。

NCCL 提供了一种集体通信，它了解 GPU 拓扑配置。通过使用各种分组和通信命令，您可以应用所需的通信任务。

一个前提是您的系统需要有多个 GPU，因为 NCCL 是一个与多个 GPU 一起工作的通信库。

以下步骤涵盖了如何调用ncclAllReduce()来测试和测量系统的 GPU 网络带宽。示例代码实现在04_nccl中：

  1. 让我们定义一个类型，它将包含、发送和接收每个 GPU 设备的缓冲区和cudaStream，如下所示：
typedef struct device
{
    float *d_send;
    float *d_recv;
    cudaStream_t stream;
} device_t;
  1. 在应用程序开始时，我们需要准备一些句柄，以便我们可以控制多个 GPU：
cudaGetDeviceCount(&num_dev);
ncclComm_t *ls_comms = new ncclComm_t[num_dev];
int *dev_ids = new int[num_dev];
for (int i = 0; i < num_dev; i++)
    dev_ids[i] = i;
  1. 然后，我们将创建一个缓冲区，假设我们有数据。对于每个设备，我们将初始化每个设备的项目，如下所示：
unsigned long long size = 512 * 1024 * 1024; // 2 GB

// allocate device buffers and initialize device handles
device_t *ls_dev = new device_t[num_dev];
for (int i = 0; i < num_dev; i++) {
    cudaSetDevice(i);
    cudaMalloc((void**)&ls_dev[i].d_send, sizeof(float) * size);
    cudaMalloc((void**)&ls_dev[i].d_recv, sizeof(float) * size);
    cudaMemset(ls_dev[i].d_send, 0, sizeof(float) * size);
    cudaMemset(ls_dev[i].d_recv, 0, sizeof(float) * size);
    cudaStreamCreate(&ls_dev[i].stream);
}
  1. 在开始 NCCL 通信之前，我们需要初始化 GPU 设备，以便它们知道它们在 GPU 组中的排名。由于我们将用单个进程测试带宽，我们可以安全地调用一个初始化所有设备的函数：
ncclCommInitAll(ls_comms, num_dev, dev_ids);

  1. 如果我们要用多个进程测试带宽，我们需要调用ncclCommInitRank()。我们需要为计算进程 ID 和 GPU 排名提供 GPU ID。

  2. 现在，我们可以使用 NCCL 完成 all-reduce 操作。以下代码是ncclAllReduce的示例实现：

ncclGroupStart();
for (int i = 0; i < num_dev; i++) {
    ncclAllReduce((const void*)ls_dev[i].d_send, 
                  (void*)ls_dev[i].d_recv,
        test_size, ncclFloat, ncclSum, 
        ls_comms[i], ls_dev[i].stream);
}
ncclGroupEnd();

对于每个设备，我们需要触发流量。为此，我们需要启动和关闭 NCCL 组通信。现在，我们已经实现了一些使用ncclAllReduce()的测试代码。让我们通过微基准测试来了解 NCCL 的工作原理。

在多 GPU 系统上测试此代码，运行以下命令：

$ nvcc -run -m64 -std=c++11 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -lnccl -o nccl ./nccl.cu

以下图表显示了在 DGX Station 中使用四个 V100 32G GPU 测得的性能。蓝线表示基于 NVLink 的带宽，而橙线表示基于 PCIe 的带宽，通过设置NCCL_P2P_DISABLE=1 ./ncd并关闭对等 GPU 来实现：

这个 NCCL 测试可能会受到系统配置的影响。这意味着结果可能会有所不同，取决于您系统的 GPU 拓扑结构。

这显示了基于 PCI Express 和基于 NVLINK 的 all-reduce 性能差异。我们可以使用nvprof来查看通信。以下屏幕截图显示了通过 NCCL 2.3.7 在 DGX Station 上的 all-reduce 通信：

NCCL 越来越快。通过引入新的 GPU 互连技术 NVLink 和 NVSwitch，我们对 NCCL 的经验正在增加，以至于我们可以实现可扩展的性能。

以下链接提供了关于 NCCL 的讨论：developer.nvidia.com/gtc/2019/video/S9656/video。

摘要

在本章中，我们介绍了多 GPU 编程的不同方法。通过示例高斯消元，我们看到了如何将单个 GPU 应用程序工作负载分割到多个 GPU 中，首先是单个节点，然后是多个节点。我们看到了系统拓扑在利用 P2P 传输和 GPUDirect RDMA 等功能方面起着重要作用。我们还看到了如何使用多个 CUDA 流来重叠多个 GPU 之间的通信和数据传输。我们还简要介绍了一些其他主题，可以帮助 CUDA 程序员优化代码，如 MPS 和使用nvprof来分析多 GPU 应用程序。

在下一章中，我们将看到大多数 HPC 应用程序中出现的常见模式以及如何在 GPU 中实现它们。

CUDA 中的并行编程模式

在本章中，我们将涵盖并行编程算法，这将帮助您了解如何并行化不同的算法并优化 CUDA。本章中我们将涵盖的技术可以应用于各种问题，例如我们在第三章中看到的并行减少问题，CUDA 线程编程，它可以用于设计神经网络操作中的高效 softmax 层。

在本章中，我们将涵盖以下主题：

  • 矩阵乘法优化

  • 图像卷积

  • 前缀和

  • 打包和拆分

  • N 体操作

  • 在 CUDA 中使用动态并行性进行快速排序

  • 基数排序

  • 直方图计算

技术要求

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

本章中的相同代码已经使用 CUDA 版本 10.1 进行开发和测试。一般来说，如果适用的话，建议使用最新的 CUDA 版本。

矩阵乘法优化

虽然我们在许多示例中使用了矩阵乘法代码，但我们并没有调查操作是否被优化。现在，让我们回顾其操作以及如何找到优化的机会。

矩阵乘法是从两个矩阵进行的一组点积运算。我们可以简单地并行化所有 CUDA 线程执行的操作，以生成元素的点积。然而，从内存使用的角度来看，这种操作效率低，因为从内存加载的数据没有被重复使用。为了确认我们的类比，让我们测量性能限制器。以下图表显示了使用 NVIDIA Nsight Compute 的 Tesla V100 卡的 GPU 利用率：

根据我们的性能限制器分析，这种利用率可以归类为内存受限。因此，我们应该审查内存利用率以减少利用率。以下截图显示了内存工作负载分析部分：

通过这个分析，我们可以看到 L2 缓存命中率低，最大带宽也低。我们可以推测这是因为原始矩阵乘法操作没有重复使用加载的数据，正如我们之前提到的。这可以通过使用共享内存来解决，即重复使用加载的数据并减少全局内存使用。现在，让我们回顾矩阵乘法以及如何优化使用具有小内存空间的共享内存。

矩阵乘法是一组点积运算，使用一些小尺寸矩阵和输出的累积。小矩阵称为瓦片，它们映射到输出矩阵上。每个瓦片将并行计算自己的输出。这个操作可以按以下步骤实现：

  1. 确定两个输入和输出矩阵的瓦片大小。

  2. 遍历输入瓦片，以及它们的方向（矩阵 A 向右移动，矩阵 B 向下移动）。

  3. 在瓦片内计算矩阵乘法。

  4. 继续第二步，直到瓦片达到末尾。

  5. 刷新输出。

以下图表显示了瓦片矩阵乘法的概念：

在上图中，我们计算矩阵乘法，C = AB。我们从矩阵 A 和矩阵 B 中计算一个较小的矩阵乘法作为瓦片（绿色）。然后，我们分别遍历输入瓦片位置。操作结果累积到先前的输出，以生成矩阵乘法的输出。

这个操作提供了一个优化机会，因为我们可以将大矩阵操作分解为小问题，并将其放置在小内存空间中。在 CUDA 编程中，我们将小矩阵放置在共享内存中，并减少全局内存访问。在我们的实现中，我们将瓦片与 CUDA 线程块匹配。瓦片的位置将由其块索引确定，这是通过tid_*变量完成的。

实现平铺方法

现在，让我们使用平铺方法实现优化的矩阵乘法。我们将重用之前在第三章中使用的矩阵乘法示例代码，即 CUDA 线程编程。优化后，我们将看看如何提高性能。按照以下步骤开始：

  1. 让我们创建一个核函数，这将是我们优化版本的矩阵乘法。我们将在sgemm操作中命名核函数为v2。这个核函数将计算，因此我们应该分别提供相关参数。我们还将使用M、N和K传递矩阵大小信息：
__global__ void sgemm_kernel_v2(const float *A, const float *B, float *C,
    int M, int N, int K, float alpha, float beta) {}
  1. 对于这个操作，我们将分别使用块索引和线程索引。正如我们之前讨论的，我们需要单独使用块索引来指定瓦片位置。我们将使用线程索引进行瓦片级矩阵乘法。因此，我们需要创建 CUDA 索引参数，如下所示：
int bid_x = blockIdx.x * blockDim.x;
int bid_y = blockIdx.y * blockDim.y;
int tid_x = threadIdx.x;
int tid_y = threadIdx.y;
  1. 之后，我们将使用共享内存作为瓦片，并使用本地寄存器保存输出值：
float element_c = 0.f;
__shared__ float s_tile_A[BLOCK_DIM][BLOCK_DIM];
__shared__ float s_tile_B[BLOCK_DIM][BLOCK_DIM];
  1. 然后，我们将编写一个控制瓦片位置的循环。以下是基于其块大小控制循环的 for 循环代码。请注意，循环大小由K决定，考虑到块应该遍历多少次：
for (int k = 0; k < K; k += BLOCK_DIM)
{
   ... {step 5 and 6 will cover } ...
}
  1. 现在，我们将编写代码，将数据输入第二个循环。正如我们之前讨论的，每个瓦片都有自己的移动方向，以及矩阵；瓦片A遍历矩阵A的列，瓦片B遍历矩阵B的行。我们根据矩阵乘法优化部分中显示的图表来放置它们。之后，我们应该在从全局内存复制数据到共享内存后放置__syncthreads()，以避免来自上一次迭代的未更新数据：
// Get sub-matrix from A
s_tile_A[tid_y][tid_x] = A[ (bid_y + tid_y) * K + tid_x + k ];
// Get sub-matrix from B 
s_tile_B[tid_y][tid_x] = B[ k * N + bid_x + tid_x ]; 

__syncthreads();
  1. 然后，我们可以从瓦片中编写矩阵乘法代码。名为element_c的本地变量将累积结果：
for (int e = 0; e < BLOCK_DIM; e++)
    element_c += s_tile_A[tid_y][e] * s_tile_B[e][tid_x];
  1. 我们将结果写入全局内存。以下操作应该放置在第二个循环完成后：
C[(bid_y + tid_y) * N + (bid_x + tid_x)] = \
 alpha * element_c + beta * C[(bid_y + tid_y) * N + (bid_x + tid_x)];
  1. 现在，让我们回顾一下这种平铺方法如何有利于矩阵乘法操作。通过在我们的平铺矩阵乘法中使用共享内存，我们可以期望通过使用输入数据减少全局内存流量，从而增强性能。我们可以轻松地通过配置文件结果来确认这一点：
$ nvcc -m64 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -o sgemm ./sgemm.cu 
$ nvprof ./sgemm 

        Type Time(%)    Time Calls      Avg      Min      Max Name
GPU activities: 47.79% 9.9691ms     1 9.9691ms 9.9691ms 9.9691ms sgemm_kernel(...)
 32.52% 6.7845ms     1 6.7845ms 6.7845ms 6.7845ms sgemm_kernel_v2(...)
  1. 由于我们设计了核心以重用输入数据，增加的块大小可能有助于性能。例如，考虑到 warp 大小和共享内存银行的数量，32 x 32 的块大小可能是最佳的，以避免银行冲突。我们可以轻松地使用配置文件获得其实验结果：
 Type Time(%)    Time Calls      Avg       Min       Max Name
GPU activities: 46.52% 8.1985ms     1 8.1985ms  8.1985ms  8.1985ms sgemm_kernel(...)
 31.24% 5.4787ms     1 5.4787ms  5.4787ms  5.4787ms sgemm_kernel_v2(...)

正如你所看到的，增加的瓦片大小有利于矩阵乘法操作的性能。现在，让我们分析其性能。

平铺方法的性能分析

之前，我们看过了平铺方法以及它如何能够实现良好的性能。让我们回顾一下平铺方法解决了什么问题，并看看接下来我们可以采取哪些步骤。总的来说，覆盖这部分是可选的，因为 NVIDIA 提供了 cuBLAS 和 CUTLASS 库，用于提供优化性能的 GEMM（通用矩阵乘法）操作。

下图显示了来自 NVIDIA Nsight Compute 的更新后的 GPU 利用率报告。较低配置文件的更新利用率输出是上配置文件的结果：

由于两个资源都利用率很高，我们应该审查每个资源的资源使用情况。首先，让我们来审查内存工作量。以下截图显示了更新后的结果：

从这个结果可以看出，全局内存访问从最大化内存带宽和减少内存吞吐量进行了优化。此外，L2 缓存命中率也得到了提高。因此，我们的平铺方法将矩阵乘法从全局内存转换为芯片级操作。

然而，这并不意味着我们已经实现了最优化的性能。从内存工作量分析中，我们可以看到内存管道太忙了。这是由于我们从共享内存进行逐元素乘法。为了解决这个问题，我们需要重新映射共享内存中的数据。我们不会在这本书中涉及到这个问题，但你可以在这篇文章中了解到：github.com/NervanaSystems/maxas/wiki/SGEMM。

正如我们之前讨论的，cuBLAS 库显示出更快的性能。我们将在第八章的cuBLAS部分中介绍其用法。然而，在这个阶段理解平铺方法是有用的，这样我们就可以理解 GPU 如何开始优化。

卷积

卷积操作（或滤波）是许多应用中常见的操作，特别是在图像和信号处理以及深度学习中。虽然这个操作是基于输入和滤波器的顺序数据的乘积，但我们对矩阵乘法有不同的方法。

CUDA 中的卷积操作

卷积操作包括源数据和滤波器。滤波器也被称为核。通过将滤波器应用于输入数据，我们可以获得修改后的结果。下图显示了二维卷积的示意图：

当我们实现卷积操作时，我们需要考虑一些概念，即核和填充。核是一组我们想要应用到源数据的系数。这也被称为滤波器。填充是源数据周围的额外虚拟空间，以便我们可以将核函数应用到边缘。当填充大小为 0 时，我们不允许滤波器移动超出源空间。然而，一般来说，填充大小是滤波器大小的一半。

为了轻松开始，我们可以考虑以下几点来设计核函数：

  • 每个 CUDA 线程生成一个滤波输出。

  • 每个 CUDA 线程将滤波器的系数应用于数据。

  • 滤波器的形状是盒状滤波器。

在满足这些条件的情况下，我们可以有一个简单的卷积操作滤波器，如下所示：

__global__ void
convolution_kernel_v1(float *d_output, float *d_input, float *d_filter, int num_row, int num_col, int filter_size)
{
    int idx_x = blockDim.x * blockIdx.x + threadIdx.x;
    int idx_y = blockDim.y * blockIdx.y + threadIdx.y;

    float result = 0.f;
    // iterates over the every value in the filter
    for (int filter_row = -filter_size / 2; 
         filter_row <= filter_size / 2; ++filter_row)
    {
        for (int filter_col = -filter_size / 2; 
             filter_col <= filter_size / 2; ++filter_col)
        {
            // Find the global position to apply the given filter
            // clamp to boundary of the source
            int image_row = min(max(idx_y + filter_row, 0), 
                                static_cast<int>(num_row - 1));
            int image_col = min(max(idx_x + filter_col, 0), 
                                static_cast<int>(num_col - 1));

            float image_value = static_cast<float>(
                                d_input[image_row * num_col + 
                                image_col]);
            float filter_value = d_filter[(filter_row + 
                                           filter_size / 2) * 
                                           filter_size 
                                           + filter_col + 
                                           filter_size / 2];

            result += image_value * filter_value;
        }
    }

    d_output[idx_y * num_col + idx_x] = result;
}

这个核函数获取输入数据和滤波器进行操作，并没有重用所有数据。考虑到内存效率带来的性能影响，我们需要设计我们的核心代码，以便可以重用加载的数据。现在，让我们编写卷积的优化版本。

优化策略

首先，卷积滤波器是一个只读矩阵，并且被所有 CUDA 线程使用。在这种情况下，我们可以使用 CUDA 的常量内存来利用其缓存操作和广播操作。

在卷积实现设计中，我们使用平铺方法，每个平铺将生成映射位置的滤波输出。我们的平铺设计有额外的空间来考虑卷积滤波器的大小，这为卷积操作提供了所需的数据。这个额外的空间被称为填充。下图显示了一个具有 6 x 6 维度和 3 x 3 大小滤波器的线程块的示例。

然后，我们需要为每个线程块在共享内存上有一个 8 x 8 大小的平铺，如下所示：

当源地址无效内存空间时，或者填充为零（零填充方法）时，填充区域可以是输入数据。通过这样做，我们可以使瓷砖替换输入全局内存而不会对边界元素产生额外影响。为了填充瓷砖，我们使用线程块大小迭代瓷砖，并通过检查输入数据的边界条件来确定应该填充哪个值。我们的实现将输入数据设置为瓷砖大小的倍数，以便边界条件与每个线程块的瓷砖的填充空间匹配。将源数据映射到瓷砖的简要图示如下：

在这个设计中，我们需要做的迭代次数来填充瓷砖是四次。然而，这应该根据滤波器大小进行更改。这样，填充瓷砖的迭代次数由瓷砖大小的上限除以线程块大小确定。其实现很简单，如下面的代码所示：

for (int row = 0; row <= tile_size / BLOCK_DIM; row++) {
    for (int col = 0; col <= tile_size / BLOCK_DIM; col++) {
        ... (filter update operation) ...
    }
}

现在，让我们使用共享内存作为盒式滤波器来实现优化的卷积操作。

使用常量内存优化滤波系数

首先，我们将学习如何优化滤波系数数据的使用。

我们将制作convolution_kernel()的修改版本。让我们复制内核代码，并将其中一个重命名为convolution_kernel_v2()：

  1. 首先，我们将创建一个常量内存空间来存储滤波系数。常量内存的大小是有限的，我们不能对内核代码进行修改。然而，我们可以使用这个常量内存，因为我们的卷积滤波器适合这种条件。我们可以这样使用常量内存：
#define MAX_FILTER_LENGTH 128
__constant__ float c_filter[MAX_FILTER_LENGTH * MAX_FILTER_LENGTH];
  1. 然后，我们可以使用cudaMemcpyToSymbol()函数将卷积滤波系数放置在常量内存中：
cudaMemcpyToSymbol(c_filter, h_filter, filter_size * filter_size * sizeof(float));
  1. 让我们切换滤波操作，这样我们就可以使用常量内存。整个内核实现如下。正如你所看到的，只有一个变量的使用发生了变化：
__global__ void
convolution_kernel_v2(float *d_output, float *d_input, float *d_filter, int num_row, int num_col, int filter_size)
{
    int idx_x = blockDim.x * blockIdx.x + threadIdx.x;
    int idx_y = blockDim.y * blockIdx.y + threadIdx.y;

    float result = 0.f;
    for (int filter_row = -filter_size / 2; 
         filter_row <= filter_size / 2; ++filter_row)
    {
        for (int filter_col = -filter_size / 2; 
             filter_col <= filter_size / 2; ++filter_col)
        {
            int image_row = idx_y + filter_row;
            int image_col = idx_x + filter_col;

            float image_value = (image_row >= 0 
                                 && image_row < num_row 
                                 && image_col >= 0
                                 && image_col < num_col) ?
                                 d_input[image_row * num_col 
                                         + image_col] : 0.f;
            float filter_value = c_filter[(filter_row 
                                          + filter_size / 2) 
                                          * filter_size 
                                          + filter_col 
                                          + filter_size / 2];

            result += image_value * filter_value;
        }
    }

    d_output[idx_y * num_col + idx_x] = result;
}
  1. 现在，我们可以通过重复使用nvprof来确认性能提升：
$ nvcc -run -m64 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -o convolution ./convolution.cu
$ nvprof ./convolution
           Type Time(%) Time Calls Avg Min Max Name
 12.85% 442.21us     1 442.21us 442.21us 442.21us convolution_kernel_v1(...)
 11.97% 412.00us     1 412.00us 412.00us 412.00us convolution_kernel_v2(...)

从这个结果中，我们可以看到减少的内核执行时间。

使用共享内存平铺输入数据

现在，我们将使用共享内存来优化输入数据的使用。为了区分我们的下一个优化步骤，让我们复制之前的卷积核函数并将其命名为convolution_kernel_v3()：

  1. 首先，我们需要预先准备共享内存空间，以便它可以存储输入数据。为了从共享内存中获得滤波操作的好处，我们需要额外的输入数据。为了创建足够的内存空间，我们需要修改内核调用，如下所示：
int shared_mem_size = (2*filter_size+BLOCK_DIM) * (2*filter_size+BLOCK_DIM) * sizeof(float);
convolution_kernel_v3<<<dimGrid, dimBlock, shared_mem_size, 0 >>>(d_output, d_input, d_filter, num_row, num_col, filter_size);
  1. 在内核代码中，我们可以声明共享内存空间如下：
extern __shared__ float s_input[];
  1. 然后，我们可以将输入数据复制到由线程块计算的共享内存中。首先，让我们声明一些帮助控制内存操作的变量：
int pad_size = filter_size / 2;
int tile_size = BLOCK_DIM + 2 * pad_size;
  1. 现在，我们可以按照之前讨论的平铺设计将加载的输入数据复制到共享内存中：
for (int row = 0; row <= tile_size / BLOCK_DIM; row++) {
    for (int col = 0; col <= tile_size / BLOCK_DIM; col++) {
        int idx_row = idx_y + BLOCK_DIM * row - pad_size; 
        // input data index row
        int idx_col = idx_x + BLOCK_DIM * col - pad_size; 
        // input data index column
        int fid_row = threadIdx.y + BLOCK_DIM * row; 
        // filter index row
        int fid_col = threadIdx.x + BLOCK_DIM * col; 
        // filter index column

        if (fid_row >= tile_size || fid_col >= tile_size) continue;

        s_input[tile_size * fid_row + fid_col] = \
            (idx_row >= 0 && idx_row < num_row && idx_col >= 0 
                && idx_col < num_col) ? 
                d_input[num_col * idx_row + idx_col] : 0.f;
    }
}

__syncthreads();
  1. 由于输入内存已更改，我们的卷积代码应该更新。我们可以将卷积代码编写如下：
float result = 0.f;
    for (int filter_row = -filter_size / 2; 
         filter_row <= filter_size / 2; ++filter_row)
    {
        for (int filter_col = -filter_size / 2; 
             filter_col <= filter_size / 2; ++filter_col)
        {
            // Find the global position to apply the given filter 
            int image_row = threadIdx.y + pad_size + filter_row;
            int image_col = threadIdx.x + pad_size + filter_col;

            float image_value = s_input[tile_size 
                                        * image_row + image_col]; 
            float filter_value = c_filter[(filter_row 
                                          + filter_size / 2) 
                                          * filter_size 
                                          + filter_col 
                                          + filter_size / 2];

            result += image_value * filter_value;
        }
    }
  1. 最后，我们可以使用nvprof来测量性能增益。从结果中，我们可以确认我们的加速速度比原始操作快了大约 35%：
$ nvcc -run -m64 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -o convolution ./convolution.cu
$ nvprof ./convolution
Processing Time (1) -> GPU: 0.48 ms
Processing Time (2) -> GPU: 0.43 ms
Processing Time (3) -> GPU: 0.30 ms
Processing Time -> Host: 4104.51 ms
... (profiler output) ...
              type Time(%)    Time Calls .    Avg      Min .    Max Name
   GPU activities: 66.85% 2.3007ms     3 766.91us 1.1840us 2.2979ms [CUDA memcpy HtoD]
                   12.85% 442.21us     1 442.21us 442.21us 442.21us convolution_kernel_v1()
                   11.97% 412.00us     1 412.00us 412.00us 412.00us convolution_kernel_v2()
                    8.33% 286.56us     1 286.56us 286.56us 286.56us convolution_kernel_v3()

现在，我们已经了解如何利用加载的数据，以便我们可以重复使用它与其他片上缓存而不是全局内存。我们将在下一节中更详细地讨论这个问题。

获得更多性能

如果滤波器是对称滤波器或可分离滤波器，我们可以将盒式滤波器分解为两个滤波器：水平滤波器和垂直滤波器。使用两个方向滤波器，我们可以在共享内存使用方面进行更多优化：内存空间和内存利用率。如果您想了解更多信息，请查看名为convolutionSeparable的 CUDA 示例，该示例位于3_Imaging/convolutionSeparable目录中。其详细说明也包含在相同目录的doc/convolutionSeparable.pdf中。

前缀和（扫描）

前缀和（扫描）用于从给定的输入数字数组中获得累积数字数组。例如，我们可以按以下方式制作前缀和序列：

输入数字123456...
前缀和136101521...

它与并行减少不同，因为减少只是从给定的输入数据生成总操作输出。另一方面，扫描从每个操作生成输出。解决这个问题的最简单方法是迭代所有输入以生成输出。但是，在 GPU 中这将花费很长时间并且效率低下。因此，温和的方法可以并行化前缀和操作，如下所示：

在这种方法中，我们可以使用多个 CUDA 核心来获得输出。但是，这种方法并不会减少迭代的总次数，因为第一个输入元素应该逐个添加到所有输出中。此外，当数组足够大时，我们无法预测输出结果，因此应该启动多个线程块。这是因为在 CUDA 架构中，并非所有计划的 CUDA 线程都同时启动，并且多个 CUDA 线程会发生冲突。为了避免这种情况，我们需要对数组采用双缓冲区方法，这是另一种低效的方法。以下代码显示了它的实现：

__global__ void
scan_v1_kernel(float *d_output, float *d_input, int length, int offset) {
    int idx = blockDim.x * blockIdx.x + threadIdx.x;

    float element = 0.f;
    for (int offset = 0; offset < length; offset++) {
        if (idx - offset >= 0)
            element += d_input[idx - offset];
    }
    d_output[idx] = element;
}

还有另一种优化方法叫做Blelloch 扫描。该方法通过指数增加和减少步长来生成前缀和输出。该方法的过程如下图所示：

基于步长控制有两个步骤。在增加步长的同时，相应地获得部分总和。然后，在减小步长的同时获得部分总和。每个步骤都有不同的操作模式，但可以根据步长大小来确定。现在，让我们来看一下 Blelloch 扫描的实现并检查更新后的性能。

Blelloch 扫描实现

以下步骤将向您展示如何实现优化的并行扫描算法：

  1. 让我们创建一个可以接受输入和输出内存以及它们的大小的内核函数：
__global__ void scan_v2_kernel(float *d_output, float *d_input, int length)
{
    ...
}
  1. 然后，我们将创建一个 CUDA 线程索引和一个全局索引来处理输入数据：
int idx = blockDim.x * blockIdx.x + threadIdx.x;
int tid = threadIdx.x;
  1. 为了加快迭代速度，我们将使用共享内存。该算法可以生成 CUDA 线程大小的两倍输出，因此我们将额外加载块大小的输入数据到共享内存中：
extern __shared__ float s_buffer[];
s_buffer[threadIdx.x] = d_input[idx];
s_buffer[threadIdx.x + BLOCK_DIM] = d_input[idx + BLOCK_DIM];
  1. 在开始迭代之前，我们将声明偏移变量，该变量计算左操作数和右操作数之间的差距：
int offset = 1;
  1. 然后，我们将添加输入数据，直到偏移量大于输入的长度为止：
while (offset < length)
{
    __syncthreads();
    int idx_a = offset * (2 * tid + 1) - 1;
    int idx_b = offset * (2 * tid + 2) - 1;
    if (idx_a >= 0 && idx_b < 2 * BLOCK_DIM) {
        s_buffer[idx_b] += s_buffer[idx_a];
    }
    offset <<= 1;
}
  1. 之后，我们将通过减小减少大小来再次迭代两次：
offset >>= 1;
while (offset > 0) {
    __syncthreads();
    int idx_a = offset * (2 * tid + 2) - 1;
    int idx_b = offset * (2 * tid + 3) - 1;
    if (idx_a >= 0 && idx_b < 2 * BLOCK_DIM) {
        s_buffer[idx_b] += s_buffer[idx_a];
    }
    offset >>= 1;
}
__syncthreads();
  1. 最后，我们将使用内核函数将输出值存储在全局内存中：
d_output[idx] = s_buffer[tid];
d_output[idx + BLOCK_DIM] = s_buffer[tid + BLOCK_DIM];
  1. 现在，我们可以调用这个扫描内核函数如下：
void scan_v2(float *d_output, float *d_input, int length)
{
    dim3 dimBlock(BLOCK_DIM);
    dim3 dimGrid(1);
    scan_v2_kernel<<<dimGrid, dimBlock, 
                     sizeof(float) * BLOCK_DIM * 2>>>
                  (d_output, d_input, length);
    cudaDeviceSynchronize();
}

您还可以使用相同的函数接口编写一个朴素扫描版本。现在，让我们回顾一下我们的新版本有多快，以及我们是否可以利用其他优化机会。

  1. 以下代码显示了朴素扫描和 Blelloch 扫描性能的分析结果：
$ nvcc -m64 -std=c++11 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -L/usr/local/cuda/lib -o scan ./scan.cu ./scan_v1.cu ./scan_v2.cu
$ nvprof ./scan
            Type Time(%)     Time Calls      Avg      Min      Max Name
 GPU activities:  68.96% 22.751us     1 22.751us 22.751us 22.751us scan_v1_kernel(float*, float*, int)
 12.71% 4.1920us    1 4.1920us 4.1920us 4.1920us scan_v2_kernel(float*, float*, int)

正如你所看到的，由于减少了开销，Blolloch 扫描比朴素扫描算法快了大约五倍。我们还可以通过比较不同实现的输出来验证操作结果：

input         :: -0.4508 -0.0210 -0.4774  0.2750 ... 0.0398 0.4869
result[cpu]   :: -0.4508 -0.4718 -0.9492 -0.6742 ... 0.3091 0.7960
result[gpu_v1]:: -0.4508 -0.4718 -0.9492 -0.6742 ... 0.3091 0.7960
SUCCESS!!
result[cpu]   :: -0.4508 -0.4718 -0.9492 -0.6742 ... 0.3091 0.7960
result[gpu_v2]:: -0.4508 -0.4718 -0.9492 -0.6742 ... 0.3091 0.7960
SUCCESS!!

到目前为止，我们已经介绍了如何设计和实现优化的单个块大小的并行前缀和操作。要在输入数据上使用前缀和操作，需要比块大小更多的数据，我们需要基于我们的块级减少代码构建一个块级前缀和操作。我们将在下一节详细讨论这个问题。

构建全局大小扫描

我们实现的前缀和操作在单个线程块内工作。由于第一步有两个输入，而我们在一个线程块中最多可以有 1,024 个 CUDA 线程，因此最大可用大小为 2,048。在不考虑其他线程块操作的情况下，线程块进行上扫描和下扫描。

然而，如果我们执行一个分块扫描操作，这个操作可以被扩大。为了做到这一点，你需要额外的步骤来收集最后一个前缀和的结果，扫描它们，并将每个线程块的结果与每个块的块级扫描值相加。这个过程可以按照以下方式实现：

追求更好的性能

我们的实现代码执行了最佳操作。然而，我们可以通过减少共享内存的银行冲突来进一步优化。在我们的实现中，CUDA 线程在某些点上访问相同的内存银行。NVIDIA 的 GPU Gem3 在第三十九章，使用 CUDA 进行并行前缀和（扫描）（developer.nvidia.com/gpugems/GPUGems3/gpugems3_ch39.html）中介绍了前缀和（扫描），并在39.2.3 避免银行冲突中指出了这个问题。你可以将解决方案调整到我们的实现，但如果这样做，你应该将NUM_BANKS更新为32，LOG_NUM_BANKS更新为5。现在，CUDA 架构有 32 个共享内存银行。

并行前缀和操作的其他应用

G.E. Blelloch 博士在 1993 年发表了一篇关于他的前缀和算法的文章前缀和及其应用（www.cs.cmu.edu/~guyb/papers/Ble93.pdf）。通过阅读他的文章，你可以了解更多关于并行前缀和算法及其应用。这些应用包括压缩、分割、分段扫描、快速排序、基数排序和归并排序。

Ahmed Sallm 博士的视频讲座，使用 CUDA 进行并行处理简介-第 4 讲第 2\3 部分（youtu.be/y2HzWKTqo3E），对此提供了很好的介绍。它提供了关于前缀和算法如何用于裁剪图形和构建稀疏矩阵的概念介绍。他还提供了关于如何使用排序算法的说明。

压缩和分割

之前，我们介绍了如何并行化顺序前缀和算法，并讨论了它如何用于其他应用。现在，让我们来介绍其中一些应用：压缩和分割。压缩操作是一种可以从数组中整合满足给定条件的值的算法。另一方面，分割操作是一种将值分配到指定位置的算法。一般来说，这些算法是顺序工作的。然而，我们将看到并行前缀和操作如何改进它的功能。

压缩操作用于将满足特定条件的特定数据收集到一个数组中。例如，如果我们想要对数组中的正元素使用压缩操作，那么操作如下：

在并行编程中，我们有一种不同的方法，可以利用并行前缀和操作使用多个核心。首先，我们标记数据以检查它是否满足条件（即谓词），然后进行前缀和操作。前缀和的输出将是标记值的索引，因此我们可以通过复制它们来获得收集的数组。下面的图表显示了压缩操作的一个示例：

由于所有这些任务都可以并行完成，我们可以在四个步骤中获得收集的数组。

另一方面，拆分意味着将数据分发到多个不同的位置。一般来说，我们会从最初的位置分发数据。下面的图表显示了它的操作示例：

这个例子显示了收集的数组元素是如何分布在它们原来的位置的。我们也可以使用前缀和并行地做到这一点。首先，我们参考谓词数组并进行前缀和操作。由于输出是每个元素的地址，我们可以很容易地分配它们。下面的图表显示了如何进行这个操作：

`

现在，让我们实现这个并讨论它们的性能限制和应用。

实现压缩

压缩操作是一个谓词、扫描、寻址和收集的序列。在这个实现中，我们将从一个随机生成的数字数组中构建一个正数数组。初始版本只能承受单个线程块操作，因为我们只会使用一个块大小的前缀和操作。然而，我们可以了解前缀和如何对其他应用有用，并将这个操作扩展到更大的数组，使用扩展的前缀和操作。

为了实现压缩操作，我们将编写几个核函数，可以为每个步骤执行所需的操作，并调用最后那些：

  1. 让我们编写一个核函数，通过检查每个元素的值是否大于零来生成一个谓词数组：
__global__ void
predicate_kernel(float *d_predicates, float *d_input, int length)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;

    if (idx >= length) return;

    d_predicates[idx] = d_input[idx] > FLT_ZERO;
}
  1. 然后，我们必须对谓词数组执行前缀和操作。我们将在这里重用之前的实现。之后，我们可以编写一个可以检测扫描数组的地址并将目标元素收集为输出的核函数：
__global__ void
pack_kernel(float *d_output, float *d_input, float *d_predicates, float *d_scanned, int length)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;

    if (idx >= length) return;

    if (d_predicates[idx] != 0.f)
    {
        // addressing
        int address = d_scanned[idx] - 1;

        // gather
        d_output[address] = d_input[idx];
    }
}
  1. 现在，让我们一起调用它们来进行压缩操作：
// predicates
predicate_kernel<<< GRID_DIM, BLOCK_DIM >>>(d_predicates, d_input, length);
// scan
scan_v2(d_scanned, d_predicates, length);
// addressing & gather (pack)
pack_kernel<<< GRID_DIM, BLOCK_DIM >>>(d_output, d_input, d_predicates, d_scanned, length);
  1. 现在，我们有了一个从随机生成的数组中收集到的正数数组：
$ nvcc -run -m64 -std=c++11 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -L/usr/local/cuda/lib -o pack_n_split ./pack_n_split.cu
input    :: -0.4508 -0.0210 -0.4774  0.2750 .... 0.0398  0.4869
pack[cpu]::  0.2750  0.3169  0.1248  0.4241 .... 0.3957  0.2958
pack[gpu]::  0.2750  0.3169  0.1248  0.4241 .... 0.3957  0.2958
SUCCESS!!

通过使用并行前缀和操作，我们可以很容易地并行实现压缩操作。我们的实现从给定数组中压缩正值，但我们可以将其切换到其他条件并且应用压缩操作而不会有困难。现在，让我们来讨论如何将这些压缩元素分发到原始数组。

实现拆分

拆分操作是一个谓词、扫描、地址和拆分的序列。在这个实现中，我们将重用在前一节中创建的地址数组。因此，我们可以跳过之前的步骤，只需从地址数组中实现拆分操作：

  1. 让我们编写拆分核函数，如下所示：
__global__ void
split_kernel(float *d_output, float *d_input, float *d_predicates, float *d_scanned, int length)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;

    if (idx >= length) return;

    if (d_predicates[idx] != 0.f)
    {
        // address
        int address = d_scanned[idx] - 1;

        // split
        d_output[idx] = d_input[address];
    }
}
  1. 现在，我们可以调用核函数，如下所示：
cudaMemcpy(d_input, d_output, sizeof(float) * length, cudaMemcpyDeviceToDevice);
    cudaMemset(d_output, 0, sizeof(float) * length);
    split_kernel<<<GRID_DIM, BLOCK_DIM>>>(d_output, d_input, d_predicates, d_scanned, length);
  1. 由于我们将使用前一步骤的扫描输出，我们将把它复制到输入并清除原始数组。总的来说，我们可以使用 CUDA 进行并行压缩和拆分。这是我们实现的输出。您可以确认它按预期运行：
$ nvcc -run -m64 -std=c++11 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -L/usr/local/cuda/lib -o pack_n_split ./pack_n_split.cu
input    :: -0.4508 -0.0210 -0.4774  0.2750 .... 0.0398  0.4869
pack[cpu]::  0.2750  0.3169  0.1248  0.4241 .... 0.3957  0.2958
pack[gpu]::  0.2750  0.3169  0.1248  0.4241 .... 0.3957  0.2958
SUCCESS!!
split[gpu]   0.0000  0.0000  0.0000  0.2750 .... 0.0398  0.4869
SUCCESS!!

在我们的实现中，我们为正值生成了一个压缩数组和一个拆分数组。由于并行前缀和，我们也可以并行地做到这一点。我们版本的一个主要限制是，它只支持少于 2,048 个元素，因为我们的实现是基于之前的并行前缀和实现的。

N-body

任何 N 体模拟都是一个在物理力的影响下演化的动力学系统的模拟。随着物体不断相互作用，进行数值近似。N 体模拟在物理学和天文学中被广泛使用，例如，科学家可以了解宇宙中粒子的动态。N 体模拟也在许多其他领域中使用，包括计算流体动力学，以便理解湍流流体流动模拟。

解决 N 体模拟的一个相对简单的方法是利用*O(N²)*复杂度的蛮力技术。这种方法在本质上是尴尬地并行的。在算法规模上有各种优化可以减少计算复杂度。可以用来确定近距离相互作用中的力，而不是将所有对应用于整个模拟。即使在这种情况下，为 CUDA 解决力量创建一个内核也是非常有用的，因为它还将提高远场组件的性能。加速一个组件将卸载其他组件的工作，因此整个应用程序都会从加速一个内核中受益。

在 GPU 上实现 N 体模拟

该算法基本上是一个计算力f[ij]的所有对算法，对于一个 NN 网格。一个物体i上的总力/加速度F[i]是该行中所有条目的总和。从并行性的角度来看，这是一个尴尬地并行的任务，复杂度为O(N²)。

从性能角度来看，该应用程序受到内存限制，并且会受到内存带宽的限制。好的一点是，许多数据可以被重复使用并存储在高带宽和低延迟的内存中，比如共享内存。在共享内存中重复使用和存储数据可以减少对全局内存的负载，从而有助于达到峰值计算性能。

以下图表显示了我们将使用的策略：

我们不再从全局内存中反复加载内存，而是利用平铺。我们已经在之前的章节中演示了在矩阵乘法中使用平铺，并在图像应用中使用了它。前面的图表显示了每一行都是并行评估的。平铺大小由可以存储在共享内存中而不影响内核占用率的最大元素数量定义。每个块将数据加载到共享内存中，然后执行同步。一旦数据加载到共享内存中，就在每个块中进行力/加速度计算。可以看到，即使单独的行是并行计算的，为了实现最佳的数据重用，每行中的相互作用是顺序进行的。

N 体模拟实现概述

让我们以伪代码格式回顾一下这个实现，然后解释它的逻辑。在这个例子中，我们使用引力势来说明所有对 N 体模拟中的基本计算形式。实现的代码可以在07_parallel_programming_pattern/05_n-body中找到。按照以下步骤开始：

  1. 用随机变量初始化 n 空间：
data[i] = 2.0f * (rand() / max) - 1.0f
  1. 在一个中间共享内存空间中声明和存储数据，以便有效地重复使用。同步以确保块内的所有线程都能看到共享内存中的更新值：
for (int tile = 0; tile < gridDim.x; tile++) {
... 
__shared__ float3 shared_position[blockDim.x];
float4 temp_position = p[tile * blockDim.x + threadIdx.x];
shared_position[threadIdx.x] = make_float3(temp_position.x, temp_position.y, temp_position.z);
__syncthreads();
...
}
  1. 通过迭代每个块来计算力：
for (int j = 0; j < BLOCK_SIZE; j++) {
    //Calculate Force
    __syncthreads();
}
  1. 最后，使用以下命令将应用程序编译为nvcc编译器：
$nvcc -run --gpu-architecture=sm_70 -o n-body n_body.cu

正如你所看到的，实现 N 体模拟是一个尴尬地并行的任务，而且非常简单。虽然我们在这里实现了基本版本的代码，但存在各种算法变体。你可以利用这个版本作为一个模板，根据对算法的更改进行改进。

直方图计算

在一个尴尬的并行作业中，理想情况下，您会将计算分配给每个线程，这些线程在独立数据上工作，从而不会发生数据竞争。到目前为止，您可能已经意识到有些模式不适合这个类别。其中一种模式是当我们计算直方图时。直方图模式显示了数据项的频率，例如，我们在每个章节中使用 CUDA 这个词的次数

章节，本章中每个字母出现的次数等。直方图采用以下形式：

在这一部分，我们将利用原子操作来串行访问数据，以便获得正确的结果。

编译和执行步骤

直方图提供了关于手头数据集的重要特征，以及有用的见解。例如，在整个图像中，只有少数区域可能存在感兴趣的区域。有时创建直方图用于找出图像中可能存在感兴趣区域的位置。在这个例子中，我们将使用在整个图像中将图像分成块来计算直方图。让我们开始吧：

  1. 准备您的 GPU 应用程序。此代码可以在07_parallel_programming_pattern/08_histogram中找到。

  2. 使用以下命令将您的应用程序编译为nvcc编译器：

$ nvcc -c scrImagePgmPpmPackage.cpp 
$ nvcc -c image_histogram.cu
$ nvcc -run -o image_histogram image_histogram.o scrImagePgmPpmPackage.o

scrImagePgmPpmPackage.cpp文件提供了我们可以用来读取和写入.pgm扩展名图像的源代码。直方图计算代码可以在image_histogram.cu中找到。

理解并行直方图

诸如直方图之类的模式需要原子操作，这意味着以串行方式更新特定地址的值，以消除多个线程之间的争用，从而更新相同的地址。这需要多个线程之间的协调。在这个七步过程中，您可能已经注意到我们使用了私有化。私有化是一种利用低延迟内存（如共享内存）来减少吞吐量和降低延迟的技术，如下图所示：

基本上，我们不是在全局内存上使用原子操作，而是在共享内存上使用原子操作。原因现在应该对您来说是相当明显的。与在共享内存/ L1 缓存上执行相同操作相比，在全局内存上执行原子操作的成本更高。从 Maxwell 架构开始，原子操作得到了硬件支持。私有化的共享内存实现应该从 Maxwell 架构开始为您提供 2 倍的性能。但是，请注意，原子操作仅限于特定的函数和数据大小。

使用 CUDA 原子函数计算直方图

主要地，我们将利用共享内存上的atomicAdd()操作来计算共享内存中每个块的直方图。按照以下步骤在内核中计算直方图：

  1. 为每个块分配与每个块的直方图大小相等的共享内存。由于这是一个 char 图像，我们期望元素在 0-255 的范围内：
__shared__ unsigned int histo_private[256];
  1. 将每个块的共享内存数组初始化为0：
if(localId <256)
    histo_private[localId] = 0;
  1. 同步这一点，以确保块内的所有线程看到初始化的数组：
__syncthreads();
  1. 从全局/纹理内存中读取图像的数据：
unsigned char imageData = tex2D<unsigned char>(texObj,(float)(tidX),(float)(tidY));
  1. 在共享内存上进行atomicAdd()操作：
atomicAdd(&(histo_private[imageData]), 1);
  1. 在写入全局内存之前，在块之间进行同步：
__syncthreads();
  1. 将每个块的直方图写入全局内存：
if(localId <256)
    imageHistogram[histStartIndex+localId] = histo_private[localId];

现在，我们已经完成了在 GPU 上实现直方图计算。

总之，使用共享原子内存很容易实现直方图。由于硬件对共享原子内存的本机支持，这种方法可以在 Maxwell 架构之后的显卡上获得高性能。

使用动态并行性在 CUDA 中进行快速排序

作为任何应用程序的基本构建块的关键算法之一是排序。有许多可用的排序算法已经得到了广泛的研究。最坏时间复杂度、最佳时间复杂度、输入数据特征（数据几乎排序好还是随机的？是键值对吗？是整数还是浮点数？）、原地或非原地内存需求等等，这些都定义了哪种算法适用于哪种应用。一些排序算法属于分治算法的范畴。这些算法适合并行处理，并适用于 GPU 等可以将要排序的数据分割进行排序的架构。其中一个这样的算法是快速排序。正如我们之前所述，快速排序属于分治范畴。它是一个三步方法，如下：

  1. 从需要排序的数组中选择一个元素。这个元素作为枢轴元素。

  2. 第二步是分区，确定所有元素的位置。所有小于枢轴的元素都移到左边，所有大于或等于枢轴的元素都移到枢轴元素的右边。这一步也被称为分区。

  3. 递归地执行步骤 1 和 2，直到所有子数组都被排序。

快速排序的最坏情况复杂度是 O(n²)，这与其他排序过程的最坏情况复杂度为 O(nlogn)相比可能不太理想（例如归并排序和堆排序）。然而，实际上快速排序被认为是有效的。枢轴元素的选择可以经过考虑，有时也可以随机选择，以使最坏情况复杂度几乎不会发生。此外，与其他排序算法相比，快速排序的内存负载和需求较少，例如归并排序需要额外的存储空间。更实际的快速排序实现使用随机化版本。随机化版本的期望时间复杂度为 O(nlogn)。最坏情况复杂度在随机化版本中也是可能的，但它不会发生在特定模式（例如排序好的数组）上，随机化快速排序在实践中表现良好。

虽然我们可以写一整章关于排序算法的特性，但我们计划只覆盖 CUDA 的特性，这将帮助您在 GPU 上高效实现快速排序。在本节中，我们将使用从 CUDA 6.0 和 GPU 架构 3.5 开始引入的动态并行性。

现在，让我们回顾一下动态并行性是如何对排序算法做出贡献的。

快速排序和 CUDA 动态并行性

快速排序算法要求递归地启动内核。到目前为止，我们所见过的算法是通过 CPU 一次调用内核。内核执行完毕后，我们返回到 CPU 线程，然后重新启动它。这样做会导致将控制权交还给 CPU，并且可能导致 CPU 和 GPU 之间的数据传输，这是一项昂贵的操作。以前在 GPU 上高效实现需要递归等特性的算法（如快速排序）曾经非常困难。从 GPU 架构 3.5 和 CUDA 5.0 开始，引入了一个名为动态并行性的新特性。

动态并行性允许内核内的线程在不将控制权返回给 CPU 的情况下从 GPU 上启动新的内核。动态一词来自于它基于运行时数据的动态性。多个内核可以同时由线程启动。以下图表简化了这个解释：

如果我们将这个概念转化为快速排序的执行方式，它会看起来像这样：

深度 0 是来自 CPU 的调用。对于每个子数组，我们启动两个内核：一个用于左数组，一个用于右数组。递归在达到内核的最大深度或元素数量小于 32（即 warp 大小）后停止。为了使内核的启动在非零流中是异步的，以便子数组内核可以独立启动，我们需要在每次内核启动之前创建一个流：

cudaStream_t s;
cudaStreamCreateWithFlags( &s, cudaStreamNonBlocking );
cdp_simple_quicksort<<< 1, 1, 0, s >>>(data, left, nright, depth+1);
cudaStreamDestroy( s );

这是一个非常重要的步骤，否则内核的启动可能会被序列化。有关流的更多细节，请参考多 GPU 内核。

CUDA 的 Quicksort

对于我们的 Quicksort 实现，我们将利用动态并行性来递归启动 GPU 内核。实现 Quicksort 的主要步骤如下：

  1. CPU 启动第一个内核：内核以一个块和一个线程启动。左元素是数组的开始，右元素是数组的最后一个元素（基本上是整个数组）：
int main(int argc, char **argv)
{ ...
    cdp_simple_quicksort<<< 1, 1 >>>(data, left, right, 0);
}
  1. 限制检查：在从内核内部启动内核之前检查两个条件。首先，检查我们是否已经达到硬件允许的最大深度限制。其次，我们需要检查子数组中要排序的元素数量是否小于 warp 大小（32）。如果其中一个条件为真，那么我们必须按顺序执行选择排序，而不是启动一个新的内核：
__global__ void cdp_simple_quicksort( unsigned int *data, int left, int right, int depth )
{ ...

if( depth >= MAX_DEPTH || right-left <= INSERTION_SORT )
 {
     selection_sort( data, left, right );
     return;
 }
  1. 分区：如果满足前面的条件，那么将数组分成两个子数组，并启动两个新的内核，一个用于左数组，另一个用于右数组。如果你仔细看下面的代码，你会发现我们是从内核内部启动内核的：
__global__ void cdp_simple_quicksort( unsigned int *data, int left, int right, int depth ) {
...
while(lptr <= rptr)
 {
     // Move the left pointer as long as the 
     // pointed element is smaller than the pivot.
     // Move the right pointer as long as the 
     // pointed element is larger than the pivot.
     // If the swap points are valid, do the swap!

     // Launch a new block to sort the left part.
     if(left < (rptr-data))
     { // Create a new stream for the eft sub array
        cdp_simple_quicksort<<< 1, 1, 0, s 
                            >>>(data, left, nright, depth+1);
     }
    // Launch a new block to sort the right part.
    if((lptr-data) < right)
     {//Create stream for the right sub array
         cdp_simple_quicksort<<< 1, 1, 0, s1 
                             >>>(data, nleft, right, depth+1);
     }
 }
  1. 执行代码：实现的代码可以在07_parallel_programming_pattern/06_quicksort中找到。使用以下命令使用nvcc编译您的应用程序：
$nvcc -o quick_sort --gpu-architecture=sm_70 -rdc=true quick_sort.cu

如你所见，我们在编译中添加了两个标志：

  • -- gpu-architecture=sm_70：这个标志告诉nvcc为 Volta GPU 编译和生成二进制/ptx。如果你没有特别添加这个标志，编译器会尝试从sm_20（即 Fermi 代）兼容的代码编译，直到新架构sm_70（即 Volta）。由于旧一代的卡不支持动态并行性，编译将失败。

  • -rdc=true：这是一个关键参数，它在 GPU 上启用动态并行性。

动态并行性指南和约束

虽然动态并行性为我们提供了在 GPU 上移植 Quicksort 等算法的机会，但需要遵循一些基本规则和指南。

编程模型规则：基本上，所有 CUDA 编程模型规则都适用：

  • 内核启动是每个线程异步的。

  • 同步只允许在块内进行。

  • 创建的流在一个块内共享。

  • 事件可用于创建流间依赖关系。

内存一致性规则：

  • 子内核在启动时看到父内核的状态。

  • 父内核只能在同步后看到子内核所做的更改。

  • 本地和共享内存通常是私有的，父内核无法传递或访问。

指南：

  • 重要的是要理解，每次内核启动都会增加延迟。从另一个内核内部启动内核的延迟随着新架构的推出逐渐减少。

  • 虽然启动吞吐量比主机高一个数量级，但最大深度可以设置限制。最新一代卡允许的最大深度是 24。

  • 从内核内部执行cudaDeviceSynchronize()是一个非常昂贵的操作，应尽量避免。

  • 在全局内存上预先分配了额外的内存，以便在启动之前存储内核。

  • 如果内核失败，错误只能从主机上看到。因此，建议您使用-lineinfo标志以及cuda-memcheck来定位错误的位置。

基数排序

另一个非常流行的排序算法是基数排序，因为它在顺序机器上非常快。基数排序的基本策略是每个元素都按位排序。让我们看一个简单的例子来解释基数排序涉及的步骤：

假设要排序的元素如下：

值71441

这些数字的等效二进制值如下：

位0111111001000001

第一步是根据第 0 位进行排序。这些数字的第 0 位如下：

0 位1001

根据第 o 位排序基本上意味着所有的零都在左边。所有的 1 都在右边，同时保持元素的顺序：

第 0 位上的排序值14471
根据第 0 位排序的位1110010001110001

第 0 位完成后，我们继续到第一位。根据第一位排序后的结果如下：

第一位上的排序值41471
根据第一位排序的位0100111001110001

然后，我们继续到下一个更高的位，直到所有的位都结束。最终结果如下：

所有位上的排序值1471
根据所有位排序的位0001010001111110

正如您所看到的，在这个例子中我们设置的上限是 4 位。对于更大的数字，比如整数，这将持续到 32 位，因为整数是 32 位的。

现在我们已经了解了这个算法，让我们看看如何在 GPU 中实现它。与本章中的其他部分相比，我们将采取两种方法来展示 CUDA 生态系统，以便我们可以实现/使用基数排序。

选项 1：我们将使用翘曲级别来对 32 个元素进行基数排序。这样做的原因是我们希望利用基数排序来向您介绍翘曲级别原语。

选项 2：我们将使用 CUDA 工具包的一部分 Thrust 库。它实现了通用基数排序。最好的实现是重用。由于 Thrust 已经提供了最好的基数排序实现之一，我们将使用它。

两种方法

为了方便您的理解，让我们从示例代码开始。在这个例子中，我们将使用翘曲级别原语和 Thrust 库来实现/使用基数排序。示例代码可以在07_parallel_programming_pattern/07_radixsort中找到。

使用以下命令使用nvcc编译器编译您的应用程序：

  • 翘曲级别原语版本：
$ nvcc -run -o radix_warp_sort radix_warp_sort.cu
  • Thrust 库版本：
$ nvcc -run -o radix_thrust_sort thrust_radix_sort.cu

这两个例子展示了 GPU 给出的排序输出。现在，让我们详细了解这些操作是如何实现的。

方法 1 - 翘曲级别原语

让我们看看 CUDA 翘曲级别原语是如何在代码中实现我们的算法的：

  1. 首先，将数据从全局内存加载到共享内存中：
__shared__ unsigned int s_data[WARP_SIZE*2];

内存的大小等于翘曲大小，*2，以便它可以实现乒乓缓冲区。

  1. 从低位到高位循环：
for (int i = MIN_BIT_POS; i <= MAX_BIT_POS; i++){ ... }
  1. 获取当前的掩码：
unsigned int bit  = data&bit_mask;
  1. 获取 1 和 0 的数量（直方图）：
unsigned int active = __activemask();
unsigned int ones = __ballot_sync(active,bit);
unsigned int zeroes = ~ones;

  1. 获取当前位数为零（0）的线程的位置（前缀和）。

  2. 获取当前位数为一（1）的线程的位置（前缀和）：

if (!bit) // threads with a zero bit
 // get my position in ping-pong buffer
 pos = __popc(zeroes&thread_mask);
 else // threads with a one bit
 // get my position in ping-pong buffer
 pos = __popc(zeroes)+__popc(ones&thread_mask);
  1. 将数据存储在乒乓共享缓冲区内存中：
 s_data[pos-1+offset] = data;

  1. 重复步骤 2-6，直到达到上限位。

  2. 从共享内存中将最终结果存储到全局内存中：

d_data[threadIdx.x] = s_data[threadIdx.x+offset];

也许对于您来说，直方图和前缀和突然出现可能不太清楚。让我们详细讨论这个实现，以便我们可以理解如何使用翘曲级别原语来实现相同的功能。

在本节的开头，我们描述了如何使用示例进行排序。然而，我们没有涵盖的是如何找出需要交换的元素的位置。基数排序可以使用基本原语（如直方图和前缀和）来实现，因此可以很容易地在 GPU 上实现。

让我们重新审视我们看过的示例，并收集其细节，包括直方图和前缀和的步骤。以下表格显示了在每个位上迭代进行的各种计算：

值71441
二进制0111111001000001
位 01001
直方图前缀和2022
偏移0011
新索引（前缀和和偏移）2013

让我们解释前面表格中显示的每一项计算，如下所示：

  1. 首先，我们为第 0 位位置的元素构建直方图，包括具有 0 和 1 的元素的数量：

直方图：零位（2 个值），一位（2 个值）

  1. 然后，我们对这些值进行排他性前缀和。前缀和可以定义为所有先前值的总和。在我们的情况下，我们分别对 0 位和 1 位进行这样的操作。

  2. 最后，我们根据前缀和的值移动元素。

我们用来找到直方图和前缀和的 warp 级原语分别是__ballot_sync()和__popc()。

__ballot_sync() API 评估 warp 的所有活动线程的谓词，并返回一个整数，其第 N 位设置为 1，当且仅当谓词对于 warp 的第 N 个线程求值为非零时。__popc()用于计算整数的数量，被设置为 1。

在 CUDA 编程模型中，我们已经看到最小的执行单元是一个 warp（32 个线程）。CUDA 提供了各种 warp 级原语，可以进行细粒度控制，在许多应用中可以实现更好的性能。我们在上一节介绍了一个这样的原语__ballot_sync()。其他重要的 warp 级原语包括shuffle指令，用于特定的 warp 级归约。shuffle指令已经在本书中介绍过。如果您已经达到了 CUDA 的忍者程序员水平，那么我们建议您查看 CUDA API 指南，以了解更多这些 warp 级原语。

这完成了使用 warp 级原语描述基数排序。现在，让我们看看基于 Thrust 库的实现。

方法 2 - 基于 Thrust 的基数排序

基于 Thrust 的基数排序是基数排序的通用实现，对于不同类型的数据（如整数、浮点数或键值对）都能很好地工作。我们想再次强调排序是一个经过深入研究的算法，因此有其并行实现。因此，我们建议在自己实现之前重用现有的库。

使用 Thrust 进行基数排序的步骤如下：

  1. 导入相关的头文件（Thrust 是一个仅包含头文件的库，类似于 STL）：
#include <thrust/device_vector.h>
#include <thrust/sort.h>
  1. 声明并初始化设备向量：
//declare a device vector of size N
thrust::device_vector<int> keys(N);
//Generate a random number generator engine
thrust::default_random_engine r(12);
//create a distribution engine which will create integer values
thrust::uniform_int_distribution<int> d(10, 99);
//Fill the array with randon values
for(size_t i = 0; i < v.size(); i++)
    v[i] = d(r);
  1. 对初始化的设备向量进行排序：
thrust::sort(keys.begin(), keys.end());

使用这个库提供了一种更简单和更健壮的方法。Thrust 提供了不同类型的排序方法，包括整数和浮点数的基数排序。或者，您可以创建一个自定义比较器来进行自定义排序，例如按照偶数后面是奇数的顺序排序，按降序排序等等。如果您想了解更多关于基于 Thrust 的排序示例，建议您查看 CUDA 提供的示例示例。

现在，我们已经看过了在 GPU 上实现基数排序的两种方法。

总结

在本章中，我们看了 CUDA 中常用算法和模式的实现。这些算法和模式是常见的。我们涵盖了矩阵乘法和卷积滤波中的基本优化技术。然后，我们扩展了讨论，介绍了如何通过使用前缀和、N 体、直方图和排序来并行化问题。为此，我们使用了专门的 GPU 知识、库和较低级别的原语。

我们所涵盖的许多算法都是在 CUDA 库中实现的。例如，矩阵乘法在 cuBLAS 库中，而卷积在 CUDNN 库中。此外，我们还涵盖了基数排序实现中的两种方法：使用 Thrust 库或 warp 级原语进行直方图计算。

现在您已经看到了这些模式如何在常用库中实现，下一个合乎逻辑的步骤是看看我们如何可以使用这些库。这就是我们将在下一章中要做的事情。

avatar
文章
阅读
粉丝