使用 oneAPI 让您的 C++ 应用程序支持 GPU
5.00/5 (4投票s)
使用 Intel 的 oneAPI 加速支持 GPU 的 C++ 文章
引言
在本文中，我们将通过将应用程序移植到使用 Intel 的 oneAPI 工具包和编译器来加速现有的 C++ 程序。我们将移植的应用程序是一个计算图像 离散余弦变换 的 C++ 程序。
此操作在图像和视频压缩编解码器（例如 H.264、WebM、JPEG、WebP 等）中非常常见，并且也非常适合并行化。首先，我们将介绍 C++ 应用程序并简要描述其工作原理。有关 DCT 底层数学原理的更全面描述，您需要查阅 维基页面，该页面不仅提供了易于理解的计算描述，还提供了计算示例。然后，我们将演示如何利用 Intel 的 oneAPI 库来并行化应用程序。
Intel 的 oneAPI 是一个开放的、基于标准的编程模型，使 C++ 开发人员能够以统一的方式定位不同的加速器架构。这与 Vulkan、CUDA 或 D3D12 等针对 GPU 的专用 API，或者针对 FPGA 的 Verilog 和 VHDL 不同。使用 oneAPI 编程时，相同的代码可以同时用于定位 CPU、GPU 和 FPGA。您只需要一个正常工作的 oneAPI 工具包 安装即可，该工具包随附 DPC++（Intel 对 Clang 编译器的扩展），用于将 C++ 代码编译为定位各种支持的加速器类型。
使用 oneAPI 的另一个好处是能够将应用程序部署到 Intel 的 DevCloud，这是一个提供计算集群访问权限的沙盒环境，这些集群配备了强大的 CPU 和 FPGA。熟悉 FPGA 工具链管理复杂性的从业者可能会乐于在一个近乎即插即用的环境中部署软件，定位 Arria 10 等 FPGA，而无需立即投资硬件。在本文中，我们将从头开始编写一个使用 oneAPI 和 DPC++ 的小型程序，然后将其部署并在 DevCloud 的硬件上运行。
为什么要移植到 GPU？
将 C++ 应用程序移植到 GPU 似乎显而易见是绝对有利的。然而，重要的是要考虑 CPU 和 GPU 性能的差异，并分析涉及的权衡。CPU 是拓扑结构简单的处理器。CPU 由多个核心组成（这些核心可能进一步细分为逻辑核心，例如 Intel 的超线程）。每个逻辑核心都有自己的指令缓存、数据缓存和命令队列。当一个线程分配给一个核心时，它对该核心的资源拥有独占访问权，并可以从内存中获取指令进行评估。
相比之下，GPU 的架构本质上是分层的，并且旨在最大化吞吐量。为了实现这一点，而不是让每个核心维护自己的缓存和命令队列，单个处理单元被分组到子组（在 Nvidia 硬件上称为 warp）中。整个子组的指令是单指令发布的，也就是说，子组内的所有线程在给定时间都操作相同的指令。此外，虽然每个线程都有自己的一组向量寄存器，但整个子组可以访问共享的快速设备内存。这种硅结构为何能带来更高的吞吐量，需要更少的资源和指令来驱动更多工作，这一点应该很容易理解。
能够清晰映射到 GPU 的算法是那些能够高度并行化的算法。特别是，如果一个任务可以分解为子任务，并且这些子任务可以相对独立地工作，那么 GPU 就是加速工作的首选。离散余弦变换就是一种易于 GPU 加速的算法，我们将在本文中进行探讨。
C++ 离散余弦变换
我们不深入探讨离散余弦变换（以下简称 DCT）的细节，只回顾一些关于 DCT 的事实：
-
与傅里叶变换一样，DCT 是从空间域到频率域的完全可逆变换。
-
与傅里叶变换不同，DCT 仅基于偶函数（余弦）。
-
DCT 常用于图像和视频压缩算法；应用 DCT 后，分配的位数会减少以存储高频分量，因为它们在感知上不太重要。
当 DCT 实际应用于图像或视频帧时，通常应用于光栅内的单个块（8x8 块是典型大小）。在此演示中，我们将执行以下操作：
-
加载图像并将其转换为灰度。
-
对于图像中的每个 8x8 块，将其映射到 8x8 频率幅度块（应用 DCT）。
-
立即反转 DCT 以恢复原始图像，并将其写入磁盘以供比较。
此应用程序不会“产生”任何结果，除了一个有望与原始图像完全相同的图像。但是，中间频率幅度值可以在典型的图像或视频编解码器中进行量化和压缩。我们在这里执行的编码和解码在很大程度上说明了该技术，解码输出是确定我们的代码是否按预期工作的有用视觉指标。虽然我们不会深入探讨 DCT 或其属性，但我们鼓励有兴趣的读者阅读 维基百科 上的精彩介绍。
闲话少说，我们开始编写代码吧！
图像加载和灰度转换
为了读取和写入图像，我们将分别使用单文件库 stb_image 和 stbimagewrite。将这些文件放在源代码树中后，创建一个 `main.cpp` 文件，并在顶部包含以下头文件：
#include <CL/sycl.hpp>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"
#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"
然后，我们可以像这样加载图像：
int channels;
int width;
int height;
uint8_t* image = stbi_load(path, &width, &height, &channels, 3);
Each individual RGB pixel can be mapped to greyscale with the following function:
float luminance(uint8_t r, uint8_t g, uint8_t b)
{
float r_lin = static_cast<float>(r) / 255;
float g_lin = static_cast<float>(g) / 255;
float b_lin = static_cast<float>(b) / 255;
// Perceptual luminance (CIE 1931)
return 0.2126f * r_lin + 0.7152 * g_lin + 0.0722 * b_lin;
}
如果您不熟悉色彩理论，只需记住绿色在感知上比红色和蓝色更亮，而蓝色是这三种颜色分量中感知上最暗的。
然后，我们可以简单地像这样将整个图像映射到灰度：
uint8_t* out = reinterpret_cast<uint8_t*>(std::malloc(width * height));
int8_t* grey = reinterpret_cast<int8_t*>(std::malloc(width * height));
// First, we convert the image to 8-bit greyscale, centered around 0
for (int y = 0; y != height; ++y)
{
for (int x = 0; x != width; ++x)
{
int offset = y * width + x;
uint8_t* pixel = &image[3 * offset];
float lum = luminance(pixel[0], pixel[1], pixel[2]);
// The DCT projects the signal to cosines (even functions) so our
// signal must be centered on zero
grey[offset] = static_cast<int>(lum * 255) - 128;
// Write out the direct greyscale conversion for visual comparison
out[offset] = static_cast<uint8_t>(lum * 255);
}
}
// Write out the greyscale image
stbi_write_png("grey.png", width, height, 1, out, width);
DCT 及其逆变换
接下来，我们将遍历图像中的每个 8x8 块，应用 DCT，反转 DCT，并将结果写入新图像。
// Some needed constants
constexpr static int block_size = 8;
constexpr static float inv_sqrt2 = 0.70710678118f;
constexpr static float pi = 3.14159265358979323846f;
// scale factor needed for the DCT to be orthonormal
// (1 / sqrt(2 * block_size))
constexpr static float scale = 0.25f;
int block_count = width / block_size;
for (int block_x = 0; block_x != block_count; ++block_x)
{
for (int block_y = 0; block_y != block_count; ++block_y)
{
// Results of the DCT will be stored here
float dct[block_size * block_size];
// Apply the discrete cosine transform to each block
for (int u = 0; u != block_size; ++u)
{
float a_u = u == 0 ? inv_sqrt2 : 1.f;
for (int v = 0; v != block_size; ++v)
{
// Compute a single amplitude per iteration
float& out = dct[v * block_size + u];
float a_v = v == 0 ? inv_sqrt2 : 1.f;
float a = scale * a_u * a_v;
// Initialize amplitude
out = 0.f;
for (int x = 0; x != block_size; ++x)
{
for (int y = 0; y != block_size; ++y)
{
int8_t value = grey[width * (block_y * block_size + y)
+ block_x * block_size + x];
float cosu = std::cosf((2 * x + 1) * u * pi / 2 / block_size);
float cosv = std::cosf((2 * y + 1) * v * pi / 2 / block_size);
out += cosu * cosv * value;
}
}
out *= a;
}
}
// Now, the dct array holds the 8x8 frequency amplitudes
// Next, invert the transform, writing the result out to the
// original image buffer
for (int x = 0; x != block_size; ++x)
{
for (int y = 0; y != block_size; ++y)
{
uint8_t& value = out[width * (block_y * block_size + y)
+ block_x * block_size + x];
// Accumulate the results in a floating point register
float acc = 0.f;
for (int u = 0; u != block_size; ++u)
{
float a_u = u == 0 ? inv_sqrt2 : 1.f;
for (int v = 0; v != block_size; ++v)
{
float a_v = v == 0 ? inv_sqrt2 : 1.f;
float a = scale * a_u * a_v;
float amplitude = dct[v * block_size + u];
float cos_u = std::cosf((2 * x + 1) * u * pi / 2 / block_size);
float cos_v = std::cosf((2 * y + 1) * v * pi / 2 / block_size);
acc += a * amplitude * cos_u * cos_v;
}
}
// Don't forget to reverse the subtraction where we
// centered the grey value on zero.
value = static_cast<uint8_t>(acc + 128);
}
}
}
}
stbi_write_png("grey_dct.png", width, height, 1, out, width);
这段代码片段乍一看可能有点长，但总体结构是我们有两组循环。最外层的双循环遍历每个 8x8 块。内层循环迭代单个 8x8 块的 64 个条目。
如果一切顺利，这段代码应该会生成一对如下所示的图像（它们应该看起来完全相同）：
这里的第二张图像经过了 DCT 和反 DCT，它们的相似度应该非常显著（您遇到的任何错误可能都与浮点误差有关）。
使用 Intel 的 DPC++ 在 GPU 上运行
这些循环看起来应该非常适合并行化。毕竟，每个块都可以独立评估，而无需了解其他块正在发生什么。此外，在一个块内，每个频率幅度都可以独立于其他幅度进行评估。一旦所有频率幅度都可用，反变换就可以再次并行进行。
我们的策略是，将每个 8x8 块作为一个独立的工作组进行分派。在工作组内，所有频率都可以并行评估。然后，我们希望工作组中的每个线程都阻塞，直到前向 DCT 完成。在此屏障之后，所有线程都将能够恢复并行执行，这次是执行反变换并写出所需的数据。
初始化 SYCL 设备和队列
为了提交硬件加速命令，我们需要构造一个 `cl::sycl::queue`。队列用于将命令和内存编排出设备。
// Use this device selector to select the Gen9 Intel GPU
class gen9_gpu_selector : public cl::sycl::device_selector
{
public:
gen9_gpu_selector()
: cl::sycl::device_selector()
{}
int operator()(cl::sycl::device const& device) const override
{
if (device.get_info<cl::sycl::info::device::name>()
== "Intel(R) Gen9 HD Graphics NEO")
{
return 100;
}
return 0;
}
};
cl::sycl::queue queue{gen9_gpu_selector{}};
在这里，我们创建了一个自定义选择器，以优先选择最近 Intel CPU 上的 Gen9 集成 GPU。您还可以选择内置选择器，例如 `cl::sycl::gpu_selector`、`cl::sycl::fpga_selector`、`cl::sycl::host_selector`，或者仅仅是 `cl::sycl::default_selector`。
灰度转换
要执行灰度转换，我们首先需要分配可由设备写入但可由主机读取的内存。为此，我们将使用 DPC++ 的统一内存抽象。请注意，此抽象不属于 SYCL 标准，而是 Intel 提供的扩展。我们还需要一个缓冲区来存放中间的 0 中心灰度数据的设备内存。
uint8_t* out = reinterpret_cast<uint8_t*>(
cl::sycl::malloc_shared(width * height, queue));
cl::sycl::buffer<int8_t, 1> greyscale_buffer{width * height};
灰度转换例程本身与我们的 C++ 代码非常相似。但是，代码不是直接在主机上执行，而是提交给之前创建的设备队列。
queue.submit([&greyscale_buffer, &image_buffer, out, width, height](
cl::sycl::handler& h) {
// Request write access to the greyscale buffer
auto data = greyscale_buffer.get_access<cl::sycl::access::mode::discard_write>(h);
// Request read access to our image data
auto image = image_buffer.get_access<cl::sycl::access::mode::read>(h);
// Parallel-for over each width*height pixel
h.parallel_for(cl::sycl::range<1>(width * height),
[image, data, out](cl::sycl::id<1> idx) {
int offset = 3 * idx[0];
float lum = luminance(image[offset],
image[offset + 1],
image[offset + 2]);
// The DCT projects the signal to cosines (even
// functions) so our signal must be centered on
// zero
data[idx[0]] = static_cast<int>(lum * 255) - 128;
// Write out the direct greyscale conversion for
// visual comparison
out[idx[0]] = static_cast<uint8_t>(lum * 255);
});
});
SYCL 缓冲区的方法 `get_access<M>` 允许我们声明提交给队列的代码将以特定方式访问内存。SYCL 运行时会相应地对队列提交以及任何必要的内存同步进行排序。此外，请注意，虽然我们使用纯 C++ 编写了亮度函数，但编译器能够将其编译为适合我们定位的设备的执行代码。
DPC++ 代码和 C++ 代码的另一个区别是循环以范围抽象的形式表示。例如，`cl::sycl::range<1>(20)` 表示一个包含 20 个项目的 1 维范围。在内核函数本身中，会传递一个 `cl::sycl::id<1>` 索引对象，以将分配给线程的范围中的项告知线程。我们将在下一节看到另一种重要的范围类型。
前向和反向 DCT 工作组提交
如前所述，我们希望并行提交 8x8 块，并在每个块内进行另一层并行化。在 SYCL 中，这通过一种不同类型的范围来表示，即 `cl::sycl::nd_range`。我们还需要一种方法来请求瞬态但跨工作组共享的内存。让我们看看我们的前向和反向 DCT 提交的结构：
// Now, we do another parallel-for, but use hierarchical parallelism to
// dispatch 8x8 blocks in parallel
int block_count = width / block_size;
queue.submit([&greyscale_buffer, out, width, height, block_count](
cl::sycl::handler& h) {
auto data = greyscale_buffer.get_access<cl::sycl::access::mode::read>(h);
// Local memory storage
cl::sycl::range<1> local_memory_range(block_size * block_size);
cl::sycl::accessor<float,
1,
cl::sycl::access::mode::read_write,
cl::sycl::access::target::local>
dct(local_memory_range, h);
cl::sycl::range<2> workgroup_count(block_count, block_count);
// An individual workgroup is 8x8
cl::sycl::range<2> workgroup_size(block_size, block_size);
// The global range spans with width and height of our image
auto global = workgroup_count * workgroup_size;
h.parallel_for(
cl::sycl::nd_range<2>(global, workgroup_size),
[data, out, width, dct](cl::sycl::nd_item<2> item) {
// Ranges from 0 to the block_count - 1
int block_x = item.get_group(0);
// Ranges from 0 to the block_count - 1
int block_y = item.get_group(1);
// Ranges from 0-7
int u = item.get_local_id(0);
// Ranges from 0-7
int v = item.get_local_id(1);
// Perform forward and inverse DCT here
// ...
});
});
在这里，我们使用 `nd_range<2>` 对象来表示一个两级分层范围。顶层分层跨越整个图像，但我们提供第二个范围（此处称为 `workgroup_size`）来指定内部分层的维度。此外，我们创建一个本地内存访问器，以分配 8x8 浮点数的空间，这些浮点数需要存储临时 DCT 数据。传递给 `parallel_for` 的 lambda 接收的 `cl::sycl::nd_item` 对象具有许多方法，用于检索调用在工作组内的位置以及工作组在全局调度上下文中的位置。可以通过工作组大小和计数推导出组和局部 ID 对应的范围，如上面的注释所示。
有了这个，并行 for 循环的内部编写起来就相对简单了：
// Within the parallel-for body:
// Scale factors
float a_u = u == 0 ? inv_sqrt2 : 1.f;
float a_v = v == 0 ? inv_sqrt2 : 1.f;
float a = scale * a_u * a_v;
float amp = 0.f;
for (int x = 0; x != block_size; ++x)
{
for (int y = 0; y != block_size; ++y)
{
int8_t value = data[width * (block_y * block_size + y)
+ block_x * block_size + x];
float cosu = cl::sycl::cos((2 * x + 1) * u * pi / 2 / block_size);
float cosv = cl::sycl::cos((2 * y + 1) * v * pi / 2 / block_size);
amp += cosu * cosv * value;
}
}
amp *= a;
// Write amplitude to local storage buffer
dct[v * block_size + u] = amp;
// Invoke a local memory barrier so all writes to dct are
// visible across the workgroup
item.mem_fence<cl::sycl::access::mode::write>(
cl::sycl::access::fence_space::local_space);
item.barrier(cl::sycl::access::fence_space::local_space);
// Next, invert the DCT
int x = item.get_local_id(0);
int y = item.get_local_id(1);
uint8_t& value = out[width * (block_y * block_size + y)
+ block_x * block_size + x];
// Accumulate the results in a floating point register
float acc = 0.f;
for (u = 0; u != 8; ++u)
{
a_u = u == 0 ? inv_sqrt2 : 1.f;
for (v = 0; v != 8; ++v)
{
a_v = v == 0 ? inv_sqrt2 : 1.f;
a = scale * a_u * a_v;
amp = dct[v * block_size + u];
float cos_u = cl::sycl::cos((2 * x + 1) * u * pi / 2 / block_size);
float cos_v = cl::sycl::cos((2 * y + 1) * v * pi / 2 / block_size);
acc += a * amp * cos_u * cos_v;
}
}
value = static_cast<uint8_t>(acc + 128);
这里有几个重要的注意事项。虽然大部分代码与我们的 C++ 程序相似，但我们也需要 `mem_fence` 和屏障。如上所示的内存围栏确保所有对本地共享内存的写入对工作组内的所有调用都是可见的，而屏障确保所有线程在每个其他线程遇到内存围栏之前都等待。
如果没有这种同步，内核后半部分中的调用可能会在写入发生或可见之前尝试读取内存。另一个注意事项是我们之前对 `std::cos` 的调用现在被替换为对 `cl::sycl::cos` 的调用。内核在程序执行期间进行即时编译，并且诸如此处所示的超越函数之类的数学函数必须来自 SYCL 运行时。
部署到 DevCloud
现在我们有了一个可运行的程序，我们可以将源代码和任何构建脚本部署到 Intel oneAPI DevCloud。注册该程序后，您应该会收到一封电子邮件，其中包含有关如何获取 SSH 凭据以登录 DevCloud 平台的说明。按照说明操作后，您应该可以使用以下命令上传您的源文件和测试图像：
scp main.cpp devcloud:~
scp stb_image.h devcloud:~
scp stb_image_write.h devcloud:~
scp peppers.png devcloud:~
这会将您的源文件上传到您分配的 DevCloud 用户帐户的主目录。然后，您可以登录以部署您的程序：
ssh devcloud
创建一个脚本来编译您的程序并像这样运行它：
#!/usr/bin/env bash
# run.sh
# Ensure the correct environment variables are set
source /opt/intel/inteloneapi/setvars.sh
# Compile our program
dpcpp main.cpp -o dct -std=c++17 -fsycl -lOpenCL
# Invoke the program with our test image
./dct peppers.png
然后，我们可以使用 `qsub`（队列提交）命令在 DevCloud 上可用的各种计算节点上调用我们的脚本。要查看可用主机的列表，请运行 `pbsnodes` 命令。这将按 ID 列出节点，并提供有关运行的处理器类型以及可用的加速器的附加信息。
例如，要将作业提交到具有 GPU 的主机，我们可以运行以下命令：
qsub -l nodes=1:gpu:ppn=2 -d . run.sh
简而言之，这些选项表示我们希望在具有 GPU 的单个节点上运行脚本，我们希望完全占用该节点（`ppn=2` 选项），我们希望工作目录是当前目录，我们希望节点运行 `run.sh` 脚本。
要查看作业的状态，您可以运行 `qstat`，它将产生类似以下的输出：
Job ID Name User Time Use S Queue
------------------------- ---------------- --------------- -------- - -----
681510.v-qsvr-1 run.sh u47956 0 R batch
作业 ID 可用作 `qdel` 命令的参数，以取消待处理的作业。
作业完成后，您将在当前目录中找到 `run.sh.o681510` 和 `run.sh.e068150` 等文件，分别对应脚本的 `stdout` 和 `stderr` 输出。如果我们的程序成功运行，您应该还会看到 `grey.png` 和 `grey_dct.png` 图像，您可以进行比较以验证其正确性。使用 `exit` 命令注销，并使用 `scp` 将图像传输回您的主机：
scp devcloud:~/grey.png .
scp devcloud:~/grey_dct.png .
结论
在本文中，我们将一个 C++ 应用程序移植到了 SYCL 运行时，并使用了 Intel DPC++ 编译器提供的附加扩展。该应用程序演示了统一编程模型如何定位不同的架构，以及 SYCL 运行时为协调内存访问和使用隐式依赖图编写并行代码提供的抽象。此外，我们展示了如何以与 GPU 架构相对应的分层方式表达代码，以提高性能。
最后，我们展示了如何将代码部署并针对 Intel DevCloud 提供的不同硬件配置文件进行测试。要详细了解 SYCL 运行时、Intel DevCloud 或 Intel DPC++ 编译器，我们鼓励您在 Intel DevZone 此处 阅读更多内容。
历史
- 2020 年 11 月 9 日：初始版本