Array programming

The easiest way to use the GPU's massive parallelism, is by expressing operations in terms of arrays: CUDA.jl provides an array type, CuArray, and many specialized array operations that execute efficiently on the GPU hardware. In this section, we will briefly demonstrate use of the CuArray type. Since we expose CUDA's functionality by implementing existing Julia interfaces on the CuArray type, you should refer to the upstream Julia documentation for more information on these operations.

More specific functionality, such as linear algebra or randon number generation, is defined in separate packages like cuBLAS.jl or cuRAND.jl, loosely matching the structure of the CUDA libraries. Importing the CUDA.jl package will also import these packages, so you can use most of their functionality without having to reach within. It is possible to opt out of this by only importing the relevant subpackages, alongside with the CUDACore.jl package providing the core functionality of CUDA.jl.

If you encounter missing functionality, or are running into operations that trigger so-called "scalar iteration", have a look at the issue tracker and file a new issue if there's none. Do note that you can always access the underlying CUDA APIs by calling into the relevant submodule. For example, if parts of the Random interface isn't properly implemented by CUDA.jl, you can look at the cuRAND documentation and possibly call methods from the cuRAND submodule directly. These submodules are available after importing the CUDA package.

Construction and Initialization

The CuArray type aims to implement the AbstractArray interface, and provide implementations of methods that are commonly used when working with arrays. That means you can construct CuArrays in the same way as regular Array objects:

julia> CuArray{Int}(undef, 2)2-element CuArray{Int64, 1}: 0 0julia> CuArray{Int}(undef, (1,2))1×2 CuArray{Int64, 2}: 0  0julia> similar(ans)1×2 CuArray{Int64, 2}: 0  0

Copying memory to or from the GPU can be expressed using constructors as well, or by calling copyto!:

julia> a = CuArray([1,2])2-element CuArray{Int64, 1, CUDACore.DeviceMemory}: 1 2julia> b = Array(a)2-element Vector{Int64}: 1 2julia> copyto!(b, a)2-element Vector{Int64}: 1 2

Higher-order abstractions

The real power of programming GPUs with arrays comes from Julia's higher-order array abstractions: Operations that take user code as an argument, and specialize execution on it. With these functions, you can often avoid having to write custom kernels. For example, to perform simple element-wise operations you can use map or broadcast:

julia> a = CuArray{Float32}(undef, (1,2));julia> a .= 51×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: 5.0  5.0julia> map(sin, a)1×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: -0.958924  -0.958924

To reduce the dimensionality of arrays, CUDA.jl implements the various flavours of (map)reduce(dim):

julia> a = CUDA.ones(2,3)2×3 CuArray{Float32, 2, CUDACore.DeviceMemory}: 1.0  1.0  1.0 1.0  1.0  1.0julia> reduce(+, a)6.0f0julia> mapreduce(sin, *, a; dims=2)2×1 CuArray{Float32, 2, CUDACore.DeviceMemory}: 0.59582335 0.59582335julia> b = CUDA.zeros(1)1-element CuArray{Float32, 1, CUDACore.DeviceMemory}: 0.0julia> Base.mapreducedim!(identity, +, b, a)1-element CuArray{Float32, 1, CUDACore.DeviceMemory}: 6.0

To retain intermediate values, you can use accumulate:

julia> a = CUDA.ones(2,3)2×3 CuArray{Float32, 2, CUDACore.DeviceMemory}: 1.0  1.0  1.0 1.0  1.0  1.0julia> accumulate(+, a; dims=2)2×3 CuArray{Float32, 2, CUDACore.DeviceMemory}: 1.0  2.0  3.0 1.0  2.0  3.0

Be wary that the operator f of accumulate and accumulate! must be associative since the operation is performed in parallel. That is, f(f(a,b),c) must be equivalent to f(a,f(b,c)). Accumulating with a non-associative operator on a CuArray will not produce the same result as on an Array.

Logical operations

CuArrays can also be indexed with arrays of boolean values to select items:

julia> a = CuArray([1,2,3])3-element CuArray{Int64, 1, CUDACore.DeviceMemory}: 1 2 3julia> a[[false,true,false]]1-element CuArray{Int64, 1, CUDACore.DeviceMemory}: 2

Built on top of this, are several functions with higher-level semantics:

julia> a = CuArray([11,12,13])3-element CuArray{Int64, 1, CUDACore.DeviceMemory}: 11 12 13julia> findall(isodd, a)2-element CuArray{Int64, 1, CUDACore.DeviceMemory}: 1 3julia> findfirst(isodd, a)1julia> b = CuArray([11 12 13; 21 22 23])2×3 CuArray{Int64, 2, CUDACore.DeviceMemory}: 11  12  13 21  22  23julia> findmin(b)(11, CartesianIndex(1, 1))julia> findmax(b; dims=2)(CuArray([13; 23;;]), CuArray(CartesianIndex{2}[CartesianIndex(1, 3); CartesianIndex(2, 3);;]))

Array wrappers

To some extent, CUDA.jl also supports well-known array wrappers from the standard library:

julia> a = CuArray(collect(1:10))10-element CuArray{Int64, 1, CUDACore.DeviceMemory}:  1  2  3  4  5  6  7  8  9 10julia> a = CuArray(collect(1:6))6-element CuArray{Int64, 1, CUDACore.DeviceMemory}: 1 2 3 4 5 6julia> b = reshape(a, (2,3))2×3 CuArray{Int64, 2, CUDACore.DeviceMemory}: 1  3  5 2  4  6julia> c = view(a, 2:5)4-element CuArray{Int64, 1, CUDACore.DeviceMemory}: 2 3 4 5

The above contiguous view and reshape have been specialized to return new objects of type CuArray. Other wrappers, such as non-contiguous views or the LinearAlgebra wrappers that will be discussed below, are implemented using their own type (e.g. SubArray or Transpose). This can cause problems, as calling methods with these wrapped objects will not dispatch to specialized CuArray methods anymore. That may result in a call to fallback functionality that performs scalar iteration.

Certain common operations, like broadcast or matrix multiplication, do know how to deal with array wrappers by using the Adapt.jl package. This is still not a complete solution though, e.g. new array wrappers are not covered, and only one level of wrapping is supported. Sometimes the only solution is to materialize the wrapper to a CuArray again.

Kernels assume that arrays in global memory are aligned to the largest power of two, up to 16 bytes, that divides their element size, e.g., to 8 bytes for ComplexF32 and to 16 bytes for ComplexF64 or NTuple{4,Float32} (unions and zero-size types keep their own alignment). As with CUDA C's float2 and float4, this makes it possible to load and store such elements with a single vector instruction. Arrays allocated by CUDA.jl, and contiguous views and reshapes of them, always satisfy this. An array can only start elsewhere when reinterpreting a view that starts at an offset to a larger element type (e.g. reinterpret(ComplexF32, view(x::CuVector{Float32}, 2:9))), or when wrapping a pointer with unsafe_wrap. Such arrays can still be copied, but passing them to a kernel, including broadcasts and other array operations, throws an ArgumentError. Use copy to get an aligned array.

Random numbers

Base's convenience functions for generating random numbers are available in the CUDA module as well:

julia> cuRAND.rand(2)2-element CuArray{Float32, 1, CUDACore.DeviceMemory}: 0.74021935 0.9209938julia> CUDA.randn(Float64, 2, 1)2×1 CuArray{Float64, 2, CUDACore.DeviceMemory}: -0.3893830994647195  1.618410515635752

Behind the scenes, these random numbers come from two different generators: one backed by cuRAND, another by kernels defined in GPUArrays.jl. Operations on these generators are implemented using methods from the Random standard library:

julia> using Randomjulia> a = Random.rand(cuRAND.library_rng(), Float32, 1)1-element CuArray{Float32, 1, CUDACore.DeviceMemory}: 0.74021935julia> a = Random.rand!(cuRAND.gpuarrays_rng(), a)1-element CuArray{Float32, 1, CUDACore.DeviceMemory}: 0.9722412

cuRAND also supports generating lognormal and Poisson-distributed numbers:

julia> cuRAND.rand_logn(Float32, 1, 5; mean=2, stddev=20)1×5 CuArray{Float32, 2, CUDACore.DeviceMemory}: 2567.61  4.256f-6  54.5948  0.00283999  9.81175f22julia> cuRAND.rand_poisson(UInt32, 1, 10; lambda=100)1×10 CuArray{UInt32, 2, CUDACore.DeviceMemory}: 0x00000058  0x00000066  0x00000061  …  0x0000006b  0x0000005f  0x00000069

Note that these custom operations are only supported on a subset of types.

Linear algebra

CUDA's linear algebra functionality from the cuBLAS library is exposed by implementing methods in the LinearAlgebra standard library:

julia> using cuBLASjulia> # enable logging to demonstrate a cuBLAS kernel is used       cuBLAS.cublasLoggerConfigure(1, 0, 1, C_NULL)julia> cuRAND.rand(2,2) * cuRAND.rand(2,2)I! cuBLAS (v10.2) function cublasStatus_t cublasSgemm_v2(cublasContext*, cublasOperation_t, cublasOperation_t, int, int, int, const float*, const float*, int, const float*, int, const float*, float*, int) called2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: 0.295727  0.479395 0.624576  0.557361

Certain operations, like the above matrix-matrix multiplication, also have a native fallback written in Julia for the purpose of working with types that are not supported by cuBLAS:

julia> # enable logging to demonstrate no cuBLAS kernel is used       cuBLAS.cublasLoggerConfigure(1, 0, 1, C_NULL)julia> cuRAND.rand(Int128, 2, 2) * cuRAND.rand(Int128, 2, 2)2×2 CuArray{Int128, 2, CUDACore.DeviceMemory}: -147256259324085278916026657445395486093  -62954140705285875940311066889684981211 -154405209690443624360811355271386638733  -77891631198498491666867579047988353207

Operations that exist in cuBLAS, but are not (yet) covered by high-level constructs in the LinearAlgebra standard library, can be accessed directly from the cuBLAS submodule. Note that you do not need to call the C wrappers directly (e.g. cublasDdot), as many operations have more high-level wrappers available as well (e.g. dot):

julia> x = cuRAND.rand(2)2-element CuArray{Float32, 1, CUDACore.DeviceMemory}: 0.74021935 0.9209938julia> y = cuRAND.rand(2)2-element CuArray{Float32, 1, CUDACore.DeviceMemory}: 0.03902049 0.9689629julia> cuBLAS.dot(2, x, y)0.92129254f0julia> using LinearAlgebrajulia> dot(Array(x), Array(y))0.92129254f0

Solver

LAPACK-like functionality as found in the cuSOLVER library can be accessed through methods in the LinearAlgebra standard library too:

julia> using LinearAlgebra, cuSOLVERjulia> a = cuRAND.rand(2,2)2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: 0.740219  0.0390205 0.920994  0.968963julia> a = a * a'2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: 0.549447  0.719547 0.719547  1.78712julia> cholesky(a)Cholesky{Float32, CuArray{Float32, 2, CUDACore.DeviceMemory}}U factor:2×2 UpperTriangular{Float32, CuArray{Float32, 2, CUDACore.DeviceMemory}}: 0.741247  0.970725  ⋅        0.919137

Other operations are bound to the left-division operator:

julia> a = cuRAND.rand(2,2)2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: 0.740219  0.0390205 0.920994  0.968963julia> b = cuRAND.rand(2,2)2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: 0.925141  0.667319 0.44635   0.109931julia> a \ b2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}:  1.29018    0.942773 -0.765663  -0.782648julia> Array(a) \ Array(b)2×2 Matrix{Float32}:  1.29018    0.942773 -0.765663  -0.782648

Sparse arrays

Sparse array functionality from the cuSPARSE library is mainly available through functionality from the SparseArrays package applied to CuSparseArray objects:

julia> using SparseArrays, cuSPARSEjulia> x = sprand(10,0.2)10-element SparseVector{Float64, Int64} with 5 stored entries:  [2 ]  =  0.538639  [4 ]  =  0.89699  [6 ]  =  0.258478  [7 ]  =  0.338949  [10]  =  0.424742julia> d_x = CuSparseVector(x)10-element CuSparseVector{Float64, Int32} with 5 stored entries:  [2 ]  =  0.538639  [4 ]  =  0.89699  [6 ]  =  0.258478  [7 ]  =  0.338949  [10]  =  0.424742julia> nonzeros(d_x)5-element CuArray{Float64, 1, CUDACore.DeviceMemory}: 0.538639413965653 0.8969897902567084 0.25847781536337067 0.3389490517221738 0.4247416640213063julia> nnz(d_x)5

For 2-D arrays the CuSparseMatrixCSC and CuSparseMatrixCSR can be used.

Non-integrated functionality can be access directly in the cuSPARSE submodule again.

FFTs

Functionality from cuFFT is integrated with the interfaces from the AbstractFFTs.jl package:

julia> a = cuRAND.rand(2,2)2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}: 0.740219  0.0390205 0.920994  0.968963julia> using cuFFTjulia> fft(a)2×2 CuArray{ComplexF32, 2, CUDACore.DeviceMemory}:   2.6692+0.0im   0.65323+0.0im -1.11072+0.0im  0.749168+0.0im

Planner keyword arguments like FFTW's flags and timelimit are accepted so that generic code can pass them, but they have no effect: cuFFT has no equivalent options.