GPU Computing

30.1 CUDA.jl (NVIDIA)

using CUDA
CUDA.functional()

a = CUDA.rand(10^7)                # array on the GPU
b = CUDA.ones(10^7)
c = a .+ 2 .* b                    # broadcast compiles to a GPU kernel
sum(c)                             # GPU reduction
A = CUDA.randn(1000, 1000); A * A  # cuBLAS
Array(c)                           # copy back to CPU
cu([1.0, 2.0])                     # move CPU data to GPU (Float32 by default)

# Custom kernel
function vadd!(c, a, b)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(c)
        @inbounds c[i] = a[i] + b[i]
    end
    return nothing
end
n = 10^6; a, b, c = CUDA.rand(n), CUDA.rand(n), CUDA.zeros(n)
@cuda threads=256 blocks=cld(n, 256) vadd!(c, a, b)
CUDA.@time f(a)

30.2 Other backends & portability

AMDGPU.jl (ROCm), Metal.jl (Apple Silicon), oneAPI.jl (Intel).

KernelAbstractions.jl — write one kernel for CPU and all GPU backends.

GPUArrays.jl generic array operations; most generic Julia code (broadcast, linear algebra, many SciML solvers, Flux/Lux models) runs on GPU arrays unchanged.

Tips: avoid scalar indexing on GPU arrays (CUDA.allowscalar(false)), use Float32, keep data on the device, batch operations.