Simon Danisch作者高璇 刘晓坤参与

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

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 上。

using CuArrays, FileIO, Colors, GPUArrays, BenchmarkTools
using CuArrays: CuArray
"""
The function calculating the Julia set
"""
function juliaset(z0, maxiter)
    c = ComplexF32(-0.5, 0.75)
    z = z0
    for i in 1:maxiter
        abs2(z) > 4f0 && return (i - 1) % UInt8
        z = z * z + c
    end
    return maxiter % UInt8 # % is used to convert without overflow check
end
range = 100:50:2^12
cutimes, jltimes = Float64[], Float64[]
function run_bench(in, out)
  # use dot syntax to apply `juliaset` to each elemt of q_converted 
  # and write the output to result
  out .= juliaset.(in, 16)
  # all calls to the GPU are scheduled asynchronous, 
  # so we need to synchronize
  GPUArrays.synchronize(out)
end
# store a reference to the last results for plotting
last_jl, last_cu = nothing, nothing
for N in range
  w, h = N, N
  q = [ComplexF32(r, i) for i=1:-(2.0/w):-1, r=-1.5:(3.0/h):1.5]
  for (times, Typ) in ((cutimes, CuArray), (jltimes, Array))
    # convert to Array or CuArray - moving the calculation to CPU/GPU
    q_converted = Typ(q)
    result = Typ(zeros(UInt8, size(q)))
    for i in 1:10 # 5 samples per size
      # benchmarking macro, all variables need to be prefixed with $
      t = Base.@elapsed begin
                run_bench(q_converted, result)
      end
      global last_jl, last_cu # we're in local scope
      if result isa CuArray
        last_cu = result
      else
          last_jl = result
      end
      push!(times, t)
    end
  end
end

cu_jl = hcat(Array(last_cu), last_jl)
cmap = colormap("Blues", 16 + 1)
color_lookup(val, cmap) = cmap[val + 1]
save("results/juliaset.png", color_lookup.(cu_jl, (cmap,)))

using Plots; plotly()
x = repeat(range, inner = 10)
speedup = jltimes ./ cutimes
Plots.scatter(
  log2.(x), [speedup, fill(1.0, length(speedup))], 
  label = ["cuda" "cpu"], markersize = 2, markerstrokewidth = 0,
  legend = :right, xlabel = "2^N", ylabel = "speedup"
)

对于大型数组,通过将计算转移到 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 的类型:

struct Test # an immutable struct
# that only contains other immutable, which makes 
# isbitstype(Test) == true
    x::Float32 
end

# the isbits property is important, since those types can be used
# without constraints on the GPU!
@assert isbitstype(Test) == true
x = (2, 2)
isa(x, Tuple{Int, Int}) # tuples are also immutable
mutable struct Test2 #-> mutable, isbits(Test2) == false
    x::Float32
end
struct Test3
    # contains a heap allocation/ reference, not isbits
    x::Vector{Float32}
    y::Test2 # Test2 is mutable and also heap allocated / a reference
end
Vector{Test} # <-  An Array with isbits elements is contigious in memory
Vector{Test2} # <- An Array with mutable elements is basically an array of heap pointers. Since it just contains cpu heap pointers, it won't work on the GPU.

"Array{Test2,1}"

所有这些 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 构造函数

using CuArrays, LinearAlgebra

# GPU Arrays can be constructed from all Julia arrays containing isbits types!
A1D = cu([1, 2, 3]) # cl for CLArrays
A1D = fill(CuArray{Int}, 0, (100,)) # CLArray for CLArrays
# Float32 array - Float32 is usually preferred and can be up to 30x faster on most GPUs than Float64
diagonal_matrix = CuArray{Float32}(I, 100, 100)
filled = fill(CuArray, 77f0, (4, 4, 4)) # 3D array filled with Float32 77
randy = rand(CuArray, Float32, 42, 42) # random numbers generated on the GPU
# The array constructor also accepts isbits iterators with a known size
# Note, that since you can also pass isbits types to a gpu kernel directly, in most cases you won't need to materialize them as an gpu array
from_iter = CuArray(1:10)
# let's create a point type to further illustrate what can be done:
struct Point
    x::Float32
    y::Float32
end
Base.convert(::Type{Point}, x::NTuple{2, Any}) = Point(x[1], x[2])
# because we defined the above convert from a tuple to a point
# [Point(2, 2)] can be written as Point[(2,2)] since all array 
# elements will get converted to Point
custom_types = cu(Point[(1, 2), (4, 3), (2, 2)])
typeof(custom_types)

"CuArray{point,1}"

数组操作

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

工作原理如下:

x = zeros(4, 4) # 4x4 array of zeros
y = zeros(4) # 4 element array
z = 2 # a scalar
# y's 1st dimension gets repeated for the 2nd dimension in x
# and the scalar z get's repeated for all dimensions
# the below is equal to `broadcast(+, broadcast(+, xx, y), z)`
x .+ y .+ z

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

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

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

using CuArrays
A = cu([1, 2, 3])
B = cu([1, 2, 3])
C = rand(CuArray, Float32, 3)
result = A .+ B .- C
test(a::T) where T = a * convert(T, 2) # convert to same type as `a`

# inplace broadcast, writes directly into `result`
result .= test.(A) # custom function work

# The cool thing is that this composes well with custom types and custom functions.
# Let's go back to our Point type and define addition for it
Base.:(+)(p1::Point, p2::Point) = Point(p1.x + p2.x, p1.y + p2.y)

# now this works:
custom_types = cu(Point[(1, 2), (4, 3), (2, 2)])

# This particular example also shows the power of broadcasting: 
# Non array types are broadcasted and repeated for the whole length
result = custom_types .+ Ref(Point(2, 2))

# So the above is equal to (minus all the allocations):
# this allocates a new array on the gpu, which we can avoid with the above broadcast
broadcasted = fill(CuArray, Point(2, 2), (3,))

result == custom_types .+ broadcasted

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:

using Flux, Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated, partition
using CuArrays

# Classify MNIST digits with a convolutional network

imgs = MNIST.images()

labels = onehotbatch(MNIST.labels(), 0:9)

# Partition into batches of size 1,000
train = [(cat(float.(imgs[i])..., dims = 4), labels[:,i])
         for i in partition(1:60_000, 1000)]

use_gpu = true # helper to easily switch between gpu/cpu

todevice(x) = use_gpu ? gpu(x) : x

train = todevice.(train)

# Prepare test set (first 1,000 images)
tX = cat(float.(MNIST.images(:test)[1:1000])..., dims = 4) |> todevice
tY = onehotbatch(MNIST.labels(:test)[1:1000], 0:9) |> todevice

m = Chain(
  Conv((2,2), 1=>16, relu),
  x -> maxpool(x, (2,2)),
  Conv((2,2), 16=>8, relu),
  x -> maxpool(x, (2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(288, 10), softmax) |> todevice

m(train[1][1])

loss(x, y) = crossentropy(m(x), y)

accuracy(x, y) = mean(onecold(m(x)) .== onecold(y))

evalcb = throttle(() -> @show(accuracy(tX, tY)), 10)
opt = ADAM(Flux.params(m));
# train
for i = 1:10
    Flux.train!(loss, train, opt, cb = evalcb)
end

using Colors, FileIO, ImageShow
N = 22
img = tX[:, :, 1:1, N:N]
println("Predicted: ", Flux.onecold(m(img)) .- 1)
Gray.(collect(tX[:, :, 1, N]))

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

另一个优势是为了有效地支持神经网络的反向传播,GPUArrays 无需明确地实现自动微分。这是因为 Julia 的自动微分库适用于任意函数,并存有可在 GPU 上高效运行的代码。这样即可利用最少的开发人员就能在 GPU 上实现 Flux,并使 Flux GPU 能够高效实现用户定义的功能。这种开箱即用的 GPUArrays + Flux 不需要协调,这是 Julia 的一大特点,详细解释如下:为什么 Numba 和 Cython 不能代替 Julia(http://www.stochasticlifestyle.com/why)。

编写 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 来实现一个简单的映射内核:

using GPUArrays, CuArrays
# Overloading the Julia Base map! function for GPUArrays
function Base.map!(f::Function, A::GPUArray, B::GPUArray)
    # our function that will run on the gpu
    function kernel(state, f, A, B)
        # If launch parameters aren't specified, linear_index gets the index
        # into the Array passed as second argument to gpu_call (`A`)
        i = linear_index(state)
            if i <= length(A)
          @inbounds A[i] = f(B[i])
        end
        return
    end
    # call kernel on the gpu
    gpu_call(kernel, A, (f, A, B))
end

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

using BenchmarkTools
function threadded_map!(f::Function, A::Array, B::Array)
    Threads.@threads for i in 1:length(A)
        A[i] = f(B[i])
    end
  A
end
x, y = rand(10^7), rand(10^7)
kernel(y) = (y / 33f0) * (732.f0/y)
# on the cpu without threads:
single_t = @belapsed map!($kernel, $x, $y)

# "on the CPU with 4 threads (2 real cores):
thread_t = @belapsed threadded_map!($kernel, $x, $y)

# on the GPU:
xgpu, ygpu = cu(x), cu(y)
gpu_t = @belapsed begin
  map!($kernel, $xgpu, $ygpu)
  GPUArrays.synchronize($xgpu)
end
times = [single_t, thread_t, gpu_t]
speedup = maximum(times) ./ times
println("speedup: $speedup")
bar(["1 core", "2 cores", "gpu"], speedup, legend = false, fillcolor = :grey, ylabel = "speedup")


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

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

using CuArrays

threads = (2, 2)
blocks = (2, 2)
T = fill(CuArray, (0, 0), (4, 4))
B = fill(CuArray, (0, 0), (4, 4))
gpu_call(T, (B, T), (blocks, threads)) do state, A, B
  # those names pretty much refer to the cuda names
    b = (blockidx_x(state), blockidx_y(state))
    bdim = (blockdim_x(state), blockdim_y(state))
    t = (threadidx_x(state), threadidx_y(state))
    idx = (bdim .* (b .- 1)) .+ t
    A[idx...] = b
    B[idx...] = t
    return
end
println("Threads index: \n", T)
println("Block index: \n", B)

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

结论

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

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


原文链接:https://nextjournal.com/sdanisch/julia-gpu-programming

工程编程语言JuliaGPU
1
相关数据
机器学习技术

机器学习是人工智能的一个分支,是一门多领域交叉学科,涉及概率论、统计学、逼近论、凸分析、计算复杂性理论等多门学科。机器学习理论主要是设计和分析一些让计算机可以自动“学习”的算法。因为学习算法中涉及了大量的统计学理论,机器学习与推断统计学联系尤为密切,也被称为统计学习理论。算法设计方面,机器学习理论关注可以实现的,行之有效的学习算法。

调度技术

调度在计算机中是分配工作所需资源的方法。资源可以指虚拟的计算资源,如线程、进程或数据流;也可以指硬件资源,如处理器、网络连接或扩展卡。 进行调度工作的程序叫做调度器。调度器通常的实现使得所有计算资源都处于忙碌状态,允许多位用户有效地同时共享系统资源,或达到指定的服务质量。 see planning for more details

Julia技术

Julia 是MIT设计的一个面向科学计算的高性能动态高级程序设计语言,项目大约于2009年中开始,2018年8月JuliaCon2018 发布会上发布Julia 1.0。据介绍,Julia 目前下载量已经达到了 200 万次,且 Julia 社区开发了超过 1900 多个扩展包。这些扩展包包含各种各样的数学库、数学运算工具和用于通用计算的库。除此之外,Julia 语言还可以轻松使用 Python、R、C/C++ 和 Java 中的库,这极大地扩展了 Julia 语言的使用范围。

基准技术

一种简单的模型或启发法,用作比较模型效果时的参考点。基准有助于模型开发者针对特定问题量化最低预期效果。

参数技术

在数学和统计学裡,参数(英语:parameter)是使用通用变量来建立函数和变量之间关系(当这种关系很难用方程来阐述时)的一个数量。

神经网络技术

(人工)神经网络是一种起源于 20 世纪 50 年代的监督式机器学习模型,那时候研究者构想了「感知器(perceptron)」的想法。这一领域的研究者通常被称为「联结主义者(Connectionist)」,因为这种模型模拟了人脑的功能。神经网络模型通常是通过反向传播算法应用梯度下降训练的。目前神经网络有两大主要类型,它们都是前馈神经网络:卷积神经网络(CNN)和循环神经网络(RNN),其中 RNN 又包含长短期记忆(LSTM)、门控循环单元(GRU)等等。深度学习是一种主要应用于神经网络帮助其取得更好结果的技术。尽管神经网络主要用于监督学习,但也有一些为无监督学习设计的变体,比如自动编码器和生成对抗网络(GAN)。

映射技术

映射指的是具有某种特殊结构的函数,或泛指类函数思想的范畴论中的态射。 逻辑和图论中也有一些不太常规的用法。其数学定义为:两个非空集合A与B间存在着对应关系f,而且对于A中的每一个元素x,B中总有有唯一的一个元素y与它对应,就这种对应为从A到B的映射,记作f:A→B。其中,y称为元素x在映射f下的象,记作:y=f(x)。x称为y关于映射f的原象*。*集合A中所有元素的象的集合称为映射f的值域,记作f(A)。同样的,在机器学习中,映射就是输入与输出之间的对应关系。

暂无评论
暂无评论~