CUDA并行计算平台的C/C++接口的简单介绍

　　本文是 CUDA C 和 C ++的一个系列，它是 CUDA 并行计算平台的 C / C ++接口。本系列文章假定您熟悉 C 语言编程。我们将针对 Fortran 程序员运行一系列关于 CUDA Fortran 的文章。这两个系列将介绍 CUDA 平台上并行计算的基本概念。从这里起，除非我另有说明，我将用“ CUDA C ”作为“ CUDA C 和 C ++”的速记。 CUDA C 本质上是 C / C ++，具有几个扩展，允许使用并行的多个线程在 GPU 上执行函数。

　　CUDA 编程模型基础

　　在我们跳转到 CUDA C 代码之前， CUDA 新手将从 CUDA 编程模型的基本描述和使用的一些术语中受益。

　　CUDA 编程模型是一个异构模型，其中使用了 CPU 和 GPU 。在 CUDA 中， host 指的是 CPU 及其存储器， device 是指 GPU 及其存储器。在主机上运行的代码可以管理主机和设备上的内存，还可以启动在设备上执行的函数 kernels 。这些内核由许多 GPU 线程并行执行。

　　鉴于 CUDA 编程模型的异构性， CUDA C 程序的典型操作序列是：

　　声明并分配主机和设备内存。

　　初始化主机数据。

　　将数据从主机传输到设备。

　　执行一个或多个内核。

　　将结果从设备传输到主机。

　　记住这个操作序列，让我们看一个 CUDA C 示例。

　　第一个 CUDA C 程序

     在最近的一篇文章中，我演示了 萨克斯比的六种方法 ，其中包括一个 CUDA C 版本。 SAXPY 代表“单精度 A * X + Y ”，是并行计算的一个很好的“ hello world ”示例。在这篇文章中，我将剖析 CUDA C SAXPY 的一个更完整的版本，详细解释它的作用和原因。完整的 SAXPY 代码是：

#include 

__global__
void saxpy(int n, float a, float *x, float *y)
{
  int i = blockIdx.x*blockDim.x + threadIdx.x;
  if (i < n) y[i] = a*x[i] + y[i];
}

int main(void)
{
  int N = 1<<20;
  float *x, *y, *d_x, *d_y;
  x = (float*)malloc(N*sizeof(float));
  y = (float*)malloc(N*sizeof(float));

  cudaMalloc(&d_x, N*sizeof(float)); 
  cudaMalloc(&d_y, N*sizeof(float));

  for (int i = 0; i < N; i++) {
    x[i] = 1.0f;
    y[i] = 2.0f;
  }

  cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

  // Perform SAXPY on 1M elements
  saxpy<<<(N+255)/256, 256>>>(N, 2.0f, d_x, d_y);

  cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

  float maxError = 0.0f;
  for (int i = 0; i < N; i++)
    maxError = max(maxError, abs(y[i]-4.0f));
  printf("Max error: %f
", maxError);

  cudaFree(d_x);
  cudaFree(d_y);
  free(x);
  free(y);
}

函数 saxpy 是在 GPU 上并行运行的内核， main 函数是宿主代码。让我们从宿主代码开始讨论这个程序。

主机代码

main 函数声明两对数组。

  float *x, *y, *d_x, *d_y;
  x = (float*)malloc(N*sizeof(float));
  y = (float*)malloc(N*sizeof(float));

  cudaMalloc(&d_x, N*sizeof(float)); 
  cudaMalloc(&d_y, N*sizeof(float));

　　指针x和y指向以典型方式使用malloc分配的主机阵列，d_x和d_y数组指向从CUDA运行时API使用cudaMalloc函数分配的设备数组。CUDA中的主机和设备有独立的内存空间，这两个空间都可以从主机代码进行管理（CUDAC内核也可以在支持它的设备上分配设备内存）。

　　然后，主机代码初始化主机数组。在这里，我们设置了一个 1 数组，以及一个 2 数组。

  for (int i = 0; i < N; i++) {
    x[i] = 1.0f;
    y[i] = 2.0f;
  }

为了初始化设备数组，我们只需使用 cudaMemcpy 将数据从 x 和 y 复制到相应的设备数组 d_x 和 d_y ，它的工作方式与标准的 C memcpy 函数一样，只是它采用了第四个参数，指定了复制的方向。在本例中，我们使用 cudaMemcpyHostToDevice 指定第一个（目标）参数是设备指针，第二个（源）参数是主机指针。

  cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

在运行内核之后，为了将结果返回到主机，我们使用 cudaMemcpy 和 cudaMemcpyDeviceToHost ，从 d_y 指向的设备数组复制到 y 指向的主机数组。

cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

启动内核

cord [EZX13 内核由以下语句启动：

saxpy<<<(N+255)/256, 256>>>(N, 2.0, d_x, d_y);

　　三个 V 形符号之间的信息是 执行配置 ，它指示有多少设备线程并行执行内核。在 CUDA 中，软件中有一个线程层次结构，它模仿线程处理器在 GPU 上的分组方式。在 CUDA 编程模型中，我们谈到启动一个 grid 为 螺纹块 的内核。执行配置中的第一个参数指定网格中线程块的数量，第二个参数指定线程块中的线程数。

　　线程块和网格可以通过为这些参数传递 dim3 （一个由 CUDA 用 x 、 y 和 z 成员定义的简单结构）值来生成一维、二维或三维的线程块和网格，但是对于这个简单的示例，我们只需要一维，所以我们只传递整数。在本例中，我们使用包含 256 个线程的线程块启动内核，并使用整数算术来确定处理数组（ （N+255）/256 ）的所有 N 元素所需的线程块数。

　　对于数组中的元素数不能被线程块大小平均整除的情况，内核代码必须检查内存访问是否越界。

清理

　　完成后，我们应该释放所有分配的内存。对于使用 cudaMalloc（） 分配的设备内存，只需调用 cudaFree（） 。对于主机内存，请像往常一样使用 free（） 。

cudaFree(d_x);
  cudaFree(d_y);
  free(x);
  free(y);

设备代码

现在我们继续讨论内核代码。

__global__
void saxpy(int n, float a, float *x, float *y)
{
  int i = blockIdx.x*blockDim.x + threadIdx.x;
  if (i < n) y[i] = a*x[i] + y[i];
}

　　在 CUDA 中，我们使用 __global__ de __global__ 说明符定义诸如 Clara 这样的内核。设备代码中定义的变量不需要指定为设备变量，因为假定它们驻留在设备上。在这种情况下， n 、 a 和 i 变量将由每个线程存储在寄存器中，指针 x 和 y 必须是指向设备内存地址空间的指针。这确实是真的，因为当我们从宿主代码启动内核时，我们将 d_x 和 d_y 传递给了内核。但是，前两个参数 n 和 a 没有在主机代码中显式传输到设备。因为函数参数在 C / C ++中是默认通过值传递的，所以 CUDA 运行时可以自动处理这些值到设备的传输。 CUDA 运行时 API 的这一特性使得在 GPU 上启动内核变得非常自然和简单——这几乎与调用 C 函数一样。

　　在我们的 saxpy 内核中只有两行。如前所述，内核由多个线程并行执行。如果我们希望每个线程处理结果数组的一个元素，那么我们需要一种区分和标识每个线程的方法。 CUDA 定义变量 blockDim 、 blockIdx 和 threadIdx 。这些预定义变量的类型为 dim3 ，类似于主机代码中的执行配置参数。预定义变量 blockDim 包含在内核启动的第二个执行配置参数中指定的每个线程块的维度。预定义变量 threadIdx 和 blockIdx 分别包含线程块中线程的索引和网格中的线程块的索引。表达式：

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

　　生成用于访问数组元素的全局索引。我们在这个例子中没有使用它，但是还有一个 gridDim ，它包含在启动的第一个执行配置参数中指定的网格维度。

　　在使用该索引访问数组元素之前，将根据元素的数量 n 检查其值，以确保没有越界内存访问。如果一个数组中的元素数不能被线程块大小平均整除，并且结果内核启动的线程数大于数组大小，则需要进行此检查。内核的第二行执行 SAXPY 的元素级工作，除了边界检查之外，它与 SAXPY 主机实现的内部循环相同。

if (i < n) y[i] = a*x[i] + y[i];

编译和运行代码

　　CUDA C 编译器 nvcc 是 NVIDIA CUDA 工具箱 的一部分。为了编译我们的 SAXPY 示例，我们将代码保存在一个扩展名为。 cu 的文件中，比如说 saxpy.cu 。然后我们可以用 nvcc 编译它。

nvcc -o saxpy saxpy.cu

然后我们可以运行代码：

% ./saxpy
Max error: 0.000000

总结与结论

　　通过对 SAXPY 的一个简单的 CUDA C 实现的演练，您现在了解了编程 CUDA C 的基本知识。将 C 代码“移植”到 CUDA C 只需要几个 C 扩展：设备内核函数的 __global__ de Clara 说明符；启动内核时使用的执行配置；内置的设备变量 blockDim 、 blockIdx 和 threadIdx 用来识别和区分并行执行内核的 GPU 线程。

　　异类 CUDA 编程模型的一个优点是，将现有代码从 C 移植到 CUDA C 可以逐步完成，一次只能移植一个内核。

　　在本系列的下一篇文章中，我们将研究一些性能度量和度量。

　　关于作者

　　Mark Harris 是 NVIDIA 杰出的工程师，致力于 RAPIDS 。 Mark 拥有超过 20 年的 GPUs 软件开发经验，从图形和游戏到基于物理的模拟，到并行算法和高性能计算。当他还是北卡罗来纳大学的博士生时，他意识到了一种新生的趋势，并为此创造了一个名字： GPGPU （图形处理单元上的通用计算）。

　　审核编辑：郭婷