本章涵盖了其他 GPU 编程方法——使用 GPU 加速库和其他语言进行编程。使用 GPU 加速库编程使我们能够使用优化的内核开发应用。此外，我们可以使用其他编程语言开发 CUDA 软件，这些语言知道 CUDA 加速。这两种方式都提高了可编程性和生产率。此外，我们不必花时间优化已经优化的常见操作。
CUDA 工具包提供了许多线性代数、图像和信号处理以及随机处理方面的 GPU 加速库。它们是 cuBLAS(基本线性代数子程序)、cuFFT(快速傅里叶变换)、cuRAND(随机数生成)、NPP(图像和信号处理)、cuSPARSE(稀疏线性代数)、nvGRAPH(图形分析)、cussolver(GPU 中的 LAPACK)、推力(CUDA 中的 STL)等等。我们也可以用 OpenCV 库编写 GPU 加速程序。我们将在本章中介绍其中的一些库。
我们也可以使用图形处理器加速使用 R，MATLAB，Octave 和 Python。如今，Python 集成流行且功能强大，因为 GPU 可以加速许多机器学习和数据科学任务。我们还将这些语言作为入门级产品。
本章将涵盖以下主题:
- 用 cuBLAS 进行线性代数运算
- 使用 cuBLAS 的混合精度操作
- 并行随机数生成的 cuRAND
- 图形处理器中快速傅里叶变换的快速傅里叶变换
- 用图形处理器进行图像和信号处理的 NPP
- 在 OpenCV 中编写 GPU 加速代码
- 编写可与 CUDA 一起使用的 Python 代码
- 针对八度音程和 R 中零编码加速的 NVBLAS
- MATLAB 中的 CUDA 加速
cuBLAS 库是一个 GPU 优化的标准实现基本线性代数子程序 ( BLAS )。使用它的应用编程接口，程序员可以将 GPU 优化的计算密集型代码写入单个或多个 GPU。cuBLAS 有三个等级。第一级执行向量-向量运算，第二级执行矩阵-向量运算，第三级执行矩阵-矩阵运算。
涵盖每一个层次都超出了本书的范围。我们只是关注如何使用 cuBLAS APIs，并为多个 GPU 扩展其性能。具体来说，这张收据将包括一次单精度浮动矩阵乘法 ( SGEMM )操作——一次三级操作。
cuBLAS 库是 CUDA 工具包的一部分，因此您可以使用 cuBLAS 而无需额外安装。此外，您可以使用cc或cpp文件扩展名，而不是.cu，因为您不需要使用特定于 CUDA 的内置关键字，如__global__或threadIdx。以下代码片段显示了 cuBLAS 函数(cubalsSgemm)的基本应用:
cublasHandle_t handle;
cublasCreate(&handle);
.. { data operation } ..
cublasSgemm(...);
.. { data operation } ..
cublasDestroy(handle);如您所见，cuBLAS APIs 与cublasHandle_t类型句柄一起工作。
GEMM 操作可由以下等式表示:
其中α和β是标量， A 、 B 和 C 是列主格式的矩阵。这与以下框中的 cuBLAS 功能界面相匹配:
cublasStatus_t cublasSgemm(cublasHandle_t handle,
cublasOperation_t transa,
cublasOperation_t transb,
int m, int n, int k,
const float *alpha,
const float *A, int lda,
const float *B, int ldb,
const float *beta,
float *C, int ldc);在使用这个 GEMM 函数之前，让我们看一下参数的细节:
transa和transb:对 cuBLAS 功能的指示:矩阵 A 和 B 是否应进行转置操作。m、n和k:矩阵的尺寸大小。alpha和beta:决定如何从源配置输出值的参数。*A、*B、*C:矩阵数据的线性缓冲区。lda:矩阵 A 的前导列维度。cuBLAS 将矩阵元素与该值对齐。ldb:矩阵 B 的前导列尺寸。cuBLAS 将矩阵元素与该值对齐。
要在主机和设备之间传输数据，可以使用 cuBLAS 的cublasSetMatrix()和cublasGetMatrix()助手功能。它们是cudaMemcpy()的包装函数，但是有矩阵的维度信息，这样有助于增强代码的可读性；当然，你可以简单地使用cudaMemcpy()来代替。
让我们使用 cuBLAS SGEMM 函数实现一个具有 GEMM 操作的应用。我们将包括cublas_v2.h来使用更新的 cuBLAS API。为了方便起见，我们将使用getMatrix()函数从给定维度获取随机生成的矩阵，并使用printMatrix()函数打印矩阵元素。这些代码在给定的示例代码中实现。在主功能中，我们将从给定的M、N和K中初始化三个矩阵— A、B和C。那么我们将如下计算cublasSgemm():
cublasHandle_t handle;
// Prepare input matrices
float *A, *B, *C;
int M, N, K;
float alpha, beta;
M = 3; N = 4; K = 7;
alpha = 1.f; beta = 0.f;
// create cuBLAS handle
cublasCreate(&handle);
srand(2019);
A = getMatrix(K, M);
B = getMatrix(N, K);
C = getMatrix(M, N);
std::cout << "A:" << std::endl;
printMatrix(A, K, M);
std::cout << "B:" << std::endl;
printMatrix(B, N, K);
std::cout << "C:" << std::endl;
printMatrix(C, M, N);
// Gemm
cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_T,
M, N, K, &alpha, A, K, B, N, &beta, C, M);
cudaDeviceSynchronize();
std::cout << "C out:" << std::endl;
printMatrix(C, M, N);
cublasDestroy(handle);
cudaFree(A); cudaFree(B); cudaFree(C);
return 0;通过链接 cuBLAS 库用nvcc编译代码:
$ nvcc -run -gencode arch=compute_70,code=sm_70 -lcublas -o cublasSgemm ./cublasSgemm.cpp下面的代码片段显示了执行的输出:
A:
0.0492 0.4790 0.0226
0.7750 0.2794 0.8169
0.3732 0.6248 0.2636
0.9241 0.5841 0.8532
0.7188 0.5052 0.5398
0.9869 0.6572 0.0520
0.6815 0.7814 0.5988
B:
0.8957 0.0481 0.7958 0.7825 0.3264 0.5189 0.5018
0.4559 0.6342 0.0759 0.5051 0.1132 0.0985 0.2802
0.3926 0.9153 0.6534 0.0174 0.1790 0.5775 0.6015
0.0322 0.2963 0.1068 0.5720 0.2832 0.7640 0.6240
C:
0.9647 0.5454 0.2229 0.8604
0.5935 0.0186 0.6430 0.9198
0.5375 0.1448 0.3757 0.1718
C out:
1.1785 2.5682 2.4854 0.6066
0.5817 0.8091 1.1724 2.0773
2.0882 1.4503 2.1331 1.8450如cublasSgemm()函数调用中所述，矩阵 A 和 B 是转置矩阵。我们将原来的引导列大小传递给cublasSgemm()函数，分别为lda、ldb和ldc，可以看到操作按预期进行。
cuBLAS 库的 cuBLAS-XT API 提供了 cuBLAS 在多个 GPU 上工作时的三级操作。有了这个 API，你的应用可以使用多 GPU 计算操作。这个片段展示了使用 cuBLAS-XT 的基本操作:
cublasXtHandle_t handle;
cublasXtCreate(&handle);
cudaGetDeviceCount(&num_of_total_devices);
devices = (int *)calloc(num_of_devices, sizeof(int));
for (int i = 0; i < num_of_devices; i++)
devices[i] = i;
cublasXtDeviceSelect(handle, num_of_devices, devices);
cublasXtSgemm( ... );
cublasXtDestroy(handle);cublasXtSgemm()界面与cublasSgemm()功能相同，可以放心使用多个 GPU 的计算性能。例如，我们可以使用带有两个图形处理器的存储库中的示例代码获得以下结果。根据您的 GPU 和系统配置，此性能可能会有所不同:
cuBLAS 库提供了很多通用的线性代数运算。所以你应该检查你的必要功能是如何在库中提供的。此外，您还需要一个如何使用该函数的示例。以下项目是文档和示例的链接。因此，当您需要实现基于 cuBLAS 的应用时，建议您经常检查这两个文档:
- 英伟达 cuBLAS 编程指南——参考指南:https://docs.nvidia.com/cuda/cublas/index.html
GPU 上的矩阵计算
: CUBLAS、CUSOLVER 和 MAGMA 举例安杰伊·chrzȩszczyk 和雅各布·安德斯:https://developer . NVIDIA . com/sites/default/files/akamai/cuda/files/misc/mygpu . pdf
cuBLAS 库支持混合精度计算。该计算意味着以不同精度操作的操作，例如，用单精度和半精度变量或用单精度和字符(INT8)进行计算。当我们需要使用较低的精度来获得更高的性能，同时在结果中获得更高的精度时，这种技术非常有用。
cuBLAS 库提供cublasGemmEx()和cublas{S/C}gemmEx()支持混合精度操作的 GEMM 操作。它们是cublas<t>gemm()的扩展，接受每个 A 、 B 和 C 矩阵的指定数据类型。下表显示了 cuBLAS 库中cublasGemmEx()和其他可替换 API 的精度支持矩阵:
| 计算类型 | a 型/ B 型 | c 型 | 可替换的 API |
| CUDA_R_16F | CUDA_R_16F | CUDA_R_16F | cublasHgemm() |
| CUDA_R_32I | CUDA_R_8I | CUDA_R_32I | 不适用的 |
| CUDA_R_32F | CUDA_R_16F | CUDA_R_16F | cublasSgemmEx() |
| CUDA_R_8I | CUDA_R_32F |
| CUDA_R_16F | CUDA_R_32F |
| CUDA_R_32F | CUDA_R_32F |
| CUDA_R_64F | CUDA_R_64F | CUDA_R_64F | cublasDgemm() |
| CUDA_C_32F | CUDA_C_8I | CUDA_C_32F | cublasCgemmEx() |
| CUDA_C_32F | CUDA_C_32F |
| CUDA_C_64F | CUDA_C_64F | CUDA_C_64F | cublasZgemm() |
可以看到cublasGemmEx()可以覆盖cublas{S/C}gemmEx()功能的操作。因此，我们将在本节中介绍cublasGemmEx()。
cublasGemmEx()函数的最后一个参数cublasGemmAlgo_t指定了矩阵-矩阵乘法的算法。有了这个参数，我们可以选择是否使用 TensorCore。CUBLAS_GEMM_DEFAULT选择 GEMM 算法，在 CUDA 核上运行。另一方面，CUBLAS_GEMM_DEFAULT_TENSOR_OP选择使用张量核的算法。如果 TensorCore 在给定条件下不可用，cuBLAS 将选择使用 CUDA 核心的算法。这种情况可能发生在没有张量核心或矩阵大小的图形处理器上，这与张量核心的工作方式不相符——四的倍数( ** 4* )。
现在，让我们使用 cuBLAS GEMM 操作来尝试混合精度。实现之后，我们将介绍矩阵大小如何影响操作。完全实现的版本在02_sgemm_mixed_precision/cublasGemmEx.cu:
- 这段代码使用了一个定制的内存管理类
CBuffer，以简化混合精度和复制的处理，但是它可以使用统一内存来代替。对于 cuBLAS 操作，我们应该在代码中包含cublas_v2.h:
#include <cublas_v2.h>
#include "helper.cuh" // for CBuffer and printMatrix()- 现在，让我们实现
main()功能。首先，我们将创建并初始化A、B和C矩阵。以下代码片段显示了如何使用CBuffer类并初始化矩阵:
int M = 4, N = 5, K = 6;
CBuffer<half> A, B;
CBuffer<float> C;
A.init(K * M, true);
B.init(N * K, true);
C.init(N * M, true);- 要指定
A、B、C的精度类型，并一起测试各种精度，需要指定一些 CUDA 数据类型参数:
cudaDataType TYPE_A, TYPE_B, TYPE_C;
if (typeid(*A.h_ptr_) == typeid(float)) {
TYPE_A = TYPE_B = CUDA_R_32F;
}
else if (typeid(*A.h_ptr_) == typeid(half)) {
TYPE_A = TYPE_B = CUDA_R_16F;
}
else if (typeid(*A.h_ptr_) == typeid(int8_t)) {
TYPE_A = TYPE_B = CUDA_R_8I;
}
else {
printf("Not supported precision\n");
return -1;
}
if (typeid(*C.h_ptr_) == typeid(float)) {
TYPE_C = CUDA_R_32F;
}
else if (typeid(*C.h_ptr_) == typeid(int)) {
TYPE_C = CUDA_R_32I;
}
else {
printf("Not supported precision\n");
return -1;
}- 对于 cuBLAS 操作，我们应该初始化
cublas_handle、alpha和beta:
float alpha = 1.f, beta = 0.f;
cublasHandle_t cublas_handle;
cublasCreate(&cublas_handle);- 然后，我们将数据复制到 GPU:
A.cuda(true);
B.cuda(true);
C.cuda(true);- 然后我们调用
cublasGemmEx()函数如下:
cublasGemmEx(cublas_handle,
CUBLAS_OP_N, CUBLAS_OP_N,
M, N, K,
&alpha, A.d_ptr_, TYPE_A, M, B.d_ptr_, TYPE_B, K,
&beta, C.d_ptr_, TYPE_C, M, TYPE_C,
CUBLAS_GEMM_DEFAULT);- 要查看矩阵值，我们可以使用
printMatrix()，它在helper.h中定义:
std::cout << "A:" << std::endl;
printMatrix(A.h_ptr_, K, M);
std::cout << "B:" << std::endl;
printMatrix(B.h_ptr_, N, K);- 使用函数覆盖方法定义
printMatrix()，以允许在其他数据类型中打印具有相同格式的半精度值。部分定义如下:
template <typename T>
void printMatrix(const T *matrix, const int ldm, const int n) {
std::cout << "[" << __FUNCTION__ << "]::
Not supported type request" << std::endl;
}
void printMatrix(const float *matrix, const int ldm, const int n) {
for (int j = 0; j < n; j++) {
for (int i = 0; i < ldm; i++)
std::cout << std::fixed << std::setw(8) <<
std::setprecision(4) <<
matrix[IDX2C(i, j, ldm)];
std::cout << std::endl;
}
}
void printMatrix(const half *matrix, const int ldm, const int n) {
for (int j = 0; j < n; j++) {
for (int i = 0; i < ldm; i++)
std::cout << std::fixed << std::setw(8) <<
std::setprecision(4) << __half2float(matrix[IDX2C(i, j, ldm)]);
std::cout << std::endl;
}
}
... ( functions for other data types ) ...- 然后，代码将对给定的
A、B和C矩阵进行 GEMM 运算。以下是当M为4、N为5、M为6时的输出示例:
$ nvcc -run -m64 -std=c++ 11 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -lcublas -o cublasGemmEx ./cublasGemmEx.cu
A:
0.0049 0.0479 0.0023 0.0775 0.0279 0.0817
0.0373 0.0625 0.0264 0.0924 0.0584 0.0853
0.0719 0.0505 0.0540 0.0987 0.0657 0.0052
0.0682 0.0781 0.0599 0.0896 0.0048 0.0796
B:
0.0624 0.0965 0.0545 0.0223 0.0861
0.0594 0.0019 0.0643 0.0920 0.0537
0.0145 0.0376 0.0172 0.0221 0.0881
0.0285 0.0319 0.0161 0.0677 0.0235
0.0814 0.0695 0.0414 0.0392 0.0296
0.0446 0.0688 0.0403 0.0018 0.0971
C:
0.0509 0.0117 0.0877 0.0445 0.0830
0.0742 0.0242 0.0136 0.0625 0.0681
0.0362 0.0046 0.0265 0.0963 0.0638
0.0070 0.0446 0.0516 0.0194 0.0089
C out:
0.0153 0.0228 0.0143 0.0292 0.0113
0.0200 0.0118 0.0214 0.0081 0.0138
0.0098 0.0168 0.0132 0.0199 0.0125
0.0269 0.0120 0.0222 0.0085 0.0228现在，让我们尝试其他数据类型，看看cublasGemmEx()如何对给定的矩阵进行操作。提供的示例还输出操作的执行时间来衡量性能:
- 如果矩阵 A 或者矩阵 B 是转置矩阵，应该修改什么？
- 有没有比手术更好的矩阵尺寸？通过更改大小来比较执行时间。
- 每种数据类型有没有更好的矩阵大小？如果你尝试
INT8精度，你会看到错误。这怎么解决？改变尺寸，看看如何在cublasGemmEx()支持INT8操作。
张量核心提供张量点运算的加速性能。支持 Volta 架构的 FP16，图灵架构的INT8和INT4。因此，我们应该使用缩减精度或混合精度来使用 TensorCore。
之前我们使用CUBLAS_GEMM_DEFAULT作为 cuBLAS GEMM 算法，在运算中使用 CUDA 核。要用 TensorCore，就要用CUBLAS_GEMM_DEFAULT_TENSOR_OP。为了利用张量核心，操作数矩阵的每个维度都应该是 4 的倍数。这是天梭内部 WMMA (简称，翘曲矩阵乘积)运营优化的单位规模。例如，用 A (8，192 × 8，192)和 B (8，192 × 8，192)进行的矩阵-矩阵乘法，与用 A (8，192 × 8，192)和 B (8，192 × 8，190)进行的运算相比，表现出高得多的性能。您也可以通过配置文件确认此操作。
以下时间线是使用矩阵 A (8，192 × 8，192)和矩阵 B (8，192 × 8，190)进行矩阵乘法的结果:
此外，该时间线图像是矩阵 A (8，192 × 8，192)和矩阵 B (8，192 × 8，192)的矩阵相乘的结果:
两个测试都在 CUDA C/C++ 中使用CUBLAS_GEMM_DEFAULT_TENSOR_OP，但是使用 TensorCore 的 GEMM 运算速度比使用 CUDA 内核快 6.7 倍。由于 TensorCore 基于矩阵大小可用，nvcc用特殊的内核函数编译代码，从volta_s884g开始。总之，如果你想获得张量核心的好处，把你的矩阵填充到 4。这可能是一个开销，但是来自 TensorCore 的性能提升可能会超过这个开销。
英伟达在其开发博客网站(https://devblogs.nvidia.com/programming-tensor-cores-cuda-9)中提供了如何使用 cuBLAS 库对 TensorCores 进行编程。本文档还介绍了其他可用的方法。但是，使用 cuBLAS 库可以为您提供最快的性能，正如橡树岭国家实验室的一篇论文所证明的那样——*NVIDIA Tensor Core professional、*Performance and Precision(https://arxiv.org/pdf/1803.04014.pdf)。
许多应用使用伪随机数进行模拟或概率分析。尽管有传统的用法，大量的随机数生成过程花费了大量的时间。一种解决方案是并行生成随机数，但是每个多线程应该有不同的随机种子，以便独立生成随机数。
库兰德库使图形处理器能够从图形处理器生成许多随机数。该库可从主机或设备代码获得。主机应用编程接口只允许使用主机代码生成随机数。因此，您可以将生成的数据直接用于其他内核函数。设备应用编程接口支持在内核代码中生成随机数，因此您可以让 CUDA 线程在执行过程中拥有自己的随机生成的数字。
首先，您需要使用curandGenerator()创建所需类型的新生成器。然后，设置所需种子和顺序的生成器选项。例如，您可以使用curandSetPseudoRandomGeneratorSeed()生成伪随机生成器:
curandGenerator_t curand_gen;
curandCreateGenerator(&curand_gen, CURAND_RNG_PSEUDO_DEFAULT);
curandSetPseudoRandomGeneratorSeed(curand_gen, 2019UL);然后，可以使用curandGenerate()生成随机数。有九种不同的生成函数。例如，您可以使用curandGenerateUnifrom()生成均匀分布的浮点值:
curandGenerateUniform(curand_gen, random_numbers_on_device_memory, length);The cuRAND programming guide provides descriptions of various kinds of generation functions: https://docs.nvidia.com/cuda/curand/host-api-overview.html#generation-functions.
使用随机数生成后，您可以使用毁灭者功能终止 cuRAND 生成器:
curandDestroyGenerator(curand_gen);现在，让我们实现一个使用几个 cuRAND APIs 生成随机数的应用。完全实现的版本是03_curand/curand_host.cpp。因此，您可以根据需要修改代码并测试其他功能。
首先，我们应该为 cuRAND 主机 API 和其他与 CPP 相关的头文件包含curand.h，如下所示:
#include <iostream>
#include <iomanip>
#include <curand.h>假设我们将创建一个用随机数初始化的矩阵。我们需要实现printMatrix()功能，以便对生成的随机数进行如下检查:
#define IDX2C(i, j, ld) (((j) * (ld)) + (i))
template <typename T>
void printMatrix(const T *matrix, const int ldm, const int n)
{
for (int j = 0; j < ldm; j++) {
for (int i = 0; i < n; i++)
std::cout << std::fixed << std::setw(12)
<< std::setprecision(4) << matrix[IDX2C(i, j, ldm)];
std::cout << std::endl;
}
}然后，我们将按如下方式分配所需的内存空间。现在，我们将实现main()函数，该函数使用printMatrix()初始化随机数并打印结果。首先，我们将如下初始化操作的 cuRAND 句柄:
// create curand generator & set random seed
curandGenerator_t curand_gen;
curandCreateGenerator(&curand_gen, CURAND_RNG_PSEUDO_DEFAULT);
curandSetPseudoRandomGeneratorSeed(curand_gen, 2019UL);您可以根据需要更改随机种子。接下来就是分配内存空间。为了便于评估操作，我们将使用统一的内存，因为 cuRAND 函数将在 GPU 上生成随机数:
size_t size = M * N;
unsigned int *np_random;
float *fp_random;
cudaMallocManaged((void**)&np_random, sizeof(*np_random) * size);
cudaMallocManaged((void**)&fp_random, sizeof(*fp_random) * size);接下来，我们将生成给定内存空间的随机数。我们将使用整数存储空间(np_random)来生成随机数，使用浮动存储空间(fp_random)来生成均匀分布的随机数，如下所示:
// random number generation
std::cout << "Generated random numbers" << std::endl;
curandGenerate(curand_gen, np_random, size);
cudaDeviceSynchronize();
printMatrix(np_random, M, N);
// uniform distributed random number generation
std::cout << "Generated uniform random numbers" << std::endl;
curandGenerateUniform(curand_gen, fp_random, size);
cudaDeviceSynchronize();
printMatrix(fp_random, M, N);因为我们使用的是统一的内存，所以我们可以允许 GPU 和主机共享同一个内存地址，我们可以通过同步它们来查看输出值。最后，我们可以终止 cuRAND 句柄和内存，如下所示:
// terminates used resources
curandDestroyGenerator(curand_gen);
cudaFree(np_random);
cudaFree(fp_random);现在，是时候编译和运行代码了。使用 cuRAND APIs 编译代码应该为nvcc编译器提供-lcurand。当M = 3和N = 5时，输出如下:
$ nvcc -run -gencode arch=compute_70,code=sm_70 -lcurand -o curand_host curand_host.cpp
Generated random numbers
3395652512 793372546 2133571103 595847267 2461872808
595847267 2461872808 500895635 498154070 2385617847
498154070 2385617847 196336856 388563169 745758309
Generated uniform random numbers
0.7134 0.0830 0.1458 0.2099 0.6066
0.2099 0.6066 0.3078 0.5122 0.8856
0.5122 0.8856 0.3530 0.8477 0.8370我们已经介绍了 CUDA 如何使用主机 API 生成随机数，但是，在某些情况下，最好设计 CUDA 内核来生成随机数。我们称之为设备 API，我们可以从每个 CUDA 线程中获取随机数。
使用设备应用编程接口，我们可以设置生成器种子，并在您的 CUDA 设备上生成随机数。
首先，我们需要准备一个curandState_t的设备存储空间，以便存储生成器种子，将随机种子并行提供给 CUDA 线程。这可以像正常的设备内存分配代码一样完成，如下所示:
cudaMalloc((void **)&devStates, length * sizeof(curandState_t));在你的内核代码中，我们需要使用curand_init()初始化随机种子。该函数需要种子、序列号和偏移量。然后，该函数设置状态。对于相同的种子，cuFFT 总是生成相同的状态。要生成随机值，使用curand()功能。和主机的生成函数一样，设备 API 也有各种生成函数。例如，均匀分布的随机数生成可以这样完成:
int idx = blockIdx.x * blockDim.x + threadIdx.x;
curand_init(2019UL, idx, 0, &state[idx]);
generated_out[idx] = curand_uniform(&state[idx]);cuRAND 库为各种数据类型和随机分布提供了各种生成函数。要找到您想要的生成函数，请查看《cuRAND 开发者指南》的设备 API 概述。随机数生成后，设备状态缓冲区应像正常内存一样终止，如下所示:
cudaFree(devStates);现在，我们将创建一个使用 cuRAND 设备 API 的应用。完全实现的代码是curand_device.cu，所以你也可以修改和测试代码。首先，我们应该将curand_kernel.h文件与其他 C++ 必需的头文件包括在内，如下所示:
#include <iostream>
#include <iomanip>
#include <curand_kernel.h>我们将编写setup_kernel()，为每个 CUDA 线程初始化一个随机种子，如下所示:
__global__ void setup_kernel(curandState_t *state)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
// Each thread gets same seed,
// a different sequence number, no offset */
curand_init(2019UL, idx, 0, &state[idx]);
}编写两个随机数生成函数:generate_kernel()和generate_uniform_kernel()。我们将使用均匀分布的随机数生成一个 32 位整数和一个浮点:
__global__ void generate_kernel(unsigned int *generated_out,
curandState_t *state)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
generated_out[idx] = curand(&state[idx]) & 0xFF;
}
__global__ void generate_uniform_kernel(float *generated_out,
curandState_t *state)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
generated_out[idx] = curand_uniform(&state[idx]);
}现在，我们将实现main()功能并初始化设备状态缓冲区:
cudaMalloc((void **)&devStates, sizeof(curandState) * size);
setup_kernel<<<(size + BLOCK_DIM - 1) / BLOCK_DIM, BLOCK_DIM>>>(devStates);然后，使用generate_kernel()生成随机数。为了方便起见，我们将使用统一内存作为空间，并验证主机的输出。之后，我们将打印出如下结果:
// random number generation
std::cout << "Generated random numbers" << std::endl;
cudaMallocManaged((void**)&np_random, sizeof(*np_random) * size);
generate_kernel<<<(size + BLOCK_DIM - 1) / BLOCK_DIM, BLOCK_DIM>>>
(np_random, const_cast<curandState_t *>(devStates));
cudaDeviceSynchronize();
printMatrix(np_random, M, N);同样，我们将使用generate_uniform_kernel()创建均匀分布的随机数，如下所示:
// uniform distributed random number generation
std::cout << "Generated uniform random numbers" << std::endl;
cudaMallocManaged((void**)&fp_random, sizeof(*fp_random) * size);
generate_uniform_kernel<<<(size + BLOCK_DIM - 1) / BLOCK_DIM, BLOCK_DIM>>>
(fp_random, const_cast<curandState_t *>(devStates));
cudaDeviceSynchronize();
printMatrix(fp_random, M, N);因为我们使用的是统一内存，所以我们可以允许 GPU 和主机共享同一个内存地址，我们可以通过同步它们来查看输出值。最后，我们可以终止 cuRAND 句柄和内存，如下所示:
// terminates used resources
curandDestroyGenerator(curand_gen);
cudaFree(np_random);
cudaFree(fp_random);现在，是时候编译和运行代码了。为了使用 cuRAND 编译代码，API 应该为nvcc编译器提供-lcurand。当M等于3、N等于5时，输出如下:
$ nvcc -run -gencode arch=compute_70,code=sm_70 -lcurand -o curand_device curand_device.cpp
Generated random numbers
3395652512 793372546 2133571103 595847267 2461872808
595847267 2461872808 500895635 498154070 2385617847
498154070 2385617847 196336856 388563169 745758309
Generated uniform random numbers
0.8064 0.2783 0.2971 0.2386 0.7491
0.2386 0.7491 0.4782 0.1060 0.2922
0.1060 0.2922 0.1823 0.6199 0.9137当您比较主机应用编程接口和设备应用编程接口的输出数字时，生成的随机数是相同的，而统一随机数不是。如果您在第二次生成随机数之前重置随机种子，就可以解决这个问题。
之前，我们使用 C++ 随机数生成器来初始化 GEMM 运算的矩阵。一般来说，当我们想要生成随机数时，这个函数很方便。但是，您可能会发现，在最后一节中，这个函数花费了很长时间来生成大随机数。在本节中，我们将介绍 cuRAND API 如何与 cuBLAS GEMM 操作协同工作。完全实现的版本是gemm_with_curand_host.cpp文件。让我们看看这是如何实现的:
- 目前，我们在 cuRAND 库中没有低精度随机数生成器。此外，我们需要将半精度数字转换为浮点数，以便评估输出。由于这些原因，我们需要在 GPU 上创建如下类型转换函数:
namespace fp16{
__global__ void float2half_kernel(half *out, float *in)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
out[idx] = __float2half(in[idx]);
}
void float2half(half *out, float *in, size_t length)
{
float2half_kernel<<< (length + BLOCK_DIM - 1) / BLOCK_DIM, BLOCK_DIM >>>(out, in);
}- 现在，我们将编写一个使用 cuRAND 主机 API 的随机数生成函数。正如我们之前所讨论的，当我们需要使用半精度数据时，我们应该将生成的随机数从浮点数转换为二分之一。该功能可以如下实现:
template <typename T>
typename std::enable_if<(std::is_same<T, float>::value), float>::type
*curand(curandGenerator_t generator, size_t length)
{
T *buffer = nullptr;
cudaMalloc((void **)&buffer, length * sizeof(float));
curandGenerateUniform(generator, buffer, length);
return buffer;
}
template <typename T>
typename std::enable_if<std::is_same<T, half>::value, half>::type
*curand(curandGenerator_t generator, size_t length)
{
T *buffer = nullptr;
float *buffer_fp32;
cudaMalloc((void **)&buffer_fp32, length * sizeof(float));
curandGenerateUniform(generator, buffer_fp32, length);
cudaMalloc((void **)&buffer, length * sizeof(T));
fp16::float2half(buffer, buffer_fp32, length);
cudaFree(buffer_fp32);
return buffer;
}- 在
main()函数中定义一些控制 GEMM 运算的局部变量:
void *d_A, *d_B, *d_C;
cudaDataType AType, BType, CType, computeType;
int M = 8192, N = 8192, K = 8192;
float alpha = 1.f, beta = 1.f;
std::string precision = "fp32";
bool tensor_core = true;在这段代码中，我们确定了 GEMM 操作的大小、数据类型和操作类型。
- 现在，让我们创建输入缓冲区数组，并设置参数以及操作精度:
if (precision == "fp32") {
auto *a = curand<float>(curand_gen, M * K);
auto *b = curand<float>(curand_gen, K * N);
auto *c = curand<float>(curand_gen, M * N);
AType = BType = CType = CUDA_R_32F;
computeType = CUDA_R_32F;
d_A = a, d_B = b, d_C = c;
}
else if (precision == "fp16") {
auto *a = curand<half>(curand_gen, M * K);
auto *b = curand<half>(curand_gen, K * N);
auto *c = curand<float>(curand_gen, M * N);
AType = BType = CUDA_R_16F, CType = CUDA_R_32F;
computeType = CUDA_R_32F;
d_A = a, d_B = b, d_C = c;
}
else {
exit(EXIT_FAILURE);
}- 按如下方式创建 cuRAND 和 cuBLAS 句柄:
cublasCreate(&cublas_handle);
curandCreateGenerator(&curand_gen, CURAND_RNG_PSEUDO_DEFAULT);
curandSetPseudoRandomGeneratorSeed(curand_gen, 2019UL);- 然后，我们应该确定使用张量核心的操作类型:
cublasGemmAlgo_t gemm_algo = (tensor_core) ?
CUBLAS_GEMM_DEFAULT_TENSOR_OP : CUBLAS_GEMM_DEFAULT;- 然后，我们可以调用
cublasGemmEx()函数，该函数提供 FP32 和 FP16 操作，如下所示:
cublasGemmEx(cublas_handle, CUBLAS_OP_N, CUBLAS_OP_N,
M, N, K,
&alpha, d_A, AType, M, d_B, BType, K,
&beta, d_C, CType, M,
computeType, gemm_algo);与之前的版本相比，GEMM 操作应该显示出类似的性能。但是，您可能会发现整个应用的速度得到了提高，因为 GPU 上的并行随机数生成比主机上的生成快得多。
《cuRAND 开发者指南》将帮助您找到其他随机数生成器、选项和分布。本文件位于https://docs.nvidia.com/pdf/CURAND_Library.pdf。
cuFFT 库为 FFT (简称，快速傅立叶变换)算法提供 GPU 加速操作。程序员可以利用图形处理器的计算能力转换真实或复杂的数据，并对转换后的信号应用图形处理器内核操作。此外，支持的功能与 FFTW 库相匹配，因此我们可以将主机项目迁移到 GPU。
为了处理快速傅立叶变换样本的尺寸信息，要求快速傅立叶变换相应地使用cufftPlan1D()、cufftPlan2D()或cufftPlan3D()创建一个平面处理。如果样本数据具有批处理和跨步布局，我们应该使用cufftPlanMany()。如果样本大小大于 4 GB，我们应该使用64作为计划函数的后缀来支持该大小。例如，cufftPlanMany64()在cufftPlanMany()功能上支持更大的样本。
cuFFT 库支持多 GPU 操作。首先，您需要使用cufftCreate()创建一个空计划。然后，我们可以使用cufftXtSetGPUs()指定将要执行操作的 GPU 列表。之后，我们可以使用正常的计划生成功能来生成计划，我们之前已经介绍过了。下表显示了计划生成功能类别:
| | 基本计划 | 多图形处理器计划 |
| 简单的计划 | cufftPlan{1d,2d,3d}() | cufftMakePlan{1d,2d,3d}() |
| 高级数据布局 | cufftPlanMany() | cufftMakePlanMany() |
| FP16 操作 | cufftXtMakePlanMany() |
然后，您可以使用cufftExec()功能对样本数据进行正向(快速傅立叶变换)和反向(IFFT)。cuFFT 库提供了三种数据转换:复杂到复杂、真实到复杂和复杂到真实。其操作数据类型可以是浮点型或双精度型:
| 变换方向 | 浮动 | 两倍 |
| 复杂到复杂 | cufftExecC2C()
cufftXtExecDescriptorC2C() | cufftExecZ2Z()
cufftXtExecDescriptorZ2Z() |
| 真实到复杂 | cufftDExecR2C()
cufftXtExecDescriptorR2C() | cufftExecD2Z()
cufftXtExecDescriptorD2Z() |
| 复杂到真实 | cufftExecC2R()
cufftXtExecDescriptorC2R() | cufftExecZ2D()
cufftXtExecDesciptorZ2D() |
| 全部 | cufftXtExec() / cufftXtExecDesciptor() |
cuFFT 操作为正向或反向，该操作应与另一方向配对。
在真实数据和复杂数据之间转换的函数，如R2C和C2R，在它们的函数名中有隐含的方向信息。此功能有助于您避免为了将真实域中的数据转换为复杂的数据类型而进行额外的操作。同时，您必须创建一个额外的计划，因为每个计划都有转型方向信息。
另一方面，要提供复到复变换的变换方向信息，如C2C、Z2Z。对于反转操作，您不必创建另一个 cuFFT 句柄，因为计划应该是相同的数据类型操作。
cufftXtExec()和cufftXtExecDescriptor()函数可以对任何给定的数据类型执行转换，因为当您创建 cuFFT 计划时，每个输入数据都应该提供它们的数据类型信息。
现在让我们尝试使用 cuFFT。完全实现的版本是04_cufft/cufft.1d.cpp文件。让我们讨论一下这是如何实现的:
- 首先从一些头文件开始:C++、CUDA、cuRAND、cuFFT:
#include <iostream>
#include <iomanip>
#include <cuda_runtime.h>
#include <cufft.h>
#include <curand.h>- 在这个快速傅立叶变换操作中，我们将有实数到复数和复数到实数的转换。因此，让我们声明一些自定义数据类型，
Real和Complex，以简化代码。这可以通过以下方式实现:
typedef cufftReal Real;
typedef cufftComplex Complex;- 现在，让我们从
main()功能开始。对于输入的样本数据，我们将使用统一内存，以便于主机和 GPU 之间的数据传输。转换后的数据只能在图形处理器上使用。因此，可以按如下方式分配内存空间:
cudaMallocManaged((void**)&p_sample, sizeof(Real) * sample_size * batch_size);
cudaMalloc((void**)&d_freq, sizeof(Complex) * sample_size * batch_size);- 然后，我们将使用 cuRAND 主机 API 初始化输入数据，如下所示:
curandGenerator_t curand_gen;
curandCreateGenerator(&curand_gen, CURAND_RNG_PSEUDO_DEFAULT);
curandSetPseudoRandomGeneratorSeed(curand_gen, 2019UL);
curandGenerateUniform(curand_gen, p_sample, sample_size * batch_size);- 并且，我们应该为正向和反向变换初始化 cuFFT 计划。由于它们具有不同的数据类型，我们应该分别为实数到复数和复数到实数的转换创建两个计划:
int rank = 1;
int stride_sample = 1, stride_freq = 1;
int dist_sample = sample_size, dist_freq = sample_size / 2 + 1;
int embed_sample[] = {0};
int embed_freq[] = {0};
cufftPlanMany(&plan_forward, rank, &sample_size,
embed_sample, stride_sample,
dist_sample,
embed_freq, stride_freq, dist_freq,
CUFFT_R2C, batch_size);
cufftPlanMany(&plan_inverse, rank, &sample_size,
embed_freq, stride_freq, dist_freq,
embed_sample, stride_sample,
dist_sample,
CUFFT_C2R, batch_size);- 现在，我们可以使用给定的 cuFFT 计划进行正向或反向变换。为了度量执行时间，我们可以使用 CUDA 事件来包含这些操作:
cufftExecR2C(plan_forward, p_sample, d_freq);
cufftExecC2R(plan_inverse, d_freq, p_sample);
- 然后，我们可以用以下命令编译代码:
$ nvcc -run -gencode arch=compute_70,code=sm_70 -lcufft -lcurand -o cufft.1d cufft.1d.cppcufft.1d命令将报告每一步的变换时间，如下所示:
FFT operation time for 1048576 elements with 512 batch..
Forward (ms): 21.5322
Inverse (ms): 21.4cuFFT 库提供扩展的 CUDA 计算功能，例如 FP16 FFT 操作。完整版本是cufft.half.cpp文件。让我们讨论一下它的实现。
在这段代码中，我们应该使用cufftXtMakePlanMany()进行计划创建，使用cufftXtExec()函数进行转换。cufftXtMakePlanMany()如果输入和输出数据类型是 FP16 或 FP32，则允许传递。此外，我们应该为正向和反向转换创建两个计划，以便涵盖实数到复数和复数到实数的转换。对于空 cuFFT 计划，cufftXtMakePlanMany()可以指定样本大小、输入数据格式和类型、批次大小等。例如，计划创建可以如下实现:
int rank = 1;
int stride_sample = 1, stride_freq = 1;
long long int dist_sample = sample_size, dist_freq = sample_size / 2 + 1;
long long embed_sample[] = {0};
long long embed_freq[] = {0};
size_t workSize = 0;
cufftCreate(&plan_forward);
cufftXtMakePlanMany(plan_forward,
rank, &sample_size,
embed_sample, stride_sample, dist_sample, CUDA_R_16F,
embed_freq, stride_freq, dist_freq, CUDA_C_16F,
batch_size, &workSize, CUDA_C_16F);
cufftCreate(&plan_inverse);
cufftXtMakePlanMany(plan_inverse,
rank, &sample_size,
embed_freq, stride_freq, dist_freq, CUDA_C_16F,
embed_sample, stride_sample, dist_sample, CUDA_R_16F,
batch_size, &workSize, CUDA_R_16F);在这个实现中，我们还必须考虑是否以半精度提供输入数据。您可以使用主机随机函数，并将它们转换为半精度数据，但是，这段代码向您展示了如何将 cuRAND 主机 API 用于此目的，如下所示:
template <typename T>
typename std::enable_if<std::is_same<T, half>::value>::type
curand(curandGenerator_t generator, T *buffer, size_t length) {
float *buffer_fp32;
cudaMalloc((void **)&buffer_fp32, length * sizeof(float));
curandGenerateUniform(generator, buffer_fp32, length);
// convert generated single floating to half floating
fp16::float2half(buffer, buffer_fp32, length);
cudaFree(buffer_fp32);
}因此，我们可以为快速傅立叶变换提供均匀分布随机数的半精度，并且我们可以使用cufftXtExec()进行正向和反向变换。转化表现如下:
FFT operation time for 1048576 elements with 512 batch..
Forward (ms): 15.3236
Inverse (ms): 15.4881快速傅立叶变换的另一种用法是使用多个图形处理器进行大型快速傅立叶变换操作。为此，我们必须使用cufftCreate()创建一个空的 cuFFT 计划，并使用cufftXtSetGPUs()提供 GPU 编号。例如，这可以按如下方式完成:
cufftHandle cufft_plan;
int n_gpu = 2, devices[2] = {0,1};
cufftCreaet(&cufft_plan); // create an empty plan
cufftXtSetGPUs(cufft_plan, n_gpu, devices); // set multi-gpu information图形处理器的总数可能因系统而异。现在，我们可以使用cufftXtMakePlanMany()指定样本信息来生成 cuFFT 计划。例如，cufftXtMakePlanMany()可以这样称呼:
size_t *work_size = (size_t*) new size_t[num_gpus];
cufftXtMakePlanMany(cufft_plan, 1 &sample_size,
nullptr, 1, 1, CUDA_C_32F,
nullptr, 1, 1, CUDA_C_32F,
batch_size, work_size, CUDA_C_32F);cuFFT 库提供cufftXtMalloc()，为目标 GPU 准备 GPU 内存空间。然后，我们可以使用cufftXtMemcpy()功能将数据复制到分配的内存中。例如，这可以如下实现:
cudaLibXtDesc *d_sample;
cufftXtMalloc(cufft_plan, &d_sample, CUFFT_XT_FORMAT_INPLACE);
cufftXtMemcpy(cufft_plan, d_sample, h_sample, CUFFT_COPY_HOST_TO_DEVICE);然后，我们可以用cufftXtExecDesciptor()功能在多 GPU 上执行 FFT:
cufftXtExecDesciptor(cufft_plan, d_sample, d_sample, CUFFT_FORWARD);使用nvidia-smi，我们可以监控跨 GPU 的分布式内存分配和执行。根据您的图形处理器和系统配置，经过的时间可能不同。
如果你想了解更多关于 cuFFT 库及其功能的知识，cuFFT 库用户指南(https://docs.nvidia.com/cuda/cufft/index.html)是你很好的参考。
CUDA 示例代码是学习如何使用 cuFFT 函数的另一个很好的参考。样本代码放在NVIDIA_CUDA-10.x_Samples/7_CUDALibraries/CUFFT*目录中。您可以学习如何使用 CUDA 内核代码，并通过使用 cuFFT 的前向/后向转换来应用过滤器操作。
NPP (简称，NVIDIA Performance Primitive)库是一个默认的 CUDA 库，拥有一组专注于成像和视频处理的 GPU 加速处理功能。虽然它支持这些领域的灵活开发，但开发人员可以节省他们的应用开发时间。
NPP 库有两个功能部分:成像处理 API 和信号处理 API。图像处理应用编程接口包括与图像过滤、压缩/解压缩、颜色转换、调整大小、颜色转换、统计操作等相关的工具。信号处理应用编程接口是过滤、转换等等。您可以访问核电厂的文档(https://docs.nvidia.com/cuda/npp)，并查看其配置和功能的完整列表。
CUDA 提供了许多基于 NPP 的样本。在本节中，我们将介绍核电厂库的基本用途，并讨论其应用。
首先，我们将介绍 NPP 库如何简化图像处理任务。在此之前，我们应该安装 FreeImage 库，以便能够轻松加载和写入 JPEG 压缩图像文件。有三个选项可用于准备库:
- 从 Ubuntu 档案安装:
$ sudo apt-get install libfreeimage-dev- 从源代码构建并安装:
$ wget http://downloads.sourceforge.net/freeimage/FreeImage3180.zip
$ unzip FreeImage3180.zip
$ cd FreeImage && make -j && sudo make install- 使用已经安装了 CUDA 工具包的库。CUDA 示例代码中的一个 NPP 示例代码
7_CUDALibraries/freeImageInteropNPP使用了 FreeImage 库。对于本示例，NPP 头文件和库文件安装在 CUDA 示例目录中的7_CUDALibrires/common/FreeImage处。如果您不想在您的机器上安装其他二进制文件，您可以使用此选项。
现在，让我们实现基于 NPP 的图像处理应用。完全实现的代码是05_npp/imageFilter.cpp。该文件以头文件开始:
#include <iostream>
#include <iomanip>
#include <cassert>
#include <cstring>
#include <cuda_runtime.h>
#include <npp.h>
#include <FreeImage.h>
#include <helper_timer.h>在该应用中，它具有ImageInfo_t结构，以便于管理图像信息和数据:
struct ImageInfo_t
{
/* image information */
FIBITMAP* dib; // FreeImage bitmap
int nHeight; // image height size
int nWidth; // image width size
int nPitch; // image pitch size
int nBPP; // Bit Per Pixel (i.e. 24 for BGR color)
int nChannel; // number of channels
BYTE* pData; // bytes from freeimage library
/* CUDA */
Npp8u *pDataCUDA; // CUDA global memory for nppi processing
int nPitchCUDA; // image pitch size on CUDA device
};编写LoadImage()函数，加载一个 JPEG 图像。FreeImage库支持任何其他图像格式，可以随意尝试其他图像。然后，我们将用加载的图像数据填充源图像信息管理结构。loadImage()功能实现如下:
void LoadImage(const char *szInputFile, ImageInfo_t &srcImage) {
FIBITMAP *pSrcImageBitmap = FreeImage_Load(FIF_JPEG, szInputFile, JPEG_DEFAULT);
if (!pSrcImageBitmap) {
std::cout << "Couldn't load " << szInputFile << std::endl;
FreeImage_DeInitialise();
exit(1);
}
srcImage.dib = pSrcImageBitmap;
srcImage.nWidth = FreeImage_GetWidth(pSrcImageBitmap);
srcImage.nHeight = FreeImage_GetHeight(pSrcImageBitmap);
srcImage.nPitch = FreeImage_GetPitch(pSrcImageBitmap);
srcImage.nBPP = FreeImage_GetBPP(pSrcImageBitmap);
srcImage.pData = FreeImage_GetBits(pSrcImageBitmap);
assert(srcImage.nBPP == (unsigned int)24); // BGR color image
srcImage.nChannel = 3;
}然后，编写一些 NPPI 助手函数，从图像结构中提供 NPPI 图像大小和 NPPI 感兴趣区域大小数据，如下所示:
NppiSize GetImageSize(ImageInfo_t imageInfo)
{
NppiSize imageSize;
imageSize.width = imageInfo.nWidth;
imageSize.height = imageInfo.nHeight;
return imageSize;
}
NppiRect GetROI(ImageInfo_t imageInfo)
{
NppiRect imageROI;
imageROI.x = 0; imageROI.y = 0;
imageROI.width = imageInfo.nWidth;
imageROI.height = imageInfo.nHeight;
return imageROI;
}然后，让我们实现基于 NPPI 的图像大小调整功能如下。在这个函数中，我们将使用nppiResize_8u_C3R()，这在开始时已经讨论过了。NPP APIs 有命名约定规则来明确说明它们的操作。根据它们的功能类别，它们的命名以nppi开始进行图像处理，以npps 开始进行信号处理。例如，NPP 图像处理函数nppiResize_8u_C3R()以nppi前缀开头，它将三个通道中无符号字符数据类型的输入数据调整到给定的 ROI(您可以在文档中了解有关该约定的更多详细信息):
int ResizeGPU(ImageInfo_t &dstImage, ImageInfo_t &srcImage,
NppiSize &dstSize, NppiRect &dstROI,
NppiSize &srcSize, NppiRect &srcROI, scale)
{
// update output image size
dstSize.width = dstROI.width = dstImage.nWidth;
dstSize.height = dstROI.height = dstImage.nHeight;
nppiResize_8u_C3R(srcImage.pDataCUDA, srcImage.nPitchCUDA,
srcSize, srcROI,
dstImage.pDataCUDA, dstImage.nPitchCUDA,
dstSize, dstROI,
NPPI_INTER_LANCZOS);
return 0;
}为了将性能与中央处理器进行比较，我们将使用 FreeImage 的函数，如下所示:
void ResizeCPU(const char* szInputFile, ImageInfo_t &dstImage) {
FreeImage_Rescale(dib, dstImage.nWidth, dstImage.nHeight, FILTER_LANCZOS3);
}现在，让我们实现main()功能。首先，我们应该初始化 FreeImage 库并加载一个映像:
FreeImage_Initialise();
ImageInfo_t srcImage, dstImage;
LoadImage(szInputFile, srcImage);然后，我们将初始化输入图像的 GPU 内存空间，如下所示。在此过程中，我们使用 NPPI 函数初始化全局内存空间，并使用cudaMemcpy2D()将加载的图像传输到全局内存中:
// copy loaded image to the device memory
srcImage.pDataCUDA =
nppiMalloc_8u_C3(srcImage.nWidth, srcImage.nHeight,
&srcImage.nPitchCUDA);
cudaMemcpy2D(srcImage.pDataCUDA, srcImage.nPitchCUDA,
srcImage.pData, srcImage.nPitch,
srcImage.nWidth * srcImage.nChannel * sizeof(Npp8u),
srcImage.nHeight,
cudaMemcpyHostToDevice);之后，我们将使用调整后的图像大小信息初始化输出内存空间，如下所示:
std::memcpy(&dstImage, &srcImage, sizeof(ImageInfo_t));
dstImage.nWidth *= scaleRatio;
srcImage.nHeight *= scaleRatio;
dstImage.pDataCUDA =
nppiMalloc_8u_C3(dstImage.nWidth, dstImage.nHeight,
&dstImage.nPitchCUDA);然后，我们调用ResizeGPU()和ResizeCPU()函数，这些函数我们已经实现了。对于每个操作，我们将使用cudaEvent来测量 GPU 上的执行时间:
RunNppResize(dstImage, srcImage, dstImageSize, dstROI, srcImageSize, srcROI, scaleRatio);
RunCpuResize(szInputFile, dstImage);为了验证，我们将结果保存到文件中。为此，我们应该创建一个 FreeImage 位图，并将调整后的图像复制到内存空间中。然后，我们可以保存一个输出图像，如下所示:
// Save resized image as file from the device
FIBITMAP *pDstImageBitmap =
FreeImage_Allocate(dstImage.nWidth, dstImage.nHeight,
dstImage.nBPP);
dstImage.nPitch = FreeImage_GetPitch(pDstImageBitmap);
dstImage.pData = FreeImage_GetBits(pDstImageBitmap);
cudaMemcpy2D(dstImage.pData, dstImage.nPitch,
dstImage.pDataCUDA, dstImage.nPitchCUDA,
dstImage.nWidth * dstImage.nChannel * sizeof(Npp8u),
dstImage.nHeight, cudaMemcpyDeviceToHost);
FreeImage_Save(FIF_JPEG, pDstImageBitmap, szOutputFile, JPEG_DEFAULT);之后，我们终于可以终止相关资源:
nppiFree(srcImage.pDataCUDA);
nppiFree(dstImage.pDataCUDA);
FreeImage_DeInitialise();使用链接的 NPP 和 FreeImage 库使用nvcc编译代码:
$ nvcc -run -m64 -std=c++ 11 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -lnppc -lnppif -lnppisu -lnppig -lnppicom -lnpps -lfreeimage -o imageFilter ./imageFilter.cpp因此，当比例因子为 0.5°f 时，图像尺寸减小如下:
$ ls -alh *.jpg
-rw-rw-r-- 1 ubuntu ubuntu 91K Nov 13 22:31 flower.jpg
-rw-rw-r-- 1 ubuntu ubuntu 23K Nov 17 02:46 output.jpg使用 V100 测量的经过时间为0.04576 ms。其时间可能因图形处理器而异:
Rescale flower.jpg in 0.5 ratio.
CPU: 23.857 ms
GPU: 0.04576 ms
Done (generated output.jpg)有关使用 NPP 进行图像处理的更多详细信息，请访问并查看链接文档:http://ON-demand . gputechconf . com/GTC/2014/presentations/HANDS-ON-LAB-S4793-图像处理-使用-npp.pdf 。
核电厂还提供信号处理功能。与图像处理 API 的主要区别在于，它们不需要与图像形状相关的信息。随着我们继续介绍 NPP 函数的基本用法，我们将了解如何从给定的数组中获得和、最小值/最大值、平均值和 L2 归一化分布值。完整编写的代码是05_npp/statisticsNPP.cpp。
首先，让我们从所需的头文件开始:
#include <iostream>
#include <cuda_runtime.h>
#include <npp.h>作为输入数据，我们将使用随机生成的数字:
void GetData(float** buffer, size_t size)
{
(*buffer) = (float*) new float[size];
for (int i = 0; i < size; i++) {
(*buffer)[i] = float(rand() % 0xFFFF) / RAND_MAX;
}
}在我们调用统计运算函数之前，我们需要一个临时的内存空间来进行它们的运算。我们可以使用与操作相关的其他 NPP 函数来获得所需的大小，并且我们可以创建一个公共的工作空间内存空间:
int GetWorkspaceSize(int signalSize)
{
int bufferSize, tempBufferSize;
nppsSumGetBufferSize_32f(signalSize, &tempBufferSize);
bufferSize = std::max(bufferSize, tempBufferSize);
nppsMinGetBufferSize_32f(signalSize, &tempBufferSize);
bufferSize = std::max(bufferSize, tempBufferSize);
nppsMaxGetBufferSize_32f(signalSize, &tempBufferSize);
bufferSize = std::max(bufferSize, tempBufferSize);
nppsMeanGetBufferSize_32f(signalSize, &tempBufferSize);
bufferSize = std::max(bufferSize, tempBufferSize);
nppsNormDiffL2GetBufferSize_32f(signalSize, &tempBufferSize);
bufferSize = std::max(bufferSize, tempBufferSize);
return bufferSize;
}让我们从main()功能开始。首先，我们将从输入数据准备开始，并了解所需的工作空间内存空间。我们将准备两种输入数据类型，并使用 NPP 比较它们的差异:
GetData(&h_input1, buf_size);
GetData(&h_input2, buf_size);
workspace_size = GetWorkspaceSize(buf_size);之后，我们将为输入/输出和工作空间分配 GPU 内存空间。我们还将按如下方式传输输入数据:
cudaMalloc((void **)&d_input1, buf_size * sizeof(float));
cudaMalloc((void **)&d_input2, buf_size * sizeof(float));
cudaMalloc((void **)&d_output, sizeof(float));
cudaMalloc((void **)&d_workspace, workspace_size * sizeof(Npp8u));现在，让我们使用 NPP 函数进行一些简单的统计操作:
nppsSum_32f(d_input1, buf_size, d_output, d_workspace);
cudaMemcpy(&h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);
std::cout << "Sum: " << h_output << std::endl;
nppsMin_32f(d_input1, buf_size, d_output, d_workspace);
cudaMemcpy(&h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);
std::cout << "Min: " << h_output << std::endl;
nppsMax_32f(d_input1, buf_size, d_output, d_workspace);
cudaMemcpy(&h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);
std::cout << "Max: " << h_output << std::endl;
nppsMean_32f(d_input1, buf_size, d_output, d_workspace);
cudaMemcpy(&h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);
std::cout << "Mean: " << h_output << std::endl;NPP 还提供报告两个输入之间差异的功能，如下所示:
nppsNormDiff_L2_32f(d_input1, d_input2, buf_size, d_output, d_workspace);
cudaMemcpy(&h_output, d_output, sizeof(float), cudaMemcpyDeviceToHost);
std::cout << "NormDiffL2: " << h_output << std::endl;然后，我们终止用过的记忆。之后，让我们用以下命令编译代码:
$ nvcc -run -m64 -std=c++ 11 -I/usr/local/cuda/samples/common/inc -gencode arch=compute_70,code=sm_70 -lnppc -lnppif -lnppisu -lnppig -lnppicom -lnpps -o statisticsNPP ./statisticsNPP.cpp结果，我们得到如下结果:
Sum: 0.00100016
Min: 1.30432e-06
Max: 3.04836e-05
Mean: 1.56275e-05
NormDiffL2: 9.46941e-05在这一节中，我们已经讨论了图像处理中的滤波和信号处理中的统计运算。虽然我们已经尝试了简单的应用，但我们可能会发现 NPP 编程比内核实现容易得多。因此，NPP 被应用于许多媒体转码过滤器、洗浴图像处理应用、计算机视觉或深度学习中的图像预处理等。
OpenCV 库是计算机视觉中相当受欢迎的库。它支持图形处理器编程，以提高计算机视觉领域的分辨率。在这一节中，我们将介绍如何使用带有 OpenGL 的 GPU。
要用 CUDA 开始 OpenCV 编程，您需要在启用 CUDA 功能的情况下编译 OpenCV 库。按照以下步骤在 Ubuntu 中启用 OpenCV:
$ sudo apt-get install -y --no-install-recommends \
cmake git libgtk2.0-dev pkg-config libavcodec-dev \
libavformat-dev libswscale-dev \
libatlas-base-dev gfortran libeigen3-dev \
libgtkglext1 libgtkglext1-dev如果您的系统可以使用 X 窗口(不是服务器)，请安装其他软件包以启用 GTK 对话框:
$ sudo apt-get install -y —no-install-recommends \
Libgtkglext1 libgtkglext1-dev下载源代码，并使用以下命令解压它们。这是用 OpenCV 测试的，在撰写本文时，OpenCV 是最新的版本:
# We are install OpenCV 4.1.1
OPENCV_VERSION=4.1.1
OPENCV_DIR=opencv
# Download OpenCV and contrib source codes
mkdir -p ${OPENCV_DIR}
wget -O ${OPENCV_DIR}/opencv-${OPENCV_VERSION}.tar.gz https://github.com/opencv/opencv/archive/${OPENCV_VERSION}.tar.gz
wget -O ${OPENCV_DIR}/opencv_contrib-${OPENCV_VERSION}.tar.gz https://github.com/opencv/opencv_contrib/archive/${OPENCV_VERSION}.tar.gz
# Untar the files
tar -C ${OPENCV_DIR} -xzf ${OPENCV_DIR}/opencv-${OPENCV_VERSION}.tar.gz
tar -C ${OPENCV_DIR} -xzf ${OPENCV_DIR}/opencv_contrib-${OPENCV_VERSION}.tar.gz现在，让我们使用以下命令编译下载的源代码。如果你愿意，你可以放其他选项。它的编译需要一段时间:
# Build the codes and install
cd ${OPENCV_DIR}/opencv-${OPENCV_VERSION}
mkdir build
cd build
cmake -D CMAKE_BUILD_TYPE=RELEASE \
-D CMAKE_INSTALL_PREFIX=/usr/local \
-D ENABLE_PRECOMPILED_HEADERS=OFF \
-D OPENCV_GENERATE_PKGCONFIG=ON \
-D WITH_CUDA=ON -D WITH_CUVID=OFF -D BUILD_opencv_cudacodec=OFF \
-D ENABLE_FAST_MATH=1 \
-D CUDA_FAST_MATH=1 \
-D OPENCV_EXTRA_MODULES_PATH=../../opencv_contrib-${OPENCV_VERSION}/modules \
-D WITH_CUBLAS=1 \
-D PYTHON_DEFAULT_EXECUTABLE=`which python3` \
-D INSTALL_PYTHON_EXAMPLES=ON \
-D BUILD_EXAMPLES=ON ..
make -j$(nproc)
sudo make install -j$(nproc) 要确认安装，请使用以下命令:
$ pkg-config —cflags opencv4 -I/usr/local/include/opencv4/opencv -I/usr/local/include/opencv4在 OpenCV 4 中，CUDA 相关的函数和类是在 CUDA 命名空间中定义的。例如，您可以使用以下命令创建 CUDA 全局内存空间:
cv::cuda::GpuMat cuda_mem = cv::cuda::GpuMat(src.rows, src.cols, CV_8UC1);然后，设备cuda_mem内存空间可以像正常的 CPU 内存类型(cv::Mat)一样处理。
现在，我们将实现一个支持 GPU 的小型 OpenCV 应用，并比较其性能。让我们从包含所需的头文件开始:
#include <iostream>
#include <string>
#include "opencv2/opencv.hpp"以下是使用 OpenCV 实现的主机模糊过滤器:
void BlurHost(std::string filename)
{
cv::Mat src = cv::imread(filename, 1);
cv::Mat dst;
cv::TickMeter tm;
tm.reset();
tm.start();
cv::bilateralFilter(src, dst, 10, 50, 50);
tm.stop();
std::cout << "CPU Time: " << tm.getTimeMilli() << " ms." << std::endl;
cv::imwrite("result_host.jpg", dst);
}这是启用 CUDA 的模糊过滤器实现:
void BlurCuda(std::string filename)
{
cv::Mat src = cv::imread(filename, 1);
cv::Mat dst;
cv::cuda::GpuMat src_cuda = cv::cuda::GpuMat(src.rows,
src.cols, CV_8UC1);
cv::cuda::GpuMat dst_cuda = cv::cuda::GpuMat(src.rows,
src.cols, CV_8UC1);
cv::TickMeter tm;
// warm-up
cv::cuda::bilateralFilter(src_cuda, dst_cuda, 10, 50, 50);
tm.reset();
tm.start();
src_cuda.upload(src);
cv::cuda::bilateralFilter(src_cuda, dst_cuda, 10, 50, 50);
dst_cuda.download(dst);
tm.stop();
std::cout << "GPU Time: " << tm.getTimeMilli()
<< " ms." << std::endl;
cv::imwrite("result_cuda.jpg", dst);
}这个收条代码展示了bilateralFilter()操作如何与主机匹配，以及 CUDA 如何与 CUDA 命名空间匹配。对于 CUDA 内存操作，设备内存使用cv::cuda::GpuMat，设备内存提供upload()和download()成员功能，如cudaMemcpy()。为了测量经过的时间，使用了cv::TickMeter。然后，main()调用两个实现，如下所示:
int main(int argc, char *argv[])
{
std::string filename("flower.JPG");
BlurHost(filename);
BlurCuda(filename);
return 0;
}现在，让我们编译代码。我们应该在你的编译选项中包含使用pkg-config --cflag opencv的 OpenCV 头文件和库。例如，编译选项可以这样编写:
$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc `pkg-config opencv4 --cflags --libs` -o blur ./blur.cpp然后，输出结果如下:
CPU Time: 57.6544 ms.
GPU Time: 2.97784 ms.根据您的系统和图形处理器，执行时间可能会有所不同。
在 OpenCV 中，CUDA 流是用cv::cuda::Stream来管理的。利用这一点，我们可以进行基于多流的流水线 GPU 操作:
- 我们知道，主机内存应该是固定内存，以便进行异步数据传输:
Mat::setDefaultAllocator(cuda::HostMem::getAllocator(cuda::HostMem::PAGE_LOCKED));- 然后，我们将创建多个流，如下所示:
const int num_stream = 4;
cuda::Stream stream[num_stream];- 并且，我们加载源图像，并基于加载的图像信息初始化 GPU 内存，如下所示:
Mat src = imread(filename, 1);
Mat dst;
cuda::GpuMat src_cuda[num_stream], dst_cuda[num_stream];
for (int i = 0; i < num_stream; i++)
src_cuda[i] = cuda::GpuMat(src);- 现在，我们将图像传输到 GPU，模糊图像，并将其与每个流一起传输回主机:
for (int i = 0; i < num_stream; i++) {
src_cuda[i].upload(src, stream[i]);
cuda::bilateralFilter(src_cuda[i], dst_cuda[i], 21, 150.f,
150.f, BORDER_DEFAULT, stream[i]);
dst_cuda[i].download(dst, stream[i]);
}- 然后，我们要同步主机和 GPU。为此，我们将使用
cv::Stream.waitForCompletion()功能，该功能可以在每个流完成向主机的数据传输后对其进行同步:
for (int i = 0; i < num_stream; i++)
stream[i].waitForCompletion();- 为了与中央处理器进行性能比较，我们也称
cv::bilateralFilter()如下:
bilateralFilter(src, dst, 21, 150, 150);其执行时间如下。GPU 执行时间是从多流执行循环到同步的测量时间的平均值:
$ nvcc -run -m64 -gencode arch=compute_70,code=sm_70 -I/usr/local/cuda/samples/common/inc `pkg-config opencv4 --cflags --libs` -o blur ./blur.cpp
CPU Time: 84.8649 ms.
GPU Time: 1.60979 ms.- 为了确认多流操作，我们可以对操作进行概要分析。下面的截图显示了这一点:
Profiling the operation
默认流上的第一个操作是向上翘曲执行，然后是四多流操作。在这里，我们可以看到 GPU 操作是重叠的。因此，平均执行时间比单流执行要短。
我们只在 OpenCV 中介绍了双边过滤。但是，OpenCV 的很多特性都支持 CUDA 加速，让你可以获得 GPU 计算的好处。它的接口与 CPU 版本一致，所以你可以轻松地将你的 CPU 版本迁移到 GPU。
作为入门，有一些来自 GTC 的有用材料:
- http://on-demand . gputechconf . com/GTC/2013/webinar/opencv-GTC-express-shalini-Gupta . pdf
- http://on-demand . gputechconf . com/GTC/2013/webinar/GTC-express-itseez-opencv-webinar . pdf
- http://developer.download.nvidia.com/GTC/PDF/1085_Fung.pdf
建议从 OpenCV 的参考指南开始:https://docs.opencv.org/4.1.1/d2/dbc/cuda_intro.html。
现在很多人用 CUDA 搭配 Python。它不仅作为二进制文件的粘合剂，还使我们能够直接编写 GPU 加速代码。作为一种粘合语言，Python 可以从 CUDA C/C++ 库中调用 API，使用pybind11(https://github.com/pybind/pybind11)或 SWIG(http://swig.org/)。但是，我们必须编写 CUDA C/C++ 代码，并将其集成到 Python 应用中。
然而，有一些 Python 包——Numba、CuPy 和 PyCUDA——支持使用 Python 进行 GPU 编程。它们为 CUDA 内核提供了本地加速的 API 和包装器。换句话说，我们不必编写 C/C++ 代码，也不必花费时间执行集成。Numba 提供了一个矢量化和 CUDA 准时制 ( jit )编译器来加速其运行。它与 NumPy 兼容，因此您可以基于 NumPy 加速您的数值计算代码。多亏了 jit 编译器，您还可以用 Python 编写灵活的 CUDA 代码。CuPy 也与 NumPy 兼容，并加速线性代数算法。它提供了 Pythonic 可编程性和透明的自定义内核编程，例如 Numba。PyCUDA 提供了一个 CUDA C/C++ 接口，这样就可以在 Python 代码中编写和使用 CUDA 内核函数。
numba(https://numba.pydata.org/)翻译 Python 函数在 GPU 上执行，无需任何 C/C++ 编程。
在 Numba 中，您可以通过将 Numba 装饰器应用于目标函数来轻松编写矢量化函数:
from numba import vectorize
@vectorize(["float32(float32, float32, float32)"], target='cuda')
def saxpy(scala, a, b):
return scala * a + b如您所见，装饰器指定参数和返回数据类型，目标指定代码将操作哪个架构。有三种目标:
|
目标
|
描述
|
推荐的数据大小和操作
|
| cuda | 瞄准 NVIDIA GPU | 大于 1 MB，计算密集型操作 |
| parallel | 针对多核中央处理器进行了优化 | 小于 1 MB，正常运行 |
| cpu | 针对单线程操作进行了优化 | 小于 1 KB，低计算密集型操作 |
如果您的函数没有返回值，请使用@guvectorize，并将参数指定为向量。
Numba 的另一个用途是用于@cuda.jit装饰器。这使您能够编写 CUDA 特定的操作，如下所示:
from numba import cuda
@cuda.jit
def matmul(d_c, d_a, d_b):
x, y = cuda.grid(2)
if (x < d_c.shape[0] and y < d_c.shape[1]):
sum = 0
for k in range(d_a.shape[1]):
sum += d_a[x, k] * d_b[k, y]
d_c[x, y] = sumcuda.grid()关键字提供网格级的 CUDA 线程索引，这样就可以用 Python 的方式编写内核代码，比如 CUDA C/C++ 代码。调用 CUDA 内核函数可以如下进行:
matmul[dimGrid, dimBlock](d_c, d_a, d_b)现在，让我们安装这个包，并尝试一些例子。
要在 Python 代码中使用 Numba，您需要安装软件包，并配置环境变量:
$ pip3 install numba
$ export NUMBAPRO_NVVM=/usr/local/cuda/nvvm/lib64/libnvvm.so
$ export NUMBAPRO_LIBDEVICE=/usr/local/cuda/nvvm/libdevice/为了便于将来使用，您需要将环境变量设置放在.bashrc或.zshrc的末尾。如果没有设置，Python 将返回以下消息:
numba.cuda.cudadrv.error.NvvmSupportError: libNVVM cannot be found. Do `conda install cudatoolkit`:
library nvvm not found我们将用一个简单的saxpy操作来测试@vectorize装饰器。这将特定功能转换为并行工作:
- 创建
numba_saxpy.py。 - 导入
numba、numpy和任何其他所需的包:
import numpy as np
from numba import vectorize
from timeit import default_timer as timer- 用
@vectorize装饰器写一个saxpy代码，用'cuda'打靶，以便在 CUDA 设备上工作:
@vectorize(["float32(float32, float32, float32)"], target='cuda')
def saxpy_cuda(scala, a, b):
return scala * a + b- 用
@vecotrize装饰器编写 saxpy 代码，用'parallel'作为目标，在多核处理器(主机)上工作:
@vectorize(["float32(float32, float32, float32)"], target='parallel')
def saxpy_host(scala, a, b):
return scala * a + b- 编写一个操作代码，用一些 NumPy 生成的输入数据调用函数:
scala = 2.0
np.random.seed(2019)
print("size \t\t CUDA \t\t CPU")
for i in range(16,20):
N = 1 << i
a = np.random.rand(N).astype(np.float32)
b = np.random.rand(N).astype(np.float32)
c = np.zeros(N, dtype=np.float32)
# warm-up
c = saxpy_cuda(scala, a, b)
# measuring execution time
start = timer()
c = saxpy_host(scala, a, b)
elapsed_time_host= (timer() - start) * 1e3
start = timer()
c = saxpy_cuda(scala, a, b)
elapsed_time_cuda = (timer() - start) * 1e3
print("[%d]: \t%.3f ms\t %.3f ms" % (N, elapsed_time_cuda, elapsed_time_host))此代码报告不同操作数大小的运行时间:
size CUDA CPU
[65536]: 1.174 ms 0.199 ms
[131072]: 1.362 ms 0.201 ms
[262144]: 2.240 ms 0.284 ms
[524288]: 2.384 ms 0.337 ms在这种情况下，CUDA 表现出比 CPU 更慢的性能，因为操作简单，但是数据传输开销很大。
我们还可以使用@cuda.jit装饰器编写复杂的操作来使用 Numba 在 GPU 上工作:
- 创建
numba_matmul.py。 - 导入
numpy、numba和任何其他所需的包:
import numpy as np
from numba import cuda
from timeit import default_timer as timer- 用
@cuda.jit装饰器写一个矩阵乘法代码:
@cuda.jit
def matmul(d_c, d_a, d_b):
x, y = cuda.grid(2)
if (x < d_c.shape[0] and y < d_c.shape[1]):
sum = 0
for k in range(d_a.shape[1]):
sum += d_a[x, k] * d_b[k, y]
d_c[x, y] = sum在这段代码中，我们使用cuda.grid(dimension_size)来指定网格之间的 CUDA 线程索引，因此，我们可以在 Python 中指定 CUDA 线程的索引。
- 将
a和b矩阵创建为数字矩阵:
N = 8192
a = np.random.rand(N, N).astype(np.float32)
b = np.random.rand(N, N).astype(np.float32)- 将 NumP 生成的数据复制到设备:
d_a = cuda.to_device(a)
d_b = cuda.to_device(b)- 创建将被放置在 CUDA 设备内存中的
c矩阵:
d_c = cuda.device_array((N, N))- 调用矩阵乘法内核函数:
start = timer()
matmul[dimGrid, dimBlock](d_c, d_a, d_b)
elapsed_time_gpu = (timer() - start) * 1e3- 将输出复制到主机:
c = d_c.copy_to_host()- 将 CUDA 操作与主机进行比较:
# matrix multiplication (cpu)
start = timer()
c_host = np.matmul(a, b)
elapsed_time_cpu = (timer() - start) * 1e3
# print elapse times
print("Elapsed Time")
print("GPU: %.3f ms" % elapsed_time_gpu)
print("CPU: %.3f ms" % elapsed_time_cpu)
if (np.allclose(c_host, c)):
print("Done.")
else:
print("GPU and host results are mismatching.")通过一个@cuda.jit装饰器和内置的cuda.grid()关键字，这个示例代码展示了在 Python 中将 Numba 实现为矩阵乘法是多么简单。此代码报告设备和主机上的运行时间:
Elapsed Time
GPU: 104.694 ms
CPU: 1539.005 ms
Done.现在，让我们讨论 CuPy，它在 CUDA 编程中支持更多的 Pythonic 编程。
CuPy(https://cupy.chainer.org)使用 Python 实现线性代数加速，使用 CUDA 库充分利用 GPU。它与 NumPy 兼容，并提供令人愉快的 Pythonic 可编程性。
让我们介绍一下它的安装、基本用法和手动内核开发。
我们可以使用pip使用以下命令安装 CuPy。然后它还安装cupy包和 CUDA 依赖项:
$ pip3 install cupy-cuda101 # CUDA 10.1
$ pip3 install cupy-cuda101 # CUDA 10.0
$ pip3 install cupy-cuda902 # CUDA 9.2现在，让我们介绍一下 CuPy 的基本用法。
我们可以编写一个 saxpy 操作，如下所示:
>>> x = cp.arange(5).astype('f')
>>> x
array([0., 1., 2., 3., 4.], dtype=float32)
>>> y = cp.arange(5).astype('f')
>>> 0.5 * x + y
array([0\. , 1.5, 3\. , 4.5, 6\. ], dtype=float32)我们也可以使用matmul()函数进行矩阵乘法，如下所示:
>>> x = cp.random.uniform(0, 1, (2, 4)).astype('float32')
>>> y = cp.random.uniform(0, 1, (4, 2)).astype('float32')
>>> cp.matmul(x, y)
array([[0.6514087, 0.826463 ],
[0.7826104, 0.2878886]], dtype=float32)正如我们之前讨论的，CuPy 与 NumPy 兼容。基本上，前一个 CuPy 的对象是 CuPy 的数组类型:
>>> x = cp.random.uniform(0, 1, (2, 4)).astype('float32')
>>> type(x)
<class 'cupy.core.core.ndarray'>但是，我们可以使用cupy.asnumpy()函数将其转换为 NumPy 数组，如下所示:
type(cp.asnumpy(x))
<class 'numpy.ndarray'>使用cupy.ascupy()功能也可以进行反向操作。因此，基于这种兼容性，我们可以执行以下操作:
>>> gpu = cp.random.uniform(0, 1, (2, 4)).astype('float32')
>>> cpu = np.random.uniform(0, 1, (2, 4)).astype('float32')
>>> gpu + cp.asarray(cpu)
array([[0.8649391 , 1.1412742 , 1.1280626 , 0.38262686],
[0.44767308, 0.738155 , 0.8397665 , 1.5165564 ]], dtype=float32)
>>> cpu + cp.asnumpy(gpu)
array([[0.8649391 , 1.1412742 , 1.1280626 , 0.38262686],
[0.44767308, 0.738155 , 0.8397665 , 1.5165564 ]], dtype=float32)如您所见，我们可以轻松切换目标计算流程，并且我们可以受益于每个平台的优势。现在，让我们用 CuPy 来介绍定制内核实现。
CuPy 提供了三种类型的定制内核函数:elementwise、reduction 和 raw 内核。elementwise 内核有助于每个元素的自动索引。因此，我们可以只写一个元素的操作。约简内核执行约简操作，同时也执行用户定义的操作。原始内核支持在 Python 代码上直接进行 CUDA C/C++ 内核编程，这样我们就可以在上面定义任何操作。在本节中，我们不会涵盖所有这些内容。但是，您可以从相关文档中了解更多信息—https://docs-cupy . chainer . org/en/stable/tutorial/kernel . html。
让我们讨论用户定义的 elementwise 内核实现。下面是 elementwise 操作的一个示例:
>>> squared_diff = cp.ElementwiseKernel(
... 'float32 x, float32 y',
... 'float32 z',
... 'z = (x - y) * (x - y)',
... 'squared_diff')然后，我们可以在没有显式索引操作的情况下执行 elementwise 操作:
>>> x = cp.random.uniform(0, 1, (2, 4)).astype('float32')
>>> y = cp.random.uniform(0, 1, (2, 4)).astype('float32')
>>> squared_diff(x, y)
array([[0.54103416, 0.01342529, 0.01425287, 0.67101586],
[0.04841561, 0.09939388, 0.46790633, 0.00203693]], dtype=float32)
>>> squared_diff(x, 0.5)
array([[0.23652133, 0.22603741, 0.08065639, 0.00647551],
[0.00029328, 0.07454127, 0.00666 , 0.18399356]], dtype=float32)如您所见，CuPy 提供了一个高度 Pythonic 化的界面，并且易于学习。内部例程很多，也兼容 NumPy—https://docs-cupy . chainer . org/en/stable/reference/routines . html。换句话说，当我们在 NumPy 中需要加速计算时，我们可以考虑使用 CuPy。
现在，我们将介绍 PyCUDA，它提供直接内核编程和隐式内存管理包装。
PyCUDA(https://documen.tician.de/pycuda/)让我们可以用 Python 代码编写 CUDA C/C++ 代码，不需要编译就可以执行。这样，您就可以编写 CUDA 特定操作的 CUDA C/C++ 代码。但是，您必须自己优化这段代码，因为 PyCUDA 不会优化您的内核函数。
这是使用 PyCUDA 生成的一段代码:
import pycuda.autoinit # initialize CUDA devices
from pycuda import driver, compiler, gpuarray
from string import Template
kernel_code_template = Template("""
__global__ void matmul_kernel(float *d_C, float *d_A, float *d_B)
{
int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
...
}
""")
mod = compiler.SourceModule(kernel_code_template.substitute(MATRIX_SIZE=N))
matmul_kernel = mod.get_function("matmul_kernel")
matmul_kernel(driver.Out(C), driver.In(A), driver.In(B), block=(dimBlock, dimBlock, 1), grid=(dimGrid, dimGrid))正如您在这段代码中看到的，我们可以使用相同的 Python 代码编写内核代码。我们还可以使用driver.In()和driver.Out()保留所需数据传输的迹象。这些表明 PyCUDA 应该在调用内核之前传输数据。数据会自动传输，我们还可以按如下方式传输数据:
d_A = driver.to_device(A) # cudaMemcpyHostToDevice
A = driver.from_device_like(d_A) # cudaMemcpyDeviceToHost现在，让我们安装 PyCUDA 并尝试一些简单的示例。
要使用 PyCUDA，还需要安装软件包。从网站下载 PyCUDA 源文件(https://pypi.org/project/pycuda/)。目前正在使用 2019.1.1 版本。
然后按如下方式安装依赖项:
$ sudo apt-get install build-essential python-dev python-setuptools libboost-python-dev libboost-thread-dev
$ tar -xzf pycuda-2019.1.2.tar.gz
$ cd pycuda-2019.1.2
$ python3 ./configure.py --cuda-root=/usr/local/cuda --cudadrv-lib-dir=/usr/lib \
--boost-inc-dir=/usr/include --boost-lib-dir=/usr/lib \
--boost-python-libname=boost_python-py36 --boost-thread-libname=boost_thread
$ python3 setup.py build
$ sudo python3 setup.py install如果你想使用 Python 2，跳过使用 Python 3 的configure.py命令。根据您的 Python 版本，配置命令可能会有所不同。
我们可以通过以下方式使用 PyCUDA 执行矩阵乘法:
- 创建
pycuda_matmul.py文件。 - 按照以下步骤导入所需的包:
import pycuda.autoinit
from pycuda import driver, compiler, gpuarray
import numpy as np
from string import Template
import timeit- 编写一个 CUDA 内核函数代码:
kernel_code_template = Template("""
__global__ void matmul_kernel(float *d_C, float *d_A, float *d_B)
{
int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
float sum = 0.f;
for (int e = 0; e < ${MATRIX_SIZE}; e++)
sum += d_A[idx_y * ${MATRIX_SIZE} + e] * d_B[e * ${MATRIX_SIZE} + idx_x];
d_C[idx_y * ${MATRIX_SIZE} + idx_x] = sum;
}
""")- 使用 NumPy 生成输入/输出矩阵:
N = 8192
np.random.seed(2019)
A = np.random.rand(N, N).astype(np.float32)
B = np.random.rand(N, N).astype(np.float32)
C = np.zeros((N, N), dtype=np.float32)- 编译内核代码:
mod = compiler.SourceModule( \
kernel_code_template.substitute(MATRIX_SIZE=N))- 从编译后的模块中获取内核函数:
matmul_kernel = mod.get_function("matmul_kernel")- 使用主机生成的输入数据创建设备存储器:
d_A = gpuarray.to_gpu(A)
d_B = gpuarray.to_gpu(B)
d_C = gpuarray.zeros((N, N), dtype=np.float32)- 配置网格和块尺寸:
dimBlock = 16
dimGrid = int((N + dimBlock - 1) / dimBlock)- 准备获取图形处理器事件:
start = driver.Event()
stop = driver.Event()- 调用内核函数:
print("Started GPU operation...")
start.record()
matmul_kernel(d_C, d_A, d_B,
block=(dimBlock, dimBlock, 1),
grid=(dimGrid, dimGrid))
stop.record()
stop.synchronize()
gpu_time = stop.time_since(start)
print("GPU Execution Time: %.3f ms" % (gpu_time))- 从主机启动矩阵乘法，并将其与设备的结果进行比较:
print("Started Host operation...")
start = timeit.default_timer()
c_host = np.matmul(A, B)
host_time = timeit.default_timer() - start
print("CPU Execution Time: %.3f ms" % (host_time * 1e3))
if (np.allclose(c_host, d_C.get())):
print("Done.")
else:
print("GPU and host results are mismatching.")此代码还报告设备和主机上的估计时间:
Started GPU operation...
GPU Execution Time: 657.547 ms
Started Host operation...
CPU Execution Time: 1531.133 ms
Done.虽然 PyCUDA 公开了 CUDA C/C++ 内核代码，但这个结果给出了一个提示，即需要手动进行内核优化，因为与 Numba 执行的操作相比，它缺乏性能。
NVBLAS 是一个 CUDA 库，用于 Octave 和 r 等其他包的 BLAS 操作，通过替换 OpenBLAS 执行的操作，Octave 或开发人员和数据科学家可以轻松享受 GPU 性能。在本章中，我们将介绍如何使用 NVBLAS 加速 Octave 和 R。
NVBLAS 是 cuBLAS 操作之上的动态库。cuBLAS 库是线性代数运算的 GPU 实现。它取代了 BLAS 库，因此我们可以轻松地用零编码工作加速任何应用。让我们看看如何从 GEMM 示例代码中做到这一点。
要在 Octave 和 R 中使用 NVBLAS，需要给 NVBLAS 提供一些工作环境变量。为此，让我们创建一个nvblas.conf文件，在那里可以找到目录，我们将使用 Octave 和 R 代码示例。nvblas.conf文件可以编写如下:
NVBLAS_CPU_BLAS_LIB libopenblas.so
NVBLAS_LOGFILE nvblas.log
NVBLAS_GPU_LIST 0
NVBLAS_AUTOPIN_MEM_ENABLED在这个文件中，我们可以看到 NVBLAS 需要知道 CPU 端的 BLAS 库。我们将在这个会话中使用 OpenBLAS，因此我们需要在 Ubuntu 中使用以下命令安装它:
$ sudo apt-get install libopenblas-base libopenblas-dev另外，我们可以通过为NVBLAS_GPU_LIST提供多个 GPU IDs 来获得多 GPU 性能。这本书提供了一个 GPU 执行结果的结果，但是如果你有多个 GPU，尽量提供多个 id。
要在 Octave 和 R 中使用 NVBLAS，我们应该为您的应用执行设置一个环境— LD_PRELOAD=libnvblas.so:
- 对于 Octave 代码，按如下方式执行代码:
$ LD_PRELOAD=libnvblas.so octave sgemm.m- 对于 R 脚本，按如下方式执行您的脚本:
$ LD_PRELOAD=libnvblas.so Rscript sgemm.R当然libnvblas.so文件应该可以从工作目录中访问。位于/usr/local/cuda/lib64/。
NVBLAS 与存档的包兼容。因此，在我们的测试中，使用 Octave-和 R-安装的命令和下面的命令很好地工作:
$ sudo apt-get install octave # for octave installation
$ sudo apt-get install r-base # for R installation现在，让我们尝试使用 NVBLAS，使用 Octave 和 R 语言。
首先，我们将使用 Octave 尝试 NVBLAS。完全实现的代码是08_nvblas/sgemm.m。具体如下:
for i = 1:5
N = 512*(2^i);
A = single(rand(N,N));
B = single(rand(N,N));
start = clock();
C = A * B;
elapsedTime = etime(clock(), start);
gFlops = 2*N*N*N/(elapsedTime * 1e+9);
printf("Elapsed Time [%d]: %.3f ms, %.3f GFlops\n", N, elapsedTime, gFlops);
end对于 GPU 操作，使用以下命令执行 Octave 脚本，并将性能与 GPU、默认情况下的 NVBLAS 环境库和 CPU 进行比较:
$ LD_PRELOAD=libnvblas.so octave sgemm.m然后，我们可以用octave sgemm.m命令启动这个。输出结果如下:
| 中央处理器 | GPU V100 | |
Elapsed Time [1024]: 0.011 ms, 188.909 GFlopsElapsed Time [2048]: 0.075 ms, 228.169 GFlopsElapsed Time [4096]: 0.212 ms, 647.022 GFlopsElapsed Time [8192]: 1.158 ms, 949.763 GFlopsElapsed Time [16384]: 7.292 ms, 1206.241 GFlops
|
Elapsed Time [1024]: 0.010 ms, 208.346 GFlopsElapsed Time [2048]: 0.024 ms, 721.731 GFlopsElapsed Time [4096]: 0.094 ms, 1465.538 GFlopsElapsed Time [8192]: 0.582 ms, 1889.193 GFlopsElapsed Time [16384]: 4.472 ms, 1967.037 GFlops
|
如您所见，随着矩阵的大小变大，GPU 显示出更高的计算吞吐量。
现在，我们将在以下步骤的帮助下，针对 R 语言尝试 NVBLAS:
- 首先，让我们编写一个
sgemm.R文件，它执行一个点操作:
set.seed(2019)
for(i in seq(1:5)) {
N = 512*(2^i)
A = matrix(rnorm(N^2, mean=0, sd=1), nrow=N)
B = matrix(rnorm(N^2, mean=0, sd=1), nrow=N)
elapsedTime = system.time({C = A %*% B})[3]
gFlops = 2*N*N*N/(elapsedTime * 1e+9);
print(sprintf("Elapsed Time [%d]: %3.3f ms, %.3f GFlops", N, elapsedTime, gFlops))
}- 使用以下命令执行 R 脚本，并比较性能:
$ LD_PRELOAD=libnvblas.so Rscript sgemm.R示例代码运行了几次，同时增加了数据大小。下表显示了前面命令的输出:
| 中央处理器 | GPU V100 | |
Elapsed Time [1024]: 0.029 ms, 74.051 GFlopsElapsed Time [2048]: 0.110 ms, 156.181 GFlopsElapsed Time [4096]: 0.471 ms, 291.802 GFlopsElapsed Time [8192]: 2.733 ms, 402.309 GFlopsElapsed Time [16384]: 18.291 ms, 480.897 GFlops
|
Elapsed Time [1024]: 0.034 ms, 63.161 GFlopsElapsed Time [2048]: 0.063 ms, 272.696 GFlopsElapsed Time [4096]: 0.286 ms, 480.556 GFlopsElapsed Time [8192]: 1.527 ms, 720.047 GFlopsElapsed Time [16384]: 9.864 ms, 891.737 GFlops
|
从结果可以看出 CPU 和 GPU 的性能差距。此外，我们能够确定，当我们增加样本量时，GPU 的性能增益会增加。
如果你对 GPU 的 R 加速感兴趣，请访问 NVIDIA 开发博客:https://devblogs.nvidia.com/accelerate-r-applications-cuda/
MATLAB 是一个高效的高级数值分析工具，具有各种工具和功能。该工具通过其并行计算工具箱从早期就支持 CUDA。本节将向我们展示如何使用该工具生成 CUDA 代码。
为了实现 GPU 加速，我们需要安装带有并行计算工具箱的 MATLAB。如果您已经有了 MATLAB，请检查您的许可证是否涵盖了并行计算工具箱。如果没有，可以试试 MATLAB 评估代码。从 MATLAB 的评估网站，你可以下载任何种类的软件包，除了控制系统。大多数包都包含并行计算工具箱，所以您可以尝试这样做。但是如果你不考虑使用 MATLAB，你可以跳过这一节。
当我们使用 MATLAB 代码在 GPU 上工作时，您需要使用 gpuArray 创建一个设备内存。就像 Numba 和 PyCUDA 将它们的主机数据发送到设备一样，MATLAB 的gpuArray()创建一个设备内存并将给定的主机数据传输到设备:
d_A = gpuArray(A);本课程将假设您已经安装了 MATLAB 和并行计算工具箱。在本节中，我们将重点介绍示例代码的实现，并比较主机和 GPU 的性能:
- 我们写一个
host.m文件，可以在 CPU 上工作。代码如下:
N = 8192;
A = single(rand(N,N));
B = single(rand(N,N));
start = clock();
C = A * B;
elapsedTime = etime(clock(), start);
gFlops = 2*N*N*N/(elapsedTime * 1e+9);
fprintf("Elapsed Time: %.3f ms, %.3f GFlops\n", elapsedTime, gFlops);现在，让我们用下面的命令来执行这两个实现。这是 MATLAB 的命令，它的输出:
$ matlab -r "run('host.m'); exit;" -nodisplay
Elapsed Time: 6.421 ms, 171.243 Gflops- 然后，我们写一个
cuda.m文件，在 GPU 上工作。我们只需将gpuArray()应用于输入矩阵，如下所示:
N = 8192;
A = single(rand(N,N));
B = single(rand(N,N));
d_A = gpuArray(A); 'GPU memory allocation
d_B = gpuArray(B); 'GPU memory allocation
start = clock();
d_C = d_A * d_B;
elapsedTime = etime(clock(), start);
gFlops = 2*N*N*N/(elapsedTime * 1e+9);
fprintf("Elapsed Time: %.3f ms, %.3f GFlops\n", elapsedTime, gFlops);这是 GPU 版本执行代码，执行结果:
$ matlab -r "run('cuda.m'); exit;" -nodisplay
Elapsed Time: 0.179 ms, 6140.739 Gflops.我们可以看到，相对于 CPU，GPU 表现出了更高的性能。
MathWorks 提供了大量用 MATLAB 进行 GPU 计算的例子。如果您想了解更多信息，请访问他们的网站:https://www . mathworks . com/examples/parallel-computing/category/GPU-computing。
在本章中，我们已经介绍了使用 CUDA 库和其他兼容语言的 CUDA 编程方法。我们还介绍了 cuBLAS 的基本用途及其混合精度操作特性。此外，我们还探索了 cuRAND、cuFFT、NPP 和 OpenCV 库。多亏了这些库，我们可以不费吹灰之力实现 GPU 应用，正如本章开头所讨论的。
我们已经使用与 CUDA 兼容的其他语言实现了一些 GPU 应用。首先，我们介绍了几个支持 Python 和 CUDA 互操作的 Python 包。它们提供 Python 编程能力以及与其他 Python 特性的兼容性。然后，我们介绍了其他科学计算语言中的 CUDA 加速，例如 Octave、R 和 MATLAB。
现在，我们还有一个 GPU 编程方法要介绍——OpenACC。有了这个，我们可以转换原始的 C/C++ 和 Fortran 主机代码，使用指令如#pragma acc kernels在 GPU 上工作。我们将在下一章讨论这个问题。