如果你正在寻找一篇详细介绍CUDA基础的文章，请花一小时阅读此文。文章链接：https://developer.nvidia.com/blog/even-easier-introduction-cuda/…

更加简单的CUDA入门指南（更新版）

来源：https://developer.nvidia.com/blog/even-easier-introduction-cuda/
注意：本文最初发表于2017年1月25日，但已进行修订以反映最新更新。

本文是一份超级简单的CUDA入门指南，CUDA是NVIDIA推出的流行并行计算平台和编程模型。我曾于2013年写过一篇《简单的CUDA入门指南》(https://developer.nvidia.com/blog/easy-introduction-cuda-c-and-c/)，多年来一直广受欢迎。但如今CUDA编程变得更加简单，GPU速度也更快了，所以是时候更新一版（甚至更简单的）入门指南了。

CUDA C++只是使用CUDA创建大规模并行应用的众多方式之一。它让你能够利用强大的C++编程语言，开发由GPU上数千个并行线程加速的高性能算法。许多开发者已通过这种方式加速了他们计算密集型和带宽密集型的应用，包括支撑当前人工智能革命（即深度学习）的库和框架（https://developer.nvidia.com/deep-learning）。

如果你已听说过CUDA，并有兴趣学习如何在你的应用中使用它，那么作为一名C++程序员，这篇博文应该能为你提供一个良好的开端。要跟着操作，你需要一台带有支持CUDA的GPU的电脑（Windows、WSL或64位Linux均可，任何NVIDIA GPU理论上都行），或者一个带有GPU的云实例（AWS、Azure、Google Colab等云服务提供商都有提供，参见https://marketplace.nvidia.com/en-us/enterprise/cloud-solutions/）。你还需要安装免费的CUDA工具包（https://developer.nvidia.com/cuda-toolkit）。

我们开始吧！

从简单开始

https://developer.nvidia.com/blog/even-easier-introduction-cuda/#starting_simple

我们先从一个简单的C++程序开始，该程序将两个各有100万个元素的数组相加。

#include <iostream>
#include <math.h>

// function to add the elements of two arrays
void add(int n, float *x, float *y)
{
    for (int i = 0; i < n; i++)
        y[i] = x[i] + y[i];
}

int main(void)
{
    int N = 1<<20; // 1M elements

    float *x = new float[N];
    float *y = new float[N];

    // initialize x and y arrays on the host
    for (int i = 0; i < N; i++) {
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    // Run kernel on 1M elements on the CPU
    add(N, x, y);

    // Check for errors (all values should be 3.0f)
    float maxError = 0.0f;
    for (int i = 0; i < N; i++)
        maxError = fmax(maxError, fabs(y[i]-3.0f));

    std::cout << "Max error: " << maxError << std::endl;

    // Free memory
    delete [] x;
    delete [] y;

    return 0;
}

首先，编译并运行这个C++程序。将上述代码保存到名为add.cpp的文件中，然后用你的C++编译器编译。我在Linux上使用g++，但你可以在Windows上使用MSVC（或WSL上的g++）。然后运行它：

> ./add
Max error: 0.000000

（在Windows上，你可能想将可执行文件命名为add.exe，并用.\add运行。）

如预期，程序打印出求和没有错误，然后退出。现在我想让这个计算在GPU的众多核心上并行运行。其实迈出第一步非常简单。首先，我只需将add函数转换成一个可以在GPU上运行的函数，在CUDA中这被称为内核。为此，我只需在函数前加上__global__限定符，它告诉CUDA C++编译器这是一个在GPU上运行、可从CPU代码调用的函数。

// Kernel function to add the elements of two arrays
__global__ void add(int n, float *sum, float *x, float *y)
{
    for (int i = 0; i < n; i++)
        sum[i] = x[i] + y[i];
}

这个__global__函数被称为CUDA内核，并在GPU上运行。在GPU上运行的代码通常称为设备代码，而在CPU上运行的代码称为主机代码。

CUDA中的内存分配

https://developer.nvidia.com/blog/even-easier-introduction-cuda/#memory-allocation

为了在GPU上进行计算，我需要分配GPU可访问的内存。CUDA中的统一内存（https://developer.nvidia.com/blog/unified-memory-in-cuda-6/）使这变得简单，它提供了一个单一的内存空间，系统中的所有GPU和CPU都可以访问。要分配统一内存，调用cudaMallocManaged()，它会返回一个指针，你可以从主机（CPU）代码或设备（GPU）代码中访问。要释放数据，只需将该指针传递给cudaFree()。我只需要将上述代码中的new调用替换为cudaMallocManaged()调用，并将delete []调用替换为cudaFree调用。

// Allocate Unified Memory -- accessible from CPU or GPU
float *x, *y, *sum;
cudaMallocManaged(&x, N*sizeof(float));
cudaMallocManaged(&y, N*sizeof(float));
...
// Free memory
cudaFree(x);
cudaFree(y);

最后，我需要启动add()内核，即在GPU上调用它。CUDA内核启动使用三尖括号语法<<<>>>指定。我只需在调用add时，在参数列表前加上它。

add<<<1, 1>>>(N, sum, x, y);

很简单！我很快会详细介绍尖括号内的内容；目前你只需要知道这行代码会启动一个GPU线程来运行add()。还有一件事：我需要让CPU等待内核完成，然后才能访问结果（因为CUDA内核启动不会阻塞调用CPU线程）。为此，我只需在CPU上进行最终错误检查之前调用cudaDeviceSynchronize()。以下是完整的代码：

#include <iostream>
#include <math.h>

// Kernel function to add the elements of two arrays
__global__ void add(int n, float *x, float *y)
{
    for (int i = 0; i < n; i++)
        y[i] = x[i] + y[i];
}

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

    // Allocate Unified Memory – accessible from CPU or GPU
    cudaMallocManaged(&x, N*sizeof(float));
    cudaMallocManaged(&y, N*sizeof(float));

    // initialize x and y arrays on the host
    for (int i = 0; i < N; i++) {
        x[i] = 1.0f;
        y[i] = 2.0f;
    }

    // Run kernel on 1M elements on the GPU
    add<<<1, 1>>>(N, x, y);

    // Wait for GPU to finish before accessing on host
    cudaDeviceSynchronize();

    // Check for errors (all values should be 3.0f)
    float maxError = 0.0f;
    for (int i = 0; i < N; i++) {
        maxError = fmax(maxError, fabs(y[i]-3.0f));
    }

    std::cout << "Max error: " << maxError << std::endl;

    // Free memory
    cudaFree(x);
    cudaFree(y);

    return 0;
}

CUDA文件的扩展名为.cu。因此，将这段代码保存到名为add.cu的文件中，并用nvcc（CUDA C++编译器）编译它。

> nvcc add.cu -o add_cuda
> ./add_cuda
Max error: 0.000000

这只是第一步，因为按照目前写法，这个内核只对单个线程正确，因为运行它的每个线程都会对整个数组进行加法操作。此外，由于多个并行线程会同时读写同一位置，存在竞争条件。注意：在Windows上，你需要确保在Microsoft Visual Studio的项目配置中将平台设置为x64。

性能分析

https://developer.nvidia.com/blog/even-easier-introduction-cuda/#profile-it

要了解内核运行时间，一个好方法是使用NSight Systems CLI工具nsys。我们只需在命令行输入nsys profile -t cuda --stats=true ./add_cuda。但这会产生冗长的统计信息，而在本文中我们只想了解内核运行所需的时间。因此，我编写了一个名为nsys_easy的包装脚本，只输出我们需要的输出，并防止nsys产生会杂乱源目录的中间文件。该脚本可在GitHub上获得（https://github.com/harrism/nsys_easy）。只需下载nsys_easy并将其放在PATH中的某个位置（甚至放在当前目录）。（注意：我已从输出中删除了一些统计信息，以更好地适应本网站宽度。）

> nsys_easy ./add_cuda
Max error: 0
Generating '/tmp/nsys-report-bb25.qdstrm'
[1/1] [========================100%]
nsys_easy.nsys-rep
Generated: /home/nfs/mharris/src/even_easier/nsys_easy.nsys-rep
Generating SQLite file nsys_easy.sqlite from nsys_easy.nsys-rep
Processing 1259 events: [======================================100%]
Processing [nsys_easy.sqlite] with [cuda_gpu_sum.py]...
** CUDA GPU Summary (Kernels/MemOps) (cuda_gpu_sum):
Time (%) Total Time (ns)  Instances  Category          Operation
------- ---------------  ---------  ----------------  --------------------------
98.5    75,403,544              1  CUDA_KERNEL       add(int, float *, float *)
1.0       768,480             48  MEMORY_OPER       [memcpy Unified H2D]
0.5       352,787             24  MEMORY_OPER       [memcpy Unified D2D]

CUDA GPU摘要表显示对add的一次调用。在NVIDIA T4 GPU上大约需要75ms。让我们通过并行化让它更快。

拿起线程

https://developer.nvidia.com/blog/even-easier-introduction-cuda/#picking-up-the-threads

既然你已经用单个线程运行了一个执行一些计算的内核，那么如何使其并行化呢？关键在于CUDA的<<<1, 1>>>语法。这称为执行配置，它告诉CUDA运行时在GPU上启动时使用多少个并行线程。这里有两个参数，但让我们从更改第二个参数开始：一个线程块中的线程数量。CUDA GPU使用大小为32的倍数的线程块来运行内核；256个线程是一个合理的选择。

add<<<1, 256>>>(N, x, y);

如果我只做这一处更改就运行代码，那么每个线程都会重复执行整个计算，而不是将计算分摊到并行线程上。为了正确地做到这一点，我需要修改内核。CUDA C++提供了一些关键字，让内核能够获取正在运行线程的索引。具体来说，threadIdx.x包含当前线程在其块内的索引，blockDim.x包含块中的线程数。我将修改循环，以并行线程的方式跨数组进行步进。

__global__ void add(int n, float *x, float *y)
{
    int index = threadIdx.x;
    int stride = blockDim.x;
    for (int i = index; i < n; i += stride)
        y[i] = x[i] + y[i];
}

add函数没有太大变化。实际上，将index设为0、stride设为1，它在语义上与第一个版本相同。将文件保存为add_block.cu，再次用nsys_easy编译并运行。在本文的剩余部分，我将只显示输出中的相关行。

Time (%) Time (ns)  Instances  Category          Operation
------- ---------  ---------  ----------------  ----------------------
79.0    4,221,011         1  CUDA_KERNEL       add(int, float *, float *)

这是一个巨大的加速（从75ms降到4ms），但并不意外，因为执行从1个线程增加到了256个线程。让我们继续前进，以获得更高的性能。

突破块的限制

https://developer.nvidia.com/blog/even-easier-introduction-cuda/#out-of-the-blocks

CUDA GPU拥有许多并行处理器，它们被分组为流多处理器（SM）。每个SM可以运行多个并发的线程块，但每个线程块只能在一个SM上运行。例如，基于Turing GPU架构的NVIDIA T4（https://www.nvidia.com/en-us/data-center/tesla-t4/）GPU具有40个SM和2560个CUDA核心，每个SM最多可支持1024个活跃线程。为了充分利用所有这些线程，我应该用多个线程块启动内核。现在你可能已经猜到，执行配置的第一个参数指定了线程块的数量。所有这些并行线程块共同构成了所谓的网格。由于我有N个元素需要处理，每个块有256个线程，我只需计算至少需要N个线程的块数。只需将N除以块大小（注意向上取整，以防N不是blockSize的倍数）。

int blockSize = 256;
int numBlocks = (N + blockSize - 1) / blockSize;
add<<<numBlocks, blockSize>>>(N, x, y);

网格块线程图像索引 CUDA
图1. CUDA内核中的网格、块和线程索引（一维）。

我还需要更新内核代码，以考虑整个网格的线程块。CUDA提供了gridDim.x（网格中的块数）和blockIdx.x（当前线程块在网格中的索引）。图1说明了在CUDA中使用blockDim.x、gridDim.x和threadIdx.x对数组（一维）进行索引的方法。其思想是，每个线程通过计算到其块开头的偏移量（块索引乘以块大小：blockIdx.x * blockDim.x）并加上线程在其块内的索引（threadIdx.x）来获得其索引。代码blockIdx.x * blockDim.x + threadIdx.x是CUDA的惯用写法。

__global__ void add(int n, float *x, float *y)
{
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;
    for (int i = index; i < n; i += stride)
        y[i] = x[i] + y[i];
}

更新后的内核还将stride设置为网格中的总线程数（blockDim.x * gridDim.x）。CUDA内核中的这种循环通常称为网格步幅循环（https://developer.nvidia.com/blog/cuda-pro-tip-write-flexible-kernels-grid-stride-loops/）。将文件保存为add_grid.cu，再次用nsys_easy编译并运行。

Time (%) Time (ns)  Instances  Category          Operation
------- ---------  ---------  ----------------  ------------------------
79.6    4,514,384         1  CUDA_KERNEL       add(int, float *, float *)

嗯，这很有趣。这一变化并没有带来加速，甚至可能略有变慢。为什么会这样？如果计算能力增加了40倍（SM数量），却没有降低总时间，那么计算就不是瓶颈。

统一内存预取

https://developer.nvidia.com/blog/even-easier-introduction-cuda/#unified_memory_prefetching

如果我们查看性能分析器的完整摘要表，就能发现瓶颈的线索：

Time (%) Time (ns)  Instances  Category          Operation
------- ---------  ---------  ----------------  ----------------------
79.6    4,514,384         1  CUDA_KERNEL       add(int, float *, float *)
14.2      807,245        64  MEMORY_OPER       [CUDA memcpy Unified H2D]
6.2       353,201        24  MEMORY_OPER       [CUDA memcpy Unified D2H]

这里我们看到有64次主机到设备（H2D）和24次设备到主机（D2H）的“统一”memcpy操作。但代码中没有显式的memcpy调用。CUDA中的统一内存是虚拟内存。单个虚拟内存页面可以驻留在系统中任何设备（GPU或CPU）的内存中，这些页面会在需要时被迁移。这个程序首先在CPU上用for循环初始化数组，然后启动内核，由GPU读取和写入这些数组。由于内核运行时所有内存页面都驻留在CPU上，因此会发生多次缺页中断，硬件会在缺页时将这些页面迁移到GPU内存。这就造成了内存瓶颈，这也是我们没有看到加速的原因。迁移成本很高，因为缺页中断是单独发生的，GPU线程在等待页面迁移时会停止运行。

由于我知道内核需要哪些内存（x和y数组），我可以使用预取来确保数据在内核需要之前已经位于GPU上。为此，我在启动内核之前使用cudaMemPrefetchAsync()函数：

// Prefetch the x and y arrays to the GPU
cudaMemPrefetchAsync(x, N*sizeof(float), 0, 0);
cudaMemPrefetchAsync(y, N*sizeof(float), 0, 0);

使用性能分析器运行此代码会产生以下输出。现在内核运行时间不到50微秒！

Time (%) Time (ns)  Instances  Category          Operation
------- ---------  ---------  ----------------  ----------------------
63.2      690,043         4  MEMORY_OPER       [CUDA memcpy Unifi