diff --git a/Project.toml b/Project.toml index 390ed55c..569f884f 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.1" 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/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/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/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 d402308e..1f098f68 100644 --- a/src/oneAPI.jl +++ b/src/oneAPI.jl @@ -71,16 +71,11 @@ end # integrations and specialized functionality include("broadcast.jl") -include("mapreduce.jl") -include("gpuarrays.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/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) 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..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) @@ -173,13 +211,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 +249,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