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)