當前位置：
首頁 > 科技 > 手把手教你如何用Julia做GPU編程（附代碼）

手把手教你如何用Julia做GPU編程（附代碼）

科技 10-22

原標題：手把手教你如何用Julia做GPU編程（附代碼）



來源：nextjournal


編輯：肖琴、三石


【新智元導讀】本文旨在快速介紹GPU的工作原理，詳細介紹當前的Julia GPU生態系統，並讓讀者了解簡單的GPU編程是多麼的容易。

GPU是如何工作的？


首先，什麼是GPU？


GPU是一個大規模並行處理器，具有幾千個並行處理單元。 例如，本文中使用的Tesla k80提供4992個並行CUDA內核。 GPU在頻率，延遲和硬體功能方面與CPU完全不同，但有點類似於擁有4992個內核的慢速CPU！



「Tesla K80」


可啟用並行線程的數量可以大幅提高GPU速度，但也讓它的使用性變得更加困難。讓我們來詳細看看在使用這種原始動力時，你會遇到哪些缺點：

  • GPU是一個獨立的硬體，具有自己的內存空間和不同的架構。 因此，從RAM到GPU存儲器（VRAM）的傳輸時間很長。 即使在GPU上啟動內核（換句話說，調度函數調用）也會帶來較大的延遲。 GPU的時間約為10us，而CPU的時間則為幾納秒。
  • 在沒有高級包裝器的情況下，設置內核會很快變得複雜
  • 較低的精度是默認值，而較高的精度計算可以輕鬆地消除所有性能增益
  • GPU函數(內核)本質上是並行的，所以編寫GPU內核至少和編寫並行CPU代碼一樣困難，但是硬體上的差異增加了相當多的複雜性
  • 與上述相關，許多演算法都不能很好地移植到GPU上。
  • 內核通常是用C/ C++編寫的，這並不是寫演算法的最佳語言。
  • CUDA和OpenCL之間存在分歧，OpenCL是用於編寫低級GPU代碼的主要框架。雖然CUDA只支持英偉達硬體，但OpenCL支持所有硬體，但有些粗糙。

Julia的誕生是個好消息！它是一種高級腳本語言，允許你在Julia本身編寫內核和周圍的代碼，同時在大多數GPU硬體上運行！


GPUArrays


大多數高度並行的演算法需要通過相當多的數據來克服所有線程和延遲開銷。因此，大多數演算法都需要數組來管理所有數據，這需要一個好的GPU數組庫（array library）作為關鍵基礎。


GPUArrays.jl是Julia的基礎。它提供了一個抽象數組實現，專門用於使用高度並行硬體的原始功能。它包含設置GPU所需的所有功能，啟動Julia GPU函數並提供一些基本的數組演算法。

抽象意味著它需要以CuArrays和CLArrays形式的具體實現。由於繼承了GPUArrays的所有功能，它們都提供完全相同的介面。唯一的區別出現在分配數組時，這會強制你決定數組是否位於CUDA或OpenCL設備上。關於這一點的更多信息，請參閱內存部分。


GPUArrays有助於減少代碼重複，因為它允許編寫獨立於硬體的GPU內核，可以通過CuArrays或CLArrays將其編譯為本機GPU代碼。因此，許多通用內核可以在繼承自GPUArrays的所有packages之間共享。


一點選擇建議：CuArrays僅適用於Nvidia GPU，而CLArrays適用於大多數可用的GPU。CuArrays比CLArrays更穩定，並且已經可以在Julia 0.7上運行。速度上差異不明顯。我建議兩者都試一下，看看哪個效果最好。


對於本文，我將選擇CuArrays，因為本文是為Julia 0.7 / 1.0而寫的，CLArrays仍然不支持。


性能


讓我們用一個簡單的互動式代碼示例來快速說明為什麼要將計算轉移到GPU上，這個示例計算julia set：


1using CuArrays, FileIO, Colors, GPUArrays, BenchmarkTools


2using CuArrays: CuArray


3"""


4The function calculating the Julia set

5"""


6function juliaset(z0, maxiter)


7c = ComplexF32(-0.5, 0.75)


8z = z0


9fori in1:maxiter


10abs2(z) > 4f0 && return(i - 1) % UInt8


11z = z * z + c


12end


13returnmaxiter % UInt8 # % is used to convert without overflow check


14end

15range = 100:50:2^12


16cutimes, jltimes = Float64[], Float64[]


17function run_bench(in, out)


18# use dot syntax to apply `juliaset` to each elemt of q_converted


19# and write the output to result


20out .= juliaset.(in, 16)


21# all calls to the GPU are scheduled asynchronous,


22# so we need to synchronize


23GPUArrays.synchronize(out)


24end

25# store a reference to the last results for plotting


26last_jl, last_cu = nothing, nothing


27forN inrange


28w, h = N, N


29q = [ComplexF32(r, i) fori=1:-(2.0/w):-1, r=-1.5:(3.0/h):1.5]


30for(times, Typ) in((cutimes, CuArray), (jltimes, Array))


31# convert to Array or CuArray - moving the calculation to CPU/GPU


32q_converted = Typ(q)


33result = Typ(zeros(UInt8, size(q)))


34fori in1:10# 5 samples per size

35# benchmarking macro, all variables need to be prefixed with $


36t = Base.@elapsed begin


37run_bench(q_converted, result)


38end


39globallast_jl, last_cu # we"re in local scope


40ifresult isa CuArray


41last_cu = result


42else


43last_jl = result


44end

45push!(times, t)


46end


47end


48end


49


50cu_jl = hcat(Array(last_cu), last_jl)


51cmap = colormap("Blues", 16+ 1)


52color_lookup(val, cmap) = cmap[val + 1]


53save("results/juliaset.png", color_lookup.(cu_jl, (cmap,)))


1using Plots; plotly()


2x = repeat(range, inner = 10)


3speedup = jltimes ./ cutimes


4Plots.scatter(


5log2.(x), [speedup, fill(1.0, length(speedup))],


6label = ["cuda""cpu"], markersize = 2, markerstrokewidth = 0,


7legend = :right, xlabel = "2^N", ylabel = "speedup"


8)



如你所見，對於大型數組，通過將計算移動到GPU可以獲得穩定的60-80倍的加速。而且非常簡單，只需將Julia array轉換為GPUArray。

有人可能認為GPU的性能受到像Julia這樣的動態語言的影響，但Julia的GPU性能應該與CUDA或OpenCL的原始性能相當。Tim Besard在集成LLVM Nvidia編譯pipeline方面做得非常出色，達到了與純CUDA C代碼相同（有時甚至更好）的性能。Tim發表了一篇非常詳細的博文，裡面進一步解釋了這一點[1]。CLArrays方法有點不同，它直接從Julia生成OpenCL C代碼，具有與OpenCL C相同的性能！


為了更好地了解性能並查看與多線程CPU代碼的比較，我收集了一些基準測試[2]。


內存（Memory）


GPU具有自己的存儲空間，包括視頻存儲器（VRAM），不同的高速緩存和寄存器。無論你做什麼，任何Julia對象都必須先轉移到GPU才能使用。並非Julia中的所有類型都可以在GPU上工作。


首先讓我們看一下Julia的類型：


1struct Test # an immutable struct


2# that only contains other immutable, which makes


3# isbitstype(Test) == true


4x::Float32


5end


6


7# the isbits property is important, since those types can be used


8# without constraints on the GPU!


9@assert isbitstype(Test) == true


10x = (2, 2)


11isa(x, Tuple{Int, Int}) # tuples are also immutable


12mutable struct Test2 #-> mutable, isbits(Test2) == false


13x::Float32


14end


15struct Test3


16# contains a heap allocation/ reference, not isbits


17x::Vector{Float32}


18y::Test2 # Test2 is mutable and also heap allocated / a reference


19end


20Vector{Test} # <- An Array with isbits elements is contigious in memory


21Vector{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的傳輸，該表還指示對象是否通過引用複製或傳遞。


Garbage Collection


使用GPU時的一個很大的區別是GPU上沒有垃圾回收（ garbage collector, GC）。這不是什麼大問題，因為為GPU編寫的高性能內核不應該一開始就創建任何GC-tracked memory。


為GPU實現GC是可能的，但請記住，每個執行的內核都是大規模並行的。在~1000 GPU線程中的每一個線程創建和跟蹤大量堆內存將很快破壞性能增益，因此這實際上是不值得的。


作為內核中堆分配數組的替代方法，你可以使用GPUArrays。GPUArray構造函數將創建GPU緩衝區並將數據傳輸到VRAM。如果調用Array(gpu_array)，數組將被轉移回RAM，表示為普通的Julia數組。這些GPU數組的Julia句柄由Julia的GC跟蹤，如果它不再使用，GPU內存將被釋放。


因此，只能在設備上使用堆棧分配，並且對其餘的預先分配的GPU緩衝區使用。由於傳輸非常昂貴的，因此在編程GPU時儘可能多地重用和預分配是很常見的。


The GPUArray Constructors


1using CuArrays, LinearAlgebra


2


3# GPU Arrays can be constructed from all Julia arrays containing isbits types!


4A1D = cu([1, 2, 3]) # cl for CLArrays


5A1D = fill(CuArray{Int}, 0, (100,)) # CLArray for CLArrays


6# Float32 array - Float32 is usually preferred and can be up to 30x faster on most GPUs than Float64


7diagonal_matrix = CuArray{Float32}(I, 100, 100)


8filled = fill(CuArray, 77f0, (4, 4, 4)) # 3D array filled with Float32 77


9randy = rand(CuArray, Float32, 42, 42) # random numbers generated on the GPU


10# The array constructor also accepts isbits iterators with a known size


11# 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


12from_iter = CuArray(1:10)


13# let"s create a point type to further illustrate what can be done:


14struct Point


15x::Float32


16y::Float32


17end


18Base.convert(::Type{Point}, x::NTuple{2, Any}) = Point(x[1], x[2])


19# because we defined the above convert from a tuple to a point


20# [Point(2, 2)] can be written as Point[(2,2)] since all array


21# elements will get converted to Point


22custom_types = cu(Point[(1, 2), (4, 3), (2, 2)])


23typeof(custom_types)


"CuArray{Point, 1}"


Array Operations


許多操作是已經定義好的。最重要的是，GPUArrays支持Julia的fusing dot broadcasting notation。這種標記法允許你將函數應用於數組的每個元素，並使用f的返回值創建一個新數組。這個功能通常稱為映射（map）。 broadcast 指的是具有不同形狀的數組被散布到相同的形狀。


它的工作方式如下：


1x = zeros(4, 4) # 4x4 array of zeros


2y = zeros(4) # 4 element array


3z = 2# a scalar


4# y"s 1st dimension gets repeated for the 2nd dimension in x


5# and the scalar z get"s repeated for all dimensions


6# the below is equal to `broadcast(+, broadcast(+, xx, y), z)`


7x .+ y .+ z


關於broadcasting如何工作的更多解釋，可以看看這個指南：


julia.guide/broadcasting


這意味著在不分配堆內存（僅創建isbits類型）的情況下運行的任何Julia函數都可以應用於GPUArray的每個元素，並且多個dot調用將融合到一個內核調用中。由於內核調用延遲很高，這種融合是一個非常重要的優化。


1using CuArrays


2A = cu([1, 2, 3])


3B = cu([1, 2, 3])


4C = rand(CuArray, Float32, 3)


5result = A .+ B .- C


6test(a::T) where T = a * convert(T, 2) # convert to same type as `a`


7


8# inplace broadcast, writes directly into `result`


9result .= test.(A) # custom function work


10


11# The cool thing is that this composes well with custom types and custom functions.


12# Let"s go back to our Point type and define addition for it


13Base.:(+)(p1::Point, p2::Point) = Point(p1.x + p2.x, p1.y + p2.y)


14


15# now this works:


16custom_types = cu(Point[(1, 2), (4, 3), (2, 2)])


17


18# This particular example also shows the power of broadcasting:


19# Non array types are broadcasted and repeated for the whole length


20result = custom_types .+ Ref(Point(2, 2))


21


22# So the above is equal to (minus all the allocations):


23# this allocates a new array on the gpu, which we can avoid with the above broadcast


24broadcasted = fill(CuArray, Point(2, 2), (3,))


25


26result == custom_types .+ broadcasted


ture


現實世界中的GPUArrays


讓我們直接看看一些很酷的用例。


如下面的視頻所示，這個GPU加速煙霧模擬是使用GPUArrays + CLArrays創建的，可在GPU或CPU上運行，GPU版本的速度提高了15倍：



還有更多的用例，包括求解微分方程，有限元模擬和求解偏微分方程。


讓我們來看一個簡單的機器學習示例，看看如何使用GPUArrays：


1using Flux, Flux.Data.MNIST, Statistics


2using Flux: onehotbatch, onecold, crossentropy, throttle


3using Base.Iterators: repeated, partition


4using CuArrays


5


6# Classify MNIST digits with a convolutional network


7


8imgs = MNIST.images()


9


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


11


12# Partition into batches of size 1,000


13train = [(cat(float.(imgs[i])..., dims = 4), labels[:,i])


14fori inpartition(1:60_000, 1000)]


15


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


17


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


19


20train = todevice.(train)


21


22# Prepare test set (first 1,000 images)


23tX = cat(float.(MNIST.images(:test)[1:1000])..., dims = 4) |> todevice


24tY = onehotbatch(MNIST.labels(:test)[1:1000], 0:9) |> todevice


25


26m = Chain(


27Conv((2,2), 1=>16, relu),


28x -> maxpool(x, (2,2)),


29Conv((2,2), 16=>8, relu),


30x -> maxpool(x, (2,2)),


31x -> reshape(x, :, size(x, 4)),


32Dense(288, 10), softmax) |> todevice


33


34m(train[1][1])


35


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


37


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


39


40evalcb = throttle(() -> @show(accuracy(tX, tY)), 10)


41opt = ADAM(Flux.params(m));


1# train


2fori = 1:10


3Flux.train!(loss, train, opt, cb = evalcb)


4end


accuracy(tX, tY) = 0.101


accuracy(tX, tY) = 0.888


accuracy(tX, tY) = 0.919


1using Colors, FileIO, ImageShow


2N = 22


3img = tX[:, :, 1:1, N:N]


4println("Predicted: ", Flux.onecold(m(img)) .- 1)


5Gray.(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上高效運行的代碼。這有助於幫助Flux以最少的開發人員在GPU上工作，並使Flux GPU能夠有效地支持用戶定義的函數。在沒有GPUArrays + Flux之間協調的情況下開箱即用是Julia的一個非常獨特的特性，詳細解釋見[3].


編寫GPU內核


只需使用GPUArrays的通用抽象數組介面，而不用編寫任何GPU內核，就可以做很多事了。但是，在某些時候，可能需要實現一個需要在GPU上運行的演算法，並且不能用通用數組演算法的組合來表示。


好的一點是，GPUArrays通過一種分層方法減少了大量的工作，這種方法允許你從高級代碼開始編寫低級內核，類似於大多數OpenCL / CUDA示例里的。它還允許你在OpenCL或CUDA設備上執行內核，從而抽象出這些框架中的任何差異。


使這成為可能的函數名為gpu_call。它可以被稱為 gpu_call(kernel, A::GPUArray, args)，並將在GPU上使用參數 (state, args...) 調用內核。State是一個後端特定對象，用於實現獲取線程索引之類的功能。GPUArray需要作為第二個參數傳遞，一遍分派到正確的後端並提供啟動參數的預設值。


讓我們使用gpu_call來實現一個簡單的map kernel：


1using GPUArrays, CuArrays


2# Overloading the Julia Base map! function for GPUArrays


3function Base.map!(f::Function, A::GPUArray, B::GPUArray)


4# our function that will run on the gpu


5function kernel(state, f, A, B)


6# If launch parameters aren"t specified, linear_index gets the index


7# into the Array passed as second argument to gpu_call (`A`)


8i = linear_index(state)


9ifi <= length(A)


10@inbounds A[i] = f(B[i])


11end


12return


13end


14# call kernel on the gpu


15gpu_call(kernel, A, (f, A, B))


16end


簡單來說，上面的代碼將在GPU上並行調用julia函數內核length(A) 次。內核的每個並行調用都有一個線程索引，我們可以使用它來安全地索引到數組A和B。如果我們計算自己的索引，而不是使用linear_index，我們需要確保沒有多個線程讀寫同一個數組位置。因此，如果我們使用線程在純Julia中編寫，其對應版本如下：


1using BenchmarkTools


2function threadded_map!(f::Function, A::Array, B::Array)


3Threads.@threads fori in1:length(A)


4A[i] = f(B[i])


5end


6A


7end


8x, y = rand(10^7), rand(10^7)


9kernel(y) = (y / 33f0) * (732.f0/y)


10# on the cpu without threads:


11single_t = @belapsed map!($kernel, $x, $y)


12


13# "on the CPU with 4 threads (2 real cores):


14thread_t = @belapsed threadded_map!($kernel, $x, $y)


15


16# on the GPU:


17xgpu, ygpu = cu(x), cu(y)


18gpu_t = @belapsed begin


19map!($kernel, $xgpu, $ygpu)


20GPUArrays.synchronize($xgpu)


21end


22times = [single_t, thread_t, gpu_t]


23speedup = maximum(times) ./ times


24println("speedup: $speedup")


25bar(["1 core", "2 cores", "gpu"], speedup, legend = false, fillcolor = :grey, ylabel = "speedup")



因為這個函數沒有做很多工作，我們看不到完美的擴展，但線程和GPU版本仍然提供了很大的加速。


GPU比線程示例展示的要複雜得多，因為硬體線程是在線程塊中布局的——gpu_call在簡單版本中抽象出來，但它也可以用於更複雜的啟動配置：


1using CuArrays


2


3threads = (2, 2)


4blocks = (2, 2)


5T = fill(CuArray, (0, 0), (4, 4))


6B = fill(CuArray, (0, 0), (4, 4))


7gpu_call(T, (B, T), (blocks, threads)) do state, A, B


8# those names pretty much refer to the cuda names


9b = (blockidx_x(state), blockidx_y(state))


10bdim = (blockdim_x(state), blockdim_y(state))


11t = (threadidx_x(state), threadidx_y(state))


12idx = (bdim .* (b .- 1)) .+ t


13A[idx...] = b


14B[idx...] = t


15return


16end


17println("Threads index: n", T)


18println("Block index: n", B)


在上面的示例中，你可以看到更複雜的啟動配置的迭代順序。確定正確的迭代+啟動配置對於達到GPU的最佳性能至關重要。


結論


在將可組合的高級編程引入高性能世界方面，Julia取得了長足的進步。現在是時候對GPU做同樣的事情了。


希望Julia降低開始在GPU上編程的標準，並且我們可以為開源GPU計算髮展可擴展的平台。第一個成功案例是通過Julia packages實現自動微分，這些軟體包甚至不是為GPU編寫，因此這給了我們很多理由相信Julia在GPU計算領域的可擴展和通用設計是成功的。


[1]https://devblogs.nvidia.com/gpu-computing-julia-programming-language/


[2]https://github.com/JuliaGPU/GPUBenchmarks.jl/blob/master/results/results.md


[3]www.stochasticlifestyle.com/why-numba-and-cython-are-not-substitutes-for-julia/


原文地址：


https://nextjournal.com/sdanisch/julia-gpu-programming


愛奇藝

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！

本站內容充實豐富，博大精深，小編精選每日熱門資訊，隨時更新，點擊「搶先收到最新資訊」瀏覽吧！


請您繼續閱讀更多來自 新智元 的精彩文章:

※離職亞馬遜加盟英偉達，明星女科學家擔綱英偉達機器學習研究
※谷歌AI超大規模圖像競賽，中國團隊獲目標檢測冠軍

TAG:新智元 |