Python GPU 编程实用指南（三）

原文：zh.annas-archive.org/md5/ef7eb3c148e0cfdfe01c331f2f01557c

译者：飞龙

协议：CC BY-NC-SA 4.0

第十章：使用已编译的 GPU 代码

在本书的过程中，我们通常依赖 PyCUDA 库自动为我们接口我们的内联 CUDA-C 代码，使用即时编译和与 Python 代码的链接。然而，我们可能还记得，有时编译过程可能需要一段时间。在第三章中，使用 PyCUDA 入门，我们甚至详细看到编译过程如何导致减速，以及内联代码何时被编译和保留可能是相当随意的。在某些情况下，这可能会给应用程序带来不便，或者在实时系统的情况下甚至是不可接受的。

为此，我们最终将看到如何从 Python 使用预编译的 GPU 代码。特别是，我们将看看三种不同的方法来做到这一点。首先，我们将看看如何通过编写一个主机端 CUDA-C 函数来间接启动 CUDA 内核。这种方法将涉及使用标准 Python Ctypes 库调用主机端函数。其次，我们将把我们的内核编译成所谓的 PTX 模块，这实际上是一个包含已编译二进制 GPU 的 DLL 文件。然后，我们可以使用 PyCUDA 加载此文件并直接启动我们的内核。最后，我们将通过查看如何编写我们自己的完整 Ctypes 接口来结束本章，以使用 CUDA Driver API。然后，我们可以使用 Driver API 中的适当函数加载我们的 PTX 文件并启动内核。

本章的学习成果如下：

使用 Ctypes 模块启动编译后（主机端）的代码
使用 Ctypes 使用主机端 CUDA C 包装器从 Python 启动内核
如何将 CUDA C 模块编译为 PTX 文件
如何将 PTX 模块加载到 PyCUDA 中以启动预编译的内核
如何编写自定义 Python 接口以使用 CUDA Driver API

使用 Ctypes 启动编译后的代码

我们现在将简要概述 Python 标准库中的 Ctypes 模块。Ctypes 用于调用来自 Linux.so（共享对象）或 Windows.DLL（动态链接库）预编译二进制文件的函数。这将使我们摆脱纯 Python 的世界，并与已用编译语言编写的库和代码进行接口，特别是 C 和 C++ - 恰好 Nvidia 只为与我们的 CUDA 设备进行接口提供这样的预编译二进制文件，因此如果我们想绕过 PyCUDA，我们将不得不使用 Ctypes。

让我们从一个非常基本的例子开始：我们将向您展示如何直接从 Ctypes 调用printf。打开一个 IPython 实例，键入import ctypes。现在我们将看看如何从 Ctypes 调用标准的printf函数。首先，我们必须导入适当的库：在 Linux 中，通过键入libc = ctypes.CDLL('libc.so.6')加载 LibC 库（在 Windows 中，将'libc.so.6'替换为'msvcrt.dll'）。现在我们可以通过在 IPython 提示符中键入libc.printf("Hello from ctypes!\n")直接调用printf。自己试试吧！

现在让我们试试其他东西：从 IPython 键入libc.printf("Pi is approximately %f.\n", 3.14)；您应该会收到一个错误。这是因为3.14没有适当地从 Python 浮点变量转换为 C 双精度变量 - 我们可以使用 Ctypes 这样做：

libc.printf("Pi is approximately %f.\n", ctypes.c_double(3.14))

输出应该如预期那样。与从 PyCUDA 启动 CUDA 内核的情况一样，我们必须同样小心地将输入转换为 Ctypes 函数。

请务必确保将任何从 Python 使用 Ctypes 调用的函数的输入适当地转换为适当的 C 数据类型（在 Ctypes 中，这些类型以 c_ 开头：c_float、c_double、c_char、c_int等）。

再次重温 Mandelbrot 集

让我们重新审视一下我们在第一章和第三章中看到的 Mandelbrot 集合，为什么使用 GPU 编程？和使用 PyCUDA 入门。首先，我们将编写一个完整的 CUDA 核函数，它将根据一组特定的参数计算 Mandelbrot 集合，以及一个适当的主机端包装函数，我们稍后可以从 Ctypes 接口调用。我们将首先将这些函数编写到一个单独的 CUDA-C.cu源文件中，然后使用 NVCC 编译成 DLL 或.so二进制文件。最后，我们将编写一些 Python 代码，以便我们可以运行我们的二进制代码并显示 Mandelbrot 集合。

我们现在将运用我们对 Ctypes 的知识，从 Python 中启动一个预编译的 CUDA 核函数，而不需要 PyCUDA 的任何帮助。这将要求我们在 CUDA-C 中编写一个主机端核函数启动器包装函数，我们可以直接调用，它本身已经编译成了一个动态库二进制文件，其中包含任何必要的 GPU 代码——即在 Windows 上的动态链接库（DLL）二进制文件，或者在 Linux 上的共享对象（so）二进制文件。

当然，我们将首先编写我们的 CUDA-C 代码，所以打开你最喜欢的文本编辑器并跟着做。我们将从标准的include语句开始：

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

我们现在将直接开始编写我们的核函数。请注意代码中的extern "C"，这将允许我们在外部链接到这个函数：

extern "C" __global__ void mandelbrot_ker(float * lattice, float * mandelbrot_graph, int max_iters, float upper_bound_squared, int lattice_size)
{

让我们思考一分钟关于这将如何工作：我们将使用一个单一的一维数组来存储实部和虚部，称为lattice，其长度为lattice_size。我们将使用这个数组来计算一个形状为(lattice_size, lattice_size)的二维 Mandelbrot 图形，存储在预先分配的数组mandelbrot_graph中。我们将指定每个点检查发散的迭代次数为max_iters，通过使用upper_bound_squared提供其平方值来指定之前的最大上限值。（我们稍后会看一下使用平方的动机。）

我们将在一维网格/块结构上启动这个核函数，每个线程对应于 Mandelbrot 集合图像中的一个点。然后我们可以确定相应点的实部/虚部 lattice 值，如下所示：

    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if ( tid < lattice_size*lattice_size )
    {
        int i = tid % lattice_size;
        int j = lattice_size - 1 - (tid / lattice_size);

        float c_re = lattice[i];
        float c_im = lattice[j];

让我们花一分钟来谈谈这个。首先，记住我们可能需要使用稍多于必要的线程，所以重要的是我们检查线程 ID 是否对应于输出图像中的某个点，使用if语句。让我们还记住，输出数组mandelbrot_graph将以按行方式存储为一维数组，表示为二维图像，我们将使用tid作为写入该数组的索引。我们将使用i和j，以及复平面上图形的x和y坐标。由于 lattice 是一系列从小到大排序的实值，我们将不得不颠倒它们的顺序以获得适当的虚值。另外，请注意，我们将在这里使用普通的浮点数，而不是某种结构或对象来表示复数值。由于每个复数中都有实部和虚部，我们将在这里使用两个浮点数来存储与该线程的 lattice 点对应的复数（c_re和c_im）。

我们将设置两个更多的变量来处理发散检查，z_re和z_im，并在检查发散之前将该线程的图上的点的初始值设置为1：

        float z_re = 0.0f;
        float z_im = 0.0f;

        mandelbrot_graph[tid] = 1;

现在我们将检查发散；如果在max_iters次迭代后发散，我们将把点设置为0。否则，它将保持为 1：

        for (int k = 0; k < max_iters; k++)
        {
            float temp;

            temp = z_re*z_re - z_im*z_im + c_re;
            z_im = 2*z_re*z_im + c_im;
            z_re = temp;

            if ( (z_re*z_re + z_im*z_im) > upper_bound_squared )
            {
                mandelbrot_graph[tid] = 0;
                break;
            }
        }

在我们继续之前，让我们谈一分钟关于这一块代码。让我们记住，曼德勃罗集的每次迭代都是通过复数乘法和加法来计算的，例如，z_new = z*z + c。由于我们不是使用将为我们处理复数值的类，前面的操作正是我们需要做的，以计算z的新的实部和虚部值。我们还需要计算绝对值并查看是否超过特定值——记住，复数的绝对值，c = x + iy，是用√(x²+y²)来计算的。在这里计算上限的平方然后将其插入内核中，实际上会节省我们一些时间，因为这样可以节省我们在每次迭代中计算z_re*z_re + z_im*z_im的平方根的时间。

我们现在基本上已经完成了这个内核——我们只需要关闭if语句并从内核返回，然后我们就完成了：

    }
    return;
}

然而，我们还没有完全完成。我们需要编写一个只有extern "C"的主机端包装函数，在 Linux 的情况下，以及在 Windows 的情况下，只有extern "C" __declspec(dllexport)。（与编译的 CUDA 内核相反，如果我们想要能够从 Ctypes 在 Windows 中访问主机端函数，这个额外的单词是必要的。）我们放入这个函数的参数将直接对应于进入内核的参数，除了这些参数将存储在主机上：

extern "C" __declspec(dllexport) void launch_mandelbrot(float * lattice,  float * mandelbrot_graph, int max_iters, float upper_bound, int lattice_size)
{

现在，我们将需要分配足够的内存来存储在 GPU 上的晶格和输出，然后使用cudaMemcpy将晶格复制到 GPU 上：

    int num_bytes_lattice = sizeof(float) * lattice_size;
    int num_bytes_graph = sizeof(float)* lattice_size*lattice_size;

    float * d_lattice;
    float * d_mandelbrot_graph;

    cudaMalloc((float **) &d_lattice, num_bytes_lattice);
    cudaMalloc((float **) &d_mandelbrot_graph, num_bytes_graph);

    cudaMemcpy(d_lattice, lattice, num_bytes_lattice, cudaMemcpyHostToDevice);

像我们的许多其他内核一样，我们将在一维网格上启动大小为 32 的一维块。我们将取输出点数除以 32 的上限值，以确定网格大小，如下所示：

    int grid_size = (int)  ceil(  ( (double) lattice_size*lattice_size ) / ( (double) 32 ) );

现在我们准备使用传统的 CUDA-C 三角形括号来启动我们的内核，以指定网格和块大小。请注意，我们在这里提前求出了上限的平方：

    mandelbrot_ker <<< grid_size, 32 >>> (d_lattice,  d_mandelbrot_graph, max_iters, upper_bound*upper_bound, lattice_size);

现在我们只需要在完成这些操作后将输出复制到主机上，然后在适当的数组上调用cudaFree。然后我们可以从这个函数返回：

    cudaMemcpy(mandelbrot_graph, d_mandelbrot_graph, num_bytes_graph, cudaMemcpyDeviceToHost);    
    cudaFree(d_lattice);
    cudaFree(d_mandelbrot_graph);
}

有了这个，我们已经完成了所有需要的 CUDA-C 代码。将其保存到名为mandelbrot.cu的文件中，然后继续下一步。

您还可以从github.com/btuomanen/handsongpuprogramming/blob/master/10/mandelbrot.cu下载此文件。

编译代码并与 Ctypes 进行接口

现在让我们将刚刚编写的代码编译成 DLL 或.so二进制文件。这实际上相当简单：如果你是 Linux 用户，请在命令行中输入以下内容将此文件编译成mandelbrot.so：

nvcc -Xcompiler -fPIC -shared -o mandelbrot.so mandelbrot.cu

如果你是 Windows 用户，请在命令行中输入以下内容将文件编译成mandelbrot.dll：

nvcc -shared -o mandelbrot.dll mandelbrot.cu

现在我们可以编写我们的 Python 接口。我们将从适当的导入语句开始，完全排除 PyCUDA，只使用 Ctypes。为了方便使用，我们将直接从 Ctypes 导入所有类和函数到默认的 Python 命名空间中，如下所示：

from __future__ import division
from time import time
import matplotlib
from matplotlib import pyplot as plt
import numpy as np
from ctypes import *

让我们使用 Ctypes 为launch_mandelbrot主机端函数设置一个接口。首先，我们将不得不加载我们编译的 DLL 或.so文件，Linux 用户当然需要将文件名更改为mandelbrot.so：

mandel_dll = CDLL('./mandelbrot.dll')

现在我们可以从库中获取对launch_mandelbrot的引用，就像这样；我们将简称它为mandel_c：

mandel_c = mandel_dll.launch_mandelbrot

现在，在使用 Ctypes 调用函数之前，我们将不得不让 Ctypes 知道输入类型是什么。让我们记住，对于launch_mandelbrot，输入是float-pointer，float-pointer，integer，float和integer。我们使用argtypes参数设置这一点，使用适当的 Ctypes 数据类型（c_float，c_int），以及 Ctypes 的POINTER类：

mandel_c.argtypes = [POINTER(c_float), POINTER(c_float), c_int, c_float, c_int]

现在让我们编写一个 Python 函数来为我们运行这个。我们将使用breadth指定正方形输出图像的宽度和高度，并且在复杂格的实部和虚部中指定最小和最大值。我们还将指定最大迭代次数以及上限：

def mandelbrot(breadth, low, high, max_iters, upper_bound):

现在，我们将使用 NumPy 的linspace函数创建我们的格点数组，就像这样：

 lattice = np.linspace(low, high, breadth, dtype=np.float32)

让我们记住，我们将不得不传递一个预先分配的浮点数组给launch_mandelbrot，以便以输出图的形式得到输出。我们可以通过调用 NumPy 的empty命令来设置一个适当形状和大小的数组，这将充当 C 的malloc调用：

    out = np.empty(shape=(lattice.size,lattice.size), dtype=np.float32)

现在，我们准备计算 Mandelbrot 图。请注意，我们可以通过使用它们的ctypes.data_as方法和相应的类型将 NumPy 数组传递给 C。在我们这样做之后，我们可以返回输出；也就是说，Mandelbrot 图以二维 NumPy 数组的形式：

 mandel_c(lattice.ctypes.data_as(POINTER(c_float)), out.ctypes.data_as(POINTER(c_float)), c_int(max_iters), c_float(upper_bound), c_int(lattice.size) ) 
 return out

现在，让我们编写我们的主函数来计算、计时和使用 Matplotlib 查看 Mandelbrot 图：

if __name__ == '__main__':
    t1 = time()
    mandel = mandelbrot(512,-2,2,256, 2)
    t2 = time()
    mandel_time = t2 - t1
    print 'It took %s seconds to calculate the Mandelbrot graph.' % mandel_time
    plt.figure(1)
    plt.imshow(mandel, extent=(-2, 2, -2, 2))
    plt.show()

我们现在将尝试运行这个。您应该会得到一个看起来与第一章的 Mandelbrot 图以及第三章的为什么使用 GPU 编程和使用 PyCUDA 入门中的图形完全相同的输出：

这个 Python 示例的代码也可以在 GitHub 存储库的mandelbrot_ctypes.py文件中找到。

编译和启动纯 PTX 代码

我们刚刚看到了如何从 Ctypes 调用纯 C 函数。在某些方面，这可能看起来有点不够优雅，因为我们的二进制文件必须包含主机代码以及编译后的 GPU 代码，这可能看起来很麻烦。我们是否可以只使用纯粹的编译后的 GPU 代码，然后适当地将其启动到 GPU 上，而不是每次都编写一个 C 包装器？幸运的是，我们可以。

NVCC 编译器将 CUDA-C 编译为PTX（Parallel Thread Execution），这是一种解释的伪汇编语言，与 NVIDIA 的各种 GPU 架构兼容。每当您使用 NVCC 将使用 CUDA 核心的程序编译为可执行的 EXE、DLL、.so或 ELF 文件时，该文件中将包含该核心的 PTX 代码。我们还可以直接编译具有 PTX 扩展名的文件，其中将仅包含从编译后的 CUDA .cu 文件中编译的 GPU 核心。幸运的是，PyCUDA 包括一个接口，可以直接从 PTX 加载 CUDA 核心，使我们摆脱了即时编译的枷锁，同时仍然可以使用 PyCUDA 的所有其他好功能。

现在让我们将刚刚编写的 Mandelbrot 代码编译成一个 PTX 文件；我们不需要对它进行任何更改。只需在 Linux 或 Windows 的命令行中输入以下内容：

nvcc -ptx -o mandelbrot.ptx mandelbrot.cu

现在让我们修改上一节的 Python 程序，以使用 PTX 代码。我们将从导入中删除ctypes并添加适当的 PyCUDA 导入：

from __future__ import division
from time import time
import matplotlib
from matplotlib import pyplot as plt
import numpy as np
import pycuda
from pycuda import gpuarray
import pycuda.autoinit

现在让我们使用 PyCUDA 的module_from_file函数加载 PTX 文件，就像这样：

mandel_mod = pycuda.driver.module_from_file('./mandelbrot.ptx')

现在我们可以使用get_function来获取对我们的核心的引用，就像我们用 PyCUDA 的SourceModule一样：

mandel_ker = mandel_mod.get_function('mandelbrot_ker')

我们现在可以重写 Mandelbrot 函数，以处理使用适当的gpuarray对象和typecast输入的核心。（我们不会逐行讨论这个，因为在这一点上它的功能应该是显而易见的。）：

def mandelbrot(breadth, low, high, max_iters, upper_bound):
    lattice = gpuarray.to_gpu(np.linspace(low, high, breadth, dtype=np.   
    out_gpu = gpuarray.empty(shape=(lattice.size,lattice.size), dtype=np.float32)
    gridsize = int(np.ceil(lattice.size**2 / 32))
    mandel_ker(lattice, out_gpu, np.int32(256), np.float32(upper_bound**2), np.int32(lattice.size), grid=(gridsize, 1, 1), block=(32,1,1))
    out = out_gpu.get()

    return out

main函数将与上一节完全相同：

if __name__ == '__main__':
    t1 = time()
    mandel = mandelbrot(512,-2,2,256,2)
    t2 = time()
    mandel_time = t2 - t1
    print 'It took %s seconds to calculate the Mandelbrot graph.' % mandel_time
    plt.figure(1)
    plt.imshow(mandel, extent=(-2, 2, -2, 2))
    plt.show()

现在，尝试运行这个来确保输出是正确的。您可能还会注意到与 Ctypes 版本相比的一些速度改进。

这段代码也可以在本书 GitHub 存储库的“10”目录下的mandelbrot_ptx.py文件中找到。

为 CUDA Driver API 编写包装器

现在我们将看看如何使用 Ctypes 编写我们自己的包装器，用于一些预打包的二进制 CUDA 库函数。特别是，我们将为 CUDA 驱动程序 API 编写包装器，这将允许我们执行所有基本 GPU 使用所需的操作，包括 GPU 初始化、内存分配/传输/释放、内核启动和上下文创建/同步/销毁。这是一个非常强大的知识；它将允许我们在不经过 PyCUDA 的情况下使用 GPU，也不需要编写任何繁琐的主机端 C 函数包装器。

现在我们将编写一个小模块，它将作为CUDA 驱动程序 API的包装库。让我们谈一分钟这对我们意味着什么。驱动程序 API 与CUDA Runtime API略有不同，技术性稍高，后者是我们在 CUDA-C 文本中一直在使用的。驱动程序 API 旨在与常规 C/C++编译器一起使用，而不是与 NVCC 一起使用，具有一些不同的约定，如使用cuLaunchKernel函数启动内核，而不是使用<<< gridsize, blocksize >>>括号表示法。这将允许我们直接访问使用 Ctypes 从 PTX 文件启动内核所需的函数。

让我们通过将所有 Ctypes 导入模块的命名空间，并导入 sys 模块来开始编写此模块。我们将通过使用sys.platform检查系统的操作系统（nvcuda.dll或libcuda.so）来加载适当的库文件，使我们的模块可以在 Windows 和 Linux 上使用。

from ctypes import *
import sys
if 'linux' in sys.platform:
 cuda = CDLL('libcuda.so')
elif 'win' in sys.platform:
 cuda = CDLL('nvcuda.dll')

我们已经成功加载了 CUDA 驱动程序 API，现在我们可以开始为基本 GPU 使用编写必要函数的包装器。随着我们的进行，我们将查看每个驱动程序 API 函数的原型，这通常是在编写 Ctypes 包装器时必要的。

鼓励读者在官方 Nvidia CUDA 驱动程序 API 文档中查找我们将在本节中使用的所有函数，该文档可在此处找到：docs.nvidia.com/cuda/cuda-driver-api/。

让我们从驱动程序 API 中最基本的函数cuInit开始，它将初始化驱动程序 API。这需要一个用于标志的无符号整数作为输入参数，并返回类型为 CUresult 的值，实际上只是一个整数值。我们可以这样编写我们的包装器：

cuInit = cuda.cuInit
cuInit.argtypes = [c_uint]
cuInit.restype = int

现在让我们开始下一个函数cuDeviceCount，它将告诉我们在计算机上安装了多少个 NVIDIA GPU。它以整数指针作为其唯一输入，实际上是通过引用返回的单个整数输出值。返回值是另一个 CUresult 整数——所有函数都将使用 CUresult，这是所有驱动程序 API 函数的错误值的标准化。例如，如果我们看到任何函数返回0，这意味着结果是CUDA_SUCCESS，而非零结果将始终意味着错误或警告：

cuDeviceGetCount = cuda.cuDeviceGetCount
cuDeviceGetCount.argtypes = [POINTER(c_int)]
cuDeviceGetCount.restype = int

现在让我们为cuDeviceGet编写一个包装器，它将通过引用在第一个输入中返回设备句柄。这将对应于第二个输入中给定的序号 GPU。第一个参数的类型是CUdevice *，实际上只是一个整数指针：

cuDeviceGet = cuda.cuDeviceGet
cuDeviceGet.argtypes = [POINTER(c_int), c_int]
cuDeviceGet.restype = int

让我们记住，每个 CUDA 会话都需要至少一个 CUDA 上下文，可以将其类比为在 CPU 上运行的进程。由于 Runtime API 会自动处理这一点，在这里我们将不得不在使用设备之前手动在设备上创建上下文（使用设备句柄），并且在 CUDA 会话结束时销毁此上下文。

我们可以使用cuCtxCreate函数创建一个 CUDA 上下文，它当然会创建一个上下文。让我们看看文档中列出的原型：

 CUresult cuCtxCreate ( CUcontext* pctx, unsigned int flags, CUdevice dev )

当然，返回值是CUresult。第一个输入是指向名为CUcontext的类型的指针，实际上它本身是 CUDA 内部使用的特定 C 结构的指针。由于我们从 Python 对CUcontext的唯一交互将是保持其值以在其他函数之间传递，我们可以将CUcontext存储为 C void *类型，用于存储任何类型的通用指针地址。由于这实际上是指向 CU 上下文的指针（再次，它本身是指向内部数据结构的指针——这是另一个按引用返回值），我们可以将类型设置为普通的void *，这是 Ctypes 中的c_void_p类型。第二个值是无符号整数，而最后一个值是要在其上创建新上下文的设备句柄——让我们记住这实际上只是一个整数。我们现在准备为cuCtxCreate创建包装器：

cuCtxCreate = cuda.cuCtxCreate
cuCtxCreate.argtypes = [c_void_p, c_uint, c_int]
cuCtxCreate.restype = int

您可以始终在 C/C++（Ctypes 中的c_void_p）中使用void *类型指向任意数据或变量，甚至结构和对象，其定义可能不可用。

下一个函数是cuModuleLoad，它将为我们加载一个 PTX 模块文件。第一个参数是一个 CUmodule 的引用（同样，我们可以在这里使用c_void_p），第二个是文件名，这将是一个典型的以空字符结尾的 C 字符串——这是一个char *，或者在 Ctypes 中是c_char_p：

cuModuleLoad = cuda.cuModuleLoad
cuModuleLoad.argtypes = [c_void_p, c_char_p]
cuModuleLoad.restype = int

下一个函数用于同步当前 CUDA 上下文中的所有启动操作，并称为cuCtxSynchronize（不带参数）：

cuCtxSynchronize = cuda.cuCtxSynchronize
cuCtxSynchronize.argtypes = []
cuCtxSynchronize.restype = int

下一个函数用于从加载的模块中检索内核函数句柄，以便我们可以将其启动到 GPU 上，这与 PyCUDA 的get_function方法完全对应，这一点我们已经看过很多次了。文档告诉我们原型是CUresult cuModuleGetFunction ( CUfunction* hfunc, CUmodule hmod, const char* name )。现在我们可以编写包装器：

cuModuleGetFunction = cuda.cuModuleGetFunction
 cuModuleGetFunction.argtypes = [c_void_p, c_void_p, c_char_p ]
 cuModuleGetFunction.restype = int

现在让我们为标准动态内存操作编写包装器；这些将是必要的，因为我们将不再有使用 PyCUDA gpuarray 对象的虚荣。这些实际上与我们之前使用过的 CUDA 运行时操作几乎相同；也就是说，cudaMalloc，cudaMemcpy和cudaFree：

cuMemAlloc = cuda.cuMemAlloc
cuMemAlloc.argtypes = [c_void_p, c_size_t]
cuMemAlloc.restype = int

cuMemcpyHtoD = cuda.cuMemcpyHtoD
cuMemcpyHtoD.argtypes = [c_void_p, c_void_p, c_size_t]
cuMemAlloc.restype = int

cuMemcpyDtoH = cuda.cuMemcpyDtoH
cuMemcpyDtoH.argtypes = [c_void_p, c_void_p, c_size_t]
cuMemcpyDtoH.restype = int

cuMemFree = cuda.cuMemFree
cuMemFree.argtypes = [c_void_p] 
cuMemFree.restype = int

现在，我们将为cuLaunchKernel函数编写一个包装器。当然，这是我们将用来在 GPU 上启动 CUDA 内核的函数，前提是我们已经初始化了 CUDA Driver API，设置了上下文，加载了一个模块，分配了内存并配置了输入，并且已经从加载的模块中提取了内核函数句柄。这个函数比其他函数复杂一些，所以我们将看一下原型：

CUresult cuLaunchKernel ( CUfunction f, unsigned int gridDimX, unsigned int gridDimY, unsigned int gridDimZ, unsigned int blockDimX, unsigned int blockDimY, unsigned int blockDimZ, unsigned int sharedMemBytes, CUstream hStream, void** kernelParams, void** extra )

第一个参数是我们要启动的内核函数的句柄，我们可以表示为c_void_p。六个gridDim和blockDim参数用于指示网格和块的维度。无符号整数sharedMemBytes用于指示在内核启动时为每个块分配多少字节的共享内存。CUstream hStream是一个可选参数，我们可以使用它来设置自定义流，或者如果希望使用默认流，则设置为 NULL（0），我们可以在 Ctypes 中表示为c_void_p。最后，kernelParams和extra参数用于设置内核的输入；这些有点复杂，所以现在只需知道我们也可以将这些表示为c_void_p：

cuLaunchKernel = cuda.cuLaunchKernel
cuLaunchKernel.argtypes = [c_void_p, c_uint, c_uint, c_uint, c_uint, c_uint, c_uint, c_uint, c_void_p, c_void_p, c_void_p]
cuLaunchKernel.restype = int

现在我们还有最后一个函数要为cuCtxDestroy编写一个包装器。我们在 CUDA 会话结束时使用它来销毁 GPU 上的上下文。唯一的输入是一个CUcontext对象，由c_void_p表示：

cuCtxDestroy = cuda.cuCtxDestroy
cuCtxDestroy.argtypes = [c_void_p]
cuCtxDestroy.restype = int

让我们把这个保存到cuda_driver.py文件中。我们现在已经完成了 Driver API 包装器模块！接下来，我们将看看如何仅使用我们的模块和 Mandelbrot PTX 加载一个 PTX 模块并启动一个内核。

这个示例也可以在本书的 GitHub 存储库中的cuda_driver.py文件中找到。

使用 CUDA Driver API

我们现在将翻译我们的小曼德布洛特生成程序，以便我们可以使用我们的包装库。让我们从适当的导入语句开始；注意我们如何将所有的包装器加载到当前命名空间中：

from __future__ import division
from time import time
import matplotlib
from matplotlib import pyplot as plt
import numpy as np
from cuda_driver import *

让我们把所有的 GPU 代码放到mandelbrot函数中，就像以前一样。我们将从使用cuInit初始化 CUDA Driver API 开始，然后检查系统上是否安装了至少一个 GPU，否则会引发异常：

def mandelbrot(breadth, low, high, max_iters, upper_bound):
 cuInit(0)
 cnt = c_int(0)
 cuDeviceGetCount(byref(cnt))
 if cnt.value == 0:
  raise Exception('No GPU device found!')

注意这里的byref：这是 Ctypes 中引用操作符(&)的等价物。我们现在将再次应用这个想法，记住设备句柄和 CUDA 上下文可以用 Ctypes 表示为c_int和c_void_p：

 cuDevice = c_int(0)
 cuDeviceGet(byref(cuDevice), 0)
 cuContext = c_void_p()
 cuCtxCreate(byref(cuContext), 0, cuDevice)

我们现在将加载我们的 PTX 模块，记得要使用c_char_p将文件名转换为 C 字符串：

 cuModule = c_void_p()
 cuModuleLoad(byref(cuModule), c_char_p('./mandelbrot.ptx'))

现在我们将在主机端设置晶格，以及一个名为graph的用于在主机端存储输出的全零 NumPy 数组。我们还将为晶格和图形输出在 GPU 上分配内存，然后使用cuMemcpyHtoD将晶格复制到 GPU 上：

 lattice = np.linspace(low, high, breadth, dtype=np.float32)
 lattice_c = lattice.ctypes.data_as(POINTER(c_float))
 lattice_gpu = c_void_p(0)
 graph = np.zeros(shape=(lattice.size, lattice.size), dtype=np.float32)
 cuMemAlloc(byref(lattice_gpu), c_size_t(lattice.size*sizeof(c_float)))
 graph_gpu = c_void_p(0)
 cuMemAlloc(byref(graph_gpu), c_size_t(lattice.size**2 * sizeof(c_float)))
 cuMemcpyHtoD(lattice_gpu, lattice_c, c_size_t(lattice.size*sizeof(c_float)))

现在我们将使用cuModuleGetFunction获取 Mandelbrot 内核的句柄，并设置一些输入：

 mandel_ker = c_void_p(0)
 cuModuleGetFunction(byref(mandel_ker), cuModule, c_char_p('mandelbrot_ker'))
 max_iters = c_int(max_iters)
 upper_bound_squared = c_float(upper_bound**2)
 lattice_size = c_int(lattice.size)

下一步有点复杂。在继续之前，我们必须了解参数是如何通过cuLaunchKernel传递到 CUDA 内核中的。让我们先看看 CUDA-C 中是如何工作的。

我们将输入参数在kernelParams中表达为一个void *值的数组，它们本身是指向我们希望插入内核的输入的指针。对于我们的曼德布洛特内核，它看起来像这样：

void * mandel_params [] = {&lattice_gpu, &graph_gpu, &max_iters, &upper_bound_squared, &lattice_size};

现在让我们看看如何在 Ctypes 中表达这一点，这并不是立即显而易见的。首先，让我们将所有的输入放入一个 Python 列表中，按正确的顺序：

mandel_args0 = [lattice_gpu, graph_gpu, max_iters, upper_bound_squared, lattice_size ]

现在我们需要每个值的指针，将其类型转换为void *类型。让我们使用 Ctypes 函数addressof来获取每个 Ctypes 变量的地址（类似于byref，只是不绑定到特定类型），然后将其转换为c_void_p。我们将这些值存储在另一个列表中：

mandel_args = [c_void_p(addressof(x)) for x in mandel_args0]

现在让我们使用 Ctypes 将这个 Python 列表转换成一个void *指针数组，就像这样：

 mandel_params = (c_void_p * len(mandel_args))(*mandel_args)

现在我们可以设置网格的大小，就像以前一样，并使用cuLaunchKernel启动我们的内核，使用这组参数。然后我们在之后同步上下文：

 gridsize = int(np.ceil(lattice.size**2 / 32))
 cuLaunchKernel(mandel_ker, gridsize, 1, 1, 32, 1, 1, 10000, None, mandel_params, None)
 cuCtxSynchronize()

我们现在将使用cuMemcpyDtoH将数据从 GPU 复制到我们的 NumPy 数组中，使用 NumPy 的array.ctypes.data成员，这是一个 C 指针，将允许我们直接从 C 中访问数组作为堆内存的一部分。我们将使用 Ctypes 的类型转换函数cast将其转换为c_void_p：

 cuMemcpyDtoH( cast(graph.ctypes.data, c_void_p), graph_gpu,  c_size_t(lattice.size**2 *sizeof(c_float)))

我们现在完成了！让我们释放在 GPU 上分配的数组，并通过销毁当前上下文来结束我们的 GPU 会话。然后我们将把图形 NumPy 数组返回给调用函数：

 cuMemFree(lattice_gpu)
 cuMemFree(graph_gpu)
 cuCtxDestroy(cuContext)
 return graph

现在我们可以像以前一样设置我们的main函数：

if __name__ == '__main__':
 t1 = time()
 mandel = mandelbrot(512,-2,2,256, 2)
 t2 = time()
 mandel_time = t2 - t1
 print 'It took %s seconds to calculate the Mandelbrot graph.' % mandel_time

 fig = plt.figure(1)
 plt.imshow(mandel, extent=(-2, 2, -2, 2))
 plt.show()

现在尝试运行这个函数，确保它产生与我们刚刚编写的其他曼德布洛特程序相同的输出。

恭喜你——你刚刚编写了一个直接接口到低级 CUDA Driver API，并成功使用它启动了一个内核！

这个程序也可以在本书的 GitHub 存储库中的目录下的mandelbrot_driver.py文件中找到。

总结

我们从简要概述 Python Ctypes 库开始了本章，该库用于直接与编译的二进制代码进行接口，特别是用 C/C++编写的动态库。然后，我们看了如何使用 CUDA-C 编写一个启动 CUDA 内核的基于 C 的包装器，然后使用这个包装器间接地从 Python 启动我们的 CUDA 内核，方法是使用 Ctypes 编写一个对这个函数的接口。然后，我们学习了如何将 CUDA 内核编译成 PTX 模块二进制文件，可以将其视为一个带有 CUDA 内核函数的 DLL，并看到如何使用 PyCUDA 加载 PTX 文件并启动预编译的内核。最后，我们编写了一系列 CUDA Driver API 的 Ctypes 包装器，并看到我们如何使用这些包装器执行基本的 GPU 操作，包括从 PTX 文件启动预编译的内核到 GPU 上。

我们现在将进入本书中可能是最技术性的一章：第十一章，《CUDA 性能优化》。在本章中，我们将学习关于 NVIDIA GPU 的一些技术细节，这将帮助我们提高应用程序的性能水平。

问题

假设您使用nvcc将包含主机和内核代码的单个.cu文件编译成 EXE 文件，还编译成 PTX 文件。哪个文件将包含主机函数，哪个文件将包含 GPU 代码？
如果我们使用 CUDA Driver API，为什么要销毁上下文？
在本章开始时，当我们首次看到如何使用 Ctypes 时，请注意在调用printf之前，我们必须将浮点值 3.14 强制转换为 Ctypes 的c_double对象。然而，在本章中我们可以看到许多不需要将类型转换为 Ctypes 的工作案例。你认为为什么printf在这里是一个例外呢？
假设您想要向我们的 Python CUDA Driver 接口模块添加功能以支持 CUDA 流。您将如何在 Ctypes 中表示单个流对象？
为什么在mandelbrot.cu中的函数要使用extern "C"？
再次查看mandelbrot_driver.py。为什么我们在 GPU 内存分配和主机/GPU 内存传输之后不使用cuCtxSynchronize函数，而只在单个内核调用之后使用？

第十一章：CUDA 中的性能优化

在这个倒数第二章中，我们将介绍一些相当高级的 CUDA 功能，可以用于低级性能优化。我们将首先学习动态并行性，它允许内核在 GPU 上启动和管理其他内核，并看看我们如何使用它直接在 GPU 上实现快速排序。我们将学习关于矢量化内存访问，可以用于增加从 GPU 全局内存读取时的内存访问加速。然后我们将看看如何使用 CUDA 原子操作，这些是线程安全的函数，可以在没有线程同步或mutex锁的情况下操作共享数据。我们将学习关于 Warp，它是 32 个或更少线程的基本块，在这些线程可以直接读取或写入彼此的变量，然后简要地涉足 PTX 汇编的世界。我们将通过直接在我们的 CUDA-C 代码中内联编写一些基本的 PTX 汇编来做到这一点，这本身将内联在我们的 Python 代码中！最后，我们将把所有这些小的低级调整结合到一个最终的例子中，我们将应用它们来制作一个快速的求和内核，并将其与 PyCUDA 的求和进行比较。

本章的学习成果如下：

CUDA 中的动态并行性
在 GPU 上使用动态并行性实现快速排序
使用矢量化类型加速设备内存访问
使用线程安全的 CUDA 原子操作
基本的 PTX 汇编
将所有这些概念应用于编写性能优化的求和内核

动态并行性

首先，我们将看一下动态并行性，这是 CUDA 中的一个功能，允许内核在 GPU 上启动和管理其他内核，而无需主机的任何交互或输入。这也使得许多通常在 GPU 上可用的主机端 CUDA-C 功能也可用于 GPU，例如设备内存分配/释放，设备到设备的内存复制，上下文范围的同步和流。

让我们从一个非常简单的例子开始。我们将创建一个小内核，覆盖N个线程，每个线程将向终端打印一条简短的消息，然后递归地启动另一个覆盖N-1个线程的内核。这个过程将持续，直到N达到 1。(当然，除了说明动态并行性如何工作之外，这个例子将毫无意义。)

让我们从 Python 中的import语句开始：

from __future__ import division
import numpy as np
from pycuda.compiler import DynamicSourceModule
import pycuda.autoinit

请注意，我们必须导入DynamicSourceModule而不是通常的SourceModule！这是因为动态并行性功能需要编译器设置特定的配置细节。否则，这将看起来和行为像通常的SourceModule操作。现在我们可以继续编写内核：

DynamicParallelismCode='''
__global__ void dynamic_hello_ker(int depth)
{
 printf("Hello from thread %d, recursion depth %d!\\n", threadIdx.x, depth);
 if (threadIdx.x == 0 && blockIdx.x == 0 && blockDim.x > 1)
 {
  printf("Launching a new kernel from depth %d .\\n", depth);
  printf("-----------------------------------------\\n");
  dynamic_hello_ker<<< 1, blockDim.x - 1 >>>(depth + 1);
 }
}'''

在这里需要注意的最重要的一点是：我们必须小心，只有一个线程启动下一个迭代的内核，使用一个良好放置的if语句来检查threadIdx和blockIdx的值。如果我们不这样做，那么每个线程将在每个深度迭代中启动远多于必要的内核实例。另外，注意我们可以以正常方式启动内核，使用通常的 CUDA-C 三重括号表示法——我们不必使用任何晦涩或低级命令来利用动态并行性。

在使用 CUDA 动态并行性功能时，一定要小心避免不必要的内核启动。这可以通过指定的线程启动下一个内核迭代来实现。

现在让我们结束这一切：

dp_mod = DynamicSourceModule(DynamicParallelismCode)
hello_ker = dp_mod.get_function('dynamic_hello_ker')
hello_ker(np.int32(0), grid=(1,1,1), block=(4,1,1))

现在我们可以运行前面的代码，这将给我们以下输出：

这个例子也可以在本书的 GitHub 存储库中的dynamic_hello.py文件中找到。

具有动态并行性的快速排序

现在让我们来看一个稍微有趣和实用的动态并行应用——快速排序算法。实际上，这是一个非常适合并行化的算法，我们将会看到。

让我们先简要回顾一下。快速排序是一种递归和原地排序算法，平均和最佳情况下的性能为O(N log N)，最坏情况下的性能为O(N²)。快速排序通过在未排序的数组中选择一个称为枢轴的任意点，然后将数组分成一个左数组（其中包含所有小于枢轴的点）、一个右数组（其中包含所有等于或大于枢轴的点），枢轴位于两个数组之间。如果一个或两个数组的长度大于 1，那么我们将在一个或两个子数组上再次递归调用快速排序，此时枢轴点已经处于最终位置。

在纯 Python 中，可以使用函数式编程一行代码实现快速排序：

qsort = lambda xs : [] if xs == [] else qsort(filter(lambda x: x < xs[-1] , xs[0:-1])) + [xs[-1]] + qsort(filter(lambda x: x >= xs[-1] , xs[0:-1]))

我们可以看到并行性将在快速排序递归调用右数组和左数组时发挥作用——我们可以看到这将从一个线程在初始大数组上操作开始，但当数组变得非常小时，应该有许多线程在它们上工作。在这里，我们实际上将通过在每个单个线程上启动所有内核来实现这一点！

让我们开始吧，并从导入语句开始。（我们将确保从标准随机模块中导入shuffle函数，以备后面的示例。）

from __future__ import division
import numpy as np
from pycuda.compiler import DynamicSourceModule
import pycuda.autoinit
from pycuda import gpuarray
from random import shuffle

现在我们将编写我们的快速排序内核。我们将为分区步骤编写一个device函数，它将接受一个整数指针、要分区的子数组的最低点和子数组的最高点。此函数还将使用此子数组的最高点作为枢轴。最终，在此函数完成后，它将返回枢轴的最终位置。

DynamicQuicksortCode='''
__device__ int partition(int * a, int lo, int hi)
{
 int i = lo;
 int pivot = a[hi];
 int temp;

 for (int k=lo; k<hi; k++)
 {
  if (a[k] < pivot)
  {
   temp = a[k];
   a[k] = a[i];
   a[i] = temp;
   i++;
  }
 }

 a[hi] = a[i];
 a[i] = pivot;

 return i;
}

现在我们可以编写实现此分区函数的内核，将其转换为并行快速排序。我们将使用 CUDA-C 约定来处理流，这是我们到目前为止还没有见过的：要在 CUDA-C 中的流s中启动内核k，我们使用k<<<grid, block, sharedMemBytesPerBlock, s>>>(...)。通过在这里使用两个流，我们可以确保它们是并行启动的。（考虑到我们不会使用共享内存，我们将把第三个启动参数设置为“0”。）流对象的创建和销毁应该是不言自明的：

__global__ void quicksort_ker(int *a, int lo, int hi)
{

 cudaStream_t s_left, s_right; 
 cudaStreamCreateWithFlags(&s_left, cudaStreamNonBlocking);
 cudaStreamCreateWithFlags(&s_right, cudaStreamNonBlocking);

 int mid = partition(a, lo, hi);

 if(mid - 1 - lo > 0)
   quicksort_ker<<< 1, 1, 0, s_left >>>(a, lo, mid - 1);
 if(hi - (mid + 1) > 0)
   quicksort_ker<<< 1, 1, 0, s_right >>>(a, mid + 1, hi);

 cudaStreamDestroy(s_left);
 cudaStreamDestroy(s_right);

}
'''

现在让我们随机洗牌一个包含 100 个整数的列表，并让我们的内核为我们进行排序。请注意我们如何在单个线程上启动内核：

qsort_mod = DynamicSourceModule(DynamicQuicksortCode)

qsort_ker = qsort_mod.get_function('quicksort_ker')

if __name__ == '__main__':
    a = range(100)
    shuffle(a)

    a = np.int32(a)

    d_a = gpuarray.to_gpu(a)

    print 'Unsorted array: %s' % a

    qsort_ker(d_a, np.int32(0), np.int32(a.size - 1), grid=(1,1,1), block=(1,1,1))

    a_sorted = list(d_a.get())

    print 'Sorted array: %s' % a_sorted

此程序也可以在本书的 GitHub 存储库中的dynamic_quicksort.py文件中找到。

矢量化数据类型和内存访问

现在我们将看一下 CUDA 的矢量化数据类型。这些是标准数据类型的矢量化版本，例如 int 或 double，它们可以存储多个值。32 位类型的矢量化版本的大小最多为 4（例如int2，int3，int4和float4），而 64 位变量只能矢量化为原始大小的两倍（例如double2和long2）。对于大小为 4 的矢量化变量，我们使用 C 的“struct”表示法访问每个单独的元素，成员为x，y，z和w，而对于 3 个成员的变量，我们使用x，y和z，对于 2 个成员的变量，我们只使用x和y。

现在这些可能看起来毫无意义，但这些数据类型可以用来提高从全局内存加载数组的性能。现在，让我们进行一个小测试，看看我们如何可以从整数数组中加载一些 int4 变量，以及从双精度数组中加载 double2 变量——我们将不得不使用 CUDA 的reinterpret_cast运算符来实现这一点：

from __future__ import division
import numpy as np
from pycuda.compiler import SourceModule
import pycuda.autoinit
from pycuda import gpuarray

VecCode='''
__global__ void vec_ker(int *ints, double *doubles) { 

 int4 f1, f2;

 f1 = *reinterpret_cast<int4*>(ints);
 f2 = *reinterpret_cast<int4*>(&ints[4]);

 printf("First int4: %d, %d, %d, %d\\n", f1.x, f1.y, f1.z, f1.w);
 printf("Second int4: %d, %d, %d, %d\\n", f2.x, f2.y, f2.z, f2.w);

 double2 d1, d2;

 d1 = *reinterpret_cast<double2*>(doubles);
 d2 = *reinterpret_cast<double2*>(&doubles[2]);

 printf("First double2: %f, %f\\n", d1.x, d1.y);
 printf("Second double2: %f, %f\\n", d2.x, d2.y);

}'''

请注意，我们必须使用dereference运算符*来设置矢量化变量，以及我们必须通过引用（&ints[4]，&doubles[2]）跳转到下一个地址来加载第二个int4和double2，使用数组上的引用运算符&：

这个例子也可以在本书的 GitHub 存储库中的vectorized_memory.py文件中找到。

线程安全的原子操作

我们现在将学习 CUDA 中的原子操作。原子操作是非常简单的、线程安全的操作，输出到单个全局数组元素或共享内存变量，否则可能会导致竞争条件。

让我们想一个例子。假设我们有一个内核，并且在某个时候我们设置了一个名为x的局部变量跨所有线程。然后我们想找到所有x中的最大值，然后将这个值设置为我们用__shared__ int x_largest声明的共享变量。我们可以通过在每个线程上调用atomicMax(&x_largest, x)来实现这一点。

让我们看一个原子操作的简单例子。我们将为两个实验编写一个小程序：

将变量设置为 0，然后为每个线程添加 1
找到所有线程中的最大线程 ID 值

让我们首先将tid整数设置为全局线程 ID，然后将全局add_out变量设置为 0。过去，我们会通过一个单独的线程使用if语句来改变变量，但现在我们可以使用atomicExch(add_out, 0)来跨所有线程进行操作。让我们导入并编写到这一点的内核：

from __future__ import division
import numpy as np
from pycuda.compiler import SourceModule
import pycuda.autoinit
from pycuda import gpuarray
import pycuda.driver as drv

AtomicCode='''
__global__ void atomic_ker(int *add_out, int *max_out) 
{

 int tid = blockIdx.x*blockDim.x + threadIdx.x;

 atomicExch(add_out, 0);

应该注意的是，虽然原子操作确实是线程安全的，但它们绝不保证所有线程将同时访问它们，它们可能会在不同的时间由不同的线程执行。这可能会有问题，因为我们将在下一步中修改add_out。这可能会导致add_out在一些线程部分修改后被重置。让我们进行块同步以防止这种情况发生：

 __syncthreads();

现在我们可以使用atomicAdd来为每个线程添加1到add_out，这将给我们总线程数：

 atomicAdd(add_out, 1);

现在让我们通过使用atomicMax来检查tid的最大值是多少。然后我们可以关闭我们的 CUDA 内核：

 atomicMax(max_out, tid);

}
'''

现在我们将添加测试代码；让我们尝试在 1 个包含 100 个线程的块上启动这个。我们这里只需要两个变量，所以我们将不得不分配一些大小为 1 的gpuarray对象。然后我们将打印输出：

atomic_mod = SourceModule(AtomicCode)
atomic_ker = atomic_mod.get_function('atomic_ker')

add_out = gpuarray.empty((1,), dtype=np.int32)
max_out = gpuarray.empty((1,), dtype=np.int32)

atomic_ker(add_out, max_out, grid=(1,1,1), block=(100,1,1))

print 'Atomic operations test:'
print 'add_out: %s' % add_out.get()[0]
print 'max_out: %s' % max_out.get()[0]

现在我们准备运行这个：

这个例子也可以在本书的 GitHub 存储库中的atomic.py文件中找到。

Warp shuffling

我们现在将看看所谓的warp shuffling。这是 CUDA 中的一个特性，允许存在于同一个 CUDA Warp 中的线程通过直接读写对方的寄存器（即它们的本地堆栈空间变量）进行通信，而无需使用shared变量或全局设备内存。Warp shuffling 实际上比其他两个选项快得多，更容易使用。这几乎听起来太好了，所以肯定有一个陷阱——确实，陷阱是这只在存在于同一个 CUDA Warp 上的线程之间起作用，这限制了对大小为 32 或更小的线程组的洗牌操作。另一个限制是我们只能使用 32 位或更小的数据类型。这意味着我们不能在 Warp 中洗牌 64 位长长整数或双精度浮点值。

只有 32 位（或更小）的数据类型可以与 CUDA Warp shuffling 一起使用！这意味着虽然我们可以使用整数、浮点数和字符，但不能使用双精度或长长整数！

在我们继续编码之前，让我们简要回顾一下 CUDA Warps。（在继续之前，您可能希望回顾第六章中名为“Warp Lockstep Property”的部分，“调试和分析您的 CUDA 代码”）。CUDA Warp 是 CUDA 中的最小执行单元，由 32 个线程或更少组成，运行在精确的 32 个 GPU 核心上。就像网格由块组成一样，块同样由一个或多个 Warps 组成，取决于块使用的线程数 - 如果一个块由 32 个线程组成，那么它将使用一个 Warp，如果它使用 96 个线程，它将由三个 Warps 组成。即使 Warp 的大小小于 32，它也被视为完整的 Warp：这意味着只有一个单个线程的块将使用 32 个核心。这也意味着 33 个线程的块将由两个 Warps 和 31 个核心组成。

要记住我们在第六章中看到的内容，“调试和分析您的 CUDA 代码”，Warp 具有所谓的Lockstep Property。这意味着 Warp 中的每个线程将完全并行地迭代每条指令，与 Warp 中的每个其他线程完全一致。也就是说，单个 Warp 中的每个线程将同时执行相同的指令，忽略任何不适用于特定线程的指令 - 这就是为什么要尽量避免单个 Warp 中线程之间的任何分歧。NVIDIA 将这种执行模型称为Single Instruction Multiple Thread，或SIMT。到目前为止，您应该明白为什么我们一直在文本中始终使用 32 个线程的块！

在我们开始之前，我们需要学习另一个术语 - Warp 中的lane是 Warp 内特定线程的唯一标识符，它将介于 0 和 31 之间。有时，这也被称为Lane ID。

让我们从一个简单的例子开始：我们将使用__shfl_xor命令在我们的 Warp 内的所有偶数和奇数编号的 Lanes（线程）之间交换特定变量的值。这实际上非常快速和容易做到，所以让我们编写我们的内核并查看一下：

from __future__ import division
import numpy as np
from pycuda.compiler import SourceModule
import pycuda.autoinit
from pycuda import gpuarray

ShflCode='''
__global__ void shfl_xor_ker(int *input, int * output) {

int temp = input[threadIdx.x];

temp = __shfl_xor (temp, 1, blockDim.x);

output[threadIdx.x] = temp;

}'''

除了__shfl_xor之外，这里的一切对我们来说都很熟悉。这是一个 CUDA 线程如何看待这个函数的方式：这个函数从当前线程接收temp的值作为输入。它对当前线程的二进制 Lane ID 执行一个XOR操作，这将是它的左邻居（如果该线程的 Lane 的最低有效位是二进制中的1）或右邻居（如果最低有效位是二进制中的“0”）。然后将当前线程的temp值发送到其邻居，同时检索邻居的 temp 值，这就是__shfl_xor。这将作为输出返回到temp中。然后我们在输出数组中设置值，这将交换我们的输入数组值。

现在让我们编写剩下的测试代码，然后检查输出：

shfl_mod = SourceModule(ShflCode)
shfl_ker = shfl_mod.get_function('shfl_xor_ker')

dinput = gpuarray.to_gpu(np.int32(range(32)))
doutout = gpuarray.empty_like(dinput)

shfl_ker(dinput, doutout, grid=(1,1,1), block=(32,1,1))

print 'input array: %s' % dinput.get()
print 'array after __shfl_xor: %s' % doutout.get()

上述代码的输出如下：

在我们继续之前，让我们做一个更多的 Warp-shuffling 示例 - 我们将实现一个操作，对 Warp 中所有线程的单个本地变量进行求和。让我们回顾一下第四章中的 Naive Parallel Sum 算法，“内核、线程、块和网格”，这个算法非常快速，但做出了一个天真的假设，即我们有与数据片段一样多的处理器 - 这是生活中为数不多的几种情况之一，我们实际上会有这么多处理器，假设我们正在处理大小为 32 或更小的数组。我们将使用__shfl_down函数在单个 Warp 中实现这一点。__shfl_down接受第一个参数中的线程变量，并通过第二个参数中指示的步数移动变量在线程之间，而第三个参数将指示 Warp 的总大小。

让我们立即实现这个。再次，如果您不熟悉 Naive Parallel Sum 或不记得为什么这应该起作用，请查看第四章，内核、线程、块和网格。我们将使用__shfl_down实现一个直接求和，然后在包括整数 0 到 31 的数组上运行这个求和。然后我们将与 NumPy 自己的sum函数进行比较，以确保正确性：

from __future__ import division
import numpy as np
from pycuda.compiler import SourceModule
import pycuda.autoinit
from pycuda import gpuarray

ShflSumCode='''
__global__ void shfl_sum_ker(int *input, int *out) {

 int temp = input[threadIdx.x];

 for (int i=1; i < 32; i *= 2)
     temp += __shfl_down (temp, i, 32);

 if (threadIdx.x == 0)
     *out = temp;

}'''

shfl_mod = SourceModule(ShflSumCode)
shfl_sum_ker = shfl_mod.get_function('shfl_sum_ker')

array_in = gpuarray.to_gpu(np.int32(range(32)))
out = gpuarray.empty((1,), dtype=np.int32)

shfl_sum_ker(array_in, out, grid=(1,1,1), block=(32,1,1))

print 'Input array: %s' % array_in.get()
print 'Summed value: %s' % out.get()[0]
print 'Does this match with Python''s sum? : %s' % (out.get()[0] == sum(array_in.get()) )

这将给我们以下输出：

本节中的示例也可在本书 GitHub 存储库的Chapter11目录下的shfl_sum.py和shfl_xor.py文件中找到。

内联 PTX 汇编

我们现在将初步了解编写 PTX（Parallel Thread eXecution）汇编语言，这是一种伪汇编语言，适用于所有 Nvidia GPU，反过来由即时（JIT）编译器编译为特定 GPU 的实际机器代码。虽然这显然不是用于日常使用，但如果必要，它将让我们在比 C 甚至更低的级别上工作。一个特定的用例是，如果没有其他源代码可用，您可以轻松地反汇编 CUDA 二进制文件（主机端可执行文件/库或 CUDA .cubin 二进制文件）并检查其 PTX 代码。这可以在 Windows 和 Linux 中使用cuobjdump.exe -ptx cuda_binary命令来完成。

如前所述，我们将只涵盖从 CUDA-C 内部使用 PTX 的一些基本用法，它具有特定的语法和用法，类似于在 GCC 中使用内联主机端汇编语言。让我们开始编写我们的 GPU 代码：

from __future__ import division
import numpy as np
from pycuda.compiler import SourceModule
import pycuda.autoinit
from pycuda import gpuarray

PtxCode='''

我们将通过将代码编写到单独的设备函数中进行几个小实验。让我们从一个简单的函数开始，将一个输入变量设置为零。（我们可以在 CUDA 中使用 C++的传址运算符&，我们将在device函数中使用它。）

__device__ void set_to_zero(int &x)
{
 asm("mov.s32 %0, 0;" : "=r"(x));
}

让我们在继续之前先分解一下。asm当然会告诉nvcc编译器我们将使用汇编，所以我们必须将代码放入引号中，以便正确处理。mov指令只是复制一个常量或其他值，并将其输入到寄存器中。（寄存器是 GPU 或 CPU 用于存储或操作值的最基本类型的芯片存储单元；这是 CUDA 中大多数本地变量的使用方式。）mov.s32中的.s32部分表示我们正在使用带符号的 32 位整数变量——PTX 汇编没有 C 中数据的类型，因此我们必须小心使用正确的特定操作。%0告诉nvcc使用与此处字符串的第0个参数对应的寄存器，并用逗号将此与mov的下一个输入分隔开，这是常量0。然后我们以分号结束汇编行，就像我们在 C 中一样，并用引号关闭这个汇编代码字符串。然后我们将使用冒号（而不是逗号！）来指示我们想要在我们的代码中使用的变量。"=r"表示两件事：=将告诉nvcc寄存器将被写入为输出，而r表示这应该被处理为 32 位整数数据类型。然后我们将要由汇编器处理的变量放在括号中，然后关闭asm，就像我们对任何 C 函数一样。

所有这些都是为了将一个变量的值设置为 0！现在，让我们编写一个小的设备函数，用于为我们添加两个浮点数：

__device__ void add_floats(float &out, float in1, float in2)
{
 asm("add.f32 %0, %1, %2 ;" : "=f"(out) : "f"(in1) , "f"(in2));
}

让我们停下来注意一些事情。首先，当然，我们使用add.f32来表示我们要将两个 32 位浮点值相加。我们还使用"=f"表示我们将写入一个寄存器，f表示我们将只从中读取。还要注意我们如何使用冒号来分隔write寄存器和only read寄存器，以供nvcc使用。

在继续之前，让我们再看一个简单的例子，即类似于 C 中的++运算符的函数，它将整数增加1：

__device__ void plusplus(int &x)
{
 asm("add.s32 %0, %0, 1;" : "+r"(x));
}

首先，请注意我们将“0th”参数用作输出和第一个输入。接下来，请注意我们使用的是+r而不是=r——+告诉nvcc这个寄存器在这个指令中将被读取和写入。

现在我们不会变得更复杂了，因为即使在汇编语言中编写一个简单的if语句也是相当复杂的。但是，让我们看一些更多的示例，这些示例在使用 CUDA Warps 时会很有用。让我们从一个小函数开始，这个函数将给出当前线程的 lane ID；这是非常有用的，实际上比使用 CUDA-C 更直接，因为 lane ID 实际上存储在一个称为%laneid的特殊寄存器中，我们无法在纯 C 中访问它。（请注意代码中我们使用了两个%符号，这将告诉nvcc直接在%laneid引用的汇编代码中使用%，而不是将其解释为asm命令的参数。）

__device__ int laneid()
{
 int id; 
 asm("mov.u32 %0, %%laneid; " : "=r"(id)); 
 return id;
}

现在让我们编写另外两个函数，这些函数对处理 CUDA Warps 将会很有用。请记住，只能使用洗牌命令在 Warp 之间传递 32 位变量。这意味着要在 Warp 上传递 64 位变量，我们必须将其拆分为两个 32 位变量，分别将这两个变量洗牌到另一个线程，然后将这两个 32 位值重新组合成原始的 64 位变量。对于将 64 位双精度拆分为两个 32 位整数，我们可以使用mov.b64命令，注意我们必须使用d来表示 64 位浮点双精度：

请注意我们在以下代码中使用volatile，这将确保这些命令在编译后按原样执行。我们这样做是因为有时编译器会对内联汇编代码进行自己的优化，但对于这样特别敏感的操作，我们希望按照原样执行。

__device__ void split64(double val, int & lo, int & hi)
{
 asm volatile("mov.b64 {%0, %1}, %2; ":"=r"(lo),"=r"(hi):"d"(val));
}

__device__ void combine64(double &val, int lo, int hi)
{
 asm volatile("mov.b64 %0, {%1, %2}; ":"=d"(val):"r"(lo),"r"(hi));
}

现在让我们编写一个简单的内核，用于测试我们编写的所有 PTX 汇编设备函数。然后我们将它启动在一个单个线程上，以便我们可以检查一切：

__global__ void ptx_test_ker() { 

 int x=123;

 printf("x is %d \\n", x);

 set_to_zero(x);

 printf("x is now %d \\n", x);

 plusplus(x);

 printf("x is now %d \\n", x);

 float f;

 add_floats(f, 1.11, 2.22 );

 printf("f is now %f \\n", f);

 printf("lane ID: %d \\n", laneid() );

 double orig = 3.1415;

 int t1, t2;

 split64(orig, t1, t2);

 double recon;

 combine64(recon, t1, t2);

 printf("Do split64 / combine64 work? : %s \\n", (orig == recon) ? "true" : "false"); 

}'''

ptx_mod = SourceModule(PtxCode)
ptx_test_ker = ptx_mod.get_function('ptx_test_ker')
ptx_test_ker(grid=(1,1,1), block=(1,1,1))

现在我们将运行前面的代码：

此示例也可在本书的 GitHub 存储库中的Chapter11目录下的ptx_assembly.py文件中找到。

性能优化的数组求和

对于本书的最后一个示例，我们将为给定的双精度数组制作一个标准的数组求和内核，但这一次我们将使用本章学到的所有技巧，使其尽可能快速。我们将检查我们求和内核的输出与 NumPy 的sum函数，然后我们将使用标准的 Python timeit函数运行一些测试，以比较我们的函数与 PyCUDA 自己的gpuarray对象的sum函数相比如何。

让我们开始导入所有必要的库，然后从一个laneid函数开始，类似于我们在上一节中使用的函数：

from __future__ import division
import numpy as np
from pycuda.compiler import SourceModule
import pycuda.autoinit
from pycuda import gpuarray
import pycuda.driver as drv
from timeit import timeit

SumCode='''
__device__ void __inline__ laneid(int & id)
{
 asm("mov.u32 %0, %%laneid; " : "=r"(id)); 
}

让我们注意一些事情——注意我们在设备函数的声明中放了一个新的内联语句。这将有效地将我们的函数变成一个宏，当我们从内核中调用这个函数时，会减少一点时间。另外，注意我们通过引用设置了id变量，而不是返回一个值——在这种情况下，实际上可能有两个整数寄存器应该被使用，并且应该有一个额外的复制命令。这样可以确保这种情况不会发生。

让我们以类似的方式编写其他设备函数。我们需要再编写两个设备函数，以便我们可以将 64 位双精度数拆分和合并为两个 32 位变量：

__device__ void __inline__ split64(double val, int & lo, int & hi)
{
 asm volatile("mov.b64 {%0, %1}, %2; ":"=r"(lo),"=r"(hi):"d"(val));
}

__device__ void __inline__ combine64(double &val, int lo, int hi)
{
 asm volatile("mov.b64 %0, {%1, %2}; ":"=d"(val):"r"(lo),"r"(hi));
}

让我们开始编写内核。我们将接收一个名为 input 的双精度数组，然后将整个总和输出到out，out应该初始化为0。我们将首先获取当前线程的 lane ID，并使用矢量化内存加载将两个值从全局内存加载到当前线程中：

__global__ void sum_ker(double *input, double *out) 
{

 int id;
 laneid(id);

 double2 vals = *reinterpret_cast<double2*> ( &input[(blockDim.x*blockIdx.x + threadIdx.x) * 2] );

现在让我们将双精度变量vals中的这些值求和到一个新的双精度变量sum_val中，它将跟踪本线程的所有求和。我们将创建两个 32 位整数s1和s2，我们将使用它们来分割这个值，并使用 Warp Shuffling 与其他线程共享一个temp变量来重构值：

 double sum_val = vals.x + vals.y;

 double temp;

 int s1, s2;

现在让我们再次在 Warp 上使用 Naive Parallel 求和，这将与在 Warp 上对 32 位整数求和相同，只是我们将在每次迭代中使用sum_val和temp上的split64和combine64PTX 函数：

 for (int i=1; i < 32; i *= 2)
 {

     // use PTX assembly to split
     split64(sum_val, s1, s2);

     // shuffle to transfer data
     s1 = __shfl_down (s1, i, 32);
     s2 = __shfl_down (s2, i, 32);

     // PTX assembly to combine
     combine64(temp, s1, s2);
     sum_val += temp;
 }

现在我们完成了，让我们让每个 Warp 的0th线程将其结束值添加到out，使用线程安全的atomicAdd：

 if (id == 0)
     atomicAdd(out, sum_val);

}'''

我们现在将编写我们的测试代码，使用timeit操作来测量我们的内核和 PyCUDA 对 10000232 个双精度数组进行 20 次迭代的平均时间：

sum_mod = SourceModule(SumCode)
sum_ker = sum_mod.get_function('sum_ker')

a = np.float64(np.random.randn(10000*2*32))
a_gpu = gpuarray.to_gpu(a)
out = gpuarray.zeros((1,), dtype=np.float64)

sum_ker(a_gpu, out, grid=(int(np.ceil(a.size/64)),1,1), block=(32,1,1))
drv.Context.synchronize()

print 'Does sum_ker produces the same value as NumPy\'s sum (according allclose)? : %s' % np.allclose(np.sum(a) , out.get()[0])

print 'Performing sum_ker / PyCUDA sum timing tests (20 each)...'

sum_ker_time = timeit('''from __main__ import sum_ker, a_gpu, out, np, drv \nsum_ker(a_gpu, out, grid=(int(np.ceil(a_gpu.size/64)),1,1), block=(32,1,1)) \ndrv.Context.synchronize()''', number=20)
pycuda_sum_time = timeit('''from __main__ import gpuarray, a_gpu, drv \ngpuarray.sum(a_gpu) \ndrv.Context.synchronize()''', number=20)

print 'sum_ker average time duration: %s, PyCUDA\'s gpuarray.sum average time duration: %s' % (sum_ker_time, pycuda_sum_time)
print '(Performance improvement of sum_ker over gpuarray.sum: %s )' % (pycuda_sum_time / sum_ker_time)

让我们从 IPython 中运行这个。确保你已经先运行了gpuarray.sum和sum_ker，以确保我们不会计算nvcc的编译时间：

因此，虽然求和通常相当无聊，但我们可以因为我们巧妙地利用硬件技巧来加速这样一个单调乏味的算法而感到兴奋。

此示例可在本书的 GitHub 存储库的Chapter11目录下的performance_sum_ker.py文件中找到。

总结

我们开始这一章时学习了动态并行性，这是一种允许我们直接在 GPU 上从其他内核启动和管理内核的范式。我们看到了我们如何可以使用这个来直接在 GPU 上实现快速排序算法。然后我们学习了 CUDA 中的矢量化数据类型，并看到了我们如何可以使用这些类型来加速从全局设备内存中读取数据。然后我们学习了 CUDA Warps，这是 GPU 上的小单位，每个 Warp 包含 32 个线程或更少，并且我们看到了单个 Warp 内的线程如何可以直接读取和写入彼此的寄存器，使用 Warp Shuffling。然后我们看了一下如何在 PTX 汇编中编写一些基本操作，包括确定 lane ID 和将 64 位变量分割为两个 32 位变量等导入操作。最后，我们通过编写一个新的性能优化求和内核来结束了这一章，该内核用于双精度数组，应用了本章学到的几乎所有技巧。我们看到，这实际上比双精度数组长度为 500,000 的标准 PyCUDA 求和更快。

我们已经完成了本书的所有技术章节！你应该为自己感到骄傲，因为现在你肯定是一个技艺高超的 GPU 程序员。现在我们将开始最后一章，在这一章中，我们将简要介绍一些不同的路径，可以帮助你应用和扩展你的 GPU 编程知识。

问题

在原子操作示例中，尝试在启动内核之前将网格大小从 1 更改为 2，同时保持总块大小为 100。如果这给出了add_out的错误输出（除了 200 以外的任何值），那么为什么是错误的，考虑到atomicExch是线程安全的呢？
在原子操作示例中，尝试移除__syncthreads，然后在原始参数的网格大小为 1，块大小为 100 的情况下运行内核。如果这给出了add_out的错误输出（除了 100 以外的任何值），那么为什么是错误的，考虑到atomicExch是线程安全的呢？
为什么我们不必使用__syncthreads来同步大小为 32 或更小的块？
我们发现sum_ker对于长度为 640,000（10000*2*32）的随机值数组比 PyCUDA 的求和操作快大约五倍。如果你尝试在这个数字的末尾加上一个零（也就是乘以 10），你会注意到性能下降到sum_ker只比 PyCUDA 的求和快大约 1.5 倍的程度。如果你在这个数字的末尾再加上一个零，你会注意到sum_ker只比 PyCUDA 的求和快 75%。你认为这是为什么？我们如何改进sum_ker以在更大的数组上更快？
哪种算法执行了更多的加法操作（计算 C +运算符的调用和将 atomicSum 视为单个操作）：sum_ker还是 PyCUDA 的sum？

第十二章：从这里往哪里走

这本书就像一次冒险的登山旅程一样……但现在，最终，我们已经到达了我们的徒步旅行的终点。我们现在站在介绍性 GPU 编程的山顶上，我们骄傲地回望我们的故乡串行编程村，微笑着想着我们旧的一维编程传统的天真，我们曾认为在 Unix 中fork一个进程就是我们对并行编程概念的全部理解。我们经历了许多险阻和危险才到达这一点，我们甚至可能犯了一些错误，比如在 Linux 中安装了一个损坏的 NVIDIA 驱动模块，或者在父母家度假时通过缓慢的 100k 连接下载了错误的 Visual Studio 版本。但这些挫折只是暂时的，留下的伤口变成了使我们更加强大对抗（GPU）自然力量的老茧。

然而，在我们的眼角，我们可以看到离我们站立的地方几米远处有两个木制标志；我们把目光从我们过去的小村庄移开，现在看着它们。第一个标志指向我们当前所面对的方向，上面只有一个词——过去。另一个指向相反的方向，也只有一个词——未来。我们转身朝着指向未来的方向走去，看到一个大而闪亮的大都市展现在我们面前，一直延伸到地平线，招手着我们。现在我们终于喘过气来，可以开始走向未来了…

在本章中，我们将介绍一些你现在可以继续学习 GPU 编程相关领域的选项。无论你是想要建立一个职业生涯，作为一个业余爱好者为了乐趣而做这个，作为一个工程学生为了课程而学习 GPU，作为一个程序员或工程师试图增强你的技术背景，还是作为一个学术科学家试图将 GPU 应用到一个研究项目中，你现在都有很多选择。就像我们比喻的大都市一样，迷失其中很容易，很难确定我们应该去哪里。我们希望在这最后一章中提供类似于简短的导游，为你提供一些下一步可以去的选项。

我们现在将在本章中看一下以下路径：

高级 CUDA 和 GPGPU 编程
图形
机器学习和计算机视觉
区块链技术

进一步学习 CUDA 和 GPGPU 编程的知识

你首先可以选择的是，当然是学习更多关于 CUDA 和通用 GPU（GPGPU）编程的知识。在这种情况下，你可能已经找到了一个很好的应用，并且想要编写更高级或优化的 CUDA 代码。你可能对它本身感兴趣，或者你想找一份 CUDA/GPU 程序员的工作。有了这本书提供的坚实的 GPU 编程基础，我们现在将看一些这个领域的高级主题，我们现在已经准备好学习了。

多 GPU 系统

首先想到的一个主要话题是学习如何为安装了多个 GPU 的系统编程。许多专业工作站和服务器都安装了多个 GPU，目的是处理需要不止一个，而是几个顶级 GPU 的数据。为此，存在一个称为多 GPU 编程的子领域。其中大部分工作集中在负载平衡上，即使用每个 GPU 的最大容量，确保没有一个 GPU 被过多的工作饱和，而另一个则未被充分利用。另一个话题是 Inter-GPU 通信，通常涉及一个 GPU 直接使用 CUDA 的 GPUDirect 点对点（P2P）内存访问将内存数组复制到另一个 GPU 或从另一个 GPU 复制。

NVIDIA 在这里提供了有关多 GPU 编程的简要介绍：www.nvidia.com/docs/IO/116711/sc11-multi-gpu.pdf。

集群计算和 MPI

另一个主题是集群计算，即编写程序，使其集体利用包含 GPU 的多台服务器。这些是服务器农场，它们在著名互联网公司（如 Facebook 和 Google）的数据处理设施中，以及政府和军方使用的科学超级计算设施中广泛存在。集群通常使用一种称为消息传递接口（MPI）的编程范式，它是与诸如 C++或 Fortran 之类的语言一起使用的接口，允许您编程连接到同一网络的许多计算机。

有关在 MPI 中使用 CUDA 的更多信息，请参阅此处：devblogs.nvidia.com/introduction-cuda-aware-mpi/。

OpenCL 和 PyOpenCL

CUDA 并不是唯一可以用来编程 GPU 的语言。CUDA 最主要的竞争对手是称为开放计算语言（Open Computing Language）或 OpenCL。CUDA 是一个封闭的专有系统，只能在 NVIDIA 硬件上运行，而 OpenCL 是一个由非营利性 Khronos Group 开发和支持的开放标准。OpenCL 不仅可以用于编程 NVIDIA GPU，还可以用于 AMD Radeon GPU 甚至 Intel HD GPU - 大多数主要技术公司都承诺在其产品中支持 OpenCL。此外，PyCUDA 的作者，UIUC 的 Andreas Kloeckner 教授，还编写了另一个出色（且免费）的 Python 库，名为 PyOpenCL，它提供了一个同样用户友好的 OpenCL 接口，几乎与 PyCUDA 具有相同的语法和概念。

有关 OpenCL 的信息由 NVIDIA 在这里提供：developer.nvidia.com/opencl。

Andreas Kloeckner 的网站上提供了有关免费 PyOpenCL 库的信息：

mathema.tician.de/software/pyopencl/。

图形

显然，GPU 中的 G 代表图形，而在本书中我们并没有看到太多关于图形的内容。尽管机器学习应用现在是 NVIDIA 的支柱产业，但一切都始于渲染出漂亮的图形。我们将在这里提供一些资源，让您开始学习，无论您是想开发视频游戏引擎、渲染 CGI 电影，还是开发 CAD 软件。CUDA 实际上可以与图形应用程序并驾齐驱，并且实际上已经用于专业软件，如 Adobe 的 Photoshop 和 After Effects，以及许多最近的视频游戏，如Mafia和Just Cause系列。我们将简要介绍一些您可能考虑从这里开始的主要 API。

OpenGL

OpenGL 是一个行业开放标准，自 90 年代初就存在。虽然在某些方面它显得有些陈旧，但它是一个稳定的 API，得到了广泛支持，如果你编写了一个使用它的程序，几乎可以保证在任何相对现代的 GPU 上都能运行。CUDA 示例文件夹实际上包含了许多示例，说明了 OpenGL 如何与 CUDA 接口（特别是在2_Graphics子目录中），因此感兴趣的读者可以考虑查看这些示例。（在 Windows 中，默认位置为C:\ProgramData\NVIDIA Corporation\CUDA Samples，在 Linux 中为/usr/local/cuda/samples。）

有关 OpenGL 的信息可以直接从 NVIDIA 这里获取：developer.nvidia.com/opengl。

PyCUDA 还提供了与 NVIDIA OpenGL 驱动程序的接口。有关信息，请访问此处：documen.tician.de/pycuda/gl.html。

DirectX 12

DirectX 12 是微软著名且得到良好支持的图形 API 的最新版本。虽然这是 Windows PC 和微软 Xbox 游戏机的专有，但这些系统显然拥有数亿用户的广泛安装基础。此外，除了 NVIDIA 显卡，Windows PC 还支持各种 GPU，并且 Visual Studio IDE 提供了很好的易用性。DirectX 12 实际上支持低级 GPGPU 编程类型的概念，并且可以利用多个 GPU。

微软的 DirectX 12 编程指南可以在这里找到：docs.microsoft.com/en-us/windows/desktop/direct3d12/directx-12-programming-guide。

Vulkan

Vulkan 可以被认为是 DirectX 12 的开放等效物，由 Khronos Group 开发，作为 OpenGL 的next-gen继任者。除了 Windows，Vulkan 也支持 macOS 和 Linux，以及索尼 PlayStation 4，任天堂 Switch 和 Xbox One 游戏机。Vulkan 具有与 DirectX 12 相同的许多功能，例如准 GPGPU 编程。Vulkan 对 DirectX 12 提供了一些严肃的竞争，例如 2016 年的《毁灭战士》重制版。

Khronos Group 在这里提供了《Vulkan 初学者指南》：www.khronos.org/blog/beginners-guide-to-vulkan。

机器学习和计算机视觉

当然，本章的重点是机器学习及其兄弟计算机视觉。不用说，机器学习（特别是深度神经网络和卷积神经网络的子领域）是如今让 NVIDIA 首席执行官黄仁勋有饭吃的东西。（好吧，我们承认这是本十年的轻描淡写……）如果你需要提醒为什么 GPU 在这个领域如此适用和有用，请再看一下第九章，实现深度神经网络。大量的并行计算和数学运算，以及用户友好的数学库，使 NVIDIA GPU 成为机器学习行业的硬件支柱。

基础知识

虽然现在你已经了解了低级 GPU 编程的许多复杂性，但你不会立即将这些知识应用到机器学习中。如果你在这个领域没有基本技能，比如如何对数据集进行基本统计分析，你真的应该停下来熟悉一下。斯坦福大学教授 Andrew Ng，谷歌 Brain 的创始人，在网上和 YouTube 上提供了许多免费的材料。Ng 教授的工作通常被认为是机器学习教育材料的金标准。

Ng 教授在这里提供了一门免费的机器学习入门课程：www.ml-class.org。

cuDNN

NVIDIA 提供了一个针对深度神经网络基元的优化 GPU 库，名为 cuDNN。这些基元包括前向传播、卷积、反向传播、激活函数（如 sigmoid、ReLU 和 tanh）和梯度下降。cuDNN 是大多数主流深度神经网络框架（如 Tensorflow）在 NVIDIA GPU 上的后端使用。这是 NVIDIA 免费提供的，但必须单独从 CUDA Toolkit 下载。

有关 cuDNN 的更多信息可以在这里找到：developer.nvidia.com/cudnn。

Tensorflow 和 Keras

Tensorflow 当然是谷歌著名的神经网络框架。这是一个免费的开源框架，可用于 Python 和 C++，自 2015 年以来一直向公众提供。

Tensorflow 的教程可以在 Google 这里找到：www.tensorflow.org/tutorials/。

Keras 是一个更高级的库，为 Tensorflow 提供了更用户友好的接口，最初由 Google Brain 的 Francois Chollet 编写。读者实际上可以考虑从 Keras 开始，然后再转向 Tensorflow。

有关 Keras 的信息在这里：keras.io/。

Chainer

Chainer 是另一个神经网络 API，由目前在日本东京大学攻读博士学位的 Seiya Tokui 开发。虽然它比 Tensorflow 更不为人知，但由于其令人难以置信的速度和效率而备受尊重。此外，读者可能会对 Chainer 特别感兴趣，因为最初是使用 PyCUDA 开发的。（后来改用了 CuPy，这是一个为了提供更类似于 NumPy 的接口而开发的 PyCUDA 分支。）

有关 Chainer 的信息在这里：chainer.org/。

OpenCV

自 2001 年以来，开源计算机视觉库（OpenCV）一直存在。该库提供了许多来自经典计算机视觉和图像处理的工具，在深度神经网络时代仍然非常有用。近年来，OpenCV 中的大多数算法都已移植到 CUDA，并且与 PyCUDA 接口非常容易。

有关 OpenCV 的信息在这里：opencv.org/。

区块链技术

最后但并非最不重要的是区块链技术。这是支持比特币和以太坊等加密货币的基础加密技术。这当然是一个非常新的领域，最初由比特币神秘创造者中本聪在 2008 年发表的白皮书中描述。在其发明后几乎立即应用了 GPU——生成货币单位归结为蛮力破解加密谜题，而 GPU 可以并行尝试蛮力破解比一般公众今天可用的任何其他硬件更多的组合。这个过程被称为挖矿。

有兴趣了解区块链技术的人建议阅读中本聪关于比特币的原始白皮书，网址在这里：bitcoin.org/bitcoin.pdf。

GUIMiner 是一个开源的基于 CUDA 的比特币矿工，网址在这里：guiminer.org/。

总结

在本章中，我们介绍了一些对于那些有兴趣进一步了解 GPU 编程背景的选项和路径，这超出了本书的范围。我们介绍的第一条路径是扩展您在纯 CUDA 和 GPGPU 编程方面的背景——您可以学习一些本书未涵盖的内容，包括使用多个 GPU 和网络集群的编程系统。我们还讨论了除 CUDA 之外的一些并行编程语言/API，如 MPI 和 OpenCL。接下来，我们讨论了一些对于有兴趣将 GPU 应用于渲染图形的人来说非常知名的 API，如 Vulkan 和 DirectX 12。然后，我们研究了机器学习，并深入了解了一些您应该具备的基本背景以及一些用于开发深度神经网络的主要框架。最后，我们简要介绍了区块链技术和基于 GPU 的加密货币挖矿。

作为作者，我想对每个人都说谢谢，感谢你们阅读到这里，到了结尾。GPU 编程是我遇到的最棘手的编程子领域之一，我希望我的文字能帮助你掌握基本要点。作为读者，你现在应该可以放纵自己，享用一块你能找到的最丰富、最高热量的巧克力蛋糕——只要知道你赚了它。（但只能吃一块！）

问题

使用谷歌或其他搜索引擎找到至少一个未在本章中介绍的 GPU 编程应用。
尝试找到至少一种编程语言或 API，可以用来编程 GPU，而这在本章中没有提到。
查找谷歌的新张量处理单元（TPU）芯片。这些与 GPU 有何不同？
你认为使用 Wi-Fi 还是有线以太网电缆将计算机连接成集群是更好的主意？

第十三章：评估

第一章，为什么要进行 GPU 编程？

前两个for循环遍历每个像素，它们的输出彼此不变；因此我们可以在这两个for循环上并行化。第三个for循环计算特定像素的最终值，这是固有的递归。
阿姆达尔定律没有考虑在 GPU 和主机之间传输内存所需的时间。
512 x 512 等于 262,144 像素。这意味着第一个 GPU 一次只能计算一半像素的输出，而第二个 GPU 可以一次计算所有像素；这意味着第二个 GPU 在这里将比第一个 GPU 快大约两倍。第三个 GPU 有足够的核心来一次计算所有像素，但正如我们在问题 1 中看到的那样，额外的核心在这里对我们没有用。因此，对于这个问题，第二个和第三个 GPU 的速度将是相同的。
将某个代码段通用地指定为与阿姆达尔定律相关的并行化存在一个问题，即假设这段代码的计算时间在处理器数量N非常大时接近 0。正如我们从上一个问题中看到的，情况并非如此。
首先，一致使用时间可能很麻烦，并且可能无法找出程序的瓶颈。其次，分析器可以告诉您从 Python 的角度来看，所有代码的确切计算时间，因此您可以确定某些库函数或操作系统的后台活动是否有问题，而不是您的代码。

第二章，设置 GPU 编程环境

不，CUDA 只支持 Nvidia GPU，不支持 Intel HD 或 AMD Radeon
本书只使用 Python 2.7 示例
设备管理器
lspci
free
.run

第三章，开始使用 PyCUDA

是的。
主机/设备之间的内存传输和编译时间。
你可以，但这将取决于你的 GPU 和 CPU 设置。
使用 C ?运算符来执行点对点和减少操作。
如果gpuarray对象超出范围，其析构函数将被调用，这将自动在 GPU 上释放它所代表的内存。
ReductionKernel可能执行多余的操作，这可能取决于底层 GPU 代码的结构。中性元素将确保没有值因这些多余的操作而改变。
我们应该将neutral设置为带符号 32 位整数的最小可能值。

第四章，内核，线程，块和网格

试试看。
并非所有线程都同时在 GPU 上运行。就像 CPU 在操作系统中在任务之间切换一样，GPU 的各个核心在内核的不同线程之间切换。
O(n/640 log n)，即 O(n log n)。
试试看。
实际上，CUDA 中没有内部网格级同步，只有块级（使用__syncthreads）。我们必须将单个块以上的任何内容与主机同步。
天真：129 次加法运算。高效工作：62 次加法运算。
同样，如果我们需要在大量块的网格上进行同步，我们不能使用__syncthreads。如果我们在主机上进行同步，我们还可以在每次迭代中启动更少的线程，从而为其他操作释放更多资源。
在天真的并行求和的情况下，我们可能只会处理一小部分数据点，这些数据点应该等于或少于 GPU 核心的总数，这些核心可能适合于块的最大大小（1032）；由于单个块可以在内部同步，我们应该这样做。只有当数据点的数量远远大于 GPU 上可用核心的数量时，我们才应该使用高效的工作算法。

第五章，流，事件，上下文和并发性

两者的性能都会提高；随着线程数量的增加，GPU 在两种情况下都达到了峰值利用率，减少了使用流所获得的收益。
是的，你可以异步启动任意数量的内核，并使用cudaDeviceSynchronize进行同步。
打开你的文本编辑器，试试看！
高标准差意味着 GPU 的使用不均匀，有时会过载 GPU，有时会过度利用它。低标准差意味着所有启动的操作通常都运行顺利。
i. 主机通常可以处理的并发线程远远少于 GPU。ii. 每个线程都需要自己的 CUDA 上下文。GPU 可能会因为上下文过多而不堪重负，因为每个上下文都有自己的内存空间，并且必须处理自己加载的可执行代码。

第六章，调试和分析 CUDA 代码

在 CUDA 中，内存分配会自动同步。
lockstep属性仅在大小为 32 或更小的单个块中有效。在这里，两个块将在没有任何lockstep的情况下正确分歧。
在这里也会发生同样的事情。这个 64 线程块实际上会被分成两个 32 线程的 warp。
Nvprof 可以计时单个内核启动、GPU 利用率和流使用；任何主机端的分析器只会看到 CUDA 主机函数的启动。
Printf 通常更容易用于规模较小、内联内核相对较短的项目。如果你编写了一个非常复杂的 CUDA 内核，有数千行代码，那么你可能会想使用 IDE 逐行步进和调试你的内核。
这告诉 CUDA 我们要使用哪个 GPU。
cudaDeviceSynchronize将确保相互依赖的内核启动和内存复制确实是同步的，并且它们不会在所有必要的操作完成之前启动。

第七章，使用 Scikit-CUDA 库

SBLAH 以 S 开头，因此该函数使用 32 位实数浮点数。ZBLEH 以 Z 开头，这意味着它使用 128 位复数浮点数。
提示：设置trans = cublas._CUBLAS_OP['T']
提示：使用 Scikit-CUDA 包装器进行点积，skcuda.cublas.cublasSdot
提示：建立在上一个问题的答案基础上。
你可以将 cuBLAS 操作放在一个 CUDA 流中，并使用该流的事件对象来精确测量 GPU 上的计算时间。
由于输入看起来对 cuFFT 来说很复杂，它将计算所有值作为 NumPy。
暗边是由于图像周围的零缓冲造成的。这可以通过在边缘上镜像图像来缓解，而不是使用零缓冲。

第八章，CUDA 设备函数库和 Thrust

试试看。(它实际上比你想象的更准确。)
一个应用：高斯分布可以用于向样本添加“白噪声”以增加机器学习中的数据集。
不，因为它们来自不同的种子，如果我们将它们连接在一起，这些列表可能会有很强的相关性。如果我们打算将它们连接在一起，我们应该使用相同种子的子序列。
试试看。
提示：记住矩阵乘法可以被看作是一系列矩阵-向量乘法，而矩阵-向量乘法可以被看作是一系列点积。
Operator()用于定义实际函数。

第九章，实现深度神经网络

一个问题可能是我们没有对训练输入进行归一化。另一个可能是训练速率太大。
使用较小的训练速率，一组权重可能会收敛得非常缓慢，或者根本不会收敛。
大的训练速率可能导致一组权重过度拟合特定的批处理值或该训练集。此外，它可能导致数值溢出/下溢，就像第一个问题一样。
Sigmoid。
Softmax。
更多更新。

第十章，使用已编译的 GPU 代码

只有 EXE 文件将包含主机函数，但 PTX 和 EXE 都将包含 GPU 代码。
cuCtxDestory。
printf带有任意输入参数。(尝试查找printf原型。)
使用 Ctypes 的c_void_p对象。
这将允许我们使用原始名称从 Ctypes 链接到函数。
CUDA 会自动同步设备内存分配和设备/主机之间的内存复制。

第十一章，CUDA 性能优化

atomicExch是线程安全的事实并不保证所有线程将同时执行此函数（因为在网格中不同的块可以在不同的时间执行）。
大小为 100 的块将在多个 warp 上执行，除非我们使用__syncthreads，否则块内部将不会同步。因此，atomicExch可能会被多次调用。
由于 warp 默认以锁步方式执行，并且大小为 32 或更小的块将使用单个 warp 执行，因此__syncthreads是不必要的。
我们在 warp 内使用天真的并行求和，但除此之外，我们使用atomicAdd进行了许多求和，就像我们使用串行求和一样。虽然 CUDA 自动并行化了许多这些atomicAdd调用，但我们可以通过实现高效的并行求和来减少所需的atomicAdd调用总数。
肯定是sum_ker。很明显 PyCUDA 的求和没有使用与我们相同的硬件技巧，因为我们在较小的数组上表现更好，但通过扩大尺寸，PyCUDA 版本更好的唯一解释是它执行的加法操作更少。

第十二章，下一步去哪里

两个例子：DNA 分析和物理模拟。
两个例子：OpenACC，Numba。
TPU 仅用于机器学习操作，缺少呈现图形所需的组件。
以太网。

标签：__,Python,编程,int,线程,CUDA,GPU,我们
From： https://www.cnblogs.com/apachecn/p/18140315

Python-GPU-编程实用指南（三）