From 7b546c862a1fc38d4a0e3181714245e93e5ad6a5 Mon Sep 17 00:00:00 2001 From: Tim Besard Date: Wed, 7 Oct 2026 21:18:53 +0200 Subject: [PATCH 1/2] Support GPUArrays 12 GPUArrays 12 implements Base's reductions on AcceleratedKernels and no longer calls the GPUArrays.mapreducedim! hook, so OpenCL.jl's reduction kernel is removed. GPUArrays.default_rng is gone too; the per-device RNG backing rand!, randn! and seed! now lives in OpenCL.jl and constructs GPUArrays' stateless RNG directly. --- Project.toml | 2 +- src/OpenCL.jl | 1 - src/gpuarrays.jl | 16 ---- src/mapreduce.jl | 182 ------------------------------------------ src/random.jl | 11 ++- test/array.jl | 6 -- test/device/random.jl | 2 +- 7 files changed, 12 insertions(+), 208 deletions(-) delete mode 100644 src/mapreduce.jl diff --git a/Project.toml b/Project.toml index c6768800..0c2328e5 100644 --- a/Project.toml +++ b/Project.toml @@ -32,7 +32,7 @@ SPIRVIntrinsics = {path = "lib/intrinsics"} [compat] Adapt = "4" -GPUArrays = "11.2.1" +GPUArrays = "12" GPUCompiler = "2.13" GPUToolbox = "3.3.2" KernelAbstractions = "0.9.38" diff --git a/src/OpenCL.jl b/src/OpenCL.jl index b0d37731..e0dae368 100644 --- a/src/OpenCL.jl +++ b/src/OpenCL.jl @@ -45,7 +45,6 @@ include("compiler/precompile.jl") # integrations and specialized functionality include("util.jl") include("broadcast.jl") -include("mapreduce.jl") include("gpuarrays.jl") include("random.jl") diff --git a/src/gpuarrays.jl b/src/gpuarrays.jl index 5bf53688..2051e538 100644 --- a/src/gpuarrays.jl +++ b/src/gpuarrays.jl @@ -5,19 +5,3 @@ function GPUArrays.derive(::Type{T}, a::CLArray, dims::Dims{N}, offset::Int) whe offset = a.offset + offset * sizeof(T) CLArray{T,N}(ref, dims; offset) end - -const GLOBAL_RNGs = Dict{cl.Device,GPUArrays.RNG}() -const global_rngs_lock = ReentrantLock() - -function GPUArrays.default_rng(::Type{<:CLArray}) - dev = cl.device() - return Base.@lock global_rngs_lock begin - get!(GLOBAL_RNGs, dev) do - N = dev.max_work_group_size - state = CLArray{NTuple{4, UInt32}}(undef, N) - rng = GPUArrays.RNG(state) - Random.seed!(rng) - rng - end - end -end diff --git a/src/mapreduce.jl b/src/mapreduce.jl deleted file mode 100644 index 37172162..00000000 --- a/src/mapreduce.jl +++ /dev/null @@ -1,182 +0,0 @@ -# TODO -# - serial version for lower latency -# - group-stride loop to delay need for second kernel launch -# - let the driver choose the local size - -# Reduce a value across a group, using local memory for communication -@inline function reduce_group(op, val::T, neutral, ::Val{maxitems}) where {T, maxitems} - items = get_local_size() - item = get_local_id() - - # local mem for a complete reduction - shared = CLLocalArray(T, (maxitems,)) - @inbounds shared[item] = val - - # perform a reduction - d = 1 - while d < items - work_group_barrier(LOCAL_MEM_FENCE) - index = 2 * d * (item-1) + 1 - @inbounds if index <= items - other_val = if index + d <= items - shared[index+d] - else - neutral - end - shared[index] = op(shared[index], other_val) - end - d *= 2 - end - - # load the final value on the first item - if item == 1 - val = @inbounds shared[item] - end - - return val -end - -Base.@propagate_inbounds _map_getindex(args::Tuple, I) = ((args[1][I]), _map_getindex(Base.tail(args), I)...) -Base.@propagate_inbounds _map_getindex(args::Tuple{Any}, I) = ((args[1][I]),) -Base.@propagate_inbounds _map_getindex(args::Tuple{}, I) = () - -# Reduce an array across the grid. All elements to be processed can be addressed by the -# product of the two iterators `Rreduce` and `Rother`, where the latter iterator will have -# singleton entries for the dimensions that should be reduced (and vice versa). -function partial_mapreduce_device(f, op, neutral, maxitems, Rreduce, Rother, R, As...) - # decompose the 1D hardware indices into separate ones for reduction (across items - # and possibly groups if it doesn't fit) and other elements (remaining groups) - localIdx_reduce = get_local_id() - localDim_reduce = get_local_size() - groupIdx_reduce, groupIdx_other = fldmod1(get_group_id(), length(Rother)) - groupDim_reduce = get_num_groups() ÷ length(Rother) - - # group-based indexing into the values outside of the reduction dimension - # (that means we can safely synchronize items within this group) - iother = groupIdx_other - @inbounds if iother <= length(Rother) - Iother = Rother[iother] - - # load the neutral value - Iout = CartesianIndex(Tuple(Iother)..., groupIdx_reduce) - neutral = if neutral === nothing - R[Iout] - else - neutral - end - - val = op(neutral, neutral) - - # reduce serially across chunks of input vector that don't fit in a group - ireduce = localIdx_reduce + (groupIdx_reduce - 1) * localDim_reduce - while ireduce <= length(Rreduce) - Ireduce = Rreduce[ireduce] - J = max(Iother, Ireduce) - val = op(val, f(_map_getindex(As, J)...)) - ireduce += localDim_reduce * groupDim_reduce - end - - val = reduce_group(op, val, neutral, maxitems) - - # write back to memory - if localIdx_reduce == 1 - R[Iout] = val - end - end - - return -end - -function GPUArrays.mapreducedim!( - f::F, op::OP, R::WrappedCLArray{T}, - A::Union{AbstractArray,Broadcast.Broadcasted}; - init=nothing) where {F, OP, T} - Base.check_reducedims(R, A) - length(A) == 0 && return R # isempty(::Broadcasted) iterates - - R_old = R - # add singleton dimensions to the output container, if needed - if ndims(R) < ndims(A) - dims = Base.fill_to_length(size(R), 1, Val(ndims(A))) - R = reshape(R, dims) - end - - # iteration domain, split in two: one part covers the dimensions that should - # be reduced, and the other covers the rest. combining both covers all values. - Rall = CartesianIndices(axes(A)) - Rother = CartesianIndices(axes(R)) - Rreduce = CartesianIndices(ifelse.(axes(A) .== axes(R), Ref(Base.OneTo(1)), axes(A))) - # NOTE: we hard-code `OneTo` (`first.(axes(A))` would work too) or we get a - # CartesianIndices object with UnitRanges that behave badly on the GPU. - @assert length(Rall) == length(Rother) * length(Rreduce) - - # allocate an additional, empty dimension to write the reduced value to. - # this does not affect the actual location in memory of the final values, - # but allows us to write a generalized kernel supporting partial reductions. - R′ = reshape(R, (size(R)..., 1)) - - # how many items do we want? - # - # items in a group work together to reduce values across the reduction dimensions; - # we want as many as possible to improve algorithm efficiency and execution occupancy. - wanted_items = length(Rreduce) - function compute_items(max_items) - if wanted_items > max_items - max_items - else - wanted_items - end - end - - # how many items can we launch? - # - # we might not be able to launch all those items to reduce each slice in one go. - # that's why each items also loops across their inputs, processing multiple values - # so that we can span the entire reduction dimension using a single item group. - - # group size is restricted by local memory - max_lmem_elements = cl.device().local_mem_size ÷ sizeof(T) - max_items = min(cl.device().max_work_group_size, - compute_items(max_lmem_elements ÷ 2)) - # TODO: dynamic local memory to avoid two compilations - - # let the driver suggest a group size - args = (f, op, init, Val(max_items), Rreduce, Rother, R′, A) - kernel_args = kernel_convert.(args) - kernel_tt = Tuple{Core.Typeof.(kernel_args)...} - kernel = clfunction(partial_mapreduce_device, kernel_tt) - wg_info = cl.work_group_info(kernel.fun, cl.device()) - reduce_items = compute_items(wg_info.size) - - # how many groups should we launch? - # - # even though we can always reduce each slice in a single item group, that may not be - # optimal as it might not saturate the GPU. we already launch some groups to process - # independent dimensions in parallel; pad that number to ensure full occupancy. - other_groups = length(Rother) - reduce_groups = cld(length(Rreduce), reduce_items) - - # determine the launch configuration - local_size = reduce_items - global_size = reduce_items*reduce_groups*other_groups - - # perform the actual reduction - if reduce_groups == 1 - # we can cover the dimensions to reduce using a single group - @opencl local_size global_size partial_mapreduce_device( - f, op, init, Val(local_size), Rreduce, Rother, R′, A) - else - # we need multiple steps to cover all values to reduce - partial = similar(R, (size(R)..., reduce_groups)) - if init === nothing - # without an explicit initializer we need to copy from the output container - partial .= R - end - @opencl local_size global_size partial_mapreduce_device( - f, op, init, Val(local_size), Rreduce, Rother, partial, A) - - GPUArrays.mapreducedim!(identity, op, R′, partial; init=init) - end - - return R_old -end diff --git a/src/random.jl b/src/random.jl index 655470d3..c5904837 100644 --- a/src/random.jl +++ b/src/random.jl @@ -1,6 +1,15 @@ using Random -gpuarrays_rng() = GPUArrays.default_rng(CLArray) +const GLOBAL_RNGs = Dict{cl.Device,GPUArrays.RNG{CLArray}}() +const global_rngs_lock = ReentrantLock() + +# one RNG per device, used by the RNG-less `rand!`/`randn!` methods and `seed!` +function gpuarrays_rng() + dev = cl.device() + return Base.@lock global_rngs_lock begin + get!(() -> GPUArrays.RNG{CLArray}(), GLOBAL_RNGs, dev) + end +end # GPUArrays in-place Random.rand!(A::WrappedCLArray) = Random.rand!(gpuarrays_rng(), A) diff --git a/test/array.jl b/test/array.jl index 1195f526..e9acd008 100644 --- a/test/array.jl +++ b/test/array.jl @@ -239,12 +239,6 @@ end @test length(b) == 1 end -@testset "mapreducedim! returning same type" begin - R = transpose(OpenCL.zeros(Float32, 2, 3)) - A = CLArray(rand(Float32, 3, 2, 10)) - @test @inferred(OpenCL.GPUArrays.mapreducedim!(identity, +, R, A)) === R -end - # finalizers run in no particular order, e.g. at exit (JuliaGPU/OpenCL.jl#279), so memory # has to remain freeable after the queue and context it was allocated with are finalized @testset "freeing after finalizing its queue and context" begin diff --git a/test/device/random.jl b/test/device/random.jl index b2c2c6b5..33e775f0 100644 --- a/test/device/random.jl +++ b/test/device/random.jl @@ -167,7 +167,7 @@ if Float16 in GPUArraysTestSuite.supported_eltypes(CLArray) end @testset "randn!(Complex{Float16}) is finite" begin - rng = OpenCL.GPUArrays.default_rng(CLArray) + rng = OpenCL.GPUArrays.RNG{CLArray}() Random.seed!(rng, 1) A = CLArray{Complex{Float16}}(undef, 4096) randn!(rng, A) From fa40095078d278de3d3debeba424e12fb564aa27 Mon Sep 17 00:00:00 2001 From: Tim Besard Date: Thu, 8 Oct 2026 10:25:46 +0200 Subject: [PATCH 2/2] Detect aliasing by allocation and byte range, like CUDA.jl Base.dataids returned the array's own pointer, which missed overlaps between a contiguous view (a CLArray at an offset) and a wrapped array of the same memory. As in CUDA.jl, identify the allocation in dataids. Read its address from the memory object instead of using pointer(A), which takes ownership of the memory, and fall back to the handle of buffers without an address. Empty arrays alias nothing, which Julia 1.10's mightalias does not check, so they return no dataids. --- src/array.jl | 31 +++++++++++++++++++++++++++---- test/array.jl | 48 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+), 4 deletions(-) diff --git a/src/array.jl b/src/array.jl index 14ca3d31..b764f038 100644 --- a/src/array.jl +++ b/src/array.jl @@ -76,14 +76,37 @@ GPUArrays.storage(a::CLArray) = a.data ## alias detection -Base.dataids(A::CLArray) = (UInt(pointer(A)),) +# The memory `A` lives in: its address, or the handle of a buffer without one (which only +# identifies that buffer), and whether it is an address. Not using `pointer(A)`, which takes +# ownership of the memory for the current queue. +function memory_base(A::CLArray) + mem = A.data[].mem + if mem isa cl.Buffer && mem.ptr === nothing + return UInt(mem.id), false + end + return UInt(pointer(mem)), true +end + +# Identify the underlying memory, not just where this array starts in it: derived arrays +# (contiguous views, reshapes, reinterprets) have an offset into their parent's memory, and +# Base compares `dataids` whenever one side is wrapped (e.g., a `SubArray`). The start +# address is included too, to match host memory wrapped from it. Empty arrays alias nothing, +# which Julia 1.10's `mightalias` does not check. +function Base.dataids(A::CLArray) + isempty(A) && return () + base, _ = memory_base(A) + return (base, base + A.offset) +end Base.unaliascopy(A::CLArray) = copy(A) function Base.mightalias(A::CLArray, B::CLArray) - rA = pointer(A):(pointer(A) + sizeof(A)) - rB = pointer(B):(pointer(B) + sizeof(B)) - return first(rA) <= first(rB) < last(rA) || first(rB) <= first(rA) < last(rB) + (baseA, addressA), (baseB, addressB) = memory_base(A), memory_base(B) + # (offsets into buffers without an address only compare within the same buffer) + addressA && addressB || baseA == baseB || return false + startA = baseA + A.offset + startB = baseB + B.offset + return startA <= startB < startA + sizeof(A) || startB <= startA < startB + sizeof(B) end diff --git a/test/array.jl b/test/array.jl index e9acd008..f68048cf 100644 --- a/test/array.jl +++ b/test/array.jl @@ -67,6 +67,54 @@ end r = reinterpret(Int64, v) # Int64 = 8 bytes; 4 is not a multiple of 8 @test Array(r) == reinterpret(Int64, @view Array(a)[2:7]) end + +@testset "aliasing" begin + x = CLArray([1, 2]) + y = view(x, 2:2) + @test Base.mightalias(x, x) + @test Base.mightalias(x, y) + z = view(x, 1:1) + @test Base.mightalias(x, z) + @test !Base.mightalias(y, z) + + a = copy(y)::typeof(x) + @test !Base.mightalias(x, a) + b = Base.unaliascopy(y)::typeof(y) + @test !Base.mightalias(x, b) + + # contiguous views are CLArrays with an offset into the parent's memory, + # which should still alias wrapped arrays (like SubArrays) of that memory + x = CLArray(1:16) + @test Base.mightalias(view(x, 2:16), view(x, 15:-1:1)) + @test Base.mightalias(view(x, 1:2:15), view(x, 2:9)) + @test Base.mightalias(view(x, 2:16), view(reinterpret(Int32, x), 1:2:31)) + @test !Base.mightalias(view(x, 2:16), view(CLArray(1:16), 15:-1:1)) + + # so in-place broadcasts between them should make a copy first + n = 2^20 + x = CLArray{Float32}(1:n) + view(x, 2:n) .= view(x, n-1:-1:1) + @test Array(x) == [1; n-1:-1:1] + + # host memory wrapped from a view's address + if OpenCL.system_memory_type() !== nothing + h = Float32.(1:16) + x = unsafe_wrap(CLArray, h) + y = view(x, 2:16) + z = unsafe_wrap(CLArray, pointer(h, 2), size(y)) + @test Base.mightalias(y, z) + @test Base.mightalias(y, view(z, 15:-1:1)) + end + + # empty arrays alias nothing, also on Julia 1.10 + @test !Base.mightalias(CLArray(Int[]), CLArray(Float32[])) + @test !Base.mightalias(view(CLArray(zeros(Float32, 2, 0)), 1:1, :), CLArray(Int[])) + + # disjoint parts of one array may be each other's source and destination + x = CLArray(collect(1:10)) + @test Array(sum!(view(x, 1:1), view(x, 2:10))) == [54] + @test Array(cumsum!(view(x, 1:5), view(x, 6:10))) == cumsum(6:10) +end # TODO: Look into how to port the @sync if cl.USMBackend() in cl.supported_memory_backends(cl.device())