From 7867774064d7a77c517444400c5a710fda0ec02d Mon Sep 17 00:00:00 2001 From: Tim Besard Date: Wed, 7 Oct 2026 07:17:33 +0200 Subject: [PATCH 1/3] Construct the GPUArrays RNG directly GPUArrays 12 removes the default_rng interface and the RNG constructor that took a state array; its RNG is now stateless and built as RNG{oneArray}(). --- src/gpuarrays.jl | 15 --------------- src/oneAPI.jl | 1 - src/random.jl | 9 ++++++++- 3 files changed, 8 insertions(+), 17 deletions(-) delete mode 100644 src/gpuarrays.jl diff --git a/src/gpuarrays.jl b/src/gpuarrays.jl deleted file mode 100644 index a0a38167..00000000 --- a/src/gpuarrays.jl +++ /dev/null @@ -1,15 +0,0 @@ -# GPUArrays.jl interface - -function GPUArrays.default_rng(::Type{<:oneArray}) - dev = device() - rngs = get!(task_local_storage(), :oneAPI_GLOBAL_RNGs) do - Dict{ZeDevice,GPUArrays.RNG}() - end - get!(rngs, dev) do - N = oneL0.compute_properties(dev).maxTotalGroupSize - state = oneArray{NTuple{4, UInt32}}(undef, N) - rng = GPUArrays.RNG(state) - Random.seed!(rng) - rng - end -end diff --git a/src/oneAPI.jl b/src/oneAPI.jl index d402308e..372287b3 100644 --- a/src/oneAPI.jl +++ b/src/oneAPI.jl @@ -72,7 +72,6 @@ end # integrations and specialized functionality include("broadcast.jl") include("mapreduce.jl") -include("gpuarrays.jl") include("random.jl") include("utils.jl") diff --git a/src/random.jl b/src/random.jl index b627c9b8..e1837730 100644 --- a/src/random.jl +++ b/src/random.jl @@ -1,6 +1,13 @@ using Random -gpuarrays_rng() = GPUArrays.default_rng(oneArray) +# one GPUArrays RNG per task and device, so that streams of random numbers are independent +function gpuarrays_rng() + dev = device() + rngs = get!(task_local_storage(), :oneAPI_GLOBAL_RNGs) do + Dict{ZeDevice,GPUArrays.RNG{oneArray}}() + end + get!(() -> GPUArrays.RNG{oneArray}(), rngs, dev) +end # GPUArrays in-place Random.rand!(A::oneWrappedArray) = Random.rand!(gpuarrays_rng(), A) From 0c2a74949efe36cabbe4f152c431212c7b2a487d Mon Sep 17 00:00:00 2001 From: Tim Besard Date: Wed, 7 Oct 2026 21:19:04 +0200 Subject: [PATCH 2/3] Use GPUArrays' reductions, scans, sorting and findall GPUArrays 12 implements these for every GPU array on AcceleratedKernels and no longer calls the mapreducedim! hook, so oneAPI's own reduction kernel, its AcceleratedKernels wrappers and its findall/logical indexing go, and so does the direct dependency on AcceleratedKernels (its oneAPI extension still loads through GPUArrays). The Intel workarounds those kernels carried (scan block size 64, widening sub-word types in local memory, densifying strided inputs on the Aurora LTS stack) are not needed by AcceleratedKernels on an Iris Xe; the tests now check the cases they covered, so CI on Aurora shows whether they still are. --- Project.toml | 4 +- src/accumulate.jl | 21 ---- src/indexing.jl | 36 ------ src/mapreduce.jl | 291 ---------------------------------------------- src/oneAPI.jl | 4 - src/sorting.jl | 3 - test/array.jl | 49 ++++++-- 7 files changed, 39 insertions(+), 369 deletions(-) delete mode 100644 src/accumulate.jl delete mode 100644 src/indexing.jl delete mode 100644 src/mapreduce.jl delete mode 100644 src/sorting.jl diff --git a/Project.toml b/Project.toml index 390ed55c..3083231c 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,6 @@ version = "2.9.2" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" -AcceleratedKernels = "6a4ca0a5-0e36-4168-a932-d9be78d558f1" Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82" ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04" @@ -34,11 +33,10 @@ oneAPI_Support_jll = "b049733a-a71d-5ed3-8eba-7d323ac00b36" [compat] AbstractFFTs = "1.5.0" -AcceleratedKernels = "0.3.1, 0.4" Adapt = "4" CEnum = "0.4, 0.5" ExprTools = "0.1" -GPUArrays = "11.5.14" +GPUArrays = "12" GPUCompiler = "2.13" GPUToolbox = "3.3.1" KernelAbstractions = "0.9.39" diff --git a/src/accumulate.jl b/src/accumulate.jl deleted file mode 100644 index d3647faa..00000000 --- a/src/accumulate.jl +++ /dev/null @@ -1,21 +0,0 @@ -import oneAPI -import oneAPI: oneArray, oneAPIBackend -import AcceleratedKernels as AK - -# Use a smaller block size on Intel GPUs to work around a scan correctness issue -# with the Blelloch parallel prefix sum at larger block sizes (>=128). -const _ACCUMULATE_BLOCK_SIZE = 64 - -# Accumulate operations using AcceleratedKernels -Base.accumulate!(op, B::oneArray, A::oneArray; init = zero(eltype(A)), - block_size = _ACCUMULATE_BLOCK_SIZE, kwargs...) = - AK.accumulate!(op, B, A, oneAPIBackend(); init, block_size, kwargs...) - -Base.accumulate(op, A::oneArray; init = zero(eltype(A)), - block_size = _ACCUMULATE_BLOCK_SIZE, kwargs...) = - AK.accumulate(op, A, oneAPIBackend(); init, block_size, kwargs...) - -Base.cumsum(src::oneArray; block_size = _ACCUMULATE_BLOCK_SIZE, kwargs...) = - AK.cumsum(src, oneAPIBackend(); block_size, kwargs...) -Base.cumprod(src::oneArray; block_size = _ACCUMULATE_BLOCK_SIZE, kwargs...) = - AK.cumprod(src, oneAPIBackend(); block_size, kwargs...) diff --git a/src/indexing.jl b/src/indexing.jl deleted file mode 100644 index 13d19471..00000000 --- a/src/indexing.jl +++ /dev/null @@ -1,36 +0,0 @@ -Base.to_index(::oneArray, I::AbstractArray{Bool}) = findall(I) - -if VERSION >= v"1.11.0-DEV.1157" - Base.to_indices(x::oneArray, I::Tuple{AbstractArray{Bool}}) = - (Base.to_index(x, I[1]),) -end - -function _ker!(ys, bools, indices) - i = get_global_id() - - @inbounds if i ≤ length(bools) && bools[i] - ii = CartesianIndices(bools)[i] - b = indices[i] # new position - ys[b] = ii - end - return -end - -function Base.findall(bools::oneArray{Bool}) - I = keytype(bools) - - indices = cumsum(reshape(bools, prod(size(bools)))) - - n = isempty(indices) ? 0 : @allowscalar indices[end] - - ys = oneArray{I}(undef, n) - - if n > 0 - kernel = @oneapi launch = false _ker!(ys, bools, indices) - group_size = launch_configuration(kernel) - kernel(ys, bools, indices; items = group_size, groups = cld(length(bools), group_size)) - end - # unsafe_free!(indices) - - return ys -end diff --git a/src/mapreduce.jl b/src/mapreduce.jl deleted file mode 100644 index e767e488..00000000 --- a/src/mapreduce.jl +++ /dev/null @@ -1,291 +0,0 @@ -# TODO -# - serial version for lower latency -# - group-stride loop to delay need for second kernel launch - -# Widen sub-word types to avoid shared memory corruption on Intel GPUs. -# Writing 1/2-byte values to local memory can clobber adjacent bytes. -# Only applies to integer/boolean types where `%` conversion is valid. -@inline _widen_type(::Type{Bool}) = Int32 -@inline _widen_type(::Type{Int8}) = Int32 -@inline _widen_type(::Type{UInt8}) = Int32 -@inline _widen_type(::Type{Int16}) = Int32 -@inline _widen_type(::Type{UInt16}) = Int32 -@inline _widen_type(::Type{T}) where T = T - -# Dispatch-based conversions so the compiler never generates `%` for non-integer types -@inline _to_wide(val, ::Type{W}) where W = val % W -@inline _to_wide(val::T, ::Type{T}) where T = val -@inline _from_wide(val, ::Type{T}) where T = val % T -@inline _from_wide(val::T, ::Type{T}) where T = val - -# 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() - - # use a wider type for shared memory to avoid sub-word corruption - W = _widen_type(T) - shared = oneLocalArray(W, (maxitems,)) - @inbounds shared[item] = _to_wide(val, W) - - # perform a reduction - d = 1 - while d < items - # Fence local memory: `barrier(0)` lowers to an OpControlBarrier without the - # WorkgroupMemory storage-class bit, which does not order the shared-local tree - # accesses across the barrier. Fence local memory so each tree step sees the - # previous step's `shared[]` writes. - barrier(SPIRVIntrinsics.LOCAL_MEM_FENCE) - index = 2 * d * (item-1) + 1 - @inbounds if index <= items - other_val = if index + d <= items - _from_wide(shared[index+d], T) - else - neutral - end - shared[index] = _to_wide(op(_from_wide(shared[index], T), other_val), W) - end - d *= 2 - end - - # load the final value on the first item - if item == 1 - val = @inbounds _from_wide(shared[item], T) - 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 - -# Coalesced reduction for when the contiguous leading dimension is NOT reduced (the reduced -# axes are strided). One work-item per output slice (Rother element), grid-strided; each -# serially reduces over Rreduce. Consecutive work-items map to consecutive output slices, -# which are consecutive in memory, so global reads are coalesced across lanes — the access -# pattern a `dims=1` reduction already uses. On the Aurora LTS stack the workgroup-per-slice -# kernel above reads a *strided* reduced dimension non-coalesced, which silently corrupts -# large reductions (e.g. `sum(A; dims=2)`); this path avoids that pattern entirely. -function coalesced_mapreduce_device(f, op, neutral, Rreduce, Rother, R, As...) - iother = (get_group_id() - 1) * get_local_size() + get_local_id() - gstride = get_num_groups() * get_local_size() - @inbounds while iother <= length(Rother) - Iother = Rother[iother] - Iout = CartesianIndex(Tuple(Iother)..., 1) - neut = neutral === nothing ? R[Iout] : neutral - val = op(neut, neut) - for ireduce in 1:length(Rreduce) - Ireduce = Rreduce[ireduce] - J = max(Iother, Ireduce) - val = op(val, f(_map_getindex(As, J)...)) - end - R[Iout] = val - iother += gstride - end - return -end - -# Aurora LTS workaround: the NEO/IGC LTS stack miscompiles *strided* (non-coalesced) global -# reads inside the reduction kernel, silently corrupting results whenever an input is read -# along a non-contiguous axis (e.g. `a == transpose(b)`, `sum(transpose(x))`, `ishermitian`). -# Elementwise copies are NOT affected. `_dense_reduce_input` returns false for any input that -# reads non-contiguous memory (a transposed/permuted/strided view, or a broadcast containing -# one), so such inputs get materialized to a dense `oneArray` before reducing. The original -# symptom was `a == transpose(b)` (the `naive_transpose` example) silently returning false. -@inline _dense_reduce_input(::oneArray) = true -@inline _dense_reduce_input(x::Base.ReshapedArray) = _dense_reduce_input(parent(x)) -@inline _dense_reduce_input(::AbstractArray) = false # Transpose/Adjoint/PermutedDims/SubArray/… -@inline _dense_reduce_input(::Any) = true # scalars/Refs/tuples carried in a broadcast -@inline _dense_reduce_input(bc::Broadcast.Broadcasted) = all(_dense_reduce_input, bc.args) - -function GPUArrays.mapreducedim!(f::F, op::OP, R::oneWrappedArray{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 - - # Aurora LTS workaround (see `_dense_reduce_input` above): materialize strided inputs to a - # dense array first so every global read in the reduction kernel is coalesced. - if oneL0.LTS[] && !_dense_reduce_input(A) - Acontig = Broadcast.materialize(Broadcast.broadcasted(f, A)) - # Only recurse if materialization actually produced a dense `oneArray`. A strided - # *device* view (Transpose/PermutedDims of a oneArray, or a broadcast over one) - # materializes to a oneArray, but a plain *host* AbstractArray materializes to - # another host array that is still not `_dense_reduce_input` — recursing on it would - # loop forever (StackOverflowError). Fall through instead so the normal path raises - # the usual kernel-conversion error for the unsupported host input. - if Acontig isa oneArray - return GPUArrays.mapreducedim!(identity, op, R, Acontig; init = init) - end - end - - 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)) - - # Aurora LTS workaround: the workgroup-per-slice kernel below reads a *strided* reduced - # dimension non-coalesced, which silently corrupts reductions on this stack (regardless of - # output count — it depends on the reduction length, not the number of slices). Whenever - # the contiguous leading dimension is NOT reduced (`size(Rreduce, 1) == 1`), use the - # coalesced one-work-item-per-slice kernel, whose lanes read consecutive memory. Few-slice - # reductions get less parallelism but stay correct; the common many-slice case is also fast. - # - # `size(Rreduce, 1) == 1` (i.e. dim 1 kept) is the correct predicate, not just a special - # case. In `partial_mapreduce_device` adjacent work-item lanes step through adjacent - # `Rreduce` entries, whose *fastest* axis is array dim 1. The miscompile bites only when - # adjacent lanes land a *full stride* apart on every lane — i.e. when the innermost reduced - # axis is itself strided, which is exactly the dim-1-kept case routed to the coalesced - # kernel above. When dim 1 is among the reduced dims, adjacent lanes instead read - # consecutive memory (stride 1) within each length-n1 run, with only an occasional jump at - # a run boundary — not the every-lane full-stride pattern — so the reads stay well-formed - # even when a strided dim (e.g. dim 3) is *also* reduced. Verified with exact Int32 - # reductions (immune to the accumulation-order rounding that makes >2^24-length Float32 - # sums differ between GPU and CPU): `dims=(1,3)`, `(1,4)` and full `(1,2,3)` reductions at - # large strided sizes match the CPU exactly, so these need no extra materialization (see - # the "strided mixed reductions" regression test in test/array.jl). - if oneL0.LTS[] && size(Rreduce, 1) == 1 - # cap the group size at what this kernel actually supports on this device (queried - # from the driver), rather than a hardcoded 256 that can exceed the kernel's max - # work-group size on some devices and fail the launch. - cargs = (f, op, init, Rreduce, Rother, R′, A) - ckernel = zefunction( - coalesced_mapreduce_device, - Tuple{Core.Typeof.(kernel_convert.(cargs))...} - ) - items = clamp(length(Rother), 1, launch_configuration(ckernel)) - groups = min(cld(length(Rother), items), 1024) - @oneapi items = items groups = groups coalesced_mapreduce_device( - f, op, init, Rreduce, Rother, R′, A - ) - return R_old - end - - # 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 (use widened type for sub-word types) - max_lmem_elements = compute_properties(device()).maxSharedLocalMemory ÷ sizeof(_widen_type(T)) - max_items = min(compute_properties(device()).maxTotalGroupSize, - 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 = zefunction(partial_mapreduce_device, kernel_tt) - reduce_items = compute_items(launch_configuration(kernel)) - - # 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 - items = reduce_items - groups = reduce_groups*other_groups - - # perform the actual reduction - if reduce_groups == 1 - # we can cover the dimensions to reduce using a single group - @oneapi items groups partial_mapreduce_device( - f, op, init, Val(items), 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 - @oneapi items groups partial_mapreduce_device( - f, op, init, Val(items), Rreduce, Rother, partial, A) - - GPUArrays.mapreducedim!(identity, op, R′, partial; init=init) - end - - return R_old -end diff --git a/src/oneAPI.jl b/src/oneAPI.jl index 372287b3..1f098f68 100644 --- a/src/oneAPI.jl +++ b/src/oneAPI.jl @@ -71,15 +71,11 @@ end # integrations and specialized functionality include("broadcast.jl") -include("mapreduce.jl") include("random.jl") include("utils.jl") include("oneAPIKernels.jl") import .oneAPIKernels: oneAPIBackend -include("accumulate.jl") -include("sorting.jl") -include("indexing.jl") export oneAPIBackend # precompilation workload (warms up the SPIR-V compilation pipeline) diff --git a/src/sorting.jl b/src/sorting.jl deleted file mode 100644 index f83ea078..00000000 --- a/src/sorting.jl +++ /dev/null @@ -1,3 +0,0 @@ -Base.sort!(x::oneArray; kwargs...) = (AK.sort!(x; kwargs...); return x) -Base.sortperm!(ix::oneArray, x::oneArray; kwargs...) = (AK.sortperm!(ix, x; kwargs...); return ix) -Base.sortperm(x::oneArray; kwargs...) = sortperm!(oneArray(1:length(x)), x; kwargs...) diff --git a/test/array.jl b/test/array.jl index f671a755..1d5e080b 100644 --- a/test/array.jl +++ b/test/array.jl @@ -173,13 +173,18 @@ end @test length(b) == 1 end -@testset "strided mixed reductions" begin - # The Aurora LTS IGC miscompiles a reduction kernel's global reads when the *innermost* - # reduced axis is strided (dim 1 kept, e.g. `dims=2`); mapreducedim! routes those to a - # coalesced kernel. Reductions that also reduce dim 1 (e.g. `dims=(1,3)`) keep a contiguous - # innermost axis and stay correct on the workgroup-per-slice kernel — including with a small - # leading dim, where the contiguous run is short. Use Int32 (exact, associative) so the - # check is immune to Float32 accumulation-order rounding. +# Reductions and scans come from GPUArrays (AcceleratedKernels). oneAPI used to work around +# Intel problems in its own kernels; these tests check the cases those workarounds covered. + +@testset "strided reductions" begin + # The Aurora LTS stack miscompiled strided global reads in oneAPI's reduction kernel, both + # for strided inputs (`a == transpose(b)`) and for reductions whose innermost reduced axis + # is strided (dim 1 kept). Int32 keeps the comparison exact. + A = rand(Int32(1):Int32(4), 64, 64) + dA = oneArray(A) + @test dA == transpose(oneArray(collect(transpose(A)))) + @test sum(transpose(dA)) == sum(A) + @test Array(sum(transpose(dA); dims=1)) == sum(transpose(A); dims=1) for (sz, dts) in ( ((2, 512, 64), ((1, 3), (2,), (3,), (2, 3), (1, 2, 3), 1)), ((3, 256, 48), ((1, 3), (2,), (1, 2, 3))), @@ -206,8 +211,30 @@ end @test Array(dR) == R end -@testset "mapreducedim! returning same type" begin - R = transpose(oneAPI.zeros(Float32, 2, 3)) - A = oneArray(rand(Float32, 3, 2, 10)) - @test @inferred(oneAPI.GPUArrays.mapreducedim!(identity, +, R, A)) === R +@testset "sub-word reductions" begin + # Writing 1- and 2-byte values to local memory clobbered adjacent bytes on some Intel GPUs + for T in (Bool, Int8, UInt8, Int16, UInt16), n in (31, 257, 4097, 1_000_003) + A = T == Bool ? rand(Bool, n) : rand(T(0):T(1), n) + dA = oneArray(A) + @test reduce(|, dA) == reduce(|, A) + @test reduce(xor, dA) == reduce(xor, A) + @test maximum(dA) == maximum(A) + @test minimum(dA) == minimum(A) + end + for T in (Bool, Int8, Int16), (sz, dims) in (((64, 1000), 1), ((64, 1000), 2), ((17, 33, 65), (1, 3))) + A = T == Bool ? rand(Bool, sz...) : rand(T(0):T(1), sz...) + @test Array(reduce(xor, oneArray(A); dims, init=zero(T))) == reduce(xor, A; dims, init=zero(T)) + @test Array(maximum(oneArray(A); dims)) == maximum(A; dims) + end +end + +@testset "scans" begin + # Scans with a block size of 128 or more were wrong on some Intel GPUs + for T in (Int32, Float32), n in (127, 128, 129, 256, 1000, 100_000, 1_000_000) + A = rand(T(0):T(3), n) + @test Array(cumsum(oneArray(A))) == cumsum(A) + end + A = rand(Int32(0):Int32(3), 300, 700) + @test Array(cumsum(oneArray(A); dims=1)) == cumsum(A; dims=1) + @test Array(cumsum(oneArray(A); dims=2)) == cumsum(A; dims=2) end From 1bb88616800a689e9bb07bcc39fca35e687da565 Mon Sep 17 00:00:00 2001 From: Tim Besard Date: Thu, 8 Oct 2026 10:23:27 +0200 Subject: [PATCH 3/3] Detect aliasing through GPUArrays' memory location hook Base.dataids returned the array's own pointer, which missed overlaps between a contiguous view (a oneArray at an offset) and a wrapped array of the same memory, and the byte-range mightalias let empty arrays alias. GPUArrays 12.1 defines `Base.dataids` and `Base.mightalias` for every GPU array from where its elements live, so implement that hook instead, reading the address from the buffer. --- Project.toml | 2 +- src/array.jl | 9 ++------- test/array.jl | 38 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 41 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index 3083231c..569f884f 100644 --- a/Project.toml +++ b/Project.toml @@ -36,7 +36,7 @@ AbstractFFTs = "1.5.0" Adapt = "4" CEnum = "0.4, 0.5" ExprTools = "0.1" -GPUArrays = "12" +GPUArrays = "12.1" GPUCompiler = "2.13" GPUToolbox = "3.3.1" KernelAbstractions = "0.9.39" diff --git a/src/array.jl b/src/array.jl index 514a4251..1a86dbb2 100644 --- a/src/array.jl +++ b/src/array.jl @@ -151,16 +151,11 @@ GPUArrays.storage(a::oneArray) = a.data ## alias detection -Base.dataids(A::oneArray) = (UInt(pointer(A)),) +# GPUArrays implements `Base.dataids` and `Base.mightalias` from where an array lives +GPUArrays.memory_location(A::oneArray) = (UInt(convert(ZePtr{Cvoid}, A.data[])), A.offset) Base.unaliascopy(A::oneArray) = copy(A) -function Base.mightalias(A::oneArray, B::oneArray) - 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) -end - ## convenience constructors diff --git a/test/array.jl b/test/array.jl index 1d5e080b..c014cf2b 100644 --- a/test/array.jl +++ b/test/array.jl @@ -87,6 +87,44 @@ end @test Array(r) == reinterpret(Int64, @view Array(a)[2:7]) end +@testset "aliasing" begin + x = oneArray([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 oneArrays with an offset into the parent's memory, + # which should still alias wrapped arrays (like SubArrays) of that memory + x = oneArray(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(oneArray(1:16), 15:-1:1)) + + # so in-place broadcasts between them should make a copy first + n = 2^20 + x = oneArray{Float32}(1:n) + view(x, 2:n) .= view(x, n-1:-1:1) + @test Array(x) == [1; n-1:-1:1] + + # empty arrays alias nothing, also on Julia 1.10 + @test !Base.mightalias(oneArray(Int[]), oneArray(Float32[])) + @test !Base.mightalias(view(oneArray(zeros(Float32, 2, 0)), 1:1, :), oneArray(Int[])) + + # disjoint parts of one array may be each other's source and destination + x = oneArray(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 + @testset "shared buffers & unsafe_wrap" begin a = oneVector{Int,oneL0.SharedBuffer}(undef, 2)