cuNumeric.jl wraps and extends the cuPyNumeric library from NVIDIA to bring distributed array computing on GPUs and CPUs to Julia. The central type is NDArray, which behaves like Julia's Array or the CuArray from CUDA.jl, but executes across multiple GPUs/CPUs. We implement array-level operations on NDArray which can be composed into larger programs without the need for explicit MPI calls or writing CUDA kernels.

cuNumeric.jl requires x86 Linux, an NVIDIA GPU, and Julia >= 1.10. If ARM support is of interest open an issue.

Quick Start

cuNumeric.jl can be installed with the Julia package manager. Activate your preferred environment and then from the Julia REPL run:

using Pkg
Pkg.add(url = "https://github.com/JuliaLegate/cuNumeric.jl", rev = "main")

The first installation can take a while because it includes several large dependencies, such as the CUDA SDK. To use a local cupynumeric build, see Build Modes.

using cuNumeric
cuNumeric.versioninfo()

Warning

Starting more than one instance of cuNumeric.jl can lead to a hard-crash. The default hardware configuration reserves all available resources.

For more details, see Hardware.

How NDArrays work

The semantics of NDArray closely mirror Julia's Array, and in most cases it is a drop-in replacement. You can use the same constructors (i.e., zeros, ones, rand), broadcasting, slicing, and linear algebra. Under the hood a few details differ from Base, and knowing them can help you write fast code.

Data may live across many devices. An NDArray is a logical array whose physical buffers can be stored across multiple GPUs and CPUs by the Legate runtime. You write ordinary array code and Legate decides where the data lives and how/when it is communicated between devices. As a result, elementwise indexing (i.e. arr[1]) is slow (and is prevented by default). Scalar indexing like this forces synchronization and blocks other tasks from executing.

Slices are views. Indexing an NDArray with ranges returns a view onto the same store, not a copy. That differs from Base Julia, where A[1:n] allocates a new Array. Mutating an NDArray slice mutates the parent and all other aliases of the underlying data.

Scalar reductions return device scalars. Full reductions such as sum(A), dot(x,y), and norm(x) return numeric wrappers (CNFloat, CNInt, CNUInt, CNBool, or CNComplex) backed by 0D NDArrays. Each subtypes the corresponding abstract numeric category. Dimension-preserving reductions still return NDArrays. Use fetch for explicit host extraction:

s = sum(A)          # CNScalar
x = fetch(s)        # native Julia scalar, e.g. Float32

Use allowfetch() do ... end or @allowfetch to permit host comparisons and numeric conversions of device scalars. Permission is task-local and disabled by default; arithmetic remains on the backend. See Device scalars.

0D arrays support scalar-shaped arithmetic. Use +, -, *, /, and ^ with two 0D NDArrays or with a 0D NDArray and a Julia number. Results stay as 0D NDArrays on the backend, following the existing promotion policy.

The Legate runtime builds a DAG asynchronously. Calling cuNumeric.zeros or A .+ B records work into a task graph rather than blocking until every GPU kernel finishes. Results are materialized when you need them [for example println, fetch, or communicating with the Julia runtime (i.e., Array(A))]. Hiding latency enables performant code.

For API details see Initialization and NDArray Reference. For common performance pitfalls, see Patterns to Avoid.

Kernel Fusion

Nested broadcast expressions fuse into a single kernel by default when on GPU. Prefer @. for multi-op elementwise code so every operator is dotted and the expression stays completely fused. Even just forgetting the . on unary negation (i.e., y .= -a .+ b .* c) will result in unfused code. Use the following pattern instead.

y .= @. -a + b * c

See Kernel Fusion and Debugging for controls and diagnostics.

The @accelerate macro

@accelerate fuses eligible GPU broadcasts within and across statements, then releases materialized temporary NDArrays after their last use on CPU or GPU. See The @accelerate Macro for usage guidance.

Benchmarks

Results and reproduction instructions live under Benchmark Results and Running the Benchmarks. Run ./scripts/get_benchmarks.sh to fetch the suite into benchmark/.

TensorOperations.jl Integration

We implement a package extension for TensorOperations.jl to enable clean syntax and optimal contraction order for tensor contractions. Simply load TensorOperations alongside cuNumeric to take advantage! Some useful macros if you are unfamiliar are @tensor, @tensoropt and @notensor.

using TensorOperations
using cuNumeric

α = randn() # prefer a Julia Number when you already have one
A = cuNumeric.randn(5, 5, 5, 5, 5, 5)
B = cuNumeric.randn(5, 5, 5)
C = cuNumeric.randn(5, 5, 5)
D = cuNumeric.zeros(5, 5, 5)

@tensor begin
    D[a, b, c] = A[a, e, f, c, f, g] * B[g, b, e] + α * C[c, a, b]
    E[a, b, c] := A[a, e, f, c, f, g] * B[g, b, e] + α * C[c, a, b]
end

Try an example

using cuNumeric

integrand(x) = exp(-x^2)

@accelerate function monte_carlo(N, x_max)
    Ω = 2 * x_max
    raw_samples = cuNumeric.rand(N)
    samples = @. Ω * raw_samples - x_max
    return (Ω / N) * sum(integrand.(samples))
end

N = 1_000_000
x_max = 10.0f0
estimate = monte_carlo(N, x_max)

println("Monte-Carlo Estimate: $(fetch(estimate))") # Should be ~ sqrt(pi)

More worked examples (initialization, Gray-Scott, …) are in the documentation sidebar under Examples.

Known Limitations & Bugs

• There is no support for Float16 or ComplexF16 or Complex{<:Integer}
• Stochastic bug in FFT can cause segfault
• Unfused broadcasting of a Vector reshaped to (1, n) will not mimic Julia semantics