如何在Julia编程中实现GPU加速

选自nextjournal

作者：Simon Danisch

机器之心编译

参与：高璇、刘晓坤

GPU 的并行线程可以大幅提升速度，但也使得代码编写变得更复杂。而 Julia 作为一种高级脚本语言，允许在其中编写内核和环境代码，并可在大多数 GPU 硬件上运行。本文旨在介绍 GPU 的工作原理，详细说明当前的 Julia GPU 环境，以及展示如何轻松运行简单 GPU 程序。

为了简化操作，可以在 nextjournal 上注册账户，点击「edit」即可直接运行文章中的简单代码了。

注册地址：https://nextjournal.com/signup

首先，什么是 GPU？

GPU 是一种大型并行处理器，有几千个并行处理单元。例如，本文使用的 Tesla k80，能提供 4992 个并行 CUDA 核。GPU 在频率、延迟和硬件性能方面与 CPU 有很大的不同，但实际上 Tesla k80 有点类似于具有 4992 核的慢速 CPU。

能够启动的并行线程可以大幅提升速度，但也令使用 GPU 变得更困难。当使用这种未加处理的能量时，会出现以下缺点：

GPU 是一种有专属内存空间和不同架构的独立硬件。因此，从 RAM 到 GPU 内存（VRAM，显存）的传输时间很长。甚至在 GPU 上启动内核（调用调度函数）也会带来很大的延迟，对于 GPU 而言是 10us 左右，而对于 CPU 只有几纳秒。

在没有高级封装的情况下，建立内核会变得复杂。

低精度是默认值，高精度的计算可以很容易地消除所有性能增益。

GPU 函数（内核）本质上是并行的，所以编写 GPU 内核不比编写并行 CPU 代码容易，而且硬件上的差异增加了一定的复杂性。

与上述情况相关的很多算法都不能很好地迁移到 GPU 上。想要了解更多的细节，请看这篇博文：https://streamhpc.com/blog/2013-06-03/the application-areas-opencl- cuda-can- used/。

内核通常是用 C/ C++语言编写的，但这并不是写算法的最好语言。

CUDA 和 OpenCL 之间有差异，OpenCL 是编写底层 GPU 代码的主要框架。虽然 CUDA 只支持英伟达硬件，OpenCL 支持所有硬件，但并不精细。要看个人需求进行选择。

而 Julia 作为一种高级脚本语言，允许在其中编写内核和环境代码，同时可在大多数 GPU 硬件上运行！

GPUArrays

大多数高度并行的算法都需要同时处理大量数据，以克服所有的多线程和延迟损耗。因此，大多数算法都需要数组来管理所有数据，这就需要一个好的 GPU 数组库作为关键的基础。

GPUArrays.jl 是 Julia 为此提供的基础。它实现了一个专门用于高度并行硬件的抽象数组。它包含了设置 GPU、启动 Julia GPU 函数、提供一些基本数组算法等所有必要功能。

抽象意味着它需要以 CuArrays 和 CLArrays 的形式实现。由于继承了 GPUArrays 的所有功能，它们提供的接口完全相同。唯一的区别出现在分配数组时，这会强制用户决定这一数组是存在于 CUDA 还是 OpenCL 设备上。关于这一点的更多信息，请参阅「内存」部分。

GPUArrays 有助于减少代码重复，因为它允许编写独立于硬件的 GPU 内核，这些内核可以通过 CuArrays 或 CLArrays 编译到本地的 GPU 代码。因此，大多通用内核可以在从 GPUArrays 继承的所有包之间共享。

选择小贴士：CuArrays 只支持 Nvidia GPU，而 CLArrays 支持大多数可用的 GPU。CuArrays 比 CLArrays 更稳定，可以在 Julia 0.7 上使用。速度上两者大同小异。我建议都试一试，看看哪种最有效。

本文中，我将选择 CuArrays，因为本文是在 Julia 0.7 / 1.0 上编写的，CLArrays 暂不支持。

性能

用一个简单的交互式代码示例来快速说明：为了计算 julia 集合（曼德勃罗集合），我们必须要将计算转移到 GPU 上。

对于大型数组，通过将计算转移到 GPU，可以稳定地将速度提高 60-80 倍。获得此加速和将 Julia 数组转换为 GPUArray 一样简单。

有人可能认为 GPU 性能会受到像 Julia 这样的动态语言影响，但 Julia 的 GPU 性能应该与 CUDA 或 OpenCL 的原始性能相当。Tim Besard 在集成 LLVM Nvidia 编译流程方面做得很好，能够实现与纯 CUDA C 语言代码相同（有时甚至更好）的性能。他在博客（https://devblogs.nvidia.com/gpu-computing-julia-programming-language/）中作了进一步解释。CLArrays 方法有点不同，它直接从 Julia 生成 OpenCL C 代码，代码性能与 OpenCL C 相同！

为了更好地了解性能并与多线程 CPU 代码进行比对，我整理了一些基准：https://github.com/JuliaGPU/GPUBenchmarks.jl/blob/master/results/results.md

内存

GPU 具有自己的存储空间，包括显存（VRAM）、不同的高速缓存和寄存器。无论做什么，运行前都要先将 Julia 对象转移到 GPU。并非 Julia 中的所有类型都可以在 GPU 上运行。

首先让我们看一下 Julia 的类型：

"Array"

所有这些 Julia 类型在传输到 GPU 或在 GPU 上创建时表现不同。下表概述了预期结果：

创建位置描述对象是在 CPU 上创建的，然后转移到 GPU 内核上，或者本身就由内核内部的 GPU 创建。该表显示创建类型的实例是否可行，对于从 CPU 到 GPU 的转移，该表还说明了对象是否能通过参照进行复制或传递。

垃圾收集

当使用 GPU 时，要注意 GPU 上没有垃圾收集器（GC）。这不会造成太大影响，因为写入 GPU 的高性能内核不应该创建任何 GC-跟踪的内存作为起始。

在 GPU 上实现 GC 不无可能，但请记住，每个执行内核都是大规模并行的。在大约 1000 个 gpu 线程中的每一个创建和跟踪大量堆内存就会马上破坏性能增益，因此实现 GC 是得不偿失的。

使用 GPUArrays 可以作为在内核中分配数组的替代方法。GPUArray 构造函数将创建 GPU 缓冲区并将数据转移到 VRAM。如果调用 Array(gpu_array)，数组将被转移回 RAM，变为普通的 Julia 数组。这些 gpu 数组的 Julia 操作由 Julia 的 GC 跟踪，如果不再使用，GPU 内存将被释放。

因此，只能在设备上使用堆栈分配，并且只能被其他的预先分配的 GPU 缓冲区使用。由于转移代价很高，因此在编写 GPU 时，往往要尽可能重用和预分配。

GPUArray 构造函数

"CuArray"

数组操作

我们已经定义了许多操作。最重要的是，GPUArrays 支持 Julia 的融合点广播表示法（fusing dot broadcasting notation）。此表示法允许你将函数应用于数组的每个元素，并使用 f 的返回值创建新数组。此功能通常称为映射（map）。broadcast 指的是形状各异的数组被 broadcast 成相同形状。

工作原理如下：

发生「融合」是因为 Julia 编译器会重写该表达式为一个传递调用树的 lazy broadcast 调用，然后可以在循环遍历数组之前将整个调用树融合到一个函数中。

如果你想要更详细的了解 broadcast，可以看该指南：julia.guide/broadcasting。

这意味着在不分配堆内存（仅创建 isbits 类型）的情况下运行的任何 Julia 函数，都可以应用于 GPUArray 的每个元素，并且多点调用会融合到一个内核调用中。由于内核调用会有很大延迟，所以这种融合是一个非常重要的优化。

true

GPUArrays 支持更多操作：

实现 GPU 数组转换为 CPU 数组和复制

多维索引和切片 (xs[1:2, 5, :])

permutedims

串联 (vcat(x, y), cat(3, xs, ys, zs))

映射，融合 broadcast(zs .= xs.^2 .+ ys .* 2)

填充 (CuArray, 0f0, dims)，填充 (gpu_array, 0)

减小尺寸 (reduce(+, xs, dims = 3), sum(x -> x^2, xs, dims = 1)

缩减为标量 (reduce(*, xs), sum(xs), prod(xs))

各种 BLAS 操作 (matrix*matrix, matrix*vector)

FFT，使用与 julia 的 FFT 相同的 API

GPUArrays 实际应用

让我们直接看一些很酷的实例。

GPU 加速烟雾模拟器是由 GPUArrays + CLArrays 创建的，可在 GPU 或 CPU 上运行，GPU 版本速度提升 15 倍：

还有更多的例子，包括求微分方程、FEM 模拟和求解偏微分方程。

演示地址：https://juliagpu.github.io/GPUShowcases.jl/latest/index.html

让我们通过一个简单的机器学习示例，看看如何使用 GPUArrays：

只需将数组转换为 GPUArrays（使用 gpu(array），就可以将整个计算移动到 GPU 并获得可观的速度提升。这要归功于 Julia 复杂的 AbstractArray 基础架构，使 GPUArray 可以无缝集成。随后，如果省略转换为 GPUArray 这一步，代码会按普通的 Julia 数组处理，但仍在 CPU 上运行。可以尝试将 use_gpu = true 改为 use_gpu = false，重新运行初始化和训练单元格。对比 GPU 和 CPU，CPU 运行时间为 975 秒，GPU 运行时间为 29 秒，速度提升约 33 倍。

编写 GPU 内核

一般情况，只使用 GPUArrays 的通用抽象数组接口即可，而不需要编写任何 GPU 内核。但是有些时候，可能需要在 GPU 上实现一个无法通过一般数组算法组合表示的算法。

好消息是，GPUArrays 通过分层法消除了大量工作，可以实现从高级代码开始，编写类似于大多数 OpenCL / CUDA 示例的低级内核。同时可以在 OpenCL 或 CUDA 设备上执行内核，从而提取出这些框架中的所有差异。

实现上述功能的函数名为 gpu_call。调用语句为 gpu_call(kernel, A::GPUArray, args)，在 GPU 上使用参数 (state, args...) 调用 kernel。State 是一个用于实现获取线程索引等功能的后端特定对象。GPUArray 需要作为第二个参数传递，以分配到正确的后端并提供启动参数的默认值。

让我们使用 gpu_call 来实现一个简单的映射内核：

简单来说，这将在 GPU 上并行调用 julia 函数 kernel length(A) 次。kernel 的每个并行调用都有一个线程索引，可以利用它索引到数组 A 和 B。如果计算索引时没有使用 linear_index，就需要确保没有多个线程读取和写入相同的数组位置。因此，如果在纯 Julia 中使用线程编写，可等效如下：

由于该函数未实现过多内容，也得不到更多的扩展，但线程化和 GPU 版本仍然有一个很好的加速。

GPU 与线程示例相比，能显示更复杂的内容，因为硬件线程是以线程块的形式分布的，gpu_call 是从简单版本中提取出来的，但它也可以用于更复杂的启动配置：

上面的示例中启动配置的迭代顺序更复杂。确定合适的迭代+启动配置对于实现最优 GPU 性能至关重要。很多关于 CUDA 和 OpenCL 的 GPU 教程都非常详细地解释了这一点，在 Julia 中编程 GPU 时这些原理是相通的。

结论

Julia 为高性能的世界带来了可组合的高级编程。现在是时候为 GPU 做同样的事了。

希望 Julia 能降低人们在 GPU 编程的门槛，我们可以为开源 GPU 计算开发可扩展的平台。第一个成功案例是通过 Julia 软件包实现自动微分解决方案，这些软件包甚至都不是为 GPU 编写的，因此可以相信 Julia 在 GPU 计算领域的扩展性和通用设计中一定会大放异彩。

本文为机器之心编译，转载请联系本公众号获得授权。

------------------------------------------------