How To Write A CUDA Program

Stephen Jones, GTC 2023

目录

• 1. 引言与背景
• 1.1 CUDA 发展时间线
• 1.2 个人背景与历史回顾
• 1.3 本次演讲的核心问题
• 1.4 不同抽象层级下的编程范式对比
• 1.5 NVIDIA 计算生态系统
• 2. CUDA 编程基础
• 2.1 环境搭建与示例获取
• 2.2 CUDA 程序结构与执行模型
• 2.3 CUDA 编程模型：网格 → 线程块 → 线程
• 3. 案例研究：GPU 驱动的递归地形生成
• 3.1 核心算法：二维递归细分
• 3.2 将问题映射到 GPU
• 3.3 快速实现：修改 OceanFFT 示例
• 3.4 挑战：并行伪随机数生成
• 3.5 实现递归地形生成
• 4. 高级技术与优化
• 4.1 地形动画与变形
• 4.2 视觉效果增强与无缝景观
• 4.3 性能优化：将数据保留在 GPU
• 4.4 动态高度与智能内存管理
• 5. 开发与调试
• 5.1 调试 GPU 程序
• 5.2 性能分析与优化
• 6. 算法优化：单遍地形生成
• 6.1 算法原理与流程
• 6.2 性能对比：单遍 vs. 多遍
• 7. 总结与展示

1. 引言与背景

本系列幻灯片旨在介绍如何编写CUDA程序。演讲从CUDA平台的历史发展和演讲者的个人经历入手，引出程序员与GPU硬件交互的核心问题。

1.1 CUDA 发展时间线

CUDA平台自2007年发布以来，经历了长期的发展和迭代。下图展示了从CUDA 1.0到CUDA 12.0的主要版本发布时间点。

1.2 个人背景与历史回顾

演讲者Stephen Jones回顾了他在CUDA 2.0（2008年8月）时期开始接触和使用CUDA的经历。当时，进行CUDA编程所使用的典型硬件是一台桌面电脑和一块早期的NVIDIA Tesla GPU。

1.3 本次演讲的核心问题

本次演讲的核心旨在解答一个根本性问题：程序员（软件）与GPU（硬件）之间是如何连接和交互的？这中间的“桥梁”——即CUDA程序，其工作原理是什么？

1.4 不同抽象层级下的编程范式对比

通过一个直观的例子，对比了两种不同抽象层级的编程方式，以引出后续关于并行计算生态系统多样性的讨论。

• 底层编程（以Visual C++为例）: 左侧展示了一段复杂的Visual C++代码片段，用于处理按钮点击事件。代码涉及底层的系统调用、内存分配、Windows句柄操作、以及与硬件设备直接通信的细节。这代表了一种需要开发者处理大量实现细节的低抽象层级编程模式。

• 高层编程（以HTML为例）: 右侧展示了一行简单的HTML代码 (<button>)，其功能是创建一个可交互的按钮。这代表了一种高抽象层级的声明式编程模式，开发者只需描述期望的结果，而无需关心其底层的实现细节。

通过对比，该页面强调了不同抽象层次在编程复杂性上的巨大差异，为介绍NVIDIA计算平台提供的多层次解决方案奠定了基础。

1.5 NVIDIA 计算生态系统

这张幻灯片详细阐述了NVIDIA计算生态系统的分层架构，展示了开发者在不同抽象级别上拥有的丰富选择，以满足多样化的应用需求。

该生态系统从上至下分为四个主要层次：

1. 开发者与应用生态系统 (Developer & Application Ecosystem):

   • 框架 (Frameworks): 提供最高级别的抽象，适用于特定领域（如AI/ML）的快速开发。例如：PyTorch, TensorFlow, Jax, Modulus, Triton。
   • 软件开发工具包 (SDKs): 为特定行业应用提供预构建的工具和库。例如：医疗设备, 能源, 自动驾驶汽车。

2. NVIDIA加速库 (NVIDIA Accelerated Libraries):

   • 主机库 (Host Libraries): 预编译的、高性能的库，可通过API调用。例如：cuBLAS, cuFFT, cuSOLVER, cuTENSOR, cuRAND。
   • 源码库 (Source Libraries): 提供CUDA C++模板的开源库，允许更深度的定制和集成。例如：libcu++, CUB, Thrust, CUTLASS。

3. 并行语言 (Parallel Languages):

   • GPU原生 (GPU-Native): 为GPU编程设计的语言，提供对硬件最直接的控制。例如：CUDA C++, OpenCL, CUDA FORTRAN。
   • 标准语言 (Standard Languages): 通过指令或扩展支持并行计算的标准编程语言。例如：OpenACC, OpenMP, C++ 标准并行。
   • 第三方 (3rd Party): 由社区或第三方开发的支持GPU加速的语言或接口。例如：Numba (Python), JuliaGPU, MATLAB。

4. 编译器目标 (Compiler Targets):

   • NVVM / LLVM IR: NVIDIA和LLVM的中间表示，是高级语言到底层指令的桥梁。
   • PTX (Parallel Thread Execution) Assembly ISA: 并行线程执行指令集架构，确保了跨代GPU的二进制兼容性。

2. CUDA 编程基础

在介绍了不同抽象层级的工具后，演讲内容转向一个具体的CUDA编程实例。第一步是设计并行算法，思考如何将问题分解，以便能够利用GPU的大量并行核心进行计算。

2.1 环境搭建与示例获取

为了开始 CUDA 编程，首先需要准备好开发环境并获取官方示例代码。

1. 下载并安装 CUDA Toolkit：通过搜索引擎访问 NVIDIA Developer 网站，下载最新版本的 CUDA Toolkit。安装程序会引导用户完成系统检查、接受许可协议，并提供快速（推荐）或自定义安装选项。
   

   

2. 获取 CUDA 示例代码：NVIDIA 官方在 GitHub 上提供了一个包含大量 CUDA 代码示例的仓库。这些示例覆盖了从基础的 GPU 计算概念到高级应用，通过克隆此仓库可以获取这些示例以供学习和参考。
   

   

2.2 CUDA 程序结构与执行模型

程序基本结构

CUDA C++ 允许在单个源文件（通常以 .cu 为后缀）中混合编写主机（CPU）代码和设备（GPU）代码，体现了其异构编程的特性。

• GPU 部分：在 GPU 上执行的代码，通常被称为核函数（Kernel）。
• CPU 部分：控制程序流程、调用 GPU 核函数的主机代码。

一个简单的 CUDA 程序 (hello.cu) 展示了这种结构。代码通过 nvcc（NVIDIA C++ 编译器）编译，然后执行。

异构编程：调用核函数

在 CUDA 中，从 CPU 调用 GPU 函数有特定的语法规则：

1. 核函数定义：GPU 核函数使用 __global__ 修饰符进行声明。这个修饰符表明该函数是从 CPU 调用，但在 GPU 上执行。
   

   

2. 核函数调用：从 CPU 调用 __global__ 函数时，必须使用 <<<...>>> 执行配置语法。这个语法用于指定核函数执行时的线程网格（Grid）和线程块（Block）维度。
   

   

CPU 与 GPU 执行模型对比

GPU 和 CPU 的执行模型有显著差异，理解这些差异对于高效的 GPU 编程至关重要。

• C++ on CPU：
  • 核心数量较少，每个核心运行 1-2 个线程。
  • 线程运行独立的任务。
  • 采用抢占式多任务（Pre-emptive multi-tasking）。
• CUDA C++ on GPU：
  • 核心（流多处理器，SM）数量非常多，每个核心可运行大量线程（例如 2048 个）。
  • 大量线程协作完成同一个任务。
  • 采用非抢占式执行（No pre-emption），线程会一直运行直到完成。

2.3 CUDA 编程模型：网格 → 线程块 → 线程

CUDA 采用一种层级化的模型来组织和管理在 GPU 上执行的大量线程，以便于将问题分解并映射到硬件上。

1. 工作网格 (Grid of work)：首先，整个计算任务被定义为一个网格（Grid）。例如，在图像处理中，整个图像可以被视为一个工作网格。
2. 线程块 (Blocks)：工作网格被划分为多个大小相同的线程块（Block）。每个线程块包含一组线程，它们可以相互协作（例如，通过共享内存）。
3. 线程 (Threads)：每个线程块由多个线程组成。每个线程执行核函数的代码，并通过内置变量（如 threadIdx, blockIdx）来获取其在块和网格中的唯一索引，从而处理数据的不同部分。

这种分层结构允许将大规模并行问题分解为更小、更易于管理的部分，并将其映射到 GPU 的硬件架构上。

3. 案例研究：GPU 驱动的递归地形生成

3.1 核心算法：二维递归细分

递归细分是一种在二维空间中生成更复杂几何形状的技术。该过程从一个简单的基元（如一个正方形网格）开始，通过递归地将其划分为更小的部分来增加细节。

1. 初始步骤：从一个基础的四边形网格开始。在每条边的中点和网格的中心点插入新的顶点，从而将一个四边形划分为四个更小的四边形。
   

   

2. 递归过程：对新生成的子四边形重复此细分过程。
   

   

3. 顶点共享：在细分过程中，新生成的顶点可能会被多个不同层级或邻近的细分单元所共享。处理这些共享顶点是确保网格连续性和避免裂缝的关键。
   

   

4. 应用实例：通过对一个平面网格进行多层递归细分，并对新生成的顶点在第三维度（例如高度）上进行位移，可以生成复杂的表面，如程序化的地形或动态的水面。
   

   

3.2 将问题映射到 GPU

递归细分的并行化

递归细分算法可以通过在 GPU 上迭代执行来并行化。每一次迭代都会增加网格的精细程度。这个过程可以映射为一系列 CUDA 核函数 (kernel) 的启动。

• 第 1 步: 从一个包含 4 个点的方块开始，由 1 个线程块处理。
• 第 2 步: 细分为 4 个方块，共 9 个点，由 4 个线程块处理。
• 第 3 步: 细分为 16 个方块，共 25 个点，由 16 个线程块处理。
• 第 4 步: 细分为 64 个方块，共 81 个点，由 64 个线程块处理。

随着细分的进行，线程块的数量呈指数级增长，从而利用 GPU 的大规模并行能力。

确定工作划分

在将一个大问题映射到 GPU 时，合理地划分工作至关重要。以一个 1024x1024 点的计算任务为例，总共有 1,048,576 (约 1M) 个点需要处理。

假设为每个点分配一个线程，我们需要选择合适的网格维度和线程块维度。下表分析了几种不同的划分方案：

经验法则 (Rule of thumb):
- 在不确定时，优先选择大小为 256 的线程块。
- 推论 (Corollary): 尽量避免使用小于 128 的线程块大小。

根据上述分析，最终选择将 1024x1024 的点阵划分为 64x64 的网格，每个线程块处理 16x16 的点（即 256 个线程/块）。

3.3 快速实现：修改 OceanFFT 示例

为了快速实现地形生成，可以复用一个现有的 CUDA 示例项目，如 oceanFFT。

oceanFFT 原始工作流使用 FFT 计算海面并渲染。我们的修改策略是替换其核心计算部分，同时保留渲染和交互逻辑。

我们将原始工作流中的数值计算步骤替换为一个新的地形生成核函数，保持 UI 和渲染相关的代码不变。

为了验证我们能够控制输出，可以编写一个简单的测试核函数，它在网格上绘制对角线。

// 一个快速验证用的测试核函数，用于在网格上绘制对角线
static __global__ void updateHeightmapKernel(float *heightMap,
                                             unsigned int width, int timestep) {
    unsigned int x = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int y = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned int i = y * width + x;

    heightMap[i] = (float)(((x + y + timestep) * 2) % 256) / 256.0f;
}

每个线程根据其唯一的 blockIdx 和 threadIdx 计算出它负责处理的像素坐标 (x, y)。

这个简单的核函数成功地在三维网格上生成了预期的对角线图案，证明了我们已成功接管了 oceanFFT 的数据生成流程。

3.4 挑战：并行伪随机数生成

问题描述

在并行环境中实现递归细分算法时，会遇到伪随机数生成的问题：

1. 随机数需求: 算法的每一步都需要使用随机数来确定新的点的高度。
2. 序列中断: 并行计算破坏了传统伪随机数生成器 (PRNG) 的顺序依赖性。
3. 边界一致性: 位于不同线程块边界上的共享点，可能会被不同的线程块计算两次。为了保证网格的无缝连接，这两次计算必须得到完全相同的结果，这意味着它们必须使用相同的随机数。

因此，需要一种能够基于点的坐标确定性地生成随机数的方法，而不是依赖于计算的顺序。

解决方案：使用 cuRAND

为了高效地在GPU上生成随机数，无需从零开始构建解决方案。NVIDIA CUDA工具包提供了一个专门的随机数生成库——cuRAND。cuRAND是NVIDIA加速库生态系统的一部分，提供了多种伪随机数和拟随机数生成器，能够高效地利用GPU的并行计算能力。

3.5 实现递归地形生成

在验证了可以控制高度图生成后，我们便可以实现递归地形生成算法。该算法通过重复启动 CUDA 网格来逐步细化地形，直至达到所需的分辨率。每一次启动的网格尺寸都比上一次更大，例如从 1x1 网格开始，依次扩展到 2x2、4x4、8x8，最终达到 64x64 的网格，从而生成高分辨率的地形。

利用并行生成的随机数，可以采用递归算法来创建复杂的地形。下图展示了通过启动一个64x64的线程网格（grid）所生成的地形。

4. 高级技术与优化

4.1 地形动画与变形

地图间的平滑变形

本节介绍如何实现地形图之间的平滑过渡。具体步骤如下：
1. 生成初始地图：首先生成一个基础地形图。
2. 生成目标地图：再生成一个形态不同的目标地形图。
3. 逐像素插值：通过逐步调整初始地图中每个像素的高度值，使其趋向于目标地图中对应像素的高度值，从而实现两个地形之间的平滑变形。

下图展示了单个像素值随时间从初始状态向目标状态进行插值的过程。

最终，通过这一系列的插值操作，可以实现两个完全不同的地形之间的无缝过渡。

无缝滚动动画

为了实现平滑连续的地形滚动效果，该技术通过组合两张地形图（“当前地图”和“下一张地图”）来生成最终的显示画面。通过一个“位移(shift)”变量，动态地从这两张地图中采样像素，从而创造出无缝滚动的动画效果。

这个过程在像素级别通过以下代码逻辑实现：

// Given (x,y), copy a pixel from one or other map
float pixel;
if (y < shift) {
    unsigned int row = (mapHeight - shift) + y;
    pixel = *(map1 + x + (row * mapWidth));
} else {
    unsigned int row = y - shift;
    pixel = *(map0 + x + (row * mapWidth));
}

随着位移值的增加，画面顶部会显示越来越多“下一张地图”的内容，而底部“当前地图”的内容则逐渐移出画面，实现了地形向上滚动的视觉效果。

4.2 视觉效果增强与无缝景观

视觉效果增强

为了让生成的地形更具真实感和美观度，可以对片元着色器（Fragment Shader）进行调整，并应用海平面效果。该过程包括以下几个步骤：
1. 生成地形图：首先生成原始的高度图地形。
2. 按高度着色：根据地形的高度值为其应用不同的颜色。
3. 应用海平面：设定一个高度阈值作为海平面。所有低于该高度的区域都被渲染为水体。

无缝景观生成

为了创建可以无限滚动的广阔世界，需要生成可无缝拼接的地形瓦片（Tile）。实现无缝拼接的关键在于处理相邻瓦片之间的边界。在生成一个新的地形瓦片时，其边缘的初始值直接复制自前一个相邻瓦片对应边缘的值，确保两个瓦片在拼接处的高度和形态完全匹配。

4.3 性能优化：将数据保留在 GPU

为了实现高性能的地形渲染，关键在于最小化CPU和GPU之间的数据传输。CPU与GPU之间的PCIe总线带宽远低于GPU内部显存的带宽。因此，最优策略是尽可能将所有计算和数据都保留在GPU上。

其工作流程如下：
1. 直接生成地形高度图：在GPU上直接运行CUDA核函数生成地形的高度图数据。
2. 映射OpenGL显示对象：将生成的高度图数据直接映射到OpenGL的顶点缓冲对象（VBO），无需将数据传回CPU。
3. 渲染三角形网格：GPU使用这些数据直接渲染地形。
4. 视点控制与用户输入：CPU仅负责处理用户输入和视点控制等逻辑。

4.4 动态高度与智能内存管理

动态高度调整

该技术应用于模拟低空飞行场景，需要根据飞行视点正下方的地形动态调整高度。

为了实现这一功能，系统需要在CPU和GPU之间协同工作。地形的高度剖面数据存储在GPU内存中，而视点控制逻辑在CPU上运行。

视点控制代码在CPU上运行，需要获取GPU中的地形高度数据来计算摄像机的新位置。

// 视点控制代码在CPU上运行
if ( autoFollow && ( bottomRow != nullptr ) ) {
    float y = bottomRow[meshWidth / 2];
    targetTransY = anim_base_y - ( y / 5.0f );
}

智能内存分配

一个关键挑战是CPU无法直接访问存储在GPU显存中的数据。解决方案是必须将所需的数据从GPU显存复制到CPU内存中。

内存分配策略:

1. 分配GPU内存：使用 cudaMalloc 在GPU上分配内存，用于存储大规模数据。
   c++ // 使用cudaMalloc分配GPU内存 cudaMalloc(&heightmap, numPoints * sizeof(float));

2. 分配锁页主机内存 (Pinned System Memory)：使用 cudaMallocHost 在CPU上分配特殊的“锁页”内存。这种内存对GPU是可见的，并且优化了与GPU之间的数据传输，可以实现高速内存拷贝，并且GPU可直接访问。

c++ // 使用cudaMallocHost分配CPU锁页内存 cudaMallocHost(&bottomRow, meshWidth * sizeof(float));

5. 开发与调试

5.1 调试 GPU 程序

NVIDIA CUDA Toolkit 提供了强大的调试工具 Nsight，它能够与主流的IDE（如Visual Studio）集成，为GPU程序提供完善的调试支持。

Nsight调试器的主要功能包括：
* 异构调试：在同一个IDE窗口内同时调试CPU和GPU代码。
* 设置断点：可以在CPU代码和GPU内核代码中设置断点。
* 监视变量：实时查看CPU和GPU上的变量值。
* 单步执行：逐行执行CPU和GPU代码。

5.2 性能分析与优化

为了对GPU算法进行深度性能调优，可以使用NVIDIA Nsight Compute工具。在使用前，需要在Windows系统上通过NVIDIA控制面板启用GPU性能计数器，以允许分析工具访问底层的硬件性能数据，然后启动NVIDIA Nsight Compute。

6. 算法优化：单遍地形生成

6.1 算法原理与流程

本节介绍一种高效的、基于GPU的单遍地形生成算法。该算法的核心思想是为网格中的每一个点启动一个线程，每个线程独立、无依赖地计算出其对应坐标点的高度值，从而避免了传统多遍算法中存在的读写依赖和多次内核启动开销。

算法对每个线程 (x,y) 执行以下步骤：
1. 生成初始角点：为整个地形区域生成初始的四个角点的高度值。
2. 计算中点：在角点之间计算中点的高度值。
3. 递归进入象限：如果当前计算区域的中心不是目标点 (x,y)，则算法会判断 (x,y) 位于当前区域的哪个象限，并以该象限的角点和新计算出的中点作为新的边界，递归地进行下一轮计算，直至精确定位到目标点 (x,y)。

为了确保所有线程在计算任何一个中间点时都能得到完全相同的值，算法使用了一个伪随机数生成器。这消除了线程间对共享内存的读写需求。

当递归过程最终收敛到目标点 (x,y) 时，计算出其最终的高度值，并将其写入到输出内存中的相应位置。整个过程的成本极低：
* 1 次网格启动：只需一次内核函数调用。
* 0 次内存读取：无需从全局内存中读取任何数据。
* 每个像素 1 次写入：每个线程只在最后写入一次其最终计算结果。
* 每个网格 0.4 毫秒：生成整个地形网格的耗时极短。

6.2 性能对比：单遍 vs. 多遍

该幻灯片将优化的单遍算法与传统的多遍算法进行了性能比较。

• 优化的单遍算法:

  • 耗时: 每个网格 0.4ms
  • 流程: 启动一个大网格，每个线程独立通过递归计算出最终点的值并写入。
  • 优势: 无需从内存中读取中间点，避免了多遍渲染中的读写延迟和同步开销。

• 多遍算法:

  • 耗时: 每个网格 0.8ms
  • 流程: 在循环中多次启动尺寸递增的网格，每一遍都从内存中读取上一遍的结果，计算新点并写出。
  • 劣势: 需要多次内核启动，并且在每一遍中都涉及大量的内存读写操作。

结论：优化的单遍算法性能是传统多遍算法的 2倍。

7. 总结与展示

通过利用 CUDA 并行计算模型，该单遍地形生成算法实现了高效、无依赖的计算，显著优于传统方法。

最终生成的地形效果图如下所示，展示了该算法能够创建出细节丰富且视觉效果出色的地形。