Skip to content

Latest commit

 

History

History
727 lines (511 loc) · 50.1 KB

04.md

File metadata and controls

727 lines (511 loc) · 50.1 KB
Raw

四、内核、线程、块和网格

在本章中，我们将看到如何编写有效的CUDA 内核*。在 GPU 编程中，内核（我们可以与CUDA 内核或内核函数等术语互换使用）是一个并行函数，可以直接从主机（CPU）启动到设备上（GPU），而设备函数是只能从内核函数或其他设备函数调用的函数。（一般来说，设备函数的外观和行为与普通串行 C/C++函数类似，只是它们在 GPU 上运行，并从内核并行调用。）*

然后，我们将了解 CUDA 如何使用线程、块和网格的概念来抽象出 GPU 的一些底层技术细节（如内核、翘曲和流式多处理器，我们将在本书后面介绍），以及我们如何使用这些概念来减轻并行编程中的认知开销。我们将学习 CUDA 中使用全局和**共享***内存**的线程同步（块级和网格级）和线程内通信。最后，我们将深入研究如何在 GPU 上实现我们自己的并行前缀类型算法的技术细节（即，我们在上一章中介绍的扫描/减少类型函数），这使我们能够将本章中学习的所有原则付诸实践。

本章的学习成果如下：

  • 理解内核和设备函数之间的区别

  • 如何在 PyCUDA 中编译和启动内核，以及如何在内核中使用设备函数

  • 在启动内核的上下文中有效地使用线程、块和网格，以及如何在内核中使用threadIdx和blockIdx

  • 如何以及为什么同步内核中的线程，使用__syncthreads()同步单个块中的所有线程，使用主机同步整个块网格中的所有线程

  • 如何使用设备全局和共享内存进行线程内通信

  • 如何使用我们新获得的关于内核的所有知识来正确实现 GPU 版本的并行前缀和

技术要求

本章要求 Linux 或 Windows 10 PC 配备现代 NVIDIA GPU（2016 年以后），并安装所有必要的 GPU 驱动程序和 CUDA 工具包（9.0 以后）。还需要使用 PyCUDA 模块安装合适的 Python 2.7（如 Anaconda Python 2.7）。

本章代码也可在 GitHub 上获得，网址为：

https://github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA

For more information about the prerequisites, check the Preface of this book; for the software and hardware requirements, check the README section in https://github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA.

果仁

与上一章一样，我们将学习如何在 Python 代码中以内联 CUDA C 的形式编写 CUDA 内核函数，并使用 PyCUDA 将它们启动到 GPU 上。在最后一章中，我们使用 PyCUDA 提供的模板来编写属于特定设计模式的内核；相比之下，我们现在将看到如何从头开始编写我们自己的内核，这样我们就可以编写各种各样的内核，而这些内核可能不属于 PyCUDA 所涵盖的任何特定设计模式，这样我们就可以对内核进行更精细的控制。当然，这些收益将以编程更复杂为代价；我们尤其需要了解线程、块和网格及其在内核中的作用，以及如何同步内核正在执行的线程，以及如何在线程之间交换数据。

让我们从简单开始，尝试重新创建上一章中看到的一些元素操作，但这次没有使用ElementwiseKernel函数；我们现在将使用SourceModule函数。这是 PyCUDA 中一个非常强大的函数，它允许我们从头开始构建内核，所以像往常一样，最好从简单开始。

PyCUDA SourceModule 函数

我们将使用 PyCUDA 中的SourceModule函数将原始的内联 CUDA C 代码编译成可用的内核，我们可以从 Python 启动这些内核。我们应该注意到，SourceModule实际上将代码编译成CUDA 模块，这就像 Python 模块或 Windows DLL，只是它包含一个编译过的 CUDA 代码集合。这意味着，在我们真正启动 PyCUDA 的get_function之前，我们必须“拉出”一个对我们想要使用的内核的引用。让我们从一个基本示例开始，介绍如何将 CUDA 内核与SourceModule一起使用。

和前面一样，我们将从使最简单的内核函数之一成为可能开始，它将向量乘以标量。我们将从导入开始：

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

现在，我们可以立即开始编写内核：

ker = SourceModule("""
__global__ void scalar_multiply_kernel(float *outvec, float scalar, float *vec)
{
 int i = threadIdx.x;
 outvec[i] = scalar*vec[i];
}
""")

那么，让我们停下来，将其与ElementwiseKernel中的做法进行对比。首先，当我们在 CUDA C 中声明一个核函数时，我们在它前面加上__global__关键字。这将把函数区分为编译器的内核。我们总是将其声明为一个void函数，因为我们总是通过传递一个指向某个空内存块的指针来获得输出值，该内存块作为参数传入。我们可以像使用任何标准 C 函数一样声明参数：首先我们有outvec，它将是我们的输出缩放向量，当然是一个浮点数组指针。接下来，我们有scalar，仅用float表示；请注意，这不是指针！如果我们希望将简单的单例输入值传递给内核，我们可以不使用指针。最后，我们有了输入向量vec，它当然是另一个浮点数组指针。

Singleton input parameters to a kernel function can be passed in directly from the host without using pointers or allocated device memory. 

在继续测试之前，让我们先看看内核。我们记得，ElementwiseKernel通过 PyCUDA 为我们设置的一个值i自动并行多个 GPU 线程；每个单独线程的标识由threadIdx值给出，我们检索如下：int i = threadIdx.x;。

threadIdx is used to tell each individual thread its identity. This is usually used to determine an index for what values should be processed on the input and output data arrays. (This can also be used for assigning particular threads different tasks than others with standard C control flow statements such as if or switch.)

现在，我们准备像以前一样并行执行标量乘法：outvec[i] = scalar*vec[i];。

现在，让我们测试一下这段代码：我们首先必须从我们刚刚用SourceModule编译的 CUDA 模块中取出对我们编译的内核函数的引用。我们可以通过 Python 的get_function获得如下内核参考：

scalar_multiply_gpu = ker.get_function("scalar_multiply_kernel")

现在，我们必须在 GPU 上放置一些数据来实际测试我们的内核。让我们设置一个包含 512 个随机值的浮点数组，然后使用gpuarray.to_gpu函数将它们复制到 GPU 全局内存中的数组中。（我们将把这个随机向量乘以 GPU 和 CPU 上的标量，看看输出是否匹配。）我们还将使用gpuarray.empty_like函数为 GPU 的全局内存分配一块空内存：

testvec = np.random.randn(512).astype(np.float32)
testvec_gpu = gpuarray.to_gpu(testvec)
outvec_gpu = gpuarray.empty_like(testvec_gpu)

我们现在准备启动我们的内核。我们将标量值设置为2。（同样，由于标量是一个单例，我们不必将此值复制到 GPU，但是我们应该注意正确地键入它。）在这里，我们必须使用block和grid参数专门将线程数设置为512。我们现在准备推出：

scalar_multiply_gpu( outvec_gpu, np.float32(2), testvec_gpu, block=(512,1,1), grid=(1,1,1))

我们现在可以使用gpuarray输出对象中的get函数检查输出是否与预期输出匹配，并将其与 NumPy 的allclose函数的正确输出进行比较：

print "Does our kernel work correctly? : {}".format(np.allclose(outvec_gpu.get() , 2*testvec) )

（此示例的代码作为simple_scalar_multiply_kernel.py文件提供，位于存储库的4下。）

现在我们开始移除在上一章中学习的 PyCUDA 内核模板的训练轮，现在我们可以直接用纯 CUDA C 编写内核，并启动它以在 GPU 上使用特定数量的线程。然而，在继续使用内核之前，我们必须进一步了解 CUDA 如何将线程构造成抽象单元集合，这些抽象单元集合称为块和网格。

线程、块和网格

到目前为止，在这本书中，我们一直认为术语线程是理所当然的。让我们退一步，看看这到底意味着什么——线程是在 GPU 的单个内核上执行的一系列指令——内核和线程不应该被视为同义词！事实上，可以启动比 GPU 上的内核使用更多线程的内核。这是因为，与 Intel 芯片可能只有四个内核，但却在 Linux 或 Windows 中运行数百个进程和数千个线程的情况类似，操作系统的调度程序可以在这些任务之间快速切换，使它们看起来同时运行。GPU 以类似的方式处理线程，允许在数万个线程上进行无缝计算。

在 GPU 上以称为块的抽象单元执行多个线程。您应该还记得我们是如何在标量乘法内核中从threadIdx.x获得线程 ID 的；结尾有一个x，因为还有threadIdx.y和threadIdx.z。这是因为您可以在三个维度上索引块，而不仅仅是一个维度。我们为什么要这样做？让我们回忆一下第 1 章*为什么要进行 GPU 编程中有关 Mandelbrot 集计算的示例？*和第三章开始学习 PyCUDA。这是在二维平面上逐点计算的。因此，对于这样的算法，在二维上对线程进行索引可能更有意义。类似地，在某些情况下，在物理模拟中使用三维可能是有意义的，我们可能需要计算三维网格中移动粒子的位置。

块进一步以称为网格的抽象批执行，最好将其视为块的*块。*与块中的线程一样，我们可以使用blockIdx.x、blockIdx.y和blockIdx.z给出的常量值，在网格中为每个块建立多达三维的索引。让我们看一个例子来帮助我们理解这些概念；为了简单起见，我们这里只使用二维。

康威的人生游戏

生命的游戏（通常简称生命是一种细胞自动机模拟，由英国数学家约翰·康威在 1970 年发明。这听起来很复杂，但其实很简单，生活是一个零玩家的游戏，由细胞组成的二维二进制晶格组成，这些细胞要么被认为是活的要么死的。晶格通过以下一组规则进行迭代更新：

  • 任何少于两个活邻居的活细胞都会死亡
  • 任何有两个或三个邻居的活细胞都能存活
  • 任何有三个以上邻居的活细胞都会死亡
  • 只有三个邻居的死囚才会复活

这四条简单的规则产生了一个复杂的模拟，它具有有趣的数学特性，当设置动画时，在美学上也很好看。然而，由于晶格中有大量的单元，它可以运行得非常慢，并且当用纯串行 Python 编程时，通常会产生断断续续的动画。但是，这是可并行的，因为很明显，晶格中的每个单元都可以由单个 CUDA 线程管理。

我们现在将把 LIFE 实现为 CUDA 内核，并使用matplotlib.animation模块对其进行动画制作。这对我们来说很有趣，因为我们将能够在这里应用我们关于块和网格的新知识。

我们将首先包括以下适当的模块：

import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
import numpy as np
import matplotlib.pyplot as plt 
import matplotlib.animation as animation

现在，让我们通过SourceModule开始编写内核。我们将首先使用 C 语言的#define指令来设置一些常量和宏，这些常量和宏将在整个内核中使用。让我们看看我们将设置的前两个，_X和_Y：

ker = SourceModule("""
#define _X  ( threadIdx.x + blockIdx.x * blockDim.x )
#define _Y  ( threadIdx.y + blockIdx.y * blockDim.y )

让我们首先记住#define在这里是如何工作的，它将在编译时用定义的值（在这里的括号中）替换_X或_Y的任何文本，也就是说，它为我们创建宏。（出于个人风格，我通常在所有 C 宏前面加下划线。）

In C and C++, #define is used for creating macros. This means that #define doesn't create any function or set up a proper constant variables—it just allows us to write things shorthand in our code by swapping text out right before compilation time.

现在，让我们来讨论一下_X和_Y的具体含义，它们将是我们用于生活的二维晶格上单个 CUDA 线程单元的笛卡尔x和y值。我们将在由二维块组成的二维网格上启动内核，该二维块将对应于整个细胞晶格。我们必须使用线程和块常量来找到晶格上的笛卡尔点。让我们看一些图表来说明这一点。驻留在二维 CUDA 块中的线程可以可视化如下：

此时，您可能想知道为什么我们不在单个块上启动内核，所以我们可以将_X设置为threadIdx.x，将_Y设置为threadIdx.y，然后就可以完成了。这是由于 CUDA 对块大小的限制，目前只支持最多 1024 个线程的块。这意味着我们最多只能制作尺寸为 32 x 32 的细胞晶格，这将导致一个相当无聊的模拟，可能最好在 CPU 上完成，因此我们必须在网格上启动多个块。（我们当前区块的尺寸将由blockDim.x和blockDim.y给出，这将帮助我们确定目标x和y坐标，我们将看到。）

类似地，如前所述，我们可以通过blockIdx.x和blockIdx.y确定我们在二维网格中所处的块：

在我们稍微考虑一下数学之后，应该很清楚，_X应该被定义为(threadIdx.x + blockIdx.x * blockDim.x)，而_Y应该被定义为( threadIdx.y + blockIdx.y * blockDim.y )。（添加括号是为了在代码中插入宏时不干扰操作顺序。）现在，让我们继续定义剩余的宏：

#define _WIDTH  ( blockDim.x * gridDim.x )
#define _HEIGHT ( blockDim.y * gridDim.y  )

#define _XM(x)  ( (x + _WIDTH) % _WIDTH )
#define _YM(y)  ( (y + _HEIGHT) % _HEIGHT )

_WIDTH和_HEIGHT宏将分别为我们提供细胞晶格的宽度和高度，从图表中可以清楚地看到。让我们讨论一下_XM和_YM宏。在我们的生命的实现中，我们将把端点“环绕”到格的另一边，例如，我们将考虑 TythT12x x To13Te-值为 Tyt5，以及一个 T1 的 T15 值，即 T6 值的值，我们也将考虑一个 Ty16 t16×x Ty17 t-。0和 ay-_HEIGHT的值为0。为什么我们需要这个？当我们计算一个给定单元的活邻居的数量时，我们可能在某个边缘，而邻居可能是外部点。定义这些宏来调节我们的点会自动覆盖这一点。请注意，在使用 C 的模运算符之前，我们必须添加宽度或高度，这是因为，unlik 在 Python 中，C 中的模运算符可以为整数返回负值。

现在我们要定义最后一个宏。我们记得 PyCUDA 将二维数组作为一维指针传递到 cudac 中；二维数组以行方式从 Python 传递到一维 C 指针。这意味着我们必须将晶格上给定单元的给定笛卡尔（x、y）点转换为对应于晶格的指针内的一维点。在这里，我们可以这样做：

#define _INDEX(x,y)  ( _XM(x)  + _YM(y) * _WIDTH )

由于我们的细胞晶格是按行存储的，因此我们必须将y值乘以偏移到相应行对应点的宽度。我们现在终于可以开始实施生活了。让我们从生命中最重要的一部分开始，计算一个给定的细胞所拥有的活邻居的数量。我们将使用 CUDA设备功能实现此功能，如下所示：

__device__ int nbrs(int x, int y, int * in)
{
     return ( in[ _INDEX(x -1, y+1) ] + in[ _INDEX(x-1, y) ] + in[ _INDEX(x-1, y-1) ] \
                   + in[ _INDEX(x, y+1)] + in[_INDEX(x, y - 1)] \
                   + in[ _INDEX(x+1, y+1) ] + in[ _INDEX(x+1, y) ] + in[ _INDEX(x+1, y-1) ] );
}

设备函数是串行编写的 C 函数，由内核中的单个 CUDA 线程调用。也就是说，这个小函数将由内核中的多个线程并行调用。我们将细胞晶格表示为 32 位整数的集合（1 表示一个活细胞，0 表示一个死细胞），因此这将适用于我们的目的；我们只需要将当前单元格周围邻居的值相加。

A CUDA device function is a serial C function that is called by an individual CUDA thread from within a kernel. While these functions are serial in themselves, they can be run in parallel by multiple GPU threads. Device functions cannot by themselves by launched by a host computer onto a GPU, only kernels.

我们现在准备编写 LIFE 的内核实现。实际上，我们已经完成了大部分艰苦的工作—检查当前线程单元的邻居数，检查当前单元是活的还是死的，然后使用适当的 switch case 语句根据生命规则确定下一次迭代的状态。我们将在此内核中使用两个整数指针数组，一个将引用上一次迭代作为输入（lattice），另一个将引用我们将计算为输出的迭代（lattice_out：

__global__ void conway_ker(int * lattice_out, int * lattice  )
{
   // x, y are the appropriate values for the cell covered by this thread
   int x = _X, y = _Y;

   // count the number of neighbors around the current cell
   int n = nbrs(x, y, lattice);

    // if the current cell is alive, then determine if it lives or dies for the next generation.
    if ( lattice[_INDEX(x,y)] == 1)
       switch(n)
       {
          // if the cell is alive: it remains alive only if it has 2 or 3 neighbors.
          case 2:
          case 3: lattice_out[_INDEX(x,y)] = 1;
                  break;
          default: lattice_out[_INDEX(x,y)] = 0;                   
       }
    else if( lattice[_INDEX(x,y)] == 0 )
         switch(n)
         {
            // a dead cell comes to life only if it has 3 neighbors that are alive.
            case 3: lattice_out[_INDEX(x,y)] = 1;
                    break;
            default: lattice_out[_INDEX(x,y)] = 0;         
         }

}
""")

conway_ker = ker.get_function("conway_ker")

我们记得用三个括号关闭内联 CUDAC 段，然后用get_function获取对 CUDAC 内核的引用。由于内核只会更新晶格一次，因此我们将在 Python 中设置一个短函数，用于支付为动画更新晶格的所有开销：

def update_gpu(frameNum, img, newLattice_gpu, lattice_gpu, N):

frameNum参数只是 Matplotlib 的动画模块为更新函数所需的一个值，我们可以忽略该值，而img将是模块所需的细胞晶格的代表性图像，该图像将被迭代显示。

让我们关注其余三个参数-newLattice_gpu和lattice_gpu将是我们将保持持久性的 PyCUDA 数组，因为我们希望尽可能避免在 GPU 上重新分配内存块。lattice_gpu将是与内核中的lattice参数相对应的当前一代单元阵列，newLattice_gpu将是下一代晶格。N将指示晶格的高度和宽度（换句话说，我们将使用N x N晶格）。

我们使用适当的参数启动内核，并按如下方式设置块和网格大小：

    conway_ker(newLattice_gpu, lattice_gpu, grid=(N/32,N/32,1), block=(32,32,1) )

我们将使用(32, 32, 1)将块大小设置为 32 x 32；因为我们的细胞晶格只使用二维，所以我们可以将z维设置为一维。请记住，块被限制为 1024 个线程-32 x 32=1024，因此这将起作用。（请记住，这里关于 32 x 32 没有什么特别之处；如果需要，我们可以使用 16 x 64 或 10 x 10 等值，只要线程总数不超过 1024 个。）

The number of threads in a CUDA block is limited to a maximum of 1,024.

现在我们看一下这里的网格值，因为我们使用的是 32 维，所以应该很清楚N（在本例中）应该可以被 32 整除。这意味着在本例中，我们仅限于 64 x 64、96 x 96、128 x 128 和 1024 x 1024 等晶格。同样，如果我们想使用不同大小的晶格，那么我们必须改变块的尺寸。（如果这没有意义，那么请查看前面的图表，并回顾我们如何在内核中定义宽度和高度宏。）

现在，我们可以通过get()功能从 GPU 内存中获取最新生成的晶格，为动画设置图像数据。最后，我们使用 PyCUDA 切片运算符[:]将新晶格数据复制到当前数据中，该运算符将复制 GPU 上先前分配的内存，这样我们就不必重新分配：

    img.set_data(newLattice_gpu.get() )    
    lattice_gpu[:] = newLattice_gpu[:]

    return img

让我们设置一个大小为 256 x 256 的晶格。现在，我们将使用numpy.random模块中的选择函数为晶格设置初始状态。我们将用 1 和 0 随机填充一个整数的NxN图；通常，如果大约 25%的点为 1，其余为 0，我们可以生成一些有趣的晶格动画，因此我们将继续：

if __name__ == '__main__':
    # set lattice size
    N = 256

    lattice = np.int32( np.random.choice([1,0], N*N, p=[0.25, 0.75]).reshape(N, N) )
    lattice_gpu = gpuarray.to_gpu(lattice)

最后，我们可以使用适当的gpuarray函数在 GPU 上设置晶格，并相应地设置 Matplotlib 动画，如下所示：

lattice_gpu = gpuarray.to_gpu(lattice)
    lattice_gpu = gpuarray.to_gpu(lattice)
    newLattice_gpu = gpuarray.empty_like(lattice_gpu) 

    fig, ax = plt.subplots()
    img = ax.imshow(lattice_gpu.get(), interpolation='nearest')
    ani = animation.FuncAnimation(fig, update_gpu, fargs=(img, newLattice_gpu, lattice_gpu, N, ) , interval=0, frames=1000, save_count=1000) 

    plt.show()

我们现在可以运行我们的程序并欣赏这个节目（代码也可以作为 GitHub 存储库中4目录下的conway_gpu.py文件获得）：

线程同步与交互

现在我们将讨论 GPU 编程中的两个重要概念：线程同步和线程交互通信。有时，我们需要确保在继续任何进一步的计算之前，每个线程都达到了代码中相同的精确行；我们称之为线程同步。同步与线程间通信密切相关，即不同线程相互传递和读取输入；在这种情况下，我们通常希望在传递任何数据之前，确保所有线程在计算中的同一步骤对齐。我们将从学习 CUDA__syncthreads设备函数开始，该函数用于同步内核中的单个块。

使用 _syncthreads（）设备函数

在康威的生命游戏的前一个例子中，我们的内核每次由主机启动时只更新晶格一次。在本例中，同步已启动内核中的所有线程没有任何问题，因为我们只需处理 lattice 上一次可用的迭代

现在让我们假设我们想做一些稍微不同的事情，我们想重新编写内核，这样它就可以在给定的细胞格上执行一定数量的迭代，而不会被主机一次又一次地重新启动。这在一开始可能看起来很简单——一个简单的解决方案是只在内联的conway_ker内核中放置一个整数参数来指示迭代次数和一个for循环，进行一些额外的简单更改，然后完成它

然而，这引发了竞赛条件的问题；这就是多个线程对同一内存地址进行读写的问题，以及由此产生的问题。我们的老conway_ker内核通过使用两个内存数组来避免这个问题，一个是严格读取的，另一个是每次迭代都严格写入的。此外，由于内核只执行一次迭代，因此我们有效地使用主机来同步线程。

我们希望在 GPU 上进行多个完全同步的生命迭代；我们还希望为晶格使用单个内存数组。我们可以通过使用一个名为__syncthreads()的 CUDA 设备函数来避免竞争条件。此函数是一个块级同步屏障-这意味着在块内执行的每个线程在到达__syncthreads()实例时将停止，并等待同一块内的每个其他线程到达__syncthreads()的相同调用在测试之前，线程继续执行后续代码行。

 __syncthreads() can only synchronize threads within a single CUDA block, not all threads within a CUDA grid! 

现在让我们创建新内核；这将是对前世内核的修改，它将执行一定数量的迭代，然后停止。这意味着我们不会将其表示为动画，而是表示为静态图像，因此我们将在开始时加载相应的 Python 模块。（此代码也可在 GitHub 存储库的conway_gpu_syncthreads.py文件中找到）：

import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
import numpy as np
import matplotlib.pyplot as plt

现在，让我们再次设置计算生命的内核：

ker = SourceModule("""

当然，我们的 CUDAC 代码将出现在这里，这与以前基本相同。我们只需对内核进行一些更改。当然，我们可以保留设备功能，nbrs。在我们的声明中，我们将只使用一个数组来表示单元晶格。我们可以这样做，因为我们将使用适当的线程同步。我们还必须用整数表示迭代次数。我们将参数设置如下：

__global__ void conway_ker(int * lattice, int iters)
{

我们将一如既往地继续，只使用一个for循环进行迭代：

 int x = _X, y = _Y; 
 for (int i = 0; i < iters; i++)
 {
     int n = nbrs(x, y, lattice); 
     int cell_value;

让我们回忆一下之前，我们直接在数组中设置新的单元格晶格值。在这里，我们将保持cell_value变量中的值，直到块中的所有线程都同步。我们继续与之前类似的操作，使用__syncthreads阻止执行，直到为当前迭代确定所有新单元值，然后才在晶格阵列内设置值：

 if ( lattice[_INDEX(x,y)] == 1)
 switch(n)
 {
 // if the cell is alive: it remains alive only if it has 2 or 3 neighbors.
 case 2:
 case 3: cell_value = 1;
 break;
 default: cell_value = 0; 
 }
 else if( lattice[_INDEX(x,y)] == 0 )
 switch(n)
 {
 // a dead cell comes to life only if it has 3 neighbors that are alive.
 case 3: cell_value = 1;
 break;
 default: cell_value = 0; 
 } 
 __syncthreads();
 lattice[_INDEX(x,y)] = cell_value; 
 __syncthreads();
 } 
}
""")

我们现在将像以前一样启动内核并显示输出，在晶格上迭代 1000000 次。请注意，我们在网格中只使用一个块，大小为 32 x 32，这是由于每个块限制 1024 个线程。（再次强调，__syncthreads只对块中的所有线程有效，而不是对网格中的所有线程有效，这就是为什么我们将自己限制在单个块中的原因）：

conway_ker = ker.get_function("conway_ker")
if __name__ == '__main__':
 # set lattice size
 N = 32
 lattice = np.int32( np.random.choice([1,0], N*N, p=[0.25, 0.75]).reshape(N, N) )
 lattice_gpu = gpuarray.to_gpu(lattice)
 conway_ker(lattice_gpu, np.int32(1000000), grid=(1,1,1), block=(32,32,1))
 fig = plt.figure(1)
 plt.imshow(lattice_gpu.get())

当我们运行程序时，我们将得到如下所示的期望输出（这就是随机生命晶格在一百万次迭代后将收敛到的结果！）：

使用共享内存

我们可以从前面的示例中看到，内核中的线程可以使用 GPU 全局内存中的数组进行相互通信；虽然大多数操作都可以使用全局内存，但我们可以通过使用共享内存来加快速度。这是一种专门用于单个 CUDA 块内线程相互通信的内存类型；与全局内存相比，使用这种内存的优点是，对于纯线程间通信来说，速度要快得多。但是，与全局内存不同的是，存储在共享内存中的内存不能由主机直接访问，必须首先由内核本身将共享内存复制回全局内存。

让我们先退一步，然后再继续思考我们的意思。让我们看看我们刚刚看到的迭代生命内核中声明的一些变量。让我们首先看一下x和y，这两个整数保持特定线程单元的笛卡尔坐标。请记住，我们正在使用_X和_Y宏设置它们的值。（尽管进行了编译器优化，但我们希望将这些值存储在变量中以减少计算量，因为直接使用_X和_Y将在每次代码中引用这些宏时重新计算x和y值）:

 int x = _X, y = _Y;

我们注意到，对于每个线程，晶格中都会有一个唯一的笛卡尔点，对应于x和y。类似地，我们使用一个变量n（用int n = nbrs(x, y, lattice);声明）来表示特定单元周围的活邻居数。这是因为，当我们通常在 CUDA 中声明变量时，默认情况下它们是每个单独线程的本地变量。请注意，即使我们在一个线程（如int a[10];）中声明了一个数组，也会有一个大小为 10 的数组，它是每个线程的本地数组。

Local thread arrays (for example, a declaration of int a[10]; within the kernel) and pointers to global GPU memory (for example, a value passed as a kernel parameter of the form int * b) may look and act similarly, but are very different. For every thread in the kernel, there will be a separate a array that the other threads cannot read, yet there is a single b that will hold the same values and be equally accessible for all of the threads.

我们准备使用共享内存。这允许我们声明单个 CUDA 块中线程之间共享的变量和数组。这种内存比使用全局内存指针要快得多（我们到目前为止一直在使用全局内存指针），并且减少了指针情况下分配内存的开销。

假设我们想要一个大小为 10 的共享整数数组。我们声明如下-__shared__ int a[10] 。注意，我们不必将自己局限于数组；我们可以按如下方式创建共享单例变量：__shared__ int x。

让我们重写几行我们在上一小节中看到的迭代版 LIFE，以利用共享内存。首先，让我们将输入指针重命名为p_lattice，这样我们就可以在共享数组中使用这个变量名，并在代码中惰性地保留对“lattice”的所有引用。由于我们将坚持使用 32 x 32 单元晶格，因此我们按照如下方式设置了新的共享lattice阵列：

__global__ void conway_ker_shared(int * p_lattice, int iters)
{
 int x = _X, y = _Y;
 __shared__ int lattice[32*32];

我们现在必须将全局内存p_lattice数组中的所有值复制到lattice中。我们将以完全相同的方式索引共享数组，因此我们可以在这里使用旧的_INDEX宏。请注意，我们确保在复制后放置__syncthreads()，以确保在继续执行 LIFE 算法之前，所有对 lattice 的内存访问都已完成：

 lattice[_INDEX(x,y)] = p_lattice[_INDEX(x,y)];
 __syncthreads();

内核的其余部分与以前完全相同，只是我们必须从共享晶格复制回 GPU 阵列。我们按照以下步骤进行操作，然后关闭内联代码：

 __syncthreads();
 p_lattice[_INDEX(x,y)] = lattice[_INDEX(x,y)];
 __syncthreads();
} """)

我们现在可以像以前一样运行它，使用相同的精确测试代码。（这个例子可以在 GitHub 存储库的conway_gpu_syncthreads_shared.py中看到。）

并行前缀算法

现在，我们将使用我们对 CUDA 内核的新知识来实现并行前缀算法，也称为扫描设计模式。我们在上一章中已经看到了 PyCUDA 的InclusiveScanKernel和ReductionKernel函数形式的简单示例。我们现在将更详细地研究这个想法。

这种设计模式的中心动机是我们有一个二元运算符，也就是说，一个函数作用于两个输入值并给出一个输出值（例如-+、、（最大值）、（最小值）），以及元素集合，我们希望从中计算效率高。此外，我们假设我们的二元运算符是关联的——这意味着，对于任何三个元素，x、y和z，我们总是有：。

我们希望保留部分结果，即n-1子计算-。并行前缀算法的目的是高效地生成这组n和。在串行操作中产生这些n和通常需要*O（n）*时间，我们希望降低时间复杂度。

When the terms "parallel prefix" or "scan" are used, it usually means an algorithm that produces all of these n results, while "reduce"/"reduction" usually means an algorithm that only yields the single final result, . (This is the case with PyCUDA.)

实际上，并行前缀算法有几种变体，我们首先从最简单（也是最古老）的版本开始，它被称为 naive 并行前缀算法。

朴素并行前缀算法

朴素并行前缀算法是该算法的原始版本；该算法是“幼稚”的，因为它假设给定n个输入元素，进一步假设n是并矢（即对于某个正整数k，我们可以在n上并行运行该算法处理器（或n线程）。显然，这将对我们可能处理的集合的基数n施加强大的限制。然而，在满足这些条件的情况下，我们得到了一个很好的结果，它的计算时间复杂度仅为O（logn）。我们可以从算法的伪代码中看到这一点。这里，我们用表示输入值，用表示输出值：

input: x0, ..., xn-1 initialize:
for k=0 to n-1:
    yk := xk begin:
parfor i=0 to n-1 :
    for j=0 to log2(n):
        if i >= 2j :
            yi := yi  yi - 2<sup class="calibre53">j</sup> else:
            continue
        end if
    end for
end parfor
end
output: y0, ..., yn-1

现在，我们可以清楚地看到，这将需要*O（logn）渐近时间，因为外环在parfor上并联，内环在log2（n）*上并联。经过几分钟的思考，应该很容易看出，*yi*值将产生我们想要的输出。

现在让我们开始实施；这里，我们的二元运算符将只是加法。因为这个例子是说明性的，所以这个内核将严格超过 1024 个线程。

让我们只需设置标题，就可以开始编写内核了：

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

naive_ker = SourceModule("""
__global__ void naive_prefix(double *vec, double *out)
{
     __shared__ double sum_buf[1024]; 
     int tid = threadIdx.x; 
     sum_buf[tid] = vec[tid];

那么，让我们看看我们拥有的：我们将输入元素表示为一个双精度 GPU 数组，即double *vec，并用double *out表示输出值。我们声明一个共享内存sum_buf数组，用于计算输出。现在，让我们看一下算法本身的实现：

 int iter = 1;
 for (int i=0; i < 10; i++)
 {
     __syncthreads();
     if (tid >= iter )
     {
         sum_buf[tid] = sum_buf[tid] + sum_buf[tid - iter]; 
     } 
     iter *= 2;
 }
 __syncthreads();

当然，tid变量上没有隐含的parfor,，表示线程编号。我们还可以省略使用log2和2i，方法是从初始化为 1 的变量开始，然后在每次迭代 i 时迭代乘以 2。（请注意，如果我们想更具技术性，我们可以使用按位移位运算符来实现这一点。）我们将i的迭代次数限制为 10，因为210=1024。现在，我们将关闭新内核，如下所示：

 __syncthreads();
 out[tid] = sum_buf[tid];
 __syncthreads();

}
""")
naive_gpu = naive_ker.get_function("naive_prefix")

现在让我们看一下内核后面的测试代码：

if __name__ == '__main__':
 testvec = np.random.randn(1024).astype(np.float64)
 testvec_gpu = gpuarray.to_gpu(testvec)

 outvec_gpu = gpuarray.empty_like(testvec_gpu)
 naive_gpu( testvec_gpu , outvec_gpu, block=(1024,1,1), grid=(1,1,1))

 total_sum = sum( testvec)
 total_sum_gpu = outvec_gpu[-1].get()

 print "Does our kernel work correctly? : {}".format(np.allclose(total_sum_gpu , total_sum) )

我们只关心输出中的最终和，我们用outvec_gpu[-1].get()检索，回忆一下“-1”索引给出了 Python 中数组的最后一个成员。这将是vec中每个元素的总和；部分和在outvec_gpu的先前值中。（这个例子可以在 GitHub 存储库中的naive_prefix.py文件中看到。）

By its nature, the parallel prefix algorithm has to run over n threads, corresponding to a size-n array, where n is dyadic (again, this means that n is some power of 2). However, we can extend this algorithm to an arbitrary non-dyadic size assuming that our operator has a identity element (or equivalently, neutral element)—that is to say, that there is some value e so that for any x value, we have—.  In the case that our operator is + , the identity element is 0; in the case that it is , it is 1; all we do then is just pad the elements  with a series of e values so that we have the a dyadic cardinality of the new set .

包含前缀与独占前缀

让我们停下来，做一个非常微妙但非常重要的区别。到目前为止，我们一直关注的是将形式为的输入作为输出，生成形式为的和的数组。产生这样的输出的前缀算法称为包含；在包含前缀算法的情况下，每个索引处的对应元素包括在输出数组的同一索引的总和中。这与独占的前缀算法形成对比。排他前缀算法的不同之处在于，它同样采用形式的n输入值，并产生长度-n输出数组

这一点很重要，因为前缀算法的一些有效变体本质上是排他性的。我们将在下一小节中看到一个示例。

Note that the exclusive algorithm yields nearly the same output as the inclusive algorithm, only it is right-shifted and omits the final value. We can therefore trivially obtain the equivalent output from either algorithm, provided we keep a copy of .

一种高效并行前缀算法

在我们继续我们的新算法之前，我们将从两个角度来看这个朴素的算法。在理想情况下，计算时间复杂度为O（logn），但这只是当我们的数据集有足够数量的处理器时；当我们的数据集的基数（元素数）n远大于处理器的数量时，我们发现这变成了一个*O（n log n）*时间算法。

让我们定义一个与二进制运算符相关的新概念—由并行算法执行的工作，这里是在执行期间，该运算符在所有线程上的调用次数。类似地，跨度是线程在内核执行期间进行的调用次数；而整个算法的跨度与每个线程中最长的跨度相同，这将告诉我们总的执行时间。

我们寻求具体减少算法在所有线程中执行的工作量，而不仅仅是关注范围。在朴素前缀的情况下，当可用处理器数量不足时，所需的额外工作会花费更多的时间；这些额外的工作只会溢出到有限数量的处理器中。

我们将提出一种新的算法，该算法工作效率高，因此更适合有限数量的处理器。这包括两个不同的部分：上扫（或减少）阶段和下扫阶段。我们还应该注意，我们将看到的算法是一个排他性前缀算法。

上扫阶段类似于单个 reduce 操作，产生 reduce 算法给出的值，即；在这种情况下，我们保留实现最终结果所需的部分总和（。然后，下扫阶段将对这些部分和进行操作，并给出最终结果。让我们看一些伪代码，从上扫阶段开始。（下一小节将立即从伪代码深入到实现中。）

工作效率并行前缀（上扫阶段）

这是上扫的伪代码。（注意j变量上的parfor，这意味着这段代码可以在j索引的线程上并行化）：

input: x0, ..., xn-1initialize:
    for i = 0 to n - 1:
        yi := xi
begin:
for k=0 to log2(n) - 1:
    parfor j=0 to n - 1: 
        if j is divisible by 2k+1:
            yj+2<sup class="calibre53">k+1</sup>-1 = yj+2<sup class="calibre53">k</sup>-1  yj +2<sup class="calibre53">k+1</sup> -1else:
            continue
end
output: y0, ..., yn-1

工作效率并行前缀（下扫阶段）

现在让我们继续下扫，它将对上扫的输出进行操作：

input: x0, ..., xn-1 initialize:
    for i = 0 to n - 2:
        yi := xi
    yn-1 := 0
begin:
for k = log2(n) - 1 to 0:
    parfor j = 0 to n - 1: 
        if j is divisible by 2k+1:
            temp := yj+2<sup class="calibre53">k</sup>-1
            yj+2<sup class="calibre53">k</sup>-1 := yj+2<sup class="calibre53">k+1</sup>-1
            yj+2<sup class="calibre53">k+1</sup>-1 := yj+2<sup class="calibre53">k+1</sup>-1  temp
        else:
            continue
end
output: y0 , y1 , ..., yn-1

高效并行前缀实现

作为本章的重点，我们将编写一个该算法的实现，它可以在 1024 个以上任意大的数组上运行。这意味着这将在网格和块上运行；因此，我们必须使用主机进行同步；此外，这将要求我们为上扫和下扫阶段实现两个单独的内核，它们将在两个阶段中充当parfor循环，以及 Python 函数，它们将充当上扫和下扫的外部for循环。

让我们从一个向上扫描内核开始。因为我们将从主机迭代地重新启动这个内核，所以我们还需要一个指示当前迭代的参数（k。我们将使用两个数组进行计算以避免竞争条件-x（对于当前迭代）和x_old（对于上一次迭代）。我们声明内核如下：

up_ker = SourceModule("""
__global__ void up_ker(double *x, double *x_old, int k)
{

现在让我们设置tid变量，它将是网格中所有**块中所有线程的当前线程标识。我们使用了与我们之前看到的康威生命游戏的原始网格级实现相同的技巧：

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

我们现在将使用 C 位移位运算符直接从k生成 2k和 2k+1。我们现在将j设置为tid乘以_2k1——这将使我们能够删除伪代码中的“如果j可以被 2k+1整除”，这样我们就可以只启动所需数量的线程：

 int _2k = 1 << k;
 int _2k1 = 1 << (k+1);

 int j = tid* _2k1;

We can easily generate dyadic (power-of-2) integers in CUDA C with the left bit-wise shift operator (<<). Recall that the integer 1 (that is 20) is represented as 0001, 2 (21) is represented as 0010, 4 (22 ) is represented as 0100, and so on. We can therefore compute 2k with the 1 << k operation.

我们现在可以用一条直线运行上扫相位，注意到j确实可以通过其构造被 2k+1整除：

 x[j + _2k1 - 1] = x_old[j + _2k -1 ] + x_old[j + _2k1 - 1];
}
""")

我们已经完成了内核的编写！当然，这并不是一个全面的向上扫描的实现。我们必须用 Python 完成其余的工作。让我们获取内核并开始实现。这主要说明它完全遵循伪代码；我们应该记得，我们是通过使用[:]从x_gpu复制来更新x_old_gpu，这将保留内存分配，并且只复制新数据，而不是重新分配。还请注意，我们如何根据需要启动的线程数量设置块和网格大小，我们试图将块大小保持为大小 32 的倍数（这是我们在本文中的经验法则，我们将在第 11 章、CUDA中的性能优化中详细介绍为什么使用 32）。我们应该将from __future__ import division放在文件的开头，因为我们将使用 Python 3 风格的除法来计算块和内核大小。

需要提及的一个问题是，我们假设x的并矢长度为 32 或更大。如果您希望通过在数组中填充零来对其他大小的数组进行操作，则可以对其进行简单修改，但是：

up_gpu = up_ker.get_function("up_ker")

def up_sweep(x):
    x = np.float64(x)
    x_gpu = gpuarray.to_gpu(np.float64(x) )
    x_old_gpu = x_gpu.copy()
    for k in range( int(np.log2(x.size) ) ) : 
        num_threads = int(np.ceil( x.size / 2**(k+1)))
        grid_size = int(np.ceil(num_threads / 32))

        if grid_size > 1:
            block_size = 32
        else:
            block_size = num_threads

        up_gpu(x_gpu, x_old_gpu, np.int32(k) , block=(block_size,1,1), grid=(grid_size,1,1))
        x_old_gpu[:] = x_gpu[:]

    x_out = x_gpu.get()
    return(x_out)

现在我们开始写下横扫。再次，让我们从内核开始，它将具有伪代码内部parfor循环的功能。与之前类似，我们将使用两个数组，因此此处不需要使用伪码中的temp变量，我们再次使用位移位运算符来获得 2k和 2k+1的值。我们计算j的方法与之前类似：

down_ker = SourceModule("""
__global__ void down_ker(double *y, double *y_old, int k)
{
 int j = blockIdx.x*blockDim.x + threadIdx.x;

 int _2k = 1 << k;
 int _2k1 = 1 << (k+1);

 int j = tid*_2k1;

 y[j + _2k - 1 ] = y_old[j + _2k1 - 1];
 y[j + _2k1 - 1] = y_old[j + _2k1 - 1] + y_old[j + _2k - 1];
}
""")

down_gpu = down_ker.get_function("down_ker")

现在我们可以编写 Python 函数来迭代地启动内核，它对应于下扫阶段的外部for循环。这类似于上扫阶段的 Python 函数。查看伪代码的一个重要区别是，我们必须从外部for循环中的最大值迭代到最小值；我们可以使用 Python 的reversed函数来实现这一点。现在我们可以实施下扫阶段：

def down_sweep(y):
    y = np.float64(y)
    y[-1] = 0
    y_gpu = gpuarray.to_gpu(y)
    y_old_gpu = y_gpu.copy()
    for k in reversed(range(int(np.log2(y.size)))):
        num_threads = int(np.ceil( y.size / 2**(k+1)))
        grid_size = int(np.ceil(num_threads / 32))

        if grid_size > 1:
            block_size = 32
        else:
            block_size = num_threads

        down_gpu(y_gpu, y_old_gpu, np.int32(k), block=(block_size,1,1), grid=(grid_size,1,1))
        y_old_gpu[:] = y_gpu[:]
    y_out = y_gpu.get()
    return(y_out)

实施了上扫和下扫阶段后，我们的最后一项任务很容易完成：

def efficient_prefix(x):
        return(down_sweep(up_sweep(x)))

我们现在已经完全实现了主机同步版本的高效并行前缀算法！（此实现可在存储库中的work-efficient_prefix.py文件中找到，以及一些测试代码。）

总结

我们从康威的生命游戏的一个实现开始，它让我们了解了 CUDA 内核的多个线程是如何组织在块网格张量类型的结构中的。然后，我们通过 CUDA 函数__syncthreads()深入研究了块级同步，以及使用共享内存进行块级线程交互；我们还看到，单个块的线程数量有限，我们可以对其进行操作，因此，当我们创建在更大的网格中使用多个块的内核时，在使用这些功能时必须小心。

我们概述了并行前缀算法的理论，最后实现了一个简单的并行前缀算法，作为单个内核，可以在 1024 个大小受限的阵列上运行（与___syncthreads同步，并在内部执行for和parfor循环），通过一个工作效率高的并行前缀算法，该算法跨两个内核和三个 Python 函数实现，可以在任意大小的数组上运行，内核充当算法的内部parfor循环，Python 函数有效地充当外部for循环循环并同步内核启动。

问题

  1. 更改simple_scalar_multiply_kernel.py中的随机向量，使其长度为 10000，并修改内核定义中的i索引，使其可以以网格的形式在多个块上使用。通过将 block 和 grid 参数设置为类似于block=(100,1,1)和grid=(100,1,1)的值，看看您现在是否可以启动超过 10000 个线程的内核。
  2. 在上一个问题中，我们启动了一个内核，它同时使用 10000 个线程；截至 2018 年，没有超过 5000 核的 NVIDIA GPU。为什么这仍然有效并产生预期的结果？
  3. 考虑到对于大小为n的数据集，我们有n或更多处理器，naive 并行前缀算法的时间复杂度为 O（log n。假设我们在一个拥有 640 个内核的 GTX 1050 GPU 上使用一个简单的并行前缀算法。在n >> 640的情况下，渐近时间复杂度会变成什么？
  4. 修改naive_prefix.py操作任意大小的数组（可能是非并矢数组），仅以 1024 为界。
  5. __syncthreads()CUDA 设备功能仅同步单个块上的线程。如何跨网格中所有块中的所有线程进行同步？
  6. 通过本练习，您可以确信第二个前缀和算法确实比朴素的前缀和算法工作效率更高。假设我们有一个大小为 32 的数据集。在这种情况下，第一个和第二个算法所需的“加法”操作的确切数目是多少？
  7. 在高效工作并行前缀的实现中，我们使用 Python 函数迭代内核并同步结果。为什么我们不能在内核中放一个for循环，而仔细使用__syncthreads()呢？
  8. 为什么在 CUDA C 中处理自身同步的单个内核中实现 naive 并行前缀比使用内核和 Python 函数实现工作效率高的并行前缀并让主机处理同步更有意义？*