Custom CUDA kernel tests

README.md

@@ -61,7 +61,7 @@ See [Kernel Fusion](https://julialegate.github.io/cuNumeric.jl/dev/perf/kernel_f
 
 ### The `@accelerate` macro
 
-`@accelerate` optimizes straight-line array code by fusing eligible CUDA broadcasts and releasing dead temporary `NDArray`s. See [The `@accelerate` Macro](https://julialegate.github.io/cuNumeric.jl/dev/perf/reduce_allocations) for usage guidance.
+`@accelerate` fuses eligible GPU broadcasts within and across statements, then releases materialized temporary `NDArray`s after their last use on CPU or GPU. See [The `@accelerate` Macro](https://julialegate.github.io/cuNumeric.jl/dev/perf/reduce_allocations) for usage guidance.
 
 ### Benchmarks
 

benchmark/benchmarks.toml

@@ -36,16 +36,40 @@ T = "Float32"
 gpus = [1, 2, 4, 8]
 cpus = 16
 fusion = [true, false]
-N = [2000, 2832, 4000, 5656]
-M = [2000, 2832, 4000, 5656]
+N = [24000, 33944, 48000, 67888]
+M = [24000, 33944, 48000, 67888]
 
-[[grayscott_accelerated]]
+[[grayscott_function_accelerated]]
 T = "Float32"
 gpus = [1, 2, 4, 8]
 cpus = 16
 fusion = [true, false]
-N = [2000, 2832, 4000, 5656]
-M = [2000, 2832, 4000, 5656]
+N = [24000, 33944, 48000, 67888]
+M = [24000, 33944, 48000, 67888]
+
+[[grayscott_begin_accelerated]]
+T = "Float32"
+gpus = [1, 2, 4, 8]
+cpus = 16
+fusion = [true, false]
+N = [24000, 33944, 48000, 67888]
+M = [24000, 33944, 48000, 67888]
+
+[[grayscott_let_accelerated]]
+T = "Float32"
+gpus = [1, 2, 4, 8]
+cpus = 16
+fusion = [true, false]
+N = [24000, 33944, 48000, 67888]
+M = [24000, 33944, 48000, 67888]
+
+[[grayscott_expression_accelerated]]
+T = "Float32"
+gpus = [1, 2, 4, 8]
+cpus = 16
+fusion = [true, false]
+N = [24000, 33944, 48000, 67888]
+M = [24000, 33944, 48000, 67888]
 
 #################################
 #   Monte-Carlo Integration     #

benchmark/src/benchmarks/grayscott.jl

@@ -92,63 +92,64 @@ function check_benchmark_correctness(
     return (u_ok && v_ok) ? "pass" : "fail"
 end
 
-# Variant description:
-# baseline: as written
-# accelerated: step wrapped in @accelerate
-let body = quote
-        # currently we don't have NDArray^x working yet. every operator is dotted
-        # so each rhs fuses into a single broadcast kernel rather than shattering
-        # into bare +/-/* binary tasks.
-        F_u = (
-            (
-                .-u[2:(end - 1), 2:(end - 1)] .*
-                (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
-            ) .+ args.f .* (1.0f0 .- u[2:(end - 1), 2:(end - 1)])
-        )
-        F_v = (
-            (
-                u[2:(end - 1), 2:(end - 1)] .*
-                (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
-            ) .- (args.f + args.k) .* v[2:(end - 1), 2:(end - 1)]
-        )
-        # 2-D Laplacian via slicing, excluding boundaries
-        u_lap = (
-            (
-                u[3:end, 2:(end - 1)] .- 2 .* u[2:(end - 1), 2:(end - 1)] .+
-                u[1:(end - 2), 2:(end - 1)]
-            ) ./ args.dx^2 .+
-            (
-                u[2:(end - 1), 3:end] .- 2 .* u[2:(end - 1), 2:(end - 1)] .+
-                u[2:(end - 1), 1:(end - 2)]
-            ) ./ args.dx^2
-        )
-        v_lap = (
-            (
-                v[3:end, 2:(end - 1)] .- 2 .* v[2:(end - 1), 2:(end - 1)] .+
-                v[1:(end - 2), 2:(end - 1)]
-            ) ./ args.dx^2 .+
-            (
-                v[2:(end - 1), 3:end] .- 2 .* v[2:(end - 1), 2:(end - 1)] .+
-                v[2:(end - 1), 1:(end - 2)]
-            ) ./ args.dx^2
-        )
-
-        # Forward-Euler step for all interior points
-        u_new[2:(end - 1), 2:(end - 1)] =
-            ((args.c_u .* u_lap) .+ F_u) .* args.dt .+ u[2:(end - 1), 2:(end - 1)]
-        v_new[2:(end - 1), 2:(end - 1)] =
-            ((args.c_v .* v_lap) .+ F_v) .* args.dt .+ v[2:(end - 1), 2:(end - 1)]
-
-        # Periodic boundary conditions
-        u_new[:, 1] = u[:, end - 1]
-        u_new[:, end] = u[:, 2]
-        u_new[1, :] = u[end - 1, :]
-        u_new[end, :] = u[2, :]
-        v_new[:, 1] = v[:, end - 1]
-        v_new[:, end] = v[:, 2]
-        v_new[1, :] = v[end - 1, :]
-        v_new[end, :] = v[2, :]
-    end
+# Shared syntax tree keeps every Gray-Scott variant on the exact same workload.
+const GRAYSCOTT_STEP_BODY = quote
+    # currently we don't have NDArray^x working yet. every operator is dotted
+    # so each rhs fuses into a single broadcast kernel rather than shattering
+    # into bare +/-/* binary tasks.
+    F_u = (
+        (
+            .-u[2:(end - 1), 2:(end - 1)] .*
+            (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
+        ) .+ args.f .* (1.0f0 .- u[2:(end - 1), 2:(end - 1)])
+    )
+    F_v = (
+        (
+            u[2:(end - 1), 2:(end - 1)] .*
+            (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
+        ) .- (args.f + args.k) .* v[2:(end - 1), 2:(end - 1)]
+    )
+    # 2-D Laplacian via slicing, excluding boundaries
+    u_lap = (
+        (
+            u[3:end, 2:(end - 1)] .- 2 .* u[2:(end - 1), 2:(end - 1)] .+
+            u[1:(end - 2), 2:(end - 1)]
+        ) ./ args.dx^2 .+
+        (
+            u[2:(end - 1), 3:end] .- 2 .* u[2:(end - 1), 2:(end - 1)] .+
+            u[2:(end - 1), 1:(end - 2)]
+        ) ./ args.dx^2
+    )
+    v_lap = (
+        (
+            v[3:end, 2:(end - 1)] .- 2 .* v[2:(end - 1), 2:(end - 1)] .+
+            v[1:(end - 2), 2:(end - 1)]
+        ) ./ args.dx^2 .+
+        (
+            v[2:(end - 1), 3:end] .- 2 .* v[2:(end - 1), 2:(end - 1)] .+
+            v[2:(end - 1), 1:(end - 2)]
+        ) ./ args.dx^2
+    )
+
+    # Forward-Euler step for all interior points
+    u_new[2:(end - 1), 2:(end - 1)] =
+        ((args.c_u .* u_lap) .+ F_u) .* args.dt .+ u[2:(end - 1), 2:(end - 1)]
+    v_new[2:(end - 1), 2:(end - 1)] =
+        ((args.c_v .* v_lap) .+ F_v) .* args.dt .+ v[2:(end - 1), 2:(end - 1)]
+
+    # Periodic boundary conditions
+    u_new[:, 1] = u[:, end - 1]
+    u_new[:, end] = u[:, 2]
+    u_new[1, :] = u[end - 1, :]
+    u_new[end, :] = u[2, :]
+    v_new[:, 1] = v[:, end - 1]
+    v_new[:, end] = v[:, 2]
+    v_new[1, :] = v[end - 1, :]
+    v_new[end, :] = v[2, :]
+end
+
+# Original baseline and recommended function-form benchmark.
+let body = deepcopy(GRAYSCOTT_STEP_BODY)
     @eval _gs_step!(b::GrayScottBaseline, u, v, u_new, v_new, args::GSParams) = $body
     @eval @accelerate function _gs_step!(
         b::GrayScottAccelerated, u, v, u_new, v_new, args::GSParams

benchmark/src/benchmarks/grayscott_accelerate_forms.jl

@@ -0,0 +1,96 @@
+# Compare the four scope contracts of `@accelerate` on one shared Gray-Scott step.
+# Each type has a distinct result name so benchmark runs produce separate CSVs.
+
+abstract type AbstractGrayScottAccelerateForm{T} <: AbstractGrayScott{T} end
+
+Base.@kwdef struct GrayScottFunctionAccelerated{T} <:
+                   AbstractGrayScottAccelerateForm{T}
+    N::Int
+    M::Int
+end
+
+Base.@kwdef struct GrayScottBeginAccelerated{T} <: AbstractGrayScottAccelerateForm{T}
+    N::Int
+    M::Int
+end
+
+Base.@kwdef struct GrayScottLetAccelerated{T} <: AbstractGrayScottAccelerateForm{T}
+    N::Int
+    M::Int
+end
+
+Base.@kwdef struct GrayScottExpressionAccelerated{T} <:
+                   AbstractGrayScottAccelerateForm{T}
+    N::Int
+    M::Int
+end
+
+name(::GrayScottFunctionAccelerated) = "grayscott_function_accelerated"
+name(::GrayScottBeginAccelerated) = "grayscott_begin_accelerated"
+name(::GrayScottLetAccelerated) = "grayscott_let_accelerated"
+name(::GrayScottExpressionAccelerated) = "grayscott_expression_accelerated"
+
+# Function form is the reusable default: arguments and the return value survive,
+# while non-returned locals may fuse across statements or die after their last use.
+let body = deepcopy(GRAYSCOTT_STEP_BODY)
+    @eval @accelerate function _gs_step!(
+        b::GrayScottFunctionAccelerated, u, v, u_new, v_new, args::GSParams
+    )
+        $body
+    end
+end
+
+# `begin` adds no scope. Every named local remains visible, so it measures the
+# multi-output/materialized path rather than eliminating named intermediates.
+let body = deepcopy(GRAYSCOTT_STEP_BODY)
+    @eval function _gs_step!(
+        b::GrayScottBeginAccelerated, u, v, u_new, v_new, args::GSParams
+    )
+        @accelerate begin
+            $body
+        end
+    end
+end
+
+# `let` is a hard one-off scope. Only its result escapes, allowing aggressive
+# inter-statement fusion and last-use cleanup for all other local temporaries.
+let body = deepcopy(GRAYSCOTT_STEP_BODY)
+    @eval function _gs_step!(
+        b::GrayScottLetAccelerated, u, v, u_new, v_new, args::GSParams
+    )
+        @accelerate let
+            $body
+        end
+    end
+end
+
+# Expression form has no multi-statement scope. Accelerating each RHS preserves
+# fusion inside that expression but deliberately materializes statement results,
+# isolating intra-expression fusion from the inter-statement rewrite cases above.
+function accelerate_grayscott_rhs(body::Expr)
+    statements = Any[]
+    for statement in body.args
+        if statement isa LineNumberNode
+            push!(statements, statement)
+        elseif statement isa Expr && statement.head === :(=)
+            lhs, rhs = statement.args
+            push!(statements, :($lhs = @accelerate $rhs))
+        else
+            error("Gray-Scott expression benchmark expected assignments; got $(repr(statement))")
+        end
+    end
+    return Expr(:block, statements...)
+end
+
+let body = accelerate_grayscott_rhs(deepcopy(GRAYSCOTT_STEP_BODY))
+    @eval function _gs_step!(
+        b::GrayScottExpressionAccelerated, u, v, u_new, v_new, args::GSParams
+    )
+        $body
+    end
+end
+
+register_benchmark("grayscott_function_accelerated", GrayScottFunctionAccelerated)
+register_benchmark("grayscott_begin_accelerated", GrayScottBeginAccelerated)
+register_benchmark("grayscott_let_accelerated", GrayScottLetAccelerated)
+register_benchmark("grayscott_expression_accelerated", GrayScottExpressionAccelerated)

docs/src/index.md

@@ -0,0 +1,96 @@
+```@raw html
+<h1>
+  <img src="./assets/logo.png" alt="cuNumeric.jl" width="50">
+  cuNumeric.jl
+</h1>
+```
+
+[![Documentation dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://julialegate.github.io/cuNumeric.jl/dev) [![codecov](https://codecov.io/github/julialegate/cuNumeric.jl/branch/main/graph/badge.svg)](https://app.codecov.io/github/JuliaLegate/cuNumeric.jl) [![License: MIT](https://img.shields.io/badge/License-MIT-green.svg)](https://opensource.org/licenses/MIT)
+
+cuNumeric.jl wraps and extends the [cuPyNumeric](https://github.com/nv-legate/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](https://github.com/juliagpu/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:
+
+```julia
+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](https://julialegate.github.io/cuNumeric.jl/dev/install).
+
+```julia
+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](https://julialegate.github.io/cuNumeric.jl/dev/configuration/hardware).
+
+### How `NDArray`s 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 partitioned over 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`. Mutations through an `NDArray` slice are visible through other aliases of the same data.
+
+**Reductions return arrays, not Julia scalars.** Reductions such as `sum(A)` produce a **0D or 1D** `NDArray` (axis reductions produce a lower-rank `NDArray`), rather than a bare `Float64` / `Float32`. That keeps the Legate task graph asynchronous instead of forcing synchronization to communicate with the Julia runtime. When you need a plain Julia number, call `unwrap`:
+
+```julia
+s = sum(A)          # NDArray{T,0}
+x = unwrap(s)       # T, e.g. Float32
+```
+
+**The Legate runtime builds a DAG asynchronously.** Calling `cuNumeric.zeros` or `A .+ B` records work into Legate's task graph rather than blocking until every GPU kernel finishes. Results are materialized when you need them (for example `println`, `unwrap`, or converting with `Array(A)`). Hiding latency enables performant code.
+
+For API details see [Initialization](https://julialegate.github.io/cuNumeric.jl/dev/api_initialization) and [NDArray Reference](https://julialegate.github.io/cuNumeric.jl/dev/api). For common performance pitfalls, see [Patterns to Avoid](https://julialegate.github.io/cuNumeric.jl/dev/perf/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.
+
+```julia
+y .= @. -a + b * c
+```
+
+See [Kernel Fusion](https://julialegate.github.io/cuNumeric.jl/dev/perf/kernel_fusion) and [Debugging](https://julialegate.github.io/cuNumeric.jl/dev/debugging) for controls and diagnostics.
+
+### The `@accelerate` macro
+
+`@accelerate` fuses eligible GPU broadcasts within and across statements, then releases materialized temporary `NDArray`s after their last use on CPU or GPU. See [The `@accelerate` Macro](https://julialegate.github.io/cuNumeric.jl/dev/perf/reduce_allocations) for usage guidance.
+
+### Benchmarks
+
+Results and reproduction instructions live under [Benchmark Results](https://julialegate.github.io/cuNumeric.jl/dev/benchmarks/results) and [How to Benchmark](https://julialegate.github.io/cuNumeric.jl/dev/benchmarks/howto).
+
+### Try an example
+
+```julia
+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: $(estimate)")
+```
+More worked examples (initialization, Gray-Scott, …) are in the documentation sidebar under **Examples**.
+
+### Known Limitations
+
+- There is no support for `Float16` or `ComplexF16`