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

Python GPU 编程实用指南（二）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第六章：调试和分析您的 CUDA 代码

在本章中，我们将最终学习如何使用多种不同的方法和工具调试和分析我们的 GPU 代码。虽然我们可以使用 Spyder 和 PyCharm 等 IDE 轻松调试纯 Python 代码，但我们无法使用这些工具来调试实际的 GPU 代码，记住 GPU 代码本身是用 CUDA-C 编写的，PyCUDA 提供了一个接口。调试 CUDA 内核的第一种最简单的方法是使用printf语句，我们实际上可以直接在 CUDA 内核中调用它来打印到标准输出。我们将看到如何在 CUDA 的上下文中使用printf以及如何有效地应用它进行调试。

接下来，我们将填补 CUDA-C 编程中的一些空白，以便我们可以直接在 NVIDIA Nsight IDE 中编写 CUDA 程序，这将允许我们为我们一直在编写的一些代码创建 CUDA-C 的测试用例。我们将看看如何使用nvcc命令行编译 CUDA-C 程序，以及如何在 Nsight IDE 中进行编译。然后，我们将看看如何在 Nsight 中进行调试，并使用 Nsight 了解 CUDA lockstep 属性。最后，我们将概述 NVIDIA 命令行和 Visual Profilers 以对我们的代码进行分析。

本章的学习成果包括以下内容：

• 有效地使用printf作为 CUDA 内核的调试工具

• 在 Python 之外编写完整的 CUDA-C 程序，特别是用于创建调试的测试用例

• 使用nvcc编译器在命令行上编译 CUDA-C 程序

• 使用 NVIDIA Nsight IDE 开发和调试 CUDA 程序

• 了解 CUDA warp lockstep 属性以及为什么我们应该避免单个 CUDA warp 内的分支分歧

• 学会有效使用 NVIDIA 命令行和 Visual Profilers 进行 GPU 代码的调试

技术要求

本章需要一台安装了现代 NVIDIA GPU（2016 年以后）的 Linux 或 Windows 10 PC，并安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（9.0 及以上）。还需要一个合适的 Python 2.7 安装（如 Anaconda Python 2.7），并安装了 PyCUDA 模块。

本章的代码也可以在 GitHub 上找到：github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA。

有关先决条件的更多信息，请查看本书的前言，有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA中的 README。

在 CUDA 内核中使用 printf

也许会让人惊讶，但我们实际上可以直接从 CUDA 内核中将文本打印到标准输出；不仅如此，每个单独的线程都可以打印自己的输出。当我们调试内核时，这将特别方便，因为我们可能需要监视特定变量或代码中特定点的计算值，这也将使我们摆脱使用调试器逐步进行调试的束缚。从 CUDA 内核中打印输出的方法是使用 C/C++编程中最基本的函数，大多数人在编写他们的第一个 C 程序Hello world时会学到的函数：printf。当然，printf是将字符串打印到标准输出的标准函数，实际上在 C 编程语言中相当于 Python 的print函数。

现在让我们简要回顾一下如何在 CUDA 中使用printf。首先要记住的是，printf总是以字符串作为其第一个参数；因此，在 C 中打印"Hello world!"是用printf("Hello world!\n");来完成的。（当然，\n表示"新行"或"返回"，它将输出在终端上移到下一行。）printf还可以在我们想要直接在 C 中打印任何常量或变量的情况下，采用可变数量的参数：如果我们想要将123整数打印到输出，我们可以使用printf("%d", 123);（其中%d表示字符串后面跟着一个整数。）

类似地，我们使用%f，%e或%g来打印浮点值（其中%f是十进制表示法，%e是科学表示法，%g是最短的表示法，无论是十进制还是科学表示法）。我们甚至可以连续打印几个值，记得按正确的顺序放置这些指示符：printf("%d is a prime number, %f is close to pi, and %d is even.\n", 17, 3.14, 4);将在终端上打印"17 is a prime number, 3.14 is close to pi, and 4 is even."。

现在，在这本书的近一半时，我们终于要开始创建我们的第一个并行Hello world程序了！我们首先导入适当的模块到 Python 中，然后编写我们的内核。我们将首先打印每个单独线程的线程和网格标识（我们只会在一维块和网格中启动这个，所以我们只需要x值）：

ker = SourceModule('''
__global__ void hello_world_ker()
{
    printf("Hello world from thread %d, in block %d!\\n", threadIdx.x, blockIdx.x);

让我们停下来注意一下，我们写的是\\n而不是\n。这是因为 Python 中的三引号本身会将\n解释为"新行"，所以我们必须使用双反斜杠来表示我们是字面意思，以便将\n直接传递给 CUDA 编译器。

我们现在将打印有关块和网格维度的一些信息，但我们希望确保它在每个线程完成其初始printf命令后打印。我们可以通过放入__syncthreads();来确保每个单独的线程在执行第一个printf函数后同步。

现在，我们只想将块和网格的维度打印到终端一次；如果我们在这里放置printf语句，每个线程都会打印相同的信息。我们可以通过只有一个指定的线程打印输出来实现这一点；让我们选择第 0 个块的第 0 个线程，这是唯一保证存在的线程，无论我们选择的块和网格的维度如何。我们可以通过 C 的if语句来实现这一点：

 if(threadIdx.x == 0 && blockIdx.x == 0)
 {

我们现在将打印出我们的块和网格的维度，并关闭if语句，这将是我们的 CUDA 内核的结束：

 printf("-------------------------------------\\n");
 printf("This kernel was launched over a grid consisting of %d blocks,\\n", gridDim.x);
 printf("where each block has %d threads.\\n", blockDim.x);
 }
}
''')

我们现在将提取内核，然后在由两个块组成的网格上启动它，每个块有五个线程：

hello_ker = ker.get_function("hello_world_ker")
hello_ker( block=(5,1,1), grid=(2,1,1) )

让我们现在运行这个程序（该程序也可以在存储库中的hello-world_gpu.py的6下找到）：

使用 printf 进行调试

让我们通过一个例子来看看如何使用printf调试 CUDA 内核，然后再继续。这种方法并没有确切的科学依据，但通过经验可以学会这种技能。我们将从一个 CUDA 内核开始，用于矩阵乘法，但其中有几个错误。（鼓励读者随着我们的步骤阅读代码，该代码可在存储库中的6目录中的broken_matrix_ker.py文件中找到。）

在继续之前，让我们简要回顾一下矩阵乘法。假设我们有两个矩阵A和B，我们将它们相乘得到另一个相同大小的矩阵C，如下所示：。我们通过迭代所有元组来做到这一点，并将的值设置为A的第i行和B的第j列的点积：。

换句话说，我们将输出矩阵C中的每个i, j元素设置如下：

假设我们已经编写了一个内核，用于执行矩阵乘法，它接受表示输入矩阵的两个数组，一个额外的预先分配的浮点数组，输出将写入其中，并且一个表示每个矩阵的高度和宽度的整数（我们将假设所有矩阵都是相同大小和正方形的）。这些矩阵都将以一维的float *数组以行优先的一维布局表示。此外，这将被实现为每个 CUDA 线程处理输出矩阵中的单个行/列元组。

我们进行了一个小的测试案例，并将其与 CUDA 中矩阵乘法的输出进行了检查，对于两个 4 x 4 矩阵，它作为断言检查失败，如下所示：

test_a = np.float32( [xrange(1,5)] * 4 )
test_b = np.float32([xrange(14,10, -1)]*4 )
output_mat = np.matmul(test_a, test_b)

test_a_gpu = gpuarray.to_gpu(test_a)
test_b_gpu = gpuarray.to_gpu(test_b)
output_mat_gpu = gpuarray.empty_like(test_a_gpu)

matrix_ker(test_a_gpu, test_b_gpu, output_mat_gpu, np.int32(4), block=(2,2,1), grid=(2,2,1))

assert( np.allclose(output_mat_gpu.get(), output_mat) )

我们现在将运行这个程序，并且不出所料地得到以下输出：

现在让我们来看一下 CUDA C 代码，其中包括一个内核和一个设备函数：

ker = SourceModule('''
// row-column dot-product for matrix multiplication
__device__ float rowcol_dot(float *matrix_a, float *matrix_b, int row, int col, int N)
{
 float val = 0;

 for (int k=0; k < N; k++)
 {
     val += matrix_a[ row + k*N ] * matrix_b[ col*N + k];
 }
 return(val);
}

// matrix multiplication kernel that is parallelized over row/column tuples.

__global__ void matrix_mult_ker(float * matrix_a, float * matrix_b, float * output_matrix, int N)
{
 int row = blockIdx.x + threadIdx.x;
 int col = blockIdx.y + threadIdx.y;

 output_matrix[col + row*N] = rowcol_dot(matrix_a, matrix_b, col, row, N);
}
''')

我们的目标是在我们的 CUDA 代码中聪明地放置printf调用，以便我们可以监视内核和设备函数中的许多适当的值和变量；我们还应该确保在每个printf调用中打印出线程和块号。

让我们从内核的入口点开始。我们看到两个变量row和col，所以我们应该立即检查这些。让我们在设置它们之后立即放上一行代码（因为这是在两个维度上并行化的，我们应该打印threadIdx和blockIdx的x和y值）：

printf("threadIdx.x,y: %d,%d blockIdx.x,y: %d,%d -- row is %d, col is %d.\\n", threadIdx.x, threadIdx.y, blockIdx.x, blockIdx.y, row, col);

再次运行代码，我们得到了这个输出：

有两件事情立即引人注目：行和列元组有重复的值（每个单独的元组应该只表示一次），而且行和列的值从未超过两，而它们都应该达到三（因为这个单元测试使用 4 x 4 的矩阵）。这应该告诉我们，我们正在错误地计算行和列的值；确实，我们忘记了将blockIdx的值乘以blockDim的值来找到目标行/列值。我们按照以下方式修复这个问题：

int row = blockIdx.x*blockDim.x + threadIdx.x;
int col = blockIdx.y*blockDim.y + threadIdx.y;

然而，如果我们再次运行程序，我们仍然会得到一个断言错误。让我们保留原始的printf调用，这样我们就可以在继续进行的过程中监视这些值。我们看到在内核中有一个对设备函数rowcol_dot的调用，所以我们决定去看一下。让我们首先确保变量被正确传递到设备函数中，通过在开始处放置这个printf调用：

printf("threadIdx.x,y: %d,%d blockIdx.x,y: %d,%d -- row is %d, col is %d, N is %d.\\n", threadIdx.x, threadIdx.y, blockIdx.x, blockIdx.y, row, col, N);

当我们运行程序时，会有更多的行输出，但我们会看到一行说—threadIdx.x,y: 0,0 blockIdx.x,y: 1,0 -- row is 2, col is 0.，还有另一行说—threadIdx.x,y: 0,0 blockIdx.x,y: 1,0 -- row is 0, col is 2, N is 4。通过threadIdx和blockIdx的值，我们看到这是同一个块中的同一个线程，但row和col的值是颠倒的。实际上，当我们查看rowcol_dot设备函数的调用时，我们看到row和col确实是与设备函数声明中的相反。我们修复了这个问题，但当我们再次运行程序时，又出现了另一个断言错误。

让我们在设备函数中的for循环内放置另一个printf调用；这当然是点积，用于在矩阵A的行与矩阵B的列之间执行点积。我们将检查我们正在相乘的矩阵的值，以及k；我们还将只查看第一个线程的值，否则我们将得到一个不连贯的输出混乱。

if(threadIdx.x == 0 && threadIdx.y == 0 && blockIdx.x == 0 && blockIdx.y == 0)
            printf("Dot-product loop: k value is %d, matrix_a value is %f, matrix_b is %f.\\n", k, matrix_a[ row + k*N ], matrix_b[ col*N + k]);

在继续之前，让我们看一下为我们的单元测试设置的A和B矩阵的值：

我们看到，当我们在列之间切换时，两个矩阵都会变化，但在行之间变化时是恒定的。因此，根据矩阵乘法的性质，矩阵A的值应该在我们的for循环中随着k的变化而变化，而B的值应该保持恒定。让我们再次运行程序并检查相关的输出：

因此，看起来我们没有以正确的方式访问矩阵的元素；记住这些矩阵是以行方式存储的，我们修改索引，以便以正确的方式访问它们的值：

val += matrix_a[ row*N + k ] * matrix_b[ col + k*N];

再次运行程序将不会产生断言错误。恭喜，您刚刚使用唯一的printf调试了一个 CUDA 内核！

用 CUDA-C 填补空白

我们现在将介绍如何编写一个完整的 CUDA-C 程序的基础知识。我们将从小处开始，将我们刚刚在上一节中调试的小矩阵乘法测试程序的修复版本翻译成纯 CUDA-C 程序，然后使用 NVIDIA 的nvcc编译器从命令行编译成本机 Windows 或 Linux 可执行文件（我们将在下一节中看到如何使用 NVIDIA 的 Nsight IDE，所以现在我们只使用文本编辑器和命令行）。同样，鼓励读者在我们进行翻译时查看我们正在翻译的 Python 代码，该代码在存储库中作为matrix_ker.py文件可用。

现在，让我们打开我们最喜欢的文本编辑器，创建一个名为matrix_ker.cu的新文件。扩展名将表明这是一个 CUDA-C 程序，可以使用nvcc编译器进行编译。

CUDA-C 程序和库源代码的文件名总是使用.cu文件扩展名。

让我们从头开始——正如 Python 在程序开头使用import关键字导入库一样，我们回忆 C 语言使用#include。在继续之前，我们需要包含一些导入库。

让我们从这些开始：

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

让我们简要地思考一下我们需要这些的原因：cuda_runtime.h是一个头文件，其中包含了我们程序所需的所有特定 CUDA 数据类型、函数和结构的声明。我们需要为我们编写的任何纯 CUDA-C 程序包含这个头文件。stdio.h当然为我们提供了主机的所有标准 I/O 函数，如printf，我们需要stdlib.h来使用主机上的malloc和free动态内存分配函数。

请记住，始终在每个纯 CUDA-C 程序的开头放置#include <cuda_runtime.h>！

现在，在我们继续之前，我们记得我们最终将不得不检查我们的内核输出与正确的已知输出，就像我们在 NumPy 的allclose函数中所做的那样。不幸的是，在 C 中，我们没有像 Python 中的 NumPy 那样的标准或易于使用的数值数学库。往往，如果是一些简单的东西，编写自己的等效函数会更容易，就像在这种情况下一样。这意味着我们现在将明确地制作我们自己的等效于 NumPy 的allclose。我们将这样做：我们将使用 C 中的#define宏来设置一个名为_EPSILON的值，它将作为一个常数来指示输出和期望输出之间的最小值，以便被认为是相同的，我们还将设置一个名为_ABS的宏，它将告诉我们两个数字之间的绝对差异。我们这样做如下：

#define _EPSILON 0.001
#define _ABS(x) ( x > 0.0f ? x : -x )

现在我们可以创建我们自己的allclose版本。这将接受两个浮点指针和一个整数值len，我们循环遍历两个数组并检查它们：如果任何点的差异超过_EPSILON，我们返回-1，否则我们返回 0，表示这两个数组确实匹配。

我们注意到一件事：由于我们使用 CUDA-C，我们在函数定义之前加上__host__，以表明这个函数是打算在 CPU 上运行而不是在 GPU 上运行的：

__host__ int allclose(float *A, float *B, int len)
{

  int returnval = 0;

  for (int i = 0; i < len; i++)
  {
    if ( _ABS(A[i] - B[i]) > _EPSILON )
    {
      returnval = -1;
      break;
    }
  }

  return(returnval);
}

现在我们可以将设备和内核函数剪切并粘贴到这里，就像它们在我们的 Python 版本中出现的那样：

__device__ float rowcol_dot(float *matrix_a, float *matrix_b, int row, int col, int N)
{
  float val = 0;

  for (int k=0; k < N; k++)
  {
        val += matrix_a[ row*N + k ] * matrix_b[ col + k*N];
  }

  return(val);
}

__global__ void matrix_mult_ker(float * matrix_a, float * matrix_b, float * output_matrix, int N)
{

    int row = blockIdx.x*blockDim.x + threadIdx.x;
    int col = blockIdx.y*blockDim.y + threadIdx.y;

  output_matrix[col + row*N] = rowcol_dot(matrix_a, matrix_b, row, col, N);
}

再次，与__host__相比，注意 CUDA 设备函数前面有__device__，而 CUDA 内核前面有__global__。

现在，就像在任何 C 程序中一样，我们需要编写main函数，它将在主机上运行，我们将在其中设置我们的测试案例，并从中显式地启动我们的 CUDA 内核到 GPU 上。再次，与普通的 C 相比，我们将明确指定这也要在 CPU 上运行，并使用__host__：

__host__ int main()
{

我们将要做的第一件事是选择和初始化我们的 GPU。我们可以使用cudaSetDevice来这样做：

cudaSetDevice(0);

cudaSetDevice(0)将选择默认的 GPU。如果您的系统中安装了多个 GPU，您可以选择并使用它们，而不是使用cudaSetDevice(1)，cudaSetDevice(2)等。

现在我们将设置N，就像在 Python 中一样，以指示矩阵的高度/宽度。由于我们的测试案例将只包括 4 x 4 矩阵，我们将其设置为4。由于我们将使用动态分配的数组和指针，我们还必须设置一个值，以指示我们的测试矩阵将需要的字节数。矩阵将由N x N浮点数组成，我们可以使用 C 中的sizeof关键字确定浮点数所需的字节数：

int N = 4;
int num_bytes = sizeof(float)*N*N;

现在我们设置我们的测试矩阵如下；这些将与我们在 Python 测试程序中看到的test_a和test_b矩阵完全对应（请注意我们如何使用h_前缀来表示这些数组存储在主机上，而不是设备上）：

 float h_A[] = { 1.0, 2.0, 3.0, 4.0, \
                 1.0, 2.0, 3.0, 4.0, \
                 1.0, 2.0, 3.0, 4.0, \
                 1.0, 2.0, 3.0, 4.0 };

 float h_B[] = { 14.0, 13.0, 12.0, 11.0, \
                 14.0, 13.0, 12.0, 11.0, \
                 14.0, 13.0, 12.0, 11.0, \
                 14.0, 13.0, 12.0, 11.0 };

现在我们设置另一个数组，它将指示先前测试矩阵的矩阵乘法的预期输出。我们将不得不明确计算这一点，并将这些值放入我们的 C 代码中。最终，我们将在程序结束时将其与 GPU 输出进行比较，但让我们先设置它并将其放在一边：

float h_AxB[] = { 140.0, 130.0, 120.0, 110.0, \
                 140.0, 130.0, 120.0, 110.0, \
                 140.0, 130.0, 120.0, 110.0, \
                 140.0, 130.0, 120.0, 110.0 };

现在我们声明一些指针，这些指针将存在于 GPU 上的数组，并且我们将复制h_A和h_B的值并指向 GPU 的输出。请注意，我们只是使用标准的浮点指针。还要注意前缀d_ - 这是另一个标准的 CUDA-C 约定，表示这些将存在于设备上：

float * d_A;
float * d_B;
float * d_output;

现在，我们将使用cudaMalloc在设备上为d_A和d_B分配一些内存，这几乎与 C 中的malloc相同；这就是 PyCUDA gpuarray函数（如empty或to_gpu）在本书中无形地调用我们分配 GPU 上的内存数组的方式：

cudaMalloc((float **) &d_A, num_bytes);
cudaMalloc((float **) &d_B, num_bytes);

让我们思考一下这是如何工作的：在 C 函数中，我们可以通过在变量前加上一个取地址运算符（&）来获取变量的地址；如果有一个整数x，我们可以用&x来获取它的地址。&x将是一个指向整数的指针，因此它的类型将是int *。我们可以使用这个来将参数的值设置到 C 函数中，而不仅仅使用纯返回值。

由于cudaMalloc通过参数设置指针而不是使用返回值（与常规的malloc相反），我们必须使用取地址运算符，它将是一个指向指针的指针，因为它是一个指向浮点指针的指针，所以我们必须使用括号显式地转换这个值，因为cudaMalloc可以分配任何类型的数组。最后，在第二个参数中，我们必须指示在 GPU 上分配多少字节；我们之前已经设置了num_bytes，以便它是我们需要保存由浮点数组成的 4 x 4 矩阵所需的字节数，所以我们将其插入并继续。

现在我们可以使用两次cudaMemcpy函数将h_A和h_B的值分别复制到d_A和d_B中：

cudaMemcpy(d_A, h_A, num_bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, num_bytes, cudaMemcpyHostToDevice);

cudaMemcpy总是以目标指针作为第一个参数，源指针作为第二个参数，要复制的字节数作为第三个参数，并且还有一个最后的参数。最后一个参数将指示我们是使用cudaMemcpyHostToDevice从主机到 GPU 进行复制，使用cudaMemcpyDeviceToHost从 GPU 到主机进行复制，还是在 GPU 上的两个数组之间进行复制cudaMemcpyDeviceToDevice。

我们现在将分配一个数组来保存我们在 GPU 上进行矩阵乘法的输出，使用cudaMalloc的另一个调用：

cudaMalloc((float **) &d_output, num_bytes);

最后，当我们想要检查内核的输出时，我们将在主机上设置一些存储 GPU 输出的内存。让我们设置一个常规的 C 浮点指针，并使用malloc分配内存，就像我们通常做的那样：

float * h_output;
h_output = (float *) malloc(num_bytes);

现在，我们几乎准备好启动我们的内核。CUDA 使用一个名为dim3的数据结构来指示内核启动的块和网格大小；我们将设置这些，因为我们想要一个 2 x 2 维度的网格和也是 2 x 2 维度的块：

dim3 block(2,2,1);
dim3 grid(2,2,1);

现在我们准备启动我们的内核；我们使用三角形括号来指示 CUDA-C 编译器内核应该启动的块和网格大小：

matrix_mult_ker <<< grid, block >>> (d_A, d_B, d_output, N);

现在，当然，在我们可以将内核的输出复制回主机之前，我们必须确保内核已经执行完毕。我们通过调用cudaDeviceSynchronize来做到这一点，这将阻止主机向 GPU 发出更多命令，直到内核执行完毕：

cudaDeviceSynchronize();

现在我们可以将内核的输出复制到我们在主机上分配的数组中：

cudaMemcpy(h_output, d_output, num_bytes, cudaMemcpyDeviceToHost);

再次，我们同步：

cudaDeviceSynchronize();

在检查输出之前，我们意识到我们不再需要在 GPU 上分配的任何数组。我们通过在每个数组上调用cudaFree来释放这些内存：

cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_output);

我们已经完成了 GPU 的工作，所以我们调用cudaDeviceReset：

cudaDeviceReset();

现在，我们最终检查我们在主机上复制的输出，使用我们在本章开头编写的allclose函数。如果实际输出与预期输出不匹配，我们打印一个错误并返回-1，否则，我们打印匹配并返回 0。然后我们在程序的main函数上放一个闭括号：

if (allclose(h_AxB, h_output, N*N) < 0)
 {
     printf("Error! Output of kernel does not match expected output.\n");
     free(h_output);
     return(-1);
 }
 else
 {
     printf("Success! Output of kernel matches expected output.\n");
     free(h_output);
     return(0);
 }
}

请注意，由于我们已经为h_output分配了内存，在这两种情况下，我们最后调用了标准的 C free 函数。

现在保存我们的文件，并从命令行编译成 Windows 或 Linux 可执行文件，使用nvcc matrix_ker.cu -o matrix_ker。这应该输出一个二进制可执行文件，matrix_ker.exe（在 Windows 中）或matrix_ker（在 Linux 中）。让我们尝试编译和运行它：

恭喜，您刚刚创建了您的第一个纯 CUDA-C 程序！（此示例在存储库中作为matrix_ker.cu在7下可用。）

使用 Nsight IDE 进行 CUDA-C 开发和调试

现在让我们学习如何使用 Nsight IDE 开发 CUDA-C 程序。我们将看到如何导入我们刚刚编写的程序，并在 Nsight 内部进行编译和调试。请注意，由于在 Windows 下它实际上是 Visual Studio IDE 的插件，在 Linux 下是 Eclipse IDE 的插件，因此 Windows 和 Linux 版本的 Nsight 之间存在差异。我们将在接下来的两个子部分中涵盖两者；如果不适用于您的操作系统，请随意跳过。

在 Windows 中使用 Visual Studio 的 Nsight

打开 Visual Studio，点击文件，然后选择新建|项目....会弹出一个窗口，您可以在其中设置项目类型：选择 NVIDIA 下拉项，然后选择 CUDA 9.2：

给项目取一个合适的名称，然后点击确定。在解决方案资源管理器窗口中应该会出现一个项目，其中包含一个简单的预制 CUDA 测试程序，由一个源文件kernel.cu组成，其中包含一个简单的并行加法内核和测试代码。如果您想查看这是否编译和运行，请点击顶部标有本地 Windows 调试器的绿色向右箭头。一个终端应该弹出，显示内核的一些文本输出，然后立即关闭。

如果您在从 Visual Studio 运行后关闭基于 Windows 终端的应用程序时遇到问题，请尝试在主函数的末尾添加getchar();，这将使终端保持打开状态，直到您按下一个键。（或者，您也可以在程序的末尾使用调试器断点。）

现在，让我们添加刚刚编写的 CUDA-C 程序。在解决方案资源管理器窗口中，右键单击kernel.cu，然后单击kernel.cu上的删除。现在，右键单击项目名称，选择添加，然后选择现有项目。现在我们可以选择一个现有文件，找到matrix_ker.cu的路径，并将其添加到项目中。点击 IDE 顶部标有本地 Windows 调试器的绿色箭头，程序应该会在 Windows 终端中再次编译和运行。这就是我们可以在 Visual Studio 中设置和编译完整的 CUDA 程序的全部步骤。

现在让我们看看如何调试我们的 CUDA 内核。让我们首先在代码的入口点matrix_mult_ker处添加一个断点，我们在那里设置了row和col的值。我们可以通过在窗口的行号左侧的灰色列上单击来添加此断点；每个我们添加的断点都应该在那里显示一个红点。（您可以忽略 Visual Studio 编辑器可能在您的代码下方放置的任何红色波浪线；这是因为 CUDA 不是 Visual Studio 的本地语言）：

现在我们可以开始调试。从顶部菜单中选择 Nsight 下拉菜单，然后选择开始 CUDA 调试。这里可能有两个选项，开始 CUDA 调试（Next-Gen）和开始 CUDA 调试（Legacy）。无论选择哪一个都可以，但是根据您的 GPU，可能会在 Next-Gen 上遇到问题；在这种情况下，请选择 Legacy。

您的程序应该启动，并且调试器应该在我们刚刚设置的内核中的断点处停止。让我们按F10跳过这一行，现在看看row变量是否被正确设置。让我们在变量资源管理器中的本地窗口中查看：

通过检查threadIdx和blockIdx的值，我们可以看到我们当前位于网格中的第一个块中的第一个线程；row设置为0，这确实对应于正确的值。现在，让我们检查一些不同线程的row值。为了做到这一点，我们必须在 IDE 中切换线程焦点；我们可以通过单击上面的 Nsight 下拉菜单，然后选择 Windows|CUDA Debug Focus...来实现这一点。应该会出现一个新菜单，允许您选择一个新的线程和块。在菜单中将线程从 0, 0, 0 更改为 1, 0, 0，然后单击确定：

当您再次检查变量时，您应该看到为此线程设置了正确的row值：

简而言之，这就是您在 Visual Studio 中使用 Nsight 进行调试的方法。我们现在已经掌握了如何在 Windows 中使用 Nsight/Visual Studio 调试 CUDA 程序的基础知识，我们可以像调试常规 Windows 程序一样使用所有常规约定（设置断点，启动调试器，继续/恢复，跳过，步入和步出）。主要的区别在于您必须知道如何在 CUDA 线程和块之间切换以检查变量，否则它基本上是一样的。

在 Linux 中使用 Nsight 与 Eclipse

现在我们将看到如何在 Linux 中使用 Nsight。您可以从桌面上选择它打开 Nsight，也可以使用nsight命令从命令行运行它。Nsight IDE 将打开。从 IDE 顶部，单击文件，然后从下拉菜单中选择新建...，然后选择新建 CUDA C/C++项目。将出现一个新窗口，在这里选择 CUDA Runtime 项目。给项目取一个合适的名字，然后点击下一步。您将被提示提供进一步的设置选项，但默认设置对我们的目的来说现在可以工作得很好。（请确保注意这里第三和第四屏幕中源文件和项目路径的位置。）您将进入最终屏幕，在这里您可以按完成来创建项目：

最后，您将在项目视图中看到您的新项目和一些占位代码；从 CUDA 9.2 开始，这将包括一个倒数内核示例。

现在我们可以导入我们的代码。您可以使用 Nsight 中的编辑器删除默认源文件中的所有代码并剪切粘贴，或者您可以手动从项目的源目录中删除文件，手动将matrix_ker.cu文件复制到源目录中，然后选择刷新 Nsight 中的源目录视图，然后按F5。现在可以使用Ctrl + B构建项目，并使用F11运行它。我们程序的输出应该出现在 IDE 的 Console 子窗口中，如下所示：

现在，我们可以在 CUDA 代码中设置断点；让我们在内核的入口点设置一个断点，那里设置了row值。我们将光标放在 Eclipse 编辑器中的该行上，然后按Ctrl + Shift + B进行设置。

现在，我们可以通过按F11（或单击 bug 图标）开始调试。程序应该在main函数的开头暂停，所以按F8继续到第一个断点。您应该在 IDE 中看到我们的 CUDA 内核中的第一行被箭头指向。让我们通过按F6跳过当前行，确保row已经设置。

现在，我们可以轻松地在 CUDA 网格中切换不同的线程和块，以检查它们当前持有的值：从 IDE 顶部，单击窗口下拉菜单，然后单击显示视图，然后选择 CUDA。应该会打开一个显示当前运行内核的窗口，从这里您可以看到此内核正在运行的所有块的列表。

点击第一个，从这里你将能够看到块内运行的所有单个线程：

现在，我们可以通过单击“变量”选项卡来查看与第一个块中的第一个线程对应的变量，这里，row 应该是 0，正如我们所期望的：

现在，我们可以通过再次转到 CUDA 选项卡，选择适当的线程并切换回来，来检查不同线程的值。让我们留在同一个块中，但这次选择线程（1,0,0），再次检查 row 的值：

我们看到 row 的值现在是 1，正如我们所期望的。

现在，我们已经掌握了如何从 Nisight/Eclipse 在 Linux 中调试 CUDA 程序的基础知识，我们可以像调试其他 IDE 中的常规 Linux 程序一样使用所有常规约定（设置断点、启动调试器、继续/恢复、步进、步入和步出）。主要的区别在于我们必须知道如何在 CUDA 线程和块之间切换以检查变量，否则，它基本上是一样的。

使用 Nisight 来理解 CUDA 中的 warp lockstep 属性

我们现在将使用 Nisight 逐步执行一些代码，以帮助我们更好地理解一些 CUDA GPU 架构，以及内核内的分支是如何处理的。这将使我们对如何编写更有效的 CUDA 内核有一些见解。通过分支，我们指的是 GPU 如何处理 CUDA 内核中的if、else或switch等控制流语句。特别是，我们对内核内的分支分歧如何处理感兴趣，即当内核中的一个线程满足成为if语句的条件时，而另一个线程不满足条件并成为else语句时会发生什么：它们是分歧的，因为它们执行不同的代码片段。

让我们写一个小的 CUDA-C 程序作为实验：我们将从一个小内核开始，如果其threadIdx.x值是偶数，则打印一个输出，如果是奇数，则打印另一个输出。然后，我们编写一个main函数，将这个内核启动到由 32 个不同线程组成的一个单一块上：

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

__global__ void divergence_test_ker()
{
    if( threadIdx.x % 2 == 0)
        printf("threadIdx.x %d : This is an even thread.\n", threadIdx.x);
    else
        printf("threadIdx.x %d : This is an odd thread.\n", threadIdx.x);
}

__host__ int main()
{
    cudaSetDevice(0);
    divergence_test_ker<<<1, 32>>>();
    cudaDeviceSynchronize();
    cudaDeviceReset();
}

（此代码也可在存储库中的divergence_test.cu中找到。）

如果我们从命令行编译和运行这个程序，我们可能天真地期望偶数和奇数线程之间会有交错的字符串序列；或者它们可能会随机交错——因为所有线程都是并发运行并且大约在同一时间分支，这是有道理的。

相反，每次我们运行这个程序，我们总是得到这个输出：

所有与偶数线程对应的字符串都先打印出来，而所有与奇数线程对应的字符串都在第二次打印出来。也许 Nisight 调试器可以解释一些问题；让我们像在上一节中那样将这个小程序导入 Nisight 项目，并在内核的第一个if语句处设置断点。然后我们将执行step over，这样调试器就会在第一个printf语句处停下来。由于 Nisight 中的默认线程是（0,0,0），这应该满足了第一个if语句，所以它会一直停在那里，直到调试器继续。

让我们切换到一个奇数线程，比如（1,0,0），看看它现在在我们的程序中的位置：

非常奇怪！线程（1,0,0）在执行中也与线程（0,0,0）处于相同的位置。实际上，如果我们在这里检查每一个其他奇数线程，它们都会停在同一个地方——在一个所有奇数线程应该跳过的printf语句处。

这是什么？这被称为warp 锁步特性。CUDA 架构中的一个warp是一个由 32 个“通道”组成的单元，在这个单元中，我们的 GPU 执行内核和网格，其中每个通道将执行一个线程。warp 的一个主要限制是，所有在一个 warp 上执行的线程必须以锁步的方式执行相同的代码；这意味着并非每个线程确实运行相同的代码，而只是忽略对它不适用的步骤。（这被称为锁步，因为就像一群士兵一起锁步行进一样——不管他们是否想要行进！）

锁步特性意味着如果一个 warp 上运行的单个线程在单个if语句中与其他 31 个线程产生分歧，所有其他 31 个线程的执行都会延迟，直到这个单独的异常线程完成并从其孤立的if分歧中返回。这是在编写内核时应该牢记的一个特性，也是为什么在 CUDA 编程中，分支分歧应该尽量减少的一般规则。

使用 NVIDIA nvprof 分析器和 Visual Profiler

最后，我们将简要概述命令行 Nvidia nvprof分析器。与 Nsight IDE 相比，我们可以自由使用我们编写的任何 Python 代码——我们不必在这里强制编写完整的、纯粹的 CUDA-C 测试函数代码。

我们可以使用nvprof program命令对二进制可执行程序进行基本分析；同样，我们可以使用python命令作为第一个参数，使用脚本作为第二个参数来对 Python 脚本进行分析，如下所示：nvprof python program.py。让我们使用nvprof matrix_ker对我们之前编写的简单矩阵乘法 CUDA-C 可执行程序进行分析：

我们看到这与我们最初使用的 Python cProfiler 模块输出非常相似，我们用它来分析第一章中的 Mandelbrot 算法——只是现在，这专门告诉我们有关执行的所有 CUDA 操作。因此，当我们专门想要在 GPU 上进行优化时，我们可以使用它，而不必关心在主机上执行的任何 Python 或其他命令。（如果我们添加--print-gpu-trace命令行选项，我们可以进一步分析每个单独的 CUDA 内核操作，包括块和网格大小的启动参数。）

让我们再看一个技巧，帮助我们可视化程序所有操作的执行时间；我们将使用nvprof来转储一个文件，然后可以由 NVIDIA Visual Profiler 读取，以图形方式显示给我们。我们将使用上一章的示例multi-kernel_streams.py（在存储库的5下可用）来做这个。让我们回忆一下，这是我们对 CUDA 流概念的介绍示例之一，它允许我们同时执行和组织多个 GPU 操作。我们将使用-o命令行选项将输出转储到一个带有.nvvp文件后缀的文件中，如下所示：nvprof -o m.nvvp python multi-kernel_streams.py。现在我们可以使用nvvp m.nvvp命令将此文件加载到 NVIDIA Visual Profiler 中。

我们应该在所有 CUDA 流上看到一个时间线（记住，此程序中使用的内核名称为mult_ker）：

不仅可以看到所有的内核启动，还可以看到内存分配、内存复制和其他操作。这对于直观和视觉上理解程序如何随着时间使用 GPU 是有用的。

总结

我们在本章开始时看到了如何在 CUDA 内核中使用printf来输出来自各个线程的数据；我们特别看到了这对于调试代码有多么有用。然后，我们涵盖了 CUDA-C 中我们知识的一些空白，以便我们可以编写完整的测试程序，将其编译成适当的可执行二进制文件：这里有很多开销在我们之前是隐藏的，我们必须非常谨慎。接下来，我们看到了如何在 Nsight IDE 中创建和编译项目以及如何使用它进行调试。我们看到了如何在 CUDA 内核中停止我们设置的任何断点，并在不同的本地变量之间切换以查看各个线程。我们还使用了 Nsight 调试器来了解 warp 锁步属性以及为什么在 CUDA 内核中避免分支分歧很重要。最后，我们对 NVIDIA 命令行nvprof分析器和 Visual Profiler 进行了非常简要的概述，用于分析我们的 GPU 代码。

问题

1. 在我们编写的第一个 CUDA-C 程序中，我们在使用cudaMalloc在 GPU 上分配内存数组之后没有使用cudaDeviceSynchronize命令。为什么这是不必要的？（提示：回顾上一章。）

2. 假设我们有一个单个内核，它在由两个块组成的网格上启动，每个块有 32 个线程。假设第一个块中的所有线程都执行一个if语句，而第二个块中的所有线程都执行相应的else语句。第二个块中的所有线程是否必须像第一个块中的线程一样“锁步”执行if语句中的命令？

3. 如果我们执行类似的代码片段，只是在由 64 个线程执行的一个单一块组成的网格上执行，其中前 32 个线程执行一个if，而后 32 个执行一个else语句，会怎么样？

4. nvprof分析器可以为我们测量哪些 Python 的 cProfiler 无法测量的内容？

5. 列举一些我们可能更喜欢使用printf来调试 CUDA 内核的情境，以及其他一些情境，其中使用 Nsight 来调试 CUDA 内核可能更容易。

6. cudaSetDevice命令在 CUDA-C 中的目的是什么？

7. 为什么我们在每次 CUDA-C 中的内核启动或内存复制后都必须使用cudaDeviceSynchronize？

第七章：使用 Scikit-CUDA 与 CUDA 库

在本章中，我们将介绍三个用于简化数值和科学计算的标准 CUDA 库。我们将首先看一下cuBLAS，这是 NVIDIA 针对 CUDA 的基本线性代数子程序（BLAS）规范的实现。（cuBLAS 是 NVIDIA 对 BLAS 的各种优化的 CPU 实现的回应，例如免费/开源的 OpenBLAS 或英特尔的专有数学核心库。）接下来我们将看一下cuFFT，它可以在 GPU 上执行几乎每种快速傅里叶变换（FFT）的变体。我们将看看如何在图像处理中使用 cuFFT 进行滤波。然后我们将看一下cuSolver，它可以执行比 cuBLAS 中更复杂的线性代数运算，例如奇异值分解（SVD）或乔列斯基分解。

到目前为止，我们主要处理了一个作为我们与 CUDA 网关的单个 Python 模块——PyCUDA。虽然 PyCUDA 是一个非常强大和多功能的 Python 库，但它的主要目的是提供一个网关来编写、编译和启动 CUDA 内核，而不是提供一个接口给 CUDA 库。幸运的是，有一个免费的 Python 模块可用，它提供了一个用户友好的包装器接口给这些库。这就是 Scikit-CUDA。

虽然您不必了解 PyCUDA 甚至理解 GPU 编程就能欣赏 Scikit-CUDA，但它与 PyCUDA 兼容，例如，Scikit-CUDA 可以轻松地与 PyCUDA 的gpuarray类一起使用，这使您可以轻松地在我们自己的 CUDA 内核例程和 Scikit-CUDA 之间传递数据。此外，大多数例程也可以与 PyCUDA 的 stream 类一起使用，这将允许我们正确地同步我们自己的自定义 CUDA 内核和 Scikit-CUDA 的包装器。

请注意，除了这三个列出的库之外，Scikit-CUDA 还为专有的 CULA 库提供了包装器，以及开源的 MAGMA 库。这两者在功能上与官方的 NVIDIA 库有很多重叠。由于这些库在标准 CUDA 安装中默认未安装，我们将选择不在本章中涵盖它们。感兴趣的读者可以分别在www.culatools.com和icl.utk.edu/magma/了解更多关于 CULA 和 MAGMA 的信息。建议读者查看 Scikit-CUDA 的官方文档，网址为：media.readthedocs.org/pdf/scikit-cuda/latest/scikit-cuda.pdf。

本章的学习成果如下：

• 了解如何安装 Scikit-CUDA

• 了解标准 CUDA 库的基本目的和区别

• 了解如何使用低级 cuBLAS 函数进行基本线性代数

• 了解如何使用 SGEMM 和 DGEMM 操作来测量 GPU 在 FLOPS 中的性能

• 了解如何使用 cuFFT 在 GPU 上执行 1D 或 2D FFT 操作

• 了解如何使用 FFT 创建 2D 卷积滤波器，并将其应用于简单的图像处理

• 了解如何使用 cuSolver 执行奇异值分解（SVD）

• 了解如何使用 cuSolver 的 SVD 算法执行基本主成分分析

技术要求

本章需要一台安装了现代 NVIDIA GPU（2016 年及以后）的 Linux 或 Windows 10 PC，并安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（9.0 及以后）。还需要一个合适的 Python 2.7 安装（如 Anaconda Python 2.7），其中包括 PyCUDA 模块。

本章的代码也可以在 GitHub 上找到，网址为：github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA。

有关先决条件的更多信息，请查看本书的前言。有关软件和硬件要求的更多信息，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA上的 README 文件。

安装 Scikit-CUDA

建议您直接从 GitHub 安装最新稳定版本的 Scikit-CUDA：github.com/lebedov/scikit-cuda。

将软件包解压缩到一个目录中，然后在此处打开命令行，并键入python setup.py install来安装模块。然后，您可以运行单元测试以确保已使用python setup.py test执行了正确的安装。（此方法建议 Windows 和 Linux 用户均使用。）或者，可以直接使用pip install scikit-cuda从 PyPI 存储库安装 Scikit-CUDA。

使用 cuBLAS 进行基本线性代数

我们将从学习如何使用 Scikit-CUDA 的 cuBLAS 包装器开始这一章。让我们花一点时间讨论 BLAS。BLAS（基本线性代数子程序）是一个基本线性代数库的规范，最早是在 1970 年代标准化的。BLAS 函数被分为几个类别，被称为级别。

Level 1 BLAS 函数包括纯粹在向量上的操作——向量-向量加法和缩放（也称为ax+y操作，或 AXPY），点积和范数。Level 2 BLAS 函数包括一般矩阵-向量操作（GEMV），例如矩阵与向量的乘法，而 Level 3 BLAS 函数包括“一般矩阵-矩阵”（GEMM）操作，例如矩阵-矩阵乘法。最初，这些库是在 1970 年代完全用 FORTRAN 编写的，因此您应该考虑到在使用和命名上可能存在一些看似过时的遗留问题，这可能对今天的新用户来说显得繁琐。

cuBLAS 是 NVIDIA 自己对 BLAS 规范的实现，当然是经过优化以充分利用 GPU 的并行性。Scikit-CUDA 提供了与 PyCUDA gpuarray对象兼容的 cuBLAS 包装器，以及与 PyCUDA 流兼容的包装器。这意味着我们可以通过 PyCUDA 将这些函数与我们自己的自定义 CUDA-C 内核耦合和接口，以及在多个流上同步这些操作。

使用 cuBLAS 进行 Level-1 AXPY

现在让我们从 cuBLAS 开始进行基本的 Level-1 ax + y（或 AXPY）操作。让我们停下来，回顾一下线性代数的一点，并思考这意味着什么。在这里，a被认为是一个标量；也就是说，一个实数，比如-10、0、1.345 或 100。x和y被认为是某个向量空间中的向量，。这意味着x和y是实数的 n 元组，因此在的情况下，这些值可以是[1,2,3]或[-0.345, 8.15, -15.867]。ax表示x的缩放乘以a，因此如果a是 10 且x是先前的第一个值，则ax是x的每个单独值乘以a；也就是[10, 20, 30]。最后，和ax + y表示我们将两个向量中每个槽的每个单独值相加以产生一个新的向量，假设y是给定的第二个向量，结果将如下所示-[9.655, 28.15, 14.133]。

现在让我们在 cuBLAS 中做这个。首先，让我们导入适当的模块：

import pycuda.autoinit
from pycuda import gpuarray
import numpy as np

现在让我们导入 cuBLAS：

from skcuda import cublas

我们现在可以设置我们的向量数组并将它们复制到 GPU。请注意，我们使用 32 位（单精度）浮点数：

a = np.float32(10)
x = np.float32([1,2,3])
y = np.float32([-.345,8.15,-15.867])
x_gpu = gpuarray.to_gpu(x)
y_gpu = gpuarray.to_gpu(y)

现在我们必须创建一个cuBLAS 上下文。这在本质上类似于 CUDA 上下文，在第五章中我们讨论过，只是这一次它被显式用于管理 cuBLAS 会话。cublasCreate函数创建一个 cuBLAS 上下文，并将其句柄作为输出。我们需要保存这个句柄，只要我们打算在此会话中使用 cuBLAS：

cublas_context_h = cublas.cublasCreate()

现在我们可以使用cublasSaxpy函数。S代表单精度，这是我们需要的，因为我们正在处理 32 位浮点数组：

cublas.cublasSaxpy(cublas_context_h, x_gpu.size, a, x_gpu.gpudata, 1, y_gpu.gpudata, 1)

让我们讨论一下我们刚刚做的事情。还要记住，这是一个直接包装到低级 C 函数的函数，因此输入可能更像 C 函数而不是真正的 Python 函数。简而言之，这执行了一个“AXPY”操作，最终将输出数据放入y_gpu数组中。让我们逐个讨论每个输入参数。

第一个输入始终是 CUDA 上下文句柄。然后我们必须指定向量的大小，因为这个函数最终将在 C 指针上操作；我们可以使用 gpuarray 的size参数来做到这一点。已经将我们的标量类型转换为 NumPy 的float32变量，我们可以将a变量直接作为标量参数传递。然后我们使用gpudata参数将x_gpu数组的底层 C 指针传递给这个函数。然后我们将第一个数组的步长设置为 1：步长指定每个输入值之间应该走多少步。 （相反，如果您正在使用来自行向矩阵的列的向量，则将步长设置为矩阵的宽度。）然后我们放入y_gpu数组的指针，并将其步长也设置为 1。

我们已经完成了计算；现在我们必须明确销毁我们的 cuBLAS 上下文：

cublas.cublasDestroy(cublas_context)

现在我们可以使用 NumPy 的allclose函数来验证这是否接近，就像这样：

print 'This is close to the NumPy approximation: %s' % np.allclose(a*x + y , y_gpu.get())

再次注意，最终输出被放入了y_gpu数组中，这也是一个输入。

始终记住，BLAS 和 CuBLAS 函数是就地操作，以节省时间和内存，而不是进行新的分配调用。这意味着输入数组也将用作输出！

我们刚刚看到如何使用cublasSaxpy函数执行AXPY操作。

让我们讨论突出的大写字母 S。正如我们之前提到的，这代表单精度，即 32 位实数浮点值（float32）。如果我们想要操作 64 位实数浮点值数组（NumPy 和 PyCUDA 中的float64），那么我们将使用cublasDaxpy函数；对于 64 位单精度复数值（complex64），我们将使用cublasCaxpy，而对于 128 位双精度复数值（complex128），我们将使用cublasZaxpy。

我们可以通过检查函数名称其余部分之前的字母来确定 BLAS 或 CuBLAS 函数操作的数据类型。使用单精度实数的函数总是以 S 开头，双精度实数以 D 开头，单精度复数以 C 开头，双精度复数以 Z 开头。

其他一级 cuBLAS 函数

让我们看看其他一些一级函数。我们不会深入介绍它们的操作，但步骤与我们刚刚介绍的类似：创建一个 cuBLAS 上下文，使用适当的数组指针调用函数（可以通过 PyCUDA 的gpuarray的gpudata参数访问），并相应地设置步长。另一件需要记住的事情是，如果函数的输出是单个值而不是数组（例如，点积函数），则函数将直接将该值输出到主机，而不是在必须从 GPU 中取出的内存数组中。（我们只会在这里介绍单精度实数版本，但其他数据类型的相应版本可以通过用适当的字母替换 S 来使用。）

我们可以对两个单精度实数gpuarray，v_gpu和w_gpu进行点积。再次，1 是为了确保我们在这个计算中使用步长 1！再次回想一下，点积是两个向量的逐点乘积的和：

dot_output = cublas.cublasSdot(cublas_context_h, v_gpu.size, v_gpu.gpudata, 1, w_gpu.gpudata, 1)

我们还可以执行向量的 L2 范数，就像这样（回想一下，对于一个向量x，这是它的 L2 范数，或长度，可以用公式来计算）：

l2_output = cublas.cublasSnrm2(cublas_context_h, v_gpu.size, v_gpu.gpudata, 1)

cuBLAS 中的二级 GEMV

让我们看看如何进行GEMV矩阵-向量乘法。对于一个m x n矩阵A，一个 n 维向量x，一个m维向量y，以及标量alpha和beta，这被定义为以下操作：

现在让我们在继续之前看一下函数的布局：

cublasSgemv(handle, trans, m, n, alpha, A, lda, x, incx, beta, y, incy)  

让我们逐个检查这些输入：

• handle指的是 cuBLAS 上下文句柄。

• trans指的是矩阵的结构——我们可以指定是否要使用原始矩阵、直接转置或共轭转置（对于复数矩阵）。这很重要要记住，因为这个函数将期望矩阵A以列主格式存储。

• m和n是我们想要使用的矩阵A的行数和列数。

• alpha是浮点值α。

• A是m x n矩阵A。

• lda指示矩阵的主维度，矩阵的总大小实际上是lda x n。这在列主格式中很重要，因为如果lda大于m，这可能会导致 cuBLAS 在尝试访问A的值时出现问题，因为该矩阵的基础结构是一个一维数组。

• 然后我们有x及其步长incx；x是被A相乘的向量的基础 C 指针。记住，x的大小必须是n；也就是说，A的列数。

• beta是浮点值β。

• 最后，我们有y及其步长incy作为最后的参数。我们应该记住，y的大小应该是m，或者A的行数。

让我们通过生成一个 10 x 100 的随机值矩阵A，和一个包含 100 个随机值的向量x来测试这个。我们将初始化y为一个包含 10 个零的矩阵。我们将 alpha 设置为 1，beta 设置为 0，以便直接进行矩阵乘法而不进行缩放：

m = 10
n = 100
alpha = 1
beta = 0
A = np.random.rand(m,n).astype('float32')
x = np.random.rand(n).astype('float32')
y = np.zeros(m).astype('float32')

我们现在必须将A转换为列主（或按列）格式。NumPy 默认将矩阵存储为行主（或按行）格式，这意味着用于存储矩阵的基础一维数组会遍历所有第一行的值，然后遍历所有第二行的值，依此类推。您应该记住，转置操作会将矩阵的列与其行交换。然而，结果将是转置矩阵的新一维数组将以列主格式表示原始矩阵。我们可以通过A.T.copy()这样做，将A的转置矩阵的副本以及x和y一起复制到 GPU 上：

A_columnwise = A.T.copy()
A_gpu = gpuarray.to_gpu(A_columnwise) 
x_gpu = gpuarray.to_gpu(x)
y_gpu = gpuarray.to_gpu(y)

由于我们现在已经正确地在 GPU 上存储了按列的矩阵，我们可以通过使用_CUBLAS_OP字典将trans变量设置为不进行转置：

trans = cublas._CUBLAS_OP['N']

由于矩阵的大小与我们想要使用的行数完全相同，我们现在将lda设置为m。x和y向量的步长再次为 1。我们现在已经设置好了所有需要的值，现在可以创建我们的 CuBLAS 上下文并存储它的句柄，就像这样：

lda = m 
incx = 1
incy = 1
handle = cublas.cublasCreate()

我们现在可以启动我们的函数。记住，A，x和y实际上是 PyCUDA gpuarray对象，所以我们必须使用gpudata参数输入到这个函数中。除了这个，这很简单：

cublas.cublasSgemv(handle, trans, m, n, alpha, A_gpu.gpudata, lda, x_gpu.gpudata, incx, beta, y_gpu.gpudata, incy)

现在我们可以销毁我们的 cuBLAS 上下文并检查返回值以确保它是正确的：

cublas.cublasDestroy(handle)
print 'cuBLAS returned the correct value: %s' % np.allclose(np.dot(A,x), y_gpu.get())

cuBLAS 中的三级 GEMM 用于测量 GPU 性能

现在我们将看看如何使用 CuBLAS 执行通用矩阵-矩阵乘法（GEMM）。实际上，我们将尝试制作一些比我们在 cuBLAS 中看到的最后几个示例更实用的东西-我们将使用这个作为我们的 GPU 性能指标，以确定它可以执行的每秒浮点运算次数（FLOPS）的数量，这将是两个单独的值：单精度和双精度的情况。使用 GEMM 是评估 FLOPS 中计算硬件性能的标准技术，因为它比使用纯时钟速度（MHz 或 GHz）更好地理解了纯计算能力。

如果您需要简要回顾，请回想一下我们在上一章中深入讨论了矩阵-矩阵乘法。如果您忘记了这是如何工作的，强烈建议您在继续本节之前复习一下这一章。

首先，让我们看看 GEMM 操作是如何定义的：

这意味着我们执行A和B的矩阵乘法，将结果缩放为alpha，然后加到我们已经通过beta缩放的C矩阵中，将最终结果放在C中。

让我们考虑一下执行实值 GEMM 操作的最终结果需要执行多少浮点运算，假设A是一个m x k（其中m是行，k是列）矩阵，B是一个k x n矩阵，C 是一个m x n矩阵。首先，让我们计算AB需要多少操作。让我们取A的一列并将其乘以B：这将导致m行中的每一行需要k次乘法和k-1次加法，这意味着这是m行上的km + (k-1)m总操作。B中有n列，因此计算AB将总共需要kmn + (k-1)mn = 2kmn - mn次操作。现在，我们使用alpha来缩放AB，这将是m**n次操作，因为这是矩阵AB的大小；类似地，通过beta缩放C是另外m**n次操作。最后，我们将这两个结果矩阵相加，这又是mn次操作。这意味着在给定的 GEMM 操作中，我们将有2kmn - mn + 3mn = 2kmn + 2mn = 2mn(k+1)次浮点运算。

现在我们唯一需要做的就是运行一个计时的 GEMM 操作，注意矩阵的不同大小，并将2kmn + 2mn除以总时间持续时间来计算我们 GPU 的 FLOPS。结果数字将非常大，因此我们将以 GFLOPS 的形式表示这一点-也就是说，每秒可以计算多少十亿（10⁹）次操作。我们可以通过将 FLOPS 值乘以 10^(-9)来计算这一点。

现在我们准备开始编写代码。让我们从我们的导入语句开始，以及time函数：

import pycuda.autoinit
from pycuda import gpuarray
import numpy as np
from skcuda import cublas
from time import time

现在我们将为我们的矩阵大小设置m，n和k变量。我们希望我们的矩阵相对较大，以便时间持续足够长，以避免除以 0 的错误。以下值对于截至 2018 年中旬或更早发布的任何 GPU 来说应该足够了；拥有更新卡的用户可能考虑增加这些值：

m = 5000
n = 10000
k = 10000

现在我们将编写一个计算单精度和双精度 GFLOPS 的函数。如果我们希望使用双精度，则将输入值设置为'D'，否则设置为'S'：

def compute_gflops(precision='S'):

if precision=='S':
    float_type = 'float32'
elif precision=='D':
    float_type = 'float64'
else:
    return -1

现在让我们生成一些随机矩阵，这些矩阵具有我们将用于计时的适当精度。GEMM 操作与我们之前看到的 GEMV 操作类似，因此我们必须在将它们复制到 GPU 之前对其进行转置。（由于我们只是在计时，这一步并不是必要的，但记住这一点是个好习惯。）

我们将为 GEMM 设置一些其他必要的变量，这些变量在这一点上应该是不言自明的（transa，lda，ldb等）：

A = np.random.randn(m, k).astype(float_type)
B = np.random.randn(k, n).astype(float_type)
C = np.random.randn(m, n).astype(float_type)
A_cm = A.T.copy()
B_cm = B.T.copy()
C_cm = C.T.copy()
A_gpu = gpuarray.to_gpu(A_cm)
B_gpu = gpuarray.to_gpu(B_cm)
C_gpu = gpuarray.to_gpu(C_cm)
alpha = np.random.randn()
beta = np.random.randn()
transa = cublas._CUBLAS_OP['N']
transb = cublas._CUBLAS_OP['N']
lda = m
ldb = k
ldc = m

现在我们可以开始计时了！首先，我们将创建一个 cuBLAS 上下文：

t = time()
handle = cublas.cublasCreate()

现在我们将启动 GEMM。请记住，对于实数情况有两个版本：cublasSgemm用于单精度，cublasDgemm用于双精度。我们可以使用一个小小的 Python 技巧执行适当的函数：我们将用cublas%sgemm和适当的参数写一个字符串，然后通过附加% precision将%s替换为 D 或 S。然后我们将使用exec函数将这个字符串作为 Python 代码执行，就像这样：

exec('cublas.cublas%sgemm(handle, transa, transb, m, n, k, alpha, A_gpu.gpudata, lda, B_gpu.gpudata, ldb, beta, C_gpu.gpudata, ldc)' % precision)

现在我们可以销毁 cuBLAS 上下文，并得到我们计算的最终时间：

cublas.cublasDestroy(handle)
t = time() - t

然后我们需要使用我们推导出的方程计算 GFLOPS，并将其作为这个函数的输出返回：

gflops = 2*m*n*(k+1)*(10**-9) / t 
return gflops

现在我们可以设置我们的主函数。我们将输出单精度和双精度情况下的 GFLOPS：

if __name__ == '__main__':
    print 'Single-precision performance: %s GFLOPS' % compute_gflops('S')
    print 'Double-precision performance: %s GFLOPS' % compute_gflops('D')

现在在运行这个程序之前，让我们做一点功课——去www.techpowerup.com搜索你的 GPU，然后注意两件事——单精度浮点性能和双精度浮点性能。我现在使用的是 GTX 1050，它的列表声称在单精度上有 1,862 GFLOPS 性能，在双精度上有 58.20 GFLOPS 性能。让我们现在运行这个程序，看看这是否符合事实：

看哪，它成功了！

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

使用 cuFFT 进行快速傅里叶变换

现在让我们看看如何使用 cuFFT 进行一些基本的快速傅里叶变换（FFT）。首先，让我们简要回顾一下傅里叶变换到底是什么。如果你上过高级微积分或分析课程，你可能已经见过傅里叶变换被定义为一个积分公式，就像这样：

这个程序将f作为x的时间域函数。这给了我们一个相应的频域函数，对应于"ξ"。这事实上是一个极其有用的工具，几乎触及到所有科学和工程的分支。

让我们记住积分可以被看作是一个求和；同样地，有一个对应的离散、有限版本的傅里叶变换，称为离散傅里叶变换（DFT）。这对长度为n的向量进行操作，并允许它们在频域中进行分析或修改。x的 DFT 被定义如下：

换句话说，我们可以将一个向量x乘以复杂的N x N矩阵 

（这里，k对应行号，n对应列号）来找到它的 DFT。我们还应该注意到逆公式，它让我们从它的 DFT 中检索x（在这里用x的 DFT 替换y，输出将是原始的x）：

通常，计算矩阵-向量操作的计算复杂度对于长度为N的向量是 O(N²)。然而，由于 DFT 矩阵中的对称性，这总是可以通过使用 FFT 减少到 O(N log N)。让我们看看如何使用 FFT 与 CuBLAS，然后我们将继续一个更有趣的例子。

一个简单的一维 FFT

让我们首先看看如何使用 cuBLAS 计算简单的一维 FFT。首先，我们将简要讨论 Scikit-CUDA 中的 cuFFT 接口。

这里有两个子模块，我们可以使用cufft和fft访问 cuFFT 库。cufft包括了一系列 cuFFT 库的低级封装，而fft提供了一个更加用户友好的接口；在本章中我们将只使用fft。

让我们从适当的导入开始，记得包括 Scikit-CUDA 的fft子模块：

import pycuda.autoinit
from pycuda import gpuarray
import numpy as np
from skcuda import fft

现在我们将设置一些随机数组并将其复制到 GPU。我们还将设置一个空的 GPU 数组，用于存储 FFT（请注意，我们使用实数 float32 数组作为输入，但输出将是 complex64 数组，因为傅立叶变换总是复值）：

x = np.asarray(np.random.rand(1000), dtype=np.float32 )
x_gpu = gpuarray.to_gpu(x)
x_hat = gpuarray.empty_like(x_gpu, dtype=np.complex64)

我们现在将为正向 FFT 变换设置一个 cuFFT 计划。这是 cuFFT 用来确定变换的形状，以及输入和输出数据类型的对象：

plan = fft.Plan(x_gpu.shape,np.float32,np.complex64)

我们还将为逆 FFT 计划对象设置一个计划。请注意，这一次我们从complex64到实数float32：

inverse_plan = fft.Plan(x.shape, in_dtype=np.complex64, out_dtype=np.float32)

现在，我们必须将x_gpu的正向 FFT 转换为x_hat，并将x_hat的逆 FFT 转换回x_gpu。请注意，在逆 FFT 中设置了scale=True；我们这样做是为了告诉 cuFFT 将逆 FFT 按 1/N 进行缩放：

fft.fft(x_gpu, x_hat, plan)
fft.ifft(x_hat, x_gpu, inverse_plan, scale=True)

我们现在将检查x_hat与x的 NumPy FFT，以及x_gpu与x本身：

y = np.fft.fft(x)
print 'cuFFT matches NumPy FFT: %s' % np.allclose(x_hat.get(), y, atol=1e-6)
print 'cuFFT inverse matches original: %s' % np.allclose(x_gpu.get(), x, atol=1e-6)

如果你运行这个程序，你会发现x_hat与y不匹配，然而，莫名其妙的是，x_gpu与x匹配。这是怎么可能的？好吧，让我们记住x是实数；如果你看一下离散傅立叶变换是如何计算的，你可以数学上证明实向量的输出在 N/2 之后会重复为它们的复共轭。虽然 NumPy FFT 会完全计算这些值，但 cuFFT 通过只计算输入为实数时的输出的前半部分来节省时间，并将其余输出设置为0。你应该通过检查前面的变量来验证这一点。

因此，如果我们将前面代码中的第一个打印语句更改为仅比较 CuFFT 和 NumPy 之间的前 N/2 个输出，那么这将返回 true：

print 'cuFFT matches NumPy FFT: %s' % np.allclose(x_hat.get()[0:N//2], y[0:N//2], atol=1e-6)

使用 FFT 进行卷积

我们现在将看看如何使用 FFT 执行卷积。首先让我们回顾一下卷积的确切定义：给定两个一维向量x和y，它们的卷积定义如下：

这对我们很有意义，因为如果x是一些长的连续信号，而y只有一小部分局部非零值，那么y将作为x的滤波器；这本身有许多应用。首先，我们可以使用滤波器平滑信号x（在数字信号处理和图像处理中很常见）。我们还可以使用它来收集信号x的样本，以表示信号或压缩它（在数据压缩或压缩感知领域很常见），或者使用滤波器为机器学习中的信号或图像识别收集特征。这个想法构成了卷积神经网络的基础。

当然，计算机无法处理无限长的向量（至少目前还不能），所以我们将考虑循环卷积。在循环卷积中，我们处理两个长度为n的向量，其索引小于 0 或大于 n-1 的部分将会环绕到另一端；也就是说，x[-1] = x[n-1]，x[-2] = x[n-2]，x[n] = x[0]，x[n+1] = x[1]，依此类推。我们定义x和y的循环卷积如下：

事实证明，我们可以很容易地使用 FFT 执行循环卷积；我们可以对x和y执行 FFT，对输出进行逐点乘法，然后对最终结果执行逆 FFT。这个结果被称为卷积定理，也可以表示如下：

我们将在两个维度上进行这个操作，因为我们希望将结果应用到信号处理上。虽然我们只看到了一维卷积和 FFT 的数学运算，但二维卷积和 FFT 的工作方式与一维类似，只是索引更复杂一些。然而，我们选择跳过这一点，以便我们可以直接进入应用。

使用 cuFFT 进行 2D 卷积

现在我们将制作一个小程序，使用基于 cuFFT 的二维卷积对图像进行高斯滤波。高斯滤波是一种使用所谓的高斯滤波器平滑粗糙图像的操作。之所以这样命名，是因为它基于统计学中的高斯（正态）分布。这是高斯滤波器在两个维度上以标准偏差σ定义的方式：

当我们用滤波器对离散图像进行卷积时，有时我们会将滤波器称为卷积核。通常，图像处理工程师会简单地称之为核，但由于我们不想将其与 CUDA 核混淆，我们将始终使用完整术语卷积核。在这里，我们将使用离散版本的高斯滤波器作为我们的卷积核。

让我们从适当的导入开始；请注意，我们将在这里使用 Scikit-CUDA 子模块linalg。这将为我们提供比 cuBLAS 更高级的接口。由于我们在这里处理图像，我们还将导入 Matplotlib 的pyplot子模块。还要注意，我们将在此处使用 Python 3 风格的除法，从第一行开始；这意味着如果我们使用/运算符除两个整数，那么返回值将是浮点数，无需类型转换（我们使用//运算符进行整数除法）：

from __future__ import division
import pycuda.autoinit
from pycuda import gpuarray
import numpy as np
from skcuda import fft
from skcuda import linalg
from matplotlib import pyplot as plt

让我们立即开始编写卷积函数。这将接受两个相同大小的 NumPy 数组x和y。我们将把它们转换为 complex64 数组，然后如果它们的大小不相同，就返回-1：

def cufft_conv(x , y):
    x = x.astype(np.complex64)
    y = y.astype(np.complex64)

    if (x.shape != y.shape):
        return -1

现在我们将设置我们的 FFT 计划和逆 FFT 计划对象：

plan = fft.Plan(x.shape, np.complex64, np.complex64)
inverse_plan = fft.Plan(x.shape, np.complex64, np.complex64)

现在我们可以将我们的数组复制到 GPU。我们还将设置一些空数组，用于保存这些数组的 FFT 的适当大小，另外一个数组将保存最终卷积的输出，out_gpu：

 x_gpu = gpuarray.to_gpu(x)
 y_gpu = gpuarray.to_gpu(y)

 x_fft = gpuarray.empty_like(x_gpu, dtype=np.complex64)
 y_fft = gpuarray.empty_like(y_gpu, dtype=np.complex64)
 out_gpu = gpuarray.empty_like(x_gpu, dtype=np.complex64)

现在我们可以进行我们的 FFT：

fft.fft(x_gpu, x_fft, plan)
fft.fft(y_gpu, y_fft, plan)

我们现在将使用linalg.multiply函数在x_fft和y_fft之间执行逐点（Hadamard）乘法。我们将设置overwrite=True，以便将最终值写入y_fft：

linalg.multiply(x_fft, y_fft, overwrite=True)

现在我们将调用逆 FFT，将最终结果输出到out_gpu。我们将这个值传输到主机并返回它：

fft.ifft(y_fft, out_gpu, inverse_plan, scale=True)
conv_out = out_gpu.get()
return conv_out

我们还没有完成。我们的卷积核将比输入图像小得多，因此我们将不得不调整我们的两个 2D 数组的大小（卷积核和图像），使它们相等，并在它们之间执行逐点乘法。我们不仅应该确保它们相等，还需要确保我们在数组上执行零填充，并适当地将卷积核居中。零填充意味着我们在图像的两侧添加零缓冲区，以防止环绕错误。如果我们使用 FFT 来执行卷积，记住它是循环卷积，因此边缘将始终环绕。当我们完成卷积后，我们可以去除图像外部的缓冲区，得到最终的输出图像。

让我们创建一个名为conv_2d的新函数，它接受卷积核ker和图像img。填充后的图像大小将是（2*ker.shape[0] + img.shape[0]，2*ker.shape[1] + img.shape[1]）。让我们首先设置填充的卷积核。我们将创建一个这个大小的零矩阵，然后将左上角子矩阵设置为我们的卷积核，如下所示：

def conv_2d(ker, img):

    padded_ker = np.zeros( (img.shape[0] + 2*ker.shape[0], img.shape[1] + 2*ker.shape[1] )).astype(np.float32)
    padded_ker[:ker.shape[0], :ker.shape[1]] = ker

现在我们需要移动卷积核，使其中心精确地位于坐标（0,0）处。我们可以使用 NumPy 的roll命令来实现这一点：

padded_ker = np.roll(padded_ker, shift=-ker.shape[0]//2, axis=0)
padded_ker = np.roll(padded_ker, shift=-ker.shape[1]//2, axis=1)

现在我们需要对输入图像进行填充：

padded_img = np.zeros_like(padded_ker).astype(np.float32)
padded_img[ker.shape[0]:-ker.shape[0], ker.shape[1]:-ker.shape[1]] = img

现在我们有两个大小相同且格式正确的数组。我们现在可以使用我们刚刚在这里编写的cufft_conv函数：

out_ = cufft_conv(padded_ker, padded_img)

现在我们可以去除图像外部的零缓冲区。然后返回结果：

output = out_[ker.shape[0]:-ker.shape[0], ker.shape[1]:-ker.shape[1]]

return output

我们还没有完成。让我们编写一些小函数来设置我们的高斯滤波器，然后我们可以继续将其应用于图像。我们可以使用 lambda 函数一行代码来编写基本滤波器本身：

gaussian_filter = lambda x, y, sigma : (1 / np.sqrt(2*np.pi*(sigma**2)) )*np.exp( -(x**2 + y**2) / (2 * (sigma**2) ))

现在我们可以编写一个使用此滤波器输出离散卷积核的函数。卷积核的高度和长度将为2*sigma + 1，这是相当标准的：

请注意，我们通过将其值求和为total_并将其除以来规范化高斯核的值。

def gaussian_ker(sigma):
    ker_ = np.zeros((2*sigma+1, 2*sigma+1))
    for i in range(2*sigma + 1):
        for j in range(2*sigma + 1):
            ker_[i,j] = gaussian_filter(i - sigma, j - sigma, sigma)
    total_ = np.sum(ker_.ravel())
    ker_ = ker_ */* total*_* return ker_

我们现在准备在图像上测试这个！作为我们的测试案例，我们将使用高斯滤波来模糊这本书编辑Akshada Iyer的彩色 JPEG 图像。（此图像在 GitHub 存储库的Chapter07目录中，文件名为akshada.jpg。）我们将使用 Matplotlib 的imread函数来读取图像；默认情况下，这将存储为 0 到 255 范围内的无符号 8 位整数数组。我们将将其强制转换为浮点数组并对其进行规范化，以便所有值的范围为 0 到 1。

对于本书印刷版的读者注意：尽管本书的印刷版是灰度图像，但这是一幅彩色图像。然后，我们将设置一个空的零数组，用于存储模糊的图像：

if __name__ == '__main__':
    akshada = np.float32(plt.imread('akshada.jpg')) / 255
    akshada_blurred = np.zeros_like(akshada)

让我们设置我们的卷积核。在这里，标准差为 15 应该足够：

ker = gaussian_ker(15)

现在我们可以模糊图像。由于这是一幅彩色图像，我们将不得不分别对每个颜色层（红色、绿色和蓝色）应用高斯滤波；这在图像数组的第三维中进行索引：

for k in range(3):
    akshada_blurred[:,:,k] = conv_2d(ker, akshada[:,:,k])

现在让我们通过使用一些 Matplotlib 技巧并排查看 Before 和 After 图像：

fig, (ax0, ax1) = plt.subplots(1,2)
fig.suptitle('Gaussian Filtering', fontsize=20)
ax0.set_title('Before')
ax0.axis('off')
ax0.imshow(akshada)
ax1.set_title('After')
ax1.axis('off')
ax1.imshow(akshada_blurred)
plt.tight_layout()
plt.subplots_adjust(top=.85)
plt.show()

现在我们可以运行程序并观察高斯滤波的效果：

该程序在此书的存储库中的Chapter07目录中的名为conv_2d.py的文件中可用。

使用 Scikit-CUDA 中的 cuSolver

我们现在将看看如何使用 Scikit-CUDA 的linalg子模块中的 cuSolver。同样，这为 cuBLAS 和 cuSolver 提供了一个高级接口，因此我们不必陷入细节。

正如我们在介绍中指出的，cuSolver 是一个用于执行比 cuBLAS 更高级的线性代数运算的库，例如奇异值分解、LU/QR/Cholesky 分解和特征值计算。由于 cuSolver、cuBLAS 和 cuFFT 一样，是另一个庞大的库，我们只会花时间来看数据科学和机器学习中最基本的操作之一——SVD。

如果您想进一步了解此库，请参考 NVIDIA 关于 cuSOLVER 的官方文档：docs.NVIDIA.com/cuda/cusolver/index.html。

奇异值分解（SVD）

SVD 接受任何m x n矩阵A，然后返回三个矩阵—U、Σ和V。在这里，U是一个m x m酉矩阵，Σ是一个m x n对角矩阵，V是一个n x n酉矩阵。这里，酉表示矩阵的列形成一个正交归一基；对角表示矩阵中的所有值都为零，除了可能沿着对角线的值。

SVD 的重要性在于它将A分解为这些矩阵，使得我们有A = UΣV^T；此外，Σ对角线上的值将全部为正值或零，并且被称为奇异值。我们很快将看到一些应用，但您应该记住，SVD 的计算复杂度为 O(mn²)——对于大矩阵来说，使用 GPU 绝对是一个好主意，因为这个算法是可并行化的。

现在让我们看看如何计算矩阵的 SVD。让我们进行适当的导入语句：

import pycuda.autoinit
from pycuda import gpuarray
import numpy as np
from skcuda import linalg

现在我们将生成一个相对较大的随机矩阵并将其传输到 GPU：

a = np.random.rand(1000,5000).astype(np.float32)
a_gpu = gpuarray.to_gpu(a)

现在我们可以执行 SVD。这将有三个输出，对应于我们刚刚描述的矩阵。第一个参数将是我们刚刚复制到 GPU 的矩阵数组。然后我们需要指定我们要使用 cuSolver 作为此操作的后端：

U_d, s_d, V_d = linalg.svd(a_gpu,  lib='cusolver')

现在让我们将这些数组从 GPU 复制到主机：

U = U_d.get()
s = s_d.get()
V = V_d.get()

s实际上存储为一维数组；我们将不得不创建一个大小为 1000 x 5000 的零矩阵，并将这些值沿对角线复制。我们可以使用 NumPy 的diag函数来做到这一点，再加上一些数组切片：

S = np.zeros((1000,5000))
S[:1000,:1000] = np.diag(s)

我们现在可以在主机上使用 NumPy 的dot函数对这些值进行矩阵相乘，以验证它们是否与我们的原始数组匹配：

print 'Can we reconstruct a from its SVD decomposition? : %s' % np.allclose(a, np.dot(U, np.dot(S, V)), atol=1e-5)

由于我们只使用 float32，并且我们的矩阵相对较大，引入了一些数值误差；我们不得不将“容差”级别（atol）设置得比平常高一点，但仍然足够小，以验证这两个数组是足够接近的。

使用 SVD 进行主成分分析（PCA）

主成分分析（PCA）是主要用于降维的工具。我们可以使用它来查看数据集，并找出哪些维度和线性子空间最显著。虽然有几种实现方法，但我们将向您展示如何使用 SVD 执行 PCA。

我们将这样做——我们将使用一个存在于 10 个维度中的数据集。我们将首先创建两个在前面有很大权重的向量，其他位置为 0：

vals = [ np.float32([10,0,0,0,0,0,0,0,0,0]) , np.float32([0,10,0,0,0,0,0,0,0,0]) ]

然后我们将添加 9000 个额外的向量：其中 6000 个将与前两个向量相同，只是加了一点随机白噪声，剩下的 3000 个将只是随机白噪声：

for i in range(3000):
    vals.append(vals[0] + 0.001*np.random.randn(10))
    vals.append(vals[1] + 0.001*np.random.randn(10))
    vals.append(0.001*np.random.randn(10))

我们现在将vals列表转换为float32的 NumPy 数组。我们对行进行平均，并从每行中减去这个值。（这是 PCA 的必要步骤。）然后我们转置这个矩阵，因为 cuSolver 要求输入矩阵的行数少于或等于列数：

vals = np.float32(vals)
vals = vals - np.mean(vals, axis=0)
v_gpu = gpuarray.to_gpu(vals.T.copy())

我们现在将运行 cuSolver，就像我们之前做的那样，并将输出值从 GPU 复制出来：

U_d, s_d, V_d = linalg.svd(v_gpu, lib='cusolver')

u = U_d.get()
s = s_d.get()
v = V_d.get()

现在我们准备开始我们的调查工作。让我们打开 IPython，仔细查看u和s。首先，让我们看看 s；它的值实际上是主要值的平方根，所以我们将它们平方然后看一看：

您会注意到，前两个主要值的数量级为 10⁵，而其余的分量的数量级为 10^(-3)。这告诉我们，实际上只有一个二维子空间与这些数据相关，这并不令人意外。这些是第一个和第二个值，它们将对应于第一个和第二个主成分，即相应的向量。让我们来看看这些向量，它们将存储在U中：

！

您会注意到这两个向量在前两个条目中有很大的权重，数量级为 10^(-1)；其余的条目都是 10^(-6)或更低，相对不相关。考虑到我们在前两个条目中使数据偏向，这正是我们应该期望的。这就是 PCA 背后的思想。

摘要

我们从查看如何使用 Scikit-CUDA 库的 cuBLAS 包装器开始了本章；在这里，我们必须记住许多细节，比如何使用列主存储，或者输入数组是否会被就地覆盖。然后，我们看了如何使用 Scikit-CUDA 的 cuFFT 执行一维和二维 FFT，以及如何创建一个简单的卷积滤波器。然后，我们向您展示了如何将其应用于图像的简单高斯模糊效果。最后，我们看了如何使用 cuSolver 在 GPU 上执行奇异值分解（SVD），这通常是一个非常计算密集的操作，但在 GPU 上可以很好地并行化。我们通过查看如何使用 SVD 进行基本的 PCA 来结束本章。

问题

1. 假设你得到了一个工作，需要将一些旧的遗留 FORTRAN BLAS 代码转换成 CUDA。你打开一个文件，看到一个名为 SBLAH 的函数，另一个名为 ZBLEH。你能在不查找的情况下告诉这两个函数使用的数据类型吗？

2. 你能修改 cuBLAS level-2 GEMV 示例，直接将矩阵A复制到 GPU，而不是在主机上进行转置以设置为列优先吗？

3. 使用 cuBLAS 32 位实数点积（cublasSdot）来实现使用一个按行矩阵和一个步幅为 1 的向量进行矩阵-向量乘法。

4. 使用cublasSdot实现矩阵-矩阵乘法。

5. 你能实现一种精确测量性能测量示例中的 GEMM 操作的方法吗？

6. 在一维 FFT 的示例中，尝试将x强制转换为complex64数组，然后将 FFT 和逆 FFT 计划都设置为complex64值。然后确认np.allclose(x, x_gpu.get())是否为真，而不检查数组的前一半。你认为为什么现在这样做会有效？

7. 请注意，在卷积示例中，模糊图像周围有一个暗边。为什么模糊图像中有这个暗边，而原始图像中没有呢？你能想到一种方法来减轻这个问题吗？

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

在上一章中，我们对通过 Scikit-CUDA 包装器模块在 CUDA 中可用的库进行了相当广泛的概述。我们现在将看一下另外一些库，我们将不得不直接从 CUDA C 中使用这些库，而不是像 Scikit-CUDA 中的包装器那样。我们将首先看一下两个标准库，其中包含我们可以从任何 CUDA C 内核中调用的设备函数 cuRAND 和 CUDA Math API。通过学习如何使用这些库，我们将了解如何在蒙特卡罗积分的上下文中使用这些库。蒙特卡罗积分是一种众所周知的随机方法，可以提供来自微积分的定积分值的估计。我们首先将看一个基本示例，演示如何使用 cuRAND 实现简单的蒙特卡罗方法来对π的值进行基本估计（如众所周知的常数，π=3.14159...），然后我们将着手进行一个更有雄心的项目，我们将构建一个 Python 类，可以对任意数学函数执行定积分，并使用 Math API 创建这样的函数。我们还将看一下如何在设计这个类时有效地使用元编程的一些思想。

然后，我们将再次使用 Thrust C++库来编写一些纯 CUDA 程序。Thrust 是一个提供 C++模板容器的库，类似于 C++标准模板库（STL）中的容器。这将使我们能够以更接近 PyCUDA 的gpuarray和 STL 的向量容器的更自然的方式从 C++中操作 CUDA C 数组。这将使我们免受在 CUDA C 中以前不断使用指针（如mallocs和frees）的困扰。

在本章中，我们将讨论以下主题：

• 理解种子在生成伪随机数列表中的作用

• 在 CUDA 内核中使用 cuRAND 设备函数生成随机数

• 理解蒙特卡罗积分的概念

• 在 Python 中使用基于字典的字符串格式化进行元编程

• 使用 CUDA Math API 设备函数库

• 理解 functor 是什么

• 在纯 CUDA C 编程时使用 Thrust 向量容器

技术要求

本章需要具有现代 NVIDIA GPU（2016 年至今）的 Linux 或 Windows 10 PC，并安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（9.0 及以上）。还需要适当的 Python 2.7 安装（如 Anaconda Python 2.7），并安装了 PyCUDA 模块。

本章的代码也可以在 GitHub 上找到，网址为github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA.

有关本章先决条件的更多信息，请查看本书的前言。有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA上的 README。

cuRAND 设备函数库

让我们从 cuRAND 开始。这是一个标准的 CUDA 库，用于在 CUDA 内核中按线程生成伪随机值，通过调用每个线程内核中的设备函数进行初始化和调用。让我们再次强调，这是一个伪随机值序列——因为数字硬件始终是确定性的，从不是随机或任意的，我们使用算法从初始种子值生成一系列表面上随机的值。通常，我们可以将种子值设置为真正的随机值（例如毫秒级的时钟时间），这将产生一系列随机值。这些生成的随机值与由相同种子生成的序列中的先前或未来值没有相关性，尽管当您组合从不同种子生成的值时，可能会出现相关性和重复。因此，您必须小心，希望彼此随机的值是由相同的种子生成的。

让我们从curand_init的函数原型开始，我们将使用适当的种子进行初始化：

__device__ void curand_init ( unsigned long long seed, unsigned long long sequence, unsigned long long offset, curandState_t *state)

在这里，所有的输入都是无符号长整型，在 C 中是无符号（非负值）的 64 位整数。首先，我们可以看到seed，这当然是种子值。一般来说，您将使用时钟值或某种变化来设置这个值。然后我们看到一个称为sequence的值，正如我们之前所述，cuRAND 生成的值只有在它们由相同的种子值生成时才是真正的数学上相互随机的。因此，如果我们有多个线程使用相同的种子值，我们使用sequence来指示当前线程使用长度为 2¹⁹⁰的随机数子序列的哪个子序列，而我们使用offset来指示从这个子序列的哪个点开始；这将在每个线程中生成所有数学上相互随机且没有相关性的值。最后，最后一个参数是指向curandState_t对象的指针；这将跟踪我们在伪随机数序列中的位置。

初始化类对象后，您将通过调用适当的设备函数从适当的随机分布生成随机值。最常见的两种分布是均匀分布和正态（高斯）分布。均匀分布（curand_uniform，在 cuRAND 中）是一个输出值在给定范围内都是等概率的函数：也就是说，对于 0 到 1 的均匀分布，值在 0 到 0.1 之间的概率是 10%，或者在 0.9 到 1 之间，或者在任何两个相距 0.1 的点之间。正态分布（curand_normal，在 cuRAND 中）具有以特定均值为中心的值，这些值将根据分布的标准差分布在众所周知的钟形曲线上。 （curand_normal的默认均值为0，标准差为 1，在 cuRAND 中，因此必须手动移位和缩放为其他值。）cuRAND 支持的另一个众所周知的分布是泊松分布（curand_poisson），用于对随机事件的发生进行建模。

在接下来的部分中，我们将主要研究如何在均匀分布的背景下使用 cuRAND，因为它们适用于蒙特卡罗积分。鼓励有兴趣学习如何使用 cuRAND 更多功能的读者查看 NVIDIA 的官方文档。

用蒙特卡洛法估算π

首先，我们将运用我们对 cuRAND 的新知识来估算众所周知的数学常数π，或圆周率，这当然是永不停止的无理数 3.14159265358979...

然而，要得到一个估计值，我们需要花点时间思考这意味着什么。让我们想想一个圆。记住，圆的半径是从圆心到圆上任意一点的长度；通常用 R 表示。直径被定义为 D = 2R，周长 C 是围绕圆的长度。然后，π 被定义为 π = C / D。我们可以使用欧几里得几何来找到圆的面积公式，结果是 A = πR²。现在，让我们考虑一个半径为 R 的圆被内切在边长为 2R 的正方形中：

因此，当然，我们知道正方形的面积是 (2R)² = 4R²。让我们考虑 R=1，这样我们就知道圆的面积恰好是 π，而正方形的面积恰好是 4。让我们进一步假设并声明，圆和正方形都以 (0,0) 为中心。现在，让我们在正方形内取一个完全随机的值 (x,y)，并查看它是否落在圆内。我们如何做到这一点？通过应用勾股定理公式：我们通过检查 x² + y² 是否小于或等于 1 来做到这一点。让我们用 iters 表示我们选择的随机点的总数，用 hits 表示命中的次数。

让我们再多想一下：在圆内选择一个点的概率应该与圆的面积与矩形的面积之比成比例；在这里，这是 π / 4。然而，如果我们选择一个非常大的随机点值，注意到我们将得到以下近似值：

这正是我们将估计 π 的方法！在我们能得出合理的 π 估计之前，我们将不得不进行非常多的迭代，但请注意这是多么好的可并行化：我们可以在不同的线程中检查“命中”，将总迭代次数分配给不同的线程。在一天结束时，我们只需将所有线程中的命中总数相加，即可得到我们的估计值。

现在，我们可以开始编写一个程序来进行蒙特卡洛估计。让我们首先导入我们在 PyCUDA 程序中需要的常规 Python 模块，再加上 SymPy 中的一个模块：

SymPy 用于在 Python 中进行完美的 符号 计算，因此当我们有非常大的整数时，我们可以使用 Rational 函数来对除法进行更准确的浮点估计。

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

现在，当我们构建内核时，我们必须做一些与正常情况不同的事情：我们需要在 SourceModule 中设置选项 no_extern_c=True。这会修改代码的编译方式，以便我们的代码可以正确地与 cuRAND 库所需的 C++ 代码链接。然后我们开始编写我们的内核并包含适当的头文件：

ker = SourceModule(no_extern_c=True, source='''
#include <curand_kernel.h>

现在，让我们包含一个用于勾股定理距离的宏。由于我们只是检查这个值是否等于或小于 1，因此可以省略平方根。我们将使用大量的无符号 64 位整数，因此让我们再定义一个宏，以免我们一遍又一遍地输入 unsigned long long：

#define _PYTHAG(a,b) (a*a + b*b)
#define ULL unsigned long long

现在，我们可以设置我们的内核。根据 PyCUDA 的性质，这将不得不编译为接口的真正的 C 函数，而不是 C++ 函数。我们可以通过 extern "C" 块来实现这一点：

extern "C" {

现在，我们可以定义我们的内核。我们将有两个参数：一个是 iters，它是每个线程的总迭代次数，另一个是一个数组，将保存每个线程的命中总数。我们将需要一个 curandState 对象：

__global__ void estimate_pi(ULL iters, ULL * hits)
{
    curandState cr_state;

让我们用一个名为 tid 的整数来保存全局线程 ID：

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

clock()是一个设备函数，输出当前时间到毫秒。我们可以将tid添加到clock()的输出中，以获得每个线程的唯一种子。我们不需要使用不同的子序列或偏移量，所以让我们都设置为 0。我们还将在这里仔细地将所有内容强制转换为 64 位无符号整数：

curand_init( (ULL) clock() + (ULL) tid, (ULL) 0, (ULL) 0, &cr_state);

让我们设置x和y值以保存矩形中的随机点：

float x, y;

然后我们将迭代iters次，看看我们在圆中得到了多少次命中。我们使用curand_uniform(&cr_state)生成这些。请注意，我们可以在 0 到 1 之间生成它们，而不是从-1 到 1，因为在_PYTHAG宏中对它们进行平方运算将消除任何负值：

for(ULL i=0; i < iters; i++)
 {
     x = curand_uniform(&cr_state);
     y = curand_uniform(&cr_state);

     if(_PYTHAG(x,y) <= 1.0f)
         hits[tid]++;
 }

我们现在可以结束并关闭我们的内核，以及extern "C"块，最后用另一个}括号结束：

return;
}
}
''')

现在，让我们用get_function将 Python 包装函数传递给我们的内核。我们还将设置块和网格大小：每个块 32 个线程，每个网格 512 个块。让我们计算总线程数，并在 GPU 上设置一个数组来保存所有的命中（当然初始化为 0）：

pi_ker = ker.get_function("estimate_pi")
threads_per_block = 32
blocks_per_grid = 512
total_threads = threads_per_block * blocks_per_grid
hits_d = gpuarray.zeros((total_threads,),dtype=np.uint64)

让我们设置每个线程的迭代总数为 2²⁴：

iters = 2**24

我们现在可以像往常一样启动内核：

pi_ker(np.uint64(iters), hits_d, grid=(blocks_per_grid,1,1), block=(threads_per_block,1,1))

现在，让我们对数组中的命中次数求和，这给我们了总命中次数。让我们还计算数组中所有线程的总迭代次数：

total_hits = np.sum( hits_d.get() )
total = np.uint64(total_threads) * np.uint64(iters)

我们现在可以用Rational进行估计，就像这样：

est_pi_symbolic =  Rational(4)*Rational(int(total_hits), int(total) )

我们现在可以将其转换为浮点值：

est_pi = np.float(est_pi_symbolic.evalf())

让我们检查我们的估计与 NumPy 的常数值numpy.pi：

print "Our Monte Carlo estimate of Pi is : %s" % est_pi
print "NumPy's Pi constant is: %s " % np.pi
print "Our estimate passes NumPy's 'allclose' : %s" % np.allclose(est_pi, np.pi)

我们现在完成了。让我们从 IPython 中运行并检查一下（这个程序也可以在本书的存储库中的Chapter08下的monte_carlo_pi.py文件中找到）。

CUDA 数学 API

现在，我们将看一下CUDA 数学 API。这是一个库，由设备函数组成，类似于标准 C math.h库中的函数，可以从内核中的单个线程调用。这里的一个区别是，单精度和双精度浮点运算被重载，因此如果我们使用sin(x)，其中x是一个浮点数，sin 函数将产生一个 32 位浮点数作为输出，而如果x是一个 64 位双精度浮点数，那么sin的输出也将是一个 64 位值（通常，这是 32 位函数的正确名称，但它在末尾有一个f，比如sinf）。还有其他内在函数。内在函数是内置到 NVIDIA CUDA 硬件中的不太准确但更快的数学函数；通常，它们的名称与原始函数相似，只是在前面加上两个下划线—因此，内在的 32 位 sin 函数是__sinf。

明确积分的简要回顾

现在，我们将在 Python 中使用一些面向对象的编程，设置一个类，我们可以使用蒙特卡洛方法来评估函数的定积分。让我们停下来，谈谈我们的意思：假设我们有一个数学函数（就像你在微积分课上可能看到的那种类型），我们称之为f(x)。当我们在笛卡尔平面上在点a和b之间绘制它时，它可能看起来像这样：

现在，让我们仔细回顾一下定积分的确切含义——让我们将这个图中的第一个灰色区域表示为I，第二个灰色区域表示为II，第三个灰色区域表示为III。请注意，这里的第二个灰色区域是小于零的。这里的f的定积分，从a到b，将是值I - II + III，我们将在数学上表示为。一般来说，从a到b的定积分就是所有在f函数和 x 轴之间的总“正”区域的总和，其中 y > 0，减去所有在f函数和 x 轴之间的总“负”区域的总和，其中 y < 0，位于a和b之间。

有许多方法可以计算或估计两点之间函数的定积分。 在微积分课程中可能见过的一种方法是找到一个封闭形式的解：找到f的反导数F，并计算F(b) - F(a)。 然而，在许多领域，我们将无法找到确切的反导数，而必须以数值方式确定定积分。 这正是蒙特卡罗积分的想法：我们在a和b之间的许多随机点上评估f，然后使用这些点来估计定积分。

使用蒙特卡罗方法计算定积分

我们现在将使用 CUDA Math API 来表示任意数学函数f，同时使用 cuRAND 库来实现蒙特卡罗积分。 我们将使用元编程来实现这一点：我们将使用 Python 从代码模板生成设备函数的代码，这将插入适当的蒙特卡罗核心以进行积分。

这里的想法是它看起来和行为类似于我们在 PyCUDA 中看到的一些元编程工具，比如ElementwiseKernel。

让我们首先将适当的模块导入到我们的新项目中：

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

我们将使用 Python 中称为基于字典的字符串格式化的技巧。 让我们在继续之前花一分钟来了解一下。 假设我们正在编写一段 CUDA C 代码，并且我们不确定是否希望特定的变量集是 float 还是 double；也许看起来像这样：code_string="float x, y; float * z;"。 我们实际上可能希望格式化代码，以便我们可以随时在浮点数和双精度之间切换。 让我们将字符串中所有对float的引用更改为%(precision)s—code_string="%(precision)s x, y; %(precision)s * z;"。 现在，我们可以设置一个适当的字典，它将用double交换%(presision)s，即code_dict = {'precision' : 'double'}，并使用code_double = code_string % code_dict获取新的双精度字符串。 让我们看一下：

现在，让我们想一想我们想要我们的新蒙特卡洛积分器如何工作。 我们还将使其接受一个使用 CUDA Math API 编写的数学方程的字符串，以定义我们想要积分的函数。 然后，我们可以使用我们刚学到的字典技巧将这个字符串嵌入到代码中，并使用它来积分任意函数。 我们还将使用模板在float和double精度之间进行切换，根据用户的自由裁量。

我们现在可以开始我们的 CUDA C 代码：

MonteCarloKernelTemplate = '''
#include <curand_kernel.h>

我们将保留以前的无符号 64 位整数宏ULL。 让我们为 x 的倒数（_R）和平方（_P2）定义一些新的宏：

#define ULL unsigned long long
#define _R(z) ( 1.0f / (z) )
#define _P2(z) ( (z) * (z) )

现在，让我们定义一个设备函数，我们的方程字符串将插入其中。 当我们必须从字典中交换文本时，我们将使用math_function值。 我们将有另一个值称为p，表示精度（将是float或double）。 我们将称这个设备函数为f。 我们将在函数的声明中放置一个inline，这将节省我们一些时间，因为从核心调用时会有一些分支：

__device__ inline %(p)s f(%(p)s x)
{
    %(p)s y;
    %(math_function)s;
    return y;
}

现在，让我们考虑一下这将如何工作——我们声明一个名为y的 32 位或 64 位浮点值，调用math_function，然后返回y。 math_function，如果它是对输入参数x进行操作并将某个值设置为y的一些代码，那么这只有意义。 让我们记住这一点，然后继续。

我们现在将开始编写我们的蒙特卡洛积分核。 让我们记住，我们必须使用extern "C"关键字使我们的 CUDA 核可从普通 C 中访问。 然后我们将设置我们的核心。

首先，我们将用iters指示内核中每个线程应该采样多少随机样本；然后我们用lo指示积分的下界（b），用hi指示上界（a），并传入一个数组ys_out来存储每个线程的部分积分集合（我们稍后将对ys_out求和，以得到从主机端的lo到hi的完整定积分值）。再次注意我们将精度称为p：

extern "C" {
__global__ void monte_carlo(int iters, %(p)s lo, %(p)s hi, %(p)s * ys_out)
{

我们将需要一个curandState对象来生成随机值。我们还需要找到全局线程 ID 和线程的总数。由于我们正在处理一维数学函数，因此在一维x中设置我们的块和网格参数是有意义的：

curandState cr_state;
int tid = blockIdx.x * blockDim.x + threadIdx.x;
int num_threads = blockDim.x * gridDim.x;

现在我们将计算单个线程将处理的lo和hi之间的面积量。我们将通过将整个积分的长度（即hi - lo）除以线程的总数来实现这一点。:

再次注意我们如何使用模板技巧，以便这个值可以是多精度的。

%(p)s t_width = (hi - lo) / ( %(p)s ) num_threads;

回想一下，我们有一个名为iters的参数；这表示每个线程将采样多少个随机值。我们需要稍后知道样本的密度是多少；也就是说，每单位距离的平均样本数。我们可以这样计算，记得将整数iters强制转换为浮点值：

%(p)s density = ( ( %(p)s ) iters ) / t_width;

回想一下，我们正在将我们正在积分的区域按线程数进行划分。这意味着每个线程将有自己的起点和终点。由于我们正在为每个线程公平地划分长度，我们可以这样计算：

%(p)s t_lo = t_width*tid + lo;
 %(p)s t_hi = t_lo + t_width;

现在我们可以像之前一样初始化 cuRAND，确保每个线程都从自己的个体种子生成随机值：

curand_init( (ULL)  clock() + (ULL) tid, (ULL) 0, (ULL) 0, &cr_state);

在开始采样之前，我们需要设置一些额外的浮点值。y将保存从t_lo到t_hi的积分估计的最终值，y_sum将保存所有采样值的总和。我们还将使用rand_val变量来保存我们生成的原始随机值，x来存储我们将要从中采样的区域的缩放随机值：

%(p)s y, y_sum = 0.0f;
%(p)s rand_val, x;

现在，让我们循环从我们的函数中取样值，将这些值加入y_sum中。需要注意的一点是curand_uniform末尾的%(p_curand)s——这个函数的 32 位浮点版本是curand_uniform，而 64 位版本是curand_uniform_double。稍后我们将根据这里使用的精度级别，用_double或空字符串来替换它。还要注意我们如何缩放rand_val，使得x落在t_lo和t_hi之间，记住 cuRAND 中的随机均匀分布只产生 0 到 1 之间的值：

for (int i=0; i < iters; i++)
{
    rand_val = curand_uniform%(p_curand)s(&cr_state);
    x = t_lo + t_width * rand_val;
    y_sum += f(x);
}

现在我们可以通过密度将y_sum除以t_hi到t_lo的子积分的值：

y = y_sum / density;

我们将这个值输出到数组中，并关闭我们的 CUDA 内核，以及extern "C"，以及最终的闭括号。我们已经写完了 CUDA C，所以我们将用三个引号结束这一部分：

ys_out[tid] = y;
}
}
'''

现在我们要做一些不同的事情——我们将设置一个类来处理我们的定积分。让我们称之为MonteCarloIntegrator。当然，我们将首先编写构造函数，也就是__init__函数。这是我们将输入对象引用self的地方。让我们将math_function的默认值设置为'y = sin(x)'，默认精度为'd'，即双精度。我们还将把lo的默认值设置为 0，hi的默认值设置为π的 NumPy 近似值。最后，我们将有每个线程将采样的随机样本数（samples_per_thread）和我们将在其上启动内核的网格大小（num_blocks）的值。

让我们从将文本字符串math_function存储在self对象中开始这个函数，以便以后使用：

def __init__(self, math_function='y = sin(x)', precision='d', lo=0, hi=np.pi, samples_per_thread=10**5, num_blocks=100):

        self.math_function = math_function

现在，让我们设置与我们选择的浮点精度相关的值，这将在以后设置我们的模板字典时特别需要，特别是为了存储lo和hi的值。让我们还在对象中存储lo和hi的值。如果用户输入无效的数据类型，或者hi实际上小于lo，让我们确保引发异常错误：

         if precision in [None, 's', 'S', 'single', np.float32]:
             self.precision = 'float'
             self.numpy_precision = np.float32
             self.p_curand = ''
         elif precision in ['d','D', 'double', np.float64]:
             self.precision = 'double'
             self.numpy_precision = np.float64
             self.p_curand = '_double'
         else:
             raise Exception('precision is invalid datatype!')

     if (hi - lo <= 0):
         raise Exception('hi - lo <= 0!')
     else:
         self.hi = hi
         self.lo = lo

现在，我们可以设置我们的代码模板字典：

MonteCarloDict = {'p' : self.precision, 'p_curand' : self.p_curand, 'math_function' : self.math_function}

现在我们可以使用基于字典的字符串格式生成实际的最终代码，并进行编译。让我们还通过在SourceModule中设置options=['-w']来关闭nvcc编译器的警告：

self.MonteCarloCode = MonteCarloKernelTemplate % MonteCarloDict

self.ker = SourceModule(no_extern_c=True , options=['-w'], source=self.MonteCarloCode)

现在，我们将在我们的对象中设置一个函数引用到我们编译的内核，使用get_function。在继续之前，让我们保存对象中的剩余两个参数：

self.f = self.ker.get_function('monte_carlo')
self.num_blocks = num_blocks
self.samples_per_thread = samples_per_thread

现在，虽然我们需要不同实例的MonteCarloIntegrator对象来评估不同数学函数或浮点精度的定积分，但我们可能希望在不同的lo和hi边界上评估相同的积分，改变线程/网格大小，或者改变我们在每个线程上取样的数量。幸运的是，这些都是容易进行的修改，并且都可以在运行时进行。

我们将设置一个特定的函数来评估给定对象的积分。我们将这些参数的默认值设置为在调用构造函数时存储的值：

def definite_integral(self, lo=None, hi=None, samples_per_thread=None, num_blocks=None):
    if lo is None or hi is None:
        lo = self.lo
        hi = self.hi
    if samples_per_thread is None:
        samples_per_thread = self.samples_per_thread
    if num_blocks is None:
        num_blocks = self.num_blocks
        grid = (num_blocks,1,1)
    else:
        grid = (num_blocks,1,1)

    block = (32,1,1)
    num_threads = 32*num_blocks

我们可以通过设置一个空数组来存储部分子积分并启动内核来完成这个函数。然后我们需要对子积分求和以获得最终值，然后返回：

self.ys = gpuarray.empty((num_threads,) , dtype=self.numpy_precision)

self.f(np.int32(samples_per_thread), self.numpy_precision(lo), self.numpy_precision(hi), self.ys, block=block, grid=grid)

self.nintegral = np.sum(self.ys.get() )

return np.sum(self.nintegral)

我们已经准备好尝试这个了。让我们只是设置一个具有默认值的类——这将从 0 到π积分y = sin(x)。如果您记得微积分，sin(x)的反导数是-cos(x)，所以我们可以这样评估定积分：

因此，我们应该得到一个接近 2 的数值。让我们看看我们得到了什么：

编写一些测试用例

现在，我们终于将看到如何使用 CUDA Math API 通过math_function参数编写一些测试用例来测试我们的类。如果您有 C/C++标准数学库的经验，这将会相当简单。同样，这些函数是重载的，这样当我们在单精度和双精度之间切换时，我们就不必更改任何名称。

我们已经看到了一个例子，即y = sin(x)。让我们尝试一些更有雄心的东西：

我们将从a=11.733 积分到b=18.472，然后检查我们的蒙特卡洛积分器的输出与另一个来源的已知值进行比较。在这里，Mathematica 指出这个定积分的值是 8.9999，所以我们将与其进行比较。

现在，让我们考虑如何表示这个函数：这里，log指的是自然对数（也称为ln），在 Math API 中就是log(x)。我们已经设置了一个宏来表示平方，所以我们可以用_P2(sin(x))来表示sin²(x)。现在我们可以用y = log(x)*_P2(sin(x))来表示整个函数。

让我们使用以下方程，从a=.9积分到b=4：

记住，_R是我们设置的倒数宏，我们可以这样用 Math API 来写这个函数：

'y = _R( 1 + sinh(2*x)*_P2(log(x)) )' 

在我们继续之前，让我们注意到 Mathematica 告诉我们这个定积分的值是.584977。

让我们再检查一个函数。让我们有点雄心勃勃地说它是这样的：

我们可以将这表示为'y = (cosh(x)*sin(x))/ sqrt( pow(x,3) + _P2(sin(x)))'；自然地，sqrt是分母中的平方根，pow允许我们取任意幂的值。当然，sin(x)是sin(x)，cosh(x)是cosh(x)。我们从a=1.85 积分到b=4.81；Mathematica 告诉我们这个积分的真实值是-3.34553。

我们现在准备检查一些测试用例，并验证我们的蒙特卡洛积分是否有效！让我们迭代一个列表，其第一个值是指示函数（使用 Math API）的字符串，第二个值指示积分的下限，第三个值指示积分的上限，最后一个值指示用 Mathematica 计算出的预期值：

if __name__ == '__main__':

    integral_tests = [('y =log(x)*_P2(sin(x))', 11.733 , 18.472, 8.9999), ('y = _R( 1 + sinh(2*x)*_P2(log(x)) )', .9, 4, .584977), ('y = (cosh(x)*sin(x))/ sqrt( pow(x,3) + _P2(sin(x)))', 1.85, 4.81, -3.34553) ]

我们现在可以迭代这个列表，看看我们的算法与 Mathematica 相比效果如何：

for f, lo, hi, expected in integral_tests:
    mci = MonteCarloIntegrator(math_function=f, precision='d', lo=lo, hi=hi)
    print 'The Monte Carlo numerical integration of the function\n \t f: x -> %s \n \t from x = %s to x = %s is : %s ' % (f, lo, hi, mci.definite_integral())
    print 'where the expected value is : %s\n' % expected

现在就运行这个：

这也可以在本书存储库的Chapter08目录下的monte_carlo_integrator.py文件中找到。

CUDA Thrust 库

我们现在将看一下 CUDA Thrust 库。这个库的核心特性是一个高级向量容器，类似于 C++自己的向量容器。虽然这听起来可能很琐碎，但这将使我们在 CUDA C 中编程时更少地依赖指针、malloc 和 free。与 C++向量容器一样，Thrust 的向量容器会自动处理元素的调整大小和连接，并且借助 C++析构函数的魔力，释放也会在 Thrust 向量对象超出范围时自动处理。

Thrust 实际上提供了两个向量容器：一个用于主机端，一个用于设备端。主机端的 Thrust 向量与 STL 向量几乎相同，主要区别在于它可以更轻松地与 GPU 交互。让我们用适当的 CUDA C 代码写一点代码，以了解它是如何工作的。

让我们从包含语句开始。我们将使用主机和设备端向量的头文件，并且还将包括 C++的iostream库，这将允许我们在终端上执行基本的 I/O 操作：

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <iostream>

让我们只使用标准的 C++命名空间（这样我们在检查输出时就不必输入std::分辨率运算符）：

using namespace std;

我们现在将制作我们的主函数，并在主机端设置一个空的 Thrust 向量。同样，这些都是 C++模板，所以我们必须在声明时使用< >括号选择数据类型。我们将把这设置为一个整数数组：

int main(void)
{
 thrust::host_vector<int> v;

现在，让我们通过使用push_back在v的末尾添加一些整数，就像我们用常规 STL 向量一样：

v.push_back(1);
v.push_back(2);
v.push_back(3);
v.push_back(4);

我们现在将迭代向量中的所有值，并输出每个值：

这里的输出应该是v[0] == 1到v[3] == 4。

for (int i = 0; i < v.size(); i++)
    cout << "v[" << i << "] == " << v[i] << endl;

到目前为止，这可能看起来很琐碎。让我们在 GPU 上设置一个 Thrust 向量，然后将内容从v复制过去：

thrust::device_vector<int> v_gpu = v;

是的，就这样了——只有一行，我们就完成了。现在主机上的v的所有内容都将被复制到设备上的v_gpu！（如果这让你感到惊讶，请再看一下第六章，调试和分析您的 CUDA 代码，想想在这之前我们需要多少行。）

让我们尝试在我们的新 GPU 向量上使用push_back，看看我们是否可以将另一个值连接到它上面：

v_gpu.push_back(5);

我们现在将检查v_gpu的内容，如下所示：

for (int i = 0; i < v_gpu.size(); i++)
    std::cout << "v_gpu[" << i << "] == " << v_gpu[i] << std::endl;

这部分应该输出v_gpu[0] == 1到v_gpu[4] == 5。

再次感谢这些对象的析构函数，我们不必进行任何清理工作，释放任何已分配内存的块。现在我们可以从程序中返回，我们完成了：

    return 0;
}

在 Thrust 中使用函数对象

让我们看看如何在 Thrust 中使用称为functors的概念。在 C++中，functor是一个看起来和行为像函数的类或结构对象；这让我们可以使用看起来和行为像函数的东西，但可以保存一些不必每次使用时都设置的参数。

让我们用适当的包含语句开始一个新的 Thrust 程序，并使用标准命名空间：

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <iostream>
using namespace std;

现在，让我们设置一个基本的函数对象。我们将使用struct来表示这个函数对象，而不是class。这将是一个加权乘法函数，我们将在一个名为w的浮点数中存储权重。我们将创建一个构造函数，用默认值1设置权重：

struct multiply_functor {
 float w;
 multiply_functor(float _w = 1) : w(_w) {}

现在，我们将使用operator()关键字设置我们的函数对象；这将告诉编译器将以下代码块视为此类型对象的default函数。请记住，这将在 GPU 上作为设备函数运行，因此我们在整个代码块之前加上__device__。我们用括号指示输入，并输出适当的值，这只是一个缩放的倍数。现在，我们可以用};关闭我们的结构的定义：

    __device__ float operator() (const float & x, const float & y) { 
        return w * x * y;
     }
};

现在，让我们使用这个来制作一个基本的点积函数；回想一下，这需要两个数组之间的逐点乘法，然后是一个reduce类型的求和。让我们首先声明我们的函数，并创建一个新的向量z，它将保存逐点乘法的值：

float dot_product(thrust::device_vector<float> &v, thrust::device_vector<float> &w ), thrust::device_vector<float> &z)
{
 thrust::device_vector<float> z(v.size());

我们现在将使用 Thrust 的transform操作，它将对v和w的输入进行逐点操作，并输出到z。请注意，我们将函数对象输入到 transform 的最后一个槽中；通过这样使用普通的闭括号，它将使用构造函数的默认值（w = 1），因此这将作为一个普通的、非加权/缩放的点积：

thrust::transform(v.begin(), v.end(), w.begin(), z.begin(), multiply_functor());

我们现在可以使用 Thrust 的 reduce 函数对z进行求和。让我们返回这个值：

return thrust::reduce(z.begin(), z.end());
}

我们完成了。现在，让我们编写一些测试代码——我们将只计算向量[1,2,3]和[1,1,1]的点积，这对我们来说很容易检查。（结果将是 6。）

让我们使用push_back设置第一个向量v：

int main(void)
{
    thrust::device_vector<float> v;
    v.push_back(1.0f);
    v.push_back(2.0f);
    v.push_back(3.0f);

现在，我们可以声明一个大小为3的向量w，并使用 Thrust 的 fill 函数将其默认值设置为1，如下所示：

thrust::device_vector<float> w(3);
thrust::fill(w.begin(), w.end(), 1.0f);

让我们进行一次检查，确保我们的值被正确设置，通过将它们的值输出到cout：

for (int i = 0; i < v.size(); i++)
 cout << "v[" << i << "] == " << v[i] << endl;

for (int i = 0; i < w.size(); i++)
 cout << "w[" << i << "] == " << w[i] << endl;

现在，我们可以检查我们的点积的输出，然后从程序中返回：

cout << "dot_product(v , w) == " << dot_product(v,w) << endl;
return 0;
}

让我们编译这个程序（在 Linux 或 Windows 的命令行中使用nvcc thrust_dot_product.cu -o thrust_dot_product）并运行它：

这个代码也可以在本书存储库的Chapter08目录中的thrust_dot_product.cu文件中找到。

摘要

在本章中，我们看了如何通过选择适当的种子在 cuRAND 中初始化随机数流。由于计算机是确定性设备，它们只能生成伪随机数列表，因此我们的种子应该是真正随机的；通常，将线程 ID 添加到毫秒级的时钟时间中就足够满足大多数目的。

然后，我们看了如何使用 cuRAND 的均匀分布来对 Pi 进行基本估计。然后，我们承担了一个更有雄心的项目，创建一个可以计算任意函数的定积分的 Python 类；我们使用了一些元编程的思想，结合 CUDA Math API 来定义这些任意函数。最后，我们简要介绍了 CUDA Thrust 库，它通常用于在 Python 之外编写纯 CUDA C 程序。Thrust 最显著的提供了一个类似于标准 C++ vector的device_vector容器。这减少了在 CUDA C 中使用指针的一些认知开销。

最后，我们简要介绍了如何使用 Thrust 和适当的函数对象来执行简单的point-wise和reduce操作，即实现一个简单的点积函数。

问题

1. 尝试重写蒙特卡洛积分示例（在monte_carlo_integrator.py的__main__函数中）以使用 CUDA 的intrinsic函数。精度与以前相比如何？

2. 在我们所有的 cuRAND 示例中，我们只使用了均匀分布。你能说出在 GPU 编程中使用正态（高斯）随机分布的一个可能的用途或应用吗？

3. 假设我们使用两个不同的种子生成一个包含 100 个伪随机数的列表。我们应该将它们连接成一个包含 200 个数字的列表吗？

4. 在最后一个示例中，在multiply_functor结构体的operator()函数定义之前添加__host__，现在，看看是否可以直接使用这个函数对象实现一个主机端的点积函数，而不需要任何进一步的修改。

5. 看一下 Thrust examples目录中的strided_range.cu文件。你能想到如何使用这个来使用 Thrust 实现通用矩阵乘法吗？

6. 在定义一个函数对象时，operator()函数的重要性是什么？

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

我们现在将利用我们对 GPU 编程的积累知识，使用 PyCUDA 实现我们自己的深度神经网络（DNN）。在过去的十年里，DNN 吸引了很多关注，因为它们为机器学习（ML）提供了一个强大而优雅的模型。DNN 也是第一个能够展示 GPU 真正强大之处的应用之一（除了渲染图形），通过利用它们的大规模并行吞吐量，最终帮助 NVIDIA 成为人工智能领域的主要参与者。

在本书的过程中，我们大多是按章节覆盖单个主题的——在这里，我们将在我们学到的许多主题的基础上构建我们自己的 DNN 实现。虽然目前有几个面向普通公众的基于 GPU 的 DNN 的开源框架，例如 Google 的 TensorFlow 和 Keras，微软的 CNTK，Facebook 的 Caffe2 和 PyTorch，但从头开始实现一个 DNN 非常有教育意义，这将使我们更深入地了解和欣赏 DNN 所需的基础技术。我们有很多材料要涵盖，所以在简要介绍一些基本概念后，我们将直奔主题。

在本章中，我们将研究以下内容：

• 理解人工神经元（AN）是什么

• 理解如何将多个 AN 组合在一起形成深度神经网络（DNN）

• 在 CUDA 和 Python 中从头开始实现 DNN

• 理解如何使用交叉熵损失来评估神经网络的输出

• 实现梯度下降来训练 NN

• 学习如何在小数据集上训练和测试 NN

技术要求

本章需要一台装有现代 NVIDIA GPU（2016 年以后）的 Linux 或 Windows 10 PC，并安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（9.0 及以上）。还需要一个合适的 Python 2.7 安装（如 Anaconda Python 2.7），并安装了 PyCUDA 模块。

本章的代码也可以在 GitHub 上找到：github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA。

有关本章的先决条件的更多信息，请查看本书的前言。有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA中的 README 文件。

人工神经元和神经网络

让我们简要回顾一些机器学习（ML）和神经网络（NNs）的基础知识。在机器学习中，我们的目标是使用具有特定标记类别或特征的数据集，并利用这些示例来训练我们的系统以预测未来数据的值。我们称根据先前训练数据预测未来数据的类别或标签的程序或函数为分类器。

有许多类型的分类器，但在这里我们将专注于 NNs。NNs 的理念是它们（据说）以类似于人脑的方式工作，通过使用一组人工神经元（ANs）来学习和分类数据，所有这些神经元连接在一起形成特定的结构。不过，让我们暂停一下，看看一个单独的 AN 是什么。在数学上，这只是从线性空间R^n到R的仿射函数，如下所示：

我们可以看到这可以被描述为一个常量权重向量w和输入向量x之间的点积，最后加上一个额外的偏置常量b。（再次强调，这个函数的唯一输入是x；其他值都是常数！）

现在，单个 AN 本身是相当无用（而且愚蠢）的，只有当它们与大量其他 AN 合作时，它们的智能才会显现出来。我们的第一步是将一系列相似的 AN 堆叠在一起，以形成我们将称之为密集层（DL）的东西。这是密集的，因为每个神经元将处理来自x的每个输入值 - 每个 AN 将接收来自R^{n**的数组或向量值，并在**R**中输出一个值。由于有*m*个神经元，这意味着它们的输出集体位于**R}m空间中。我们将注意到，如果我们堆叠我们层中每个神经元的权重，以形成一个m x n的权重矩阵，然后我们可以通过矩阵乘法计算每个神经元的输出，然后加上适当的偏差：

现在，假设我们想要构建一个能够对k个不同类别进行分类的 NN 分类器；我们可以创建一个新的附加密集层，该层接收来自先前密集层的m个值，并输出k个值。假设我们对每一层都有适当的权重和偏差值（这显然不容易找到），并且在每一层之后也有适当的激活函数设置（我们稍后会定义），这将作为我们k个不同类别之间的分类器，根据最终层的输出给出x落入每个类别的概率。当然，我们在这里走得太远了，但这就是 NN 的工作原理。

现在，似乎我们可以将密集层连接到长链中以实现分类。这就是所谓的 DNN。当我们有一层不直接连接到输入或输出时，这就是一个隐藏层。DNN 的优势在于额外的层允许 NN 捕捉浅层 NN 无法捕捉到的数据的抽象和细微差别。

实现人工神经元的密集层

现在，让我们实现 NN 最重要的构建模块，密集层。让我们从声明 CUDA 内核开始，就像这样：

__global__ void dense_eval(int num_outputs, int num_inputs, int relu, int sigmoid, float * w, float * b, float * x, float *y, int batch_size, int w_t, int b_t, float delta)

让我们逐个检查输入。num_outputs当然表示这一层的总输出数量；这正好是该层的神经元数量。num_inputs告诉我们输入数据的大小。为relu和sigmoid设置正值将表示我们应该在该层的输出上使用相应的激活函数，我们稍后会定义。w和b是包含该层权重和偏差的数组，而x和y将作为我们的输入和输出。通常，我们希望一次对多个数据进行分类。我们可以通过将batch_size设置为我们希望预测的点数来指示这一点。最后，w_t、b_t和delta将在训练过程中使用，通过梯度下降来确定该层的适当权重和偏差。（我们将在后面的部分中更多地了解梯度下降。）

现在，让我们开始编写我们的内核。我们将在每个输出上并行计算，因此我们将设置一个整数i作为全局线程 ID，然后使用适当的if语句让任何不必要的额外线程不执行任何操作：

{
 int i = blockDim.x*blockIdx.x + threadIdx.x;

 if (i < num_outputs)
 {

现在，让我们用适当的for循环迭代批处理中的每个数据点：

for(int k=0; k < batch_size; k++)
 { 

我们将从权重和输入中的 32 位浮点数中进行乘法和累加，得到 64 位双精度temp，然后加上适当的偏差点。然后我们将这个值强制转换回 32 位浮点数并将值放入输出数组中，然后关闭对k的循环：

double temp = 0.0f;
 for (int j = 0; j < num_inputs; j++)
 {
   temp += ((double) w[(num_inputs)*i + j ] ) * ( (double) x[k*num_inputs + j]);
 }
 temp += (double) b[i];
 y[k * num_outputs + i] = (float) temp;  
}

乘法和累加类型的操作通常会导致严重的数值精度损失。这可以通过使用更高精度的临时变量来存储操作过程中的值，然后在操作完成后将此变量强制转换回原始精度来减轻。

要训练一个 NN，我们最终将不得不计算我们的 NN 对每个权重和偏差的导数（来自微积分），这是针对特定输入批次的。请记住，数学函数f在值x处的导数可以估计为f**(x + δ) - f(x) / δ，其中 delta（δ）是某个足够小的正值。我们将使用输入值w_t和b_t来指示内核是否要计算相对于特定权重或偏差的导数；否则，我们将将这些输入值设置为负值，仅对此层进行评估。我们还将设置 delta 为适当小的值，用于计算导数，并使用它来增加适当偏差或权重的值：

if( w_t >= 0 && i == (w_t / num_inputs))
 {
 int j = w_t % num_inputs;
 for(int k=0; k < batch_size; k++)
  y[k*num_outputs + i] += delta*x[k*num_inputs+j];
}
if( b_t >= 0 && i == b_t )
 {
  for(int k=0; k < batch_size; k++)
  y[k*num_outputs + i] += delta;
 }

现在，我们将添加一些代码，用于所谓的修正线性单元（或ReLU）和Sigmoid 激活函数。这些用于处理密集神经层的直接输出。ReLU 只将所有负值设为 0，同时作为正输入的恒等式，而 sigmoid 只计算每个值上的sigmoid函数的值（1 / (1 + e^(-x))）。ReLU（或任何其他激活函数）用于 NN 中隐藏层之间，作为使整个 NN 成为非线性函数的手段；否则，整个 NN 将构成一个微不足道（且计算效率低下）的矩阵操作。（虽然可以在层之间使用许多其他非线性激活函数，但发现 ReLU 对训练来说是一个特别有效的函数。）Sigmoid 用作 NN 中用于标签的最终层，即可能为给定输入分配多个标签的层，而不是将输入分配给单个类别。

让我们在文件中稍微上移一点，甚至在我们开始定义这个 CUDA 内核之前，将这些操作定义为 C 宏。我们还要记得在我们写完 CUDA-C 代码后将其放入其中：

DenseEvalCode = '''
#define _RELU(x) ( ((x) > 0.0f) ? (x) : 0.0f )
#define _SIGMOID(x) ( 1.0f / (1.0f + expf(-(x)) ))

现在，我们将使用内核输入relu和sigmoid来指示我们是否应该使用这些额外的层；我们将从这些中获取积极的输入，以指示它们应该分别使用。我们可以添加这个，关闭我们的内核，并将其编译成可用的 Python 函数：

if(relu > 0 || sigmoid > 0)
for(int k=0; k < batch_size; k++)
 { 
   float temp = y[k * num_outputs + i];
   if (relu > 0)
    temp = _RELU(temp);
   if (sigmoid > 0)
    temp = _SIGMOID(temp);
   y[k * num_outputs + i] = temp; 
  }
 }
 return;
}
'''
eval_mod = SourceModule(DenseEvalCode)
eval_ker = eval_mod.get_function('dense_eval')

现在，让我们回到文件的开头，并设置适当的导入语句。请注意，我们将包括csv模块，该模块将用于处理测试和训练的数据输入：

from __future__ import division
import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
from pycuda.elementwise import ElementwiseKernel
import numpy as np
from Queue import Queue
import csv
import time

现在，让我们继续设置我们的密集层；我们将希望将其包装在一个 Python 类中以便使用，这将使我们在开始将这些密集层连接成一个完整的 NN 时更加轻松。我们将称之为class DenseLayer并开始编写构造函数。这里的大部分输入和设置应该是不言自明的：我们绝对应该添加一个选项，从预训练的网络中加载权重和偏差，并且我们还将包括指定默认delta值以及默认流的选项。（如果没有给出权重或偏差，权重将被初始化为随机值，而所有偏差将设置为 0。）我们还将在这里指定是否使用 ReLU 或 sigmoid 层。最后，请注意我们如何设置块和网格大小：

class DenseLayer:
    def __init__(self, num_inputs=None, num_outputs=None, weights=None, b=None, stream=None, relu=False, sigmoid=False, delta=None):
        self.stream = stream

        if delta is None:
            self.delta = np.float32(0.001)
        else:
            self.delta = np.float32(delta)

        if weights is None:
            weights = np.random.rand(num_outputs, num_inputs) - .5
            self.num_inputs = np.int32(num_inputs)
        self.num_outputs = np.int32(num_outputs) 

        if type(weights) != pycuda.gpuarray.GPUArray:
            self.weights = gpuarray.to_gpu_async(np.array(weights, 
            dtype=np.float32) , stream = self.stream)
        else:
            self.weights = weights

        if num_inputs is None or num_outputs is None:
            self.num_inputs = np.int32(self.weights.shape[1])
            self.num_outputs = np.int32(self.weights.shape[0])

        else:
            self.num_inputs = np.int32(num_inputs)
            self.num_outputs = np.int32(num_outputs)

        if b is None:
            b = gpuarray.zeros((self.num_outputs,),dtype=np.float32)

        if type(b) != pycuda.gpuarray.GPUArray:
            self.b = gpuarray.to_gpu_async(np.array(b, 
            dtype=np.float32) , stream = self.stream)
        else:
            self.b = b 

        self.relu = np.int32(relu)
        self.sigmoid = np.int32(sigmoid)

        self.block = (32,1,1)
        self.grid = (int(np.ceil(self.num_outputs / 32)), 1,1)

现在，我们将在这个类中设置一个函数来评估来自这一层的输入；我们将仔细检查输入（x）以确定它是否已在 GPU 上（如果没有，则将其转移到gpuarray），并让用户指定预分配的gpuarray用于输出（y），如果没有指定，则手动分配一个输出数组。我们还将检查训练时的 delta 和w_t/b_t值，以及batch_size。然后我们将在x输入上运行核函数，输出到y，最后将y作为输出值返回：

def eval_(self, x, y=None, batch_size=None, stream=None, delta=None, w_t = None, b_t = None):

if stream is None:
    stream = self.stream

if type(x) != pycuda.gpuarray.GPUArray:
    x = gpuarray.to_gpu_async(np.array(x,dtype=np.float32), stream=self.stream)

if batch_size is None:
    if len(x.shape) == 2:
        batch_size = np.int32(x.shape[0])
    else:
        batch_size = np.int32(1)

if delta is None:
    delta = self.delta

delta = np.float32(delta)

if w_t is None:
    w_t = np.int32(-1)

if b_t is None:
    b_t = np.int32(-1)

if y is None:
    if batch_size == 1:
        y = gpuarray.empty((self.num_outputs,), dtype=np.float32)
    else:
        y = gpuarray.empty((batch_size, self.num_outputs), dtype=np.float32)

    eval_ker(self.num_outputs, self.num_inputs, self.relu, self.sigmoid, self.weights, self.b, x, y, np.int32(batch_size), w_t, b_t, delta , block=self.block, grid=self.grid , stream=stream)

 return y

现在，我们已经完全实现了一个密集层！

实现 softmax 层

现在，让我们看看如何实现softmax 层。正如我们已经讨论过的，sigmoid 层用于为类分配标签——也就是说，如果您想要从输入中推断出多个非排他性特征，您应该使用 sigmoid 层。softmax 层用于仅通过推断为样本分配单个类别——这是通过计算每个可能类别的概率来实现的（当然，所有类别的概率总和为 100%）。然后我们可以选择具有最高概率的类别来给出最终分类。

现在，让我们看看 softmax 层到底做了什么——给定一组N个实数(c[0], ..., c[N-1])，我们首先计算每个数字的指数函数的总和，然后计算每个数字除以这个总和的指数，得到 softmax：

让我们从我们的实现开始。我们将首先编写两个非常简短的 CUDA 内核：一个用于计算每个输入的指数，另一个用于计算所有点的平均值：

SoftmaxExpCode='''
__global__ void softmax_exp( int num, float *x, float *y, int batch_size)
{
 int i = blockIdx.x * blockDim.x + threadIdx.x;

 if (i < num)
 {
  for (int k=0; k < batch_size; k++)
  {
   y[num*k + i] = expf(x[num*k+i]);
  }
 }
}
'''
exp_mod = SourceModule(SoftmaxExpCode)
exp_ker = exp_mod.get_function('softmax_exp')

SoftmaxMeanCode='''
__global__ void softmax_mean( int num, float *x, float *y, int batch_size)
{
 int i = blockDim.x*blockIdx.x + threadIdx.x;

 if (i < batch_size)
 {
  float temp = 0.0f;

  for(int k=0; k < num; k++)
   temp += x[i*num + k];

  for(int k=0; k < num; k++)
   y[i*num+k] = x[i*num+k] / temp;
 }

 return;
}'''

mean_mod = SourceModule(SoftmaxMeanCode)
mean_ker = mean_mod.get_function('softmax_mean')

现在，让我们编写一个 Python 包装类，就像我们之前做的那样。首先，我们将从构造函数开始，并使用num指示输入和输出的数量。如果需要，我们还可以指定默认流：

class SoftmaxLayer:
    def __init__(self, num=None, stream=None):
     self.num = np.int32(num)
     self.stream = stream

现在，让我们编写一个类似于密集层的eval_函数：

def eval_(self, x, y=None, batch_size=None, stream=None):
 if stream is None:
 stream = self.stream

 if type(x) != pycuda.gpuarray.GPUArray:
  temp = np.array(x,dtype=np.float32)
  x = gpuarray.to_gpu_async( temp , stream=stream)

 if batch_size==None:
  if len(x.shape) == 2:
   batch_size = np.int32(x.shape[0])
  else:
   batch_size = np.int32(1)
 else:
  batch_size = np.int32(batch_size)

 if y is None:
  if batch_size == 1:
   y = gpuarray.empty((self.num,), dtype=np.float32)
 else:
  y = gpuarray.empty((batch_size, self.num), dtype=np.float32)

 exp_ker(self.num, x, y, batch_size, block=(32,1,1), grid=(int( np.ceil( self.num / 32) ), 1, 1), stream=stream)

 mean_ker(self.num, y, y, batch_size, block=(32,1,1), grid=(int( np.ceil( batch_size / 32)), 1,1), stream=stream)

 return y

实现交叉熵损失

现在，让我们实现所谓的交叉熵损失函数。这用于在训练过程中测量神经网络在数据子集上的准确性；损失函数输出的值越大，我们的神经网络在正确分类给定数据方面就越不准确。我们通过计算期望输出和实际输出之间的标准平均对数熵差来实现这一点。为了数值稳定性，我们将限制输出值为1：

MAX_ENTROPY = 1

def cross_entropy(predictions=None, ground_truth=None):

 if predictions is None or ground_truth is None:
  raise Exception("Error! Both predictions and ground truth must be float32 arrays")

 p = np.array(predictions).copy()
 y = np.array(ground_truth).copy()

 if p.shape != y.shape:
  raise Exception("Error! Both predictions and ground_truth must have same shape.")

 if len(p.shape) != 2:
  raise Exception("Error! Both predictions and ground_truth must be 2D arrays.")

 total_entropy = 0

 for i in range(p.shape[0]):
  for j in range(p.shape[1]):
   if y[i,j] == 1: 
    total_entropy += min( np.abs( np.nan_to_num( np.log( p[i,j] ) ) ) , MAX_ENTROPY) 
   else: 
    total_entropy += min( np.abs( np.nan_to_num( np.log( 1 - p[i,j] ) ) ), MAX_ENTROPY)

 return total_entropy / p.size

实现一个顺序网络

现在，让我们实现最后一个类，将多个密集层和 softmax 层对象组合成一个连贯的前馈顺序神经网络。这将作为另一个类来实现，它将包含其他类。让我们首先从编写构造函数开始——我们将能够在这里设置最大批处理大小，这将影响为此网络分配多少内存——我们将在列表变量network_mem中存储一些用于权重和每个层的输入/输出的分配内存。我们还将在列表 network 中存储 DenseLayer 和 SoftmaxLayer 对象，并在 network_summary 中存储有关 NN 中每个层的信息。请注意，我们还可以在这里设置一些训练参数，包括 delta，用于梯度下降的流的数量（稍后我们将看到），以及训练时期的数量。

我们还可以看到开始时的另一个输入称为 layers。在这里，我们可以通过描述每个层来指示 NN 的构造，构造函数将通过迭代 layers 的每个元素并调用add_layer方法来创建它，接下来我们将实现这个方法：

class SequentialNetwork:
 def __init__(self, layers=None, delta=None, stream = None, max_batch_size=32, max_streams=10, epochs = 10):

 self.network = []
 self.network_summary = []
 self.network_mem = []

 if stream is not None:
  self.stream = stream
 else:
  self.stream = drv.Stream()

 if delta is None:
  delta = 0.0001

 self.delta = delta
 self.max_batch_size=max_batch_size
 self.max_streams = max_streams
 self.epochs = epochs

 if layers is not None:
  for layer in layers:
   add_layer(self, layer)

现在，让我们实现add_layer方法。我们将使用字典数据类型将有关该层的所有相关信息传递给顺序网络，包括层的类型（密集、softmax 等）、输入/输出的数量、权重和偏差。这将向对象的网络和network_summary列表变量追加适当的对象和信息，并适当地为network_mem列表分配gpuarray对象：

def add_layer(self, layer):
 if layer['type'] == 'dense':
  if len(self.network) == 0:
   num_inputs = layer['num_inputs']
  else:
   num_inputs = self.network_summary[-1][2]

  num_outputs = layer['num_outputs']
  sigmoid = layer['sigmoid']
  relu = layer['relu']
  weights = layer['weights']
  b = layer['bias']

  self.network.append(DenseLayer(num_inputs=num_inputs, num_outputs=num_outputs, sigmoid=sigmoid, relu=relu, weights=weights, b=b))
  self.network_summary.append( ('dense', num_inputs, num_outputs))

  if self.max_batch_size > 1:
   if len(self.network_mem) == 0:
self.network_mem.append(gpuarray.empty((self.max_batch_size, self.network_summary[-1][1]), dtype=np.float32))
 self.network_mem.append(gpuarray.empty((self.max_batch_size, self.network_summary[-1][2] ), dtype=np.float32 ) ) 
 else:
 if len(self.network_mem) == 0:
 self.network_mem.append( gpuarray.empty( (self.network_summary[-1][1], ), dtype=np.float32 ) )
 self.network_mem.append( gpuarray.empty((self.network_summary[-1][2], ), dtype=np.float32 ) ) 

 elif layer['type'] == 'softmax':

  if len(self.network) == 0:
   raise Exception("Error! Softmax layer can't be first!")

  if self.network_summary[-1][0] != 'dense':
   raise Exception("Error! Need a dense layer before a softmax layer!")

  num = self.network_summary[-1][2]
  self.network.append(SoftmaxLayer(num=num))
  self.network_summary.append(('softmax', num, num))

  if self.max_batch_size > 1:
   self.network_mem.append(gpuarray.empty((self.max_batch_size, self.network_summary[-1][2] ), dtype=np.float32)) 
  else:
   self.network_mem.append( gpuarray.empty((self.network_summary[-1][2], ), dtype=np.float32))

推断方法的实现

我们现在将在我们的SequentialNetwork类中添加两种推断方法——即，根据特定输入预测输出的方法。我们将首先称之为predict，这将由最终用户使用。在训练过程中，我们将不得不基于仅部分层的部分结果进行预测，我们将为此制作另一种方法，称为partial_predict。

让我们从实现predict开始。这将接受两个输入——以一维或二维 NumPy 数组形式的样本集，可能还有用户定义的 CUDA 流。我们将从样本（这里称为x）进行一些类型检查和格式化，记住样本将以行方式存储：

def predict(self, x, stream=None):

 if stream is None:
  stream = self.stream

 if type(x) != np.ndarray:
  temp = np.array(x, dtype = np.float32)
  x = temp

 if(x.size == self.network_mem[0].size):
  self.network_mem[0].set_async(x, stream=stream)
 else:

  if x.size > self.network_mem[0].size:
   raise Exception("Error: batch size too large for input.")

  x0 = np.zeros((self.network_mem[0].size,), dtype=np.float32)
  x0[0:x.size] = x.ravel()
  self.network_mem[0].set_async(x0.reshape( self.network_mem[0].shape), stream=stream)

 if(len(x.shape) == 2):
  batch_size = x.shape[0]
 else:
  batch_size = 1

现在，让我们执行实际的推断步骤。我们只需遍历整个神经网络，在每一层上执行eval_：

for i in xrange(len(self.network)):
 self.network[i].eval_(x=self.network_mem[i], y= self.network_mem[i+1], batch_size=batch_size, stream=stream)

现在，我们将提取 NN 的最终输出，GPU，并将其返回给用户。如果x中的样本数量实际上小于最大批处理大小，则在返回之前我们将适当地切片输出数组：

y = self.network_mem[-1].get_async(stream=stream)

if len(y.shape) == 2:
 y = y[0:batch_size, :]

return y

现在，完成了这一步，让我们实现partial_predict。让我们简要讨论一下这个想法。当我们处于训练过程中时，我们将评估一组样本，然后看看如何通过逐个添加delta到每个权重和偏差上会影响输出。为了节省时间，我们可以计算每一层的输出并将它们存储给定的样本集，然后只重新计算我们改变权重的层的输出，以及所有后续层的输出。我们很快就会更深入地了解这个想法，但现在，我们可以这样实现：

def partial_predict(self, layer_index=None, w_t=None, b_t=None, partial_mem=None, stream=None, batch_size=None, delta=None):

 self.network[layer_index].eval_(x=self.network_mem[layer_index], y = partial_mem[layer_index+1], batch_size=batch_size, stream = stream, w_t=w_t, b_t=b_t, delta=delta)

 for i in xrange(layer_index+1, len(self.network)):
  self.network[i].eval_(x=partial_mem[i], y =partial_mem[i+1], batch_size=batch_size, stream = stream)

梯度下降

我们现在将以批量随机梯度下降（BSGD）的形式对我们的 NN 进行训练方法的完整实现。让我们逐字逐句地思考这意味着什么。批量意味着这个训练算法将一次操作一组训练样本，而不是同时处理所有样本，而随机表示每个批次是随机选择的。梯度意味着我们将使用微积分中的梯度，这里是每个权重和偏差对损失函数的导数集合。最后，下降意味着我们试图减少损失函数——我们通过迭代地对权重和偏差进行微小的更改来实现这一点，通过减去梯度。

从微积分中我们知道，一个点的梯度总是指向最大增加的方向，其相反方向是最大减少的方向。因为我们想要减少，所以我们减去梯度。

我们现在将在我们的SequentialNetwork类中实现 BSGD 作为bsgd方法。让我们逐一讨论bsgd的输入参数：

• training将是一个二维 NumPy 数组的训练样本

• labels将是 NN 最终层的期望输出，对应于每个训练样本

• delta将指示我们在计算导数时应该增加权重多少

• max_streams将指示 BSGD 将在其上执行计算的最大并发 CUDA 流的数量

• batch_size将指示我们希望对每次更新权重计算损失函数的批次有多大

• epochs将指示我们对当前样本集的顺序进行多少次洗牌，分成一系列批次，然后进行 BSGD

• training_rate将指示我们使用梯度计算更新权重和偏差的速率

我们将像往常一样开始这个方法，并进行一些检查和类型转换，将 CUDA 流对象的集合设置为 Python 列表，并在另一个列表中分配一些额外需要的 GPU 内存：

def bsgd(self, training=None, labels=None, delta=None, max_streams = None, batch_size = None, epochs = 1, training_rate=0.01):

 training_rate = np.float32(training_rate)

 training = np.float32(training)
 labels = np.float32(labels)

 if( training.shape[0] != labels.shape[0] ):
  raise Exception("Number of training data points should be same as labels!")

 if max_streams is None:
  max_streams = self.max_streams

 if epochs is None:
 epochs = self.epochs

 if delta is None:
 delta = self.delta

 streams = []
 bgd_mem = []

 # create the streams needed for training
 for _ in xrange(max_streams):
  streams.append(drv.Stream())
  bgd_mem.append([])

 # allocate memory for each stream
 for i in xrange(len(bgd_mem)):
  for mem_bank in self.network_mem:
   bgd_mem[i].append( gpuarray.empty_like(mem_bank) )

现在，我们可以开始训练。我们将从每个epoch开始执行整个 BSGD 的迭代，对每个 epoch 对整个数据集进行随机洗牌。我们还将在终端打印一些信息，以便用户在训练过程中获得一些状态更新：

num_points = training.shape[0]

if batch_size is None:
 batch_size = self.max_batch_size

index = range(training.shape[0])

for k in xrange(epochs): 

 print '-----------------------------------------------------------'
 print 'Starting training epoch: %s' % k
 print 'Batch size: %s , Total number of training samples: %s' % (batch_size, num_points)
 print '-----------------------------------------------------------'

 all_grad = []

 np.random.shuffle(index)

现在，我们将循环遍历洗牌数据集中的每个批次。我们首先计算当前批次的熵，然后将其打印出来。如果用户看到熵的减少，那么他们将知道梯度下降在这里起作用：

for r in xrange(int(np.floor(training.shape[0]/batch_size))):

 batch_index = index[r*batch_size:(r+1)*batch_size] 

 batch_training = training[batch_index, :]
 batch_labels = labels[batch_index, :]

 batch_predictions = self.predict(batch_training)

 cur_entropy = cross_entropy(predictions=batch_predictions, ground_truth=batch_labels)

 print 'entropy: %s' % cur_entropy

我们现在将迭代我们的 NN 的每个密集层，计算整套权重和偏差的梯度。我们将把这些导数存储在扁平化（一维）数组中，这将对应于我们的 CUDA 内核中的w_t和b_t索引，它们也是扁平化的。由于我们将有多个流处理不同权重的不同输出，我们将使用 Python 队列容器来存储尚未处理的这一批权重和偏差：然后我们可以从该容器的顶部弹出值到下一个可用流（我们将这些存储为元组，第一个元素指示这是权重还是偏差）：

for i in xrange(len(self.network)):

 if self.network_summary[i][0] != 'dense':
  continue

 all_weights = Queue()

 grad_w = np.zeros((self.network[i].weights.size,), dtype=np.float32)
 grad_b = np.zeros((self.network[i].b.size,), dtype=np.float32)

 for w in xrange( self.network[i].weights.size ):
  all_weights.put( ('w', np.int32(w) ) )

 for b in xrange( self.network[i].b.size ):
  all_weights.put(('b', np.int32(b) ) )

现在，我们需要迭代每一个权重和偏差，我们可以使用while循环来检查我们刚刚设置的queue对象是否为空。我们将设置另一个队列stream_weights，这将帮助我们组织每个流处理的权重和偏差。适当设置权重和偏差输入后，我们现在可以使用partial_predict，使用当前流和相应的 GPU 内存数组：

请注意，我们已经对这批样本执行了predict来计算熵，所以我们现在可以对这批样本执行partial_predict，只要我们小心使用内存和层。

while not all_weights.empty():

 stream_weights = Queue()

 for j in xrange(max_streams):

  if all_weights.empty():
    break

  wb = all_weights.get()

  if wb[0] == 'w':
   w_t = wb[1]
   b_t = None
  elif wb[0] == 'b':
   b_t = wb[1]
   w_t = None

  stream_weights.put( wb )

  self.partial_predict(layer_index=i, w_t=w_t, b_t=b_t, partial_mem=bgd_mem[j], stream=streams[j], batch_size=batch_size, delta=delta)

我们只计算了一小部分权重和偏差的输出预测。我们将不得不为每个计算熵，然后将导数值存储在扁平化的数组中：

for j in xrange(max_streams):

 if stream_weights.empty():
  break

 wb = stream_weights.get()

 w_predictions = bgd_mem[j][-1].get_async(stream=streams[j])

 w_entropy = cross_entropy(predictions=w_predictions[ :batch_size,:], ground_truth=batch_labels)

 if wb[0] == 'w':
  w_t = wb[1]
  grad_w[w_t] = -(w_entropy - cur_entropy) / delta

 elif wb[0] == 'b':
  b_t = wb[1]
  grad_b[b_t] = -(w_entropy - cur_entropy) / delta

我们现在已经完成了while循环。一旦我们到达外部，我们将知道我们已经计算了这个特定层的所有权重和偏差的导数。在迭代到下一层之前，我们将把当前权重和偏差的梯度计算值附加到all_grad列表中。我们还将把扁平化的权重列表重新整形成原始形状：

all_grad.append([np.reshape(grad_w,self.network[i].weights.shape) , grad_b])

在迭代每一层之后，我们可以对这一批的 NN 的权重和偏差进行优化。请注意，如果training_rate变量远小于1，这将减少权重更新的速度：

for i in xrange(len(self.network)):
 if self.network_summary[i][0] == 'dense':
  new_weights = self.network[i].weights.get()
  new_weights += training_rate*all_grad[i][0]
  new_bias = self.network[i].b.get()
  new_bias += training_rate*all_grad[i][1]
  self.network[i].weights.set(new_weights)
  self.network[i].b.set(new_bias)

我们已经完全实现了（非常简单的）基于 GPU 的 DNN！

数据的调整和归一化

在我们继续训练和测试全新的 NN 之前，我们需要退后一步，谈谈数据调整和数据归一化。NN 对数值误差非常敏感，特别是当输入的规模差异很大时。这可以通过正确调整我们的训练数据来减轻，这意味着对于输入样本中的每个点，我们将计算所有样本中每个点的平均值和方差，然后在输入到 NN 进行训练或推断（预测）之前，对每个样本中的每个点减去平均值并除以标准差。这种方法称为归一化。让我们组合一个小的 Python 函数来为我们做这个：

def condition_data(data, means=None, stds=None):

 if means is None:
  means = np.mean(data, axis=0)

 if stds is None:
  stds = np.std(data, axis = 0)

 conditioned_data = data.copy()
 conditioned_data -= means
 conditioned_data /= stds

 return (conditioned_data, means, stds)

鸢尾花数据集

我们现在将为一个真实问题构建我们自己的 DNN：根据花瓣的测量来对花的类型进行分类。我们将使用众所周知的鸢尾花数据集。该数据集存储为逗号分隔值（CSV）文本文件，每行包含四个不同的数值（花瓣测量），然后是花的类型（这里有三个类别——山鸢尾、变色鸢尾和维吉尼亚鸢尾）。我们现在将设计一个小型 DNN，根据这个数据集对鸢尾花的类型进行分类。

在我们继续之前，请下载鸢尾花数据集并将其放入您的工作目录。这可以从 UC Irvine 机器学习存储库中获取，网址为：archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data。

我们将首先将此文件处理成适当的数据数组，以便用于训练和验证我们的 DNN。让我们从打开我们的主函数开始；我们需要将花的名称转换为 DNN 可以输出的实际类别，因此让我们创建一个小字典，为每个类别提供相应的标签。我们还将设置一些空列表来存储我们的训练数据和标签：

if __name__ == '__main__':
 to_class = { 'Iris-setosa' : [1,0,0] , 'Iris-versicolor' : [0,1,0], 'Iris-virginica' : [0,0,1]}

 iris_data = []
 iris_labels = []

现在，让我们从 CSV 文件中读取。我们将使用 Python 的csv模块中的reader函数，这是我们之前导入的：

with open('C:/Users/btuom/examples/9/iris.data', 'rb') as csvfile:
 csvreader = csv.reader(csvfile, delimiter=',')
 for row in csvreader:
  newrow = []
  if len(row) != 5:
   break
  for i in range(4):
   newrow.append(row[i])
  iris_data.append(newrow)
  iris_labels.append(to_class[row[4]])

我们现在将随机洗牌数据，并使用三分之二的样本作为训练数据。剩下的三分之一将用于测试（验证）数据：

iris_len = len(iris_data)
shuffled_index = list(range(iris_len))
np.random.shuffle(shuffled_index)

iris_data = np.float32(iris_data)
iris_labels = np.float32(iris_labels)
iris_data = iris_data[shuffled_index, :]
iris_labels = iris_labels[shuffled_index,:]

t_len = (2*iris_len) // 3

iris_train = iris_data[:t_len, :]
label_train = iris_labels[:t_len, :]

iris_test = iris_data[t_len:,:]
label_test = iris_labels[t_len:, :]

现在，最后，我们可以开始构建我们的 DNN！首先，让我们创建一个SequentialNetwork对象。我们将max_batch_size设置为32：

sn = SequentialNetwork( max_batch_size=32 )

现在，让我们创建我们的 NN。这将包括四个密集层（两个隐藏层）和一个 softmax 层。我们将逐层增加神经元的数量，直到最后一层，最后一层只有三个输出（每个类别一个）。每层神经元数量的增加允许我们捕捉数据的一些细微差别：

sn.add_layer({'type' : 'dense', 'num_inputs' : 4, 'num_outputs' : 10, 'relu': True, 'sigmoid': False, 'weights' : None, 'bias' : None} ) 
sn.add_layer({'type' : 'dense', 'num_inputs' : 10, 'num_outputs' : 15, 'relu': True, 'sigmoid': False, 'weights': None, 'bias' : None} ) 
sn.add_layer({'type' : 'dense', 'num_inputs' : 15, 'num_outputs' : 20, 'relu': True, 'sigmoid': False, 'weights': None, 'bias' : None} ) 
sn.add_layer({'type' : 'dense', 'num_inputs' : 20, 'num_outputs' : 3, 'relu': True, 'sigmoid': False, 'weights': None , 'bias': None } ) 
sn.add_layer({'type' : 'softmax'})

我们现在将调整我们的训练数据，并使用我们刚刚实现的 BSGD 方法进行训练。我们将使用batch_size设置为16，max_streams设置为10，epochs的数量设置为 100，delta设置为 0.0001，training_rate设置为 1——这些将是几乎任何现代 GPU 可接受的参数。我们还将计时训练过程，这可能会非常耗时。

ctrain, means, stds = condition_data(iris_train)

t1 = time()
sn.bsgd(training=ctrain, labels=label_train, batch_size=16, max_streams=20, epochs=100 , delta=0.0001, training_rate=1)
training_time = time() - t1

现在，我们的 DNN 已经完全训练好了。我们准备开始验证过程！让我们设置一个名为hits的 Python 变量来计算正确分类的总数。我们还需要对验证/测试数据进行条件设置。还有一件事——我们通过 DNN 的 softmax 层的最大值对应的索引来确定类别。我们可以使用 NumPy 的argmax函数来检查这是否给出了正确的分类，就像这样：

hits = 0
ctest, _, _ = condition_data(iris_test, means=means, stds=stds)
for i in range(ctest.shape[0]):
 if np.argmax(sn.predict(ctest[i,:])) == np.argmax( label_test[i,:]):
  hits += 1

现在，我们准备检查我们的 DNN 实际工作得有多好。让我们输出准确率以及总训练时间：

print 'Percentage Correct Classifications: %s' % (float(hits ) / ctest.shape[0])
print 'Total Training Time: %s' % training_time

现在，我们完成了。我们现在可以完全用 Python 和 CUDA 实现一个 DNN！一般来说，您可以期望在这个特定问题上的准确率在 80%-97%之间，使用任何 Pascal 级别的 GPU 的训练时间为 10-20 分钟。

本章的代码可在本书的 GitHub 存储库的适当目录下的deep_neural_network.py文件中找到。

总结

在本章中，我们首先给出了人工神经网络的定义，并向您展示了如何将单个 AN 组合成密集层，然后再将其组合成完整的深度神经网络。然后，我们在 CUDA-C 中实现了一个密集层，并制作了一个相应的 Python 包装类。我们还包括了在密集层的输出上添加 ReLU 和 sigmoid 层的功能。我们看到了使用 softmax 层的定义和动机，这用于分类问题，然后在 CUDA-C 和 Python 中实现了这一功能。最后，我们实现了一个 Python 类，以便我们可以从先前的类构建一个顺序前馈 DNN；我们实现了一个交叉熵损失函数，然后在我们的梯度下降实现中使用这个损失函数来训练我们 DNN 中的权重和偏差。最后，我们使用我们的实现在真实数据集上构建、训练和测试了一个 DNN。

现在我们对我们的 CUDA 编程能力有了很大的自信，因为我们可以编写自己基于 GPU 的 DNN！我们现在将在接下来的两章中学习一些非常高级的内容，我们将看看如何编写我们自己的接口到编译后的 CUDA 代码，以及一些关于 NVIDIA GPU 非常技术性的细节。

问题

1. 假设您构建了一个 DNN，并在训练后，它只产生垃圾。经过检查，您发现所有的权重和偏差要么是巨大的数字，要么是 NaN。问题可能是什么？

2. 小training_rate值可能存在的一个问题是什么？

3. 大training_rate值可能存在的一个问题是什么？

4. 假设我们想要训练一个 DNN，将多个标签分配给动物的图像（“黏滑的”，“毛茸茸的”，“红色的”，“棕色的”等等）。在 DNN 的末尾应该使用 sigmoid 还是 softmax 层？

5. 假设我们想要将单个动物的图像分类为猫或狗。我们应该使用 sigmoid 还是 softmax？

6. 如果我们减小批量大小，梯度下降训练过程中权重和偏差的更新会更多还是更少？