Prepare the GPUArrays 12 release

.buildkite/pipeline.yml

@@ -29,10 +29,9 @@ steps:
                 println("--- :julia: Installing CUDA.jl")
                 withenv("JULIA_PKG_PRECOMPILE_AUTO" => 0,
                         "JULIA_PKG_DEVDIR" => devdir) do
-                  Pkg.develop("CUDA")
-                  Pkg.activate(package)
-
                   try
+                    Pkg.develop("CUDA")
+                    Pkg.activate(package)
                     Pkg.develop([
                       PackageSpec(path=gpuarrays), PackageSpec(path=gpuarrayscore),
                       PackageSpec(path=joinpath(package, "CUDACore")),
@@ -89,10 +88,9 @@ steps:
                 println("--- :julia: Installing CUDA.jl")
                 withenv("JULIA_PKG_PRECOMPILE_AUTO" => 0,
                         "JULIA_PKG_DEVDIR" => devdir) do
-                  Pkg.develop("CUDA")
-                  Pkg.activate(package)
-
                   try
+                    Pkg.develop("CUDA")
+                    Pkg.activate(package)
                     Pkg.develop([PackageSpec(path=gpuarrays), PackageSpec(path=gpuarrayscore)])
                     TestEnv.activate()
                   catch err
@@ -141,10 +139,9 @@ steps:
                 println("--- :julia: Installing AMDGPU.jl")
                 withenv("JULIA_PKG_PRECOMPILE_AUTO" => 0,
                         "JULIA_PKG_DEVDIR" => devdir) do
-                  Pkg.develop("AMDGPU")
-                  Pkg.activate(package)
-
                   try
+                    Pkg.develop("AMDGPU")
+                    Pkg.activate(package)
                     Pkg.develop([PackageSpec(path=gpuarrays), PackageSpec(path=gpuarrayscore)])
                     TestEnv.activate()
                   catch err
@@ -195,11 +192,10 @@ steps:
                 println("--- :julia: Installing oneAPI.jl")
                 withenv("JULIA_PKG_PRECOMPILE_AUTO" => 0,
                         "JULIA_PKG_DEVDIR" => devdir) do
-                  Pkg.develop("oneAPI")
-                  include(joinpath(package, "deps", "build_ci.jl"))
-                  Pkg.activate(package)
-
                   try
+                    Pkg.develop("oneAPI")
+                    include(joinpath(package, "deps", "build_ci.jl"))
+                    Pkg.activate(package)
                     Pkg.develop([PackageSpec(path=gpuarrays), PackageSpec(path=gpuarrayscore)])
                     TestEnv.activate()
                   catch err
@@ -249,10 +245,9 @@ steps:
                 println("--- :julia: Installing Metal.jl")
                 withenv("JULIA_PKG_PRECOMPILE_AUTO" => 0,
                         "JULIA_PKG_DEVDIR" => devdir) do
-                  Pkg.develop("Metal")
-                  Pkg.activate(package)
-
                   try
+                    Pkg.develop("Metal")
+                    Pkg.activate(package)
                     Pkg.develop([PackageSpec(path=gpuarrays), PackageSpec(path=gpuarrayscore)])
                     TestEnv.activate()
                   catch err

.github/workflows/Test.yml

@@ -85,14 +85,17 @@ jobs:
                     println("--- :julia: Installing OpenCL.jl")
                     withenv("JULIA_PKG_PRECOMPILE_AUTO" => 0,
                             "JULIA_PKG_DEVDIR" => devdir) do
-                      Pkg.develop("OpenCL")
-                      Pkg.activate(package)
-
                       try
+                        Pkg.develop("OpenCL")
+                        Pkg.activate(package)
                         Pkg.develop([PackageSpec(path=gpuarrays), PackageSpec(path=gpuarrayscore)])
                         TestEnv.activate()
                       catch err
+                        # OpenCL.jl may not support this GPUArrays version yet (e.g. after a
+                        # breaking release): skip, like Buildkite soft-fails, but annotate
+                        # the run so the skip is visible
                         @error "Could not install OpenCL.jl" exception=(err,catch_backtrace())
+                        println("::warning title=OpenCL.jl not tested::Could not install OpenCL.jl with this version of GPUArrays")
                         exit(0)
                       finally
                         Pkg.activate(package)

NEWS.md

@@ -0,0 +1,47 @@
+# GPUArrays.jl release history
+
+This document lists the noteworthy changes in every breaking release of GPUArrays.jl,
+newest first, aimed at the authors of back-end packages. For the complete list of merged
+pull requests, see the [GitHub releases](https://github.com/JuliaGPU/GPUArrays.jl/releases).
+
+
+## v12.0
+
+Base's reductions, sorting, scans, `findall`, logical indexing and `reverse` are now
+implemented once for every `AnyGPUArray`, on top of
+[AcceleratedKernels.jl](https://github.com/JuliaGPU/AcceleratedKernels.jl) (AK) 0.5. A
+back-end gets all of them by defining nothing beyond the array interface and a
+KernelAbstractions back-end, and they follow Base's semantics: stable sorting by default,
+`init` applied once, no neutral element needed, Base's result types, empty results and
+errors.
+
+*Breaking changes for back-ends*:
+
+- The `GPUArrays.mapreducedim!` hook is removed and GPUArrays no longer picks neutral
+  elements. Every Base reduction goes to AK directly. Back-ends must delete their methods
+  for the hook, together with the reduction kernels only it used.
+- `GPUArrays.default_rng` and the deprecated `RNG(state::AbstractGPUArray)` and
+  `seed!(rng, ::Vector{UInt32})` methods are removed. Construct the RNG with `RNG{AT}()`
+  and seed it with an integer.
+- KernelAbstractions 0.9.43 or later and Adapt 4.7.3 or later are required.
+
+*What back-ends should delete*:
+
+Methods narrower than GPUArrays' keep winning dispatch, so a back-end that keeps its own
+`findall`, `to_index`/`to_indices`, `_accumulate!`/`accumulate`, `sort!`/`sortperm!`/
+`partialsort!`, or `reverse!` bypasses GPUArrays' Base conformance, and some of those
+methods are ambiguous with GPUArrays' (e.g. `findall(f, ::MyArray)` against
+`findall(::Fix2{typeof(in)}, ::AnyGPUArray)`). Delete them.
+
+A vendor sort or scan (e.g. one built on MPSGraph) can stay as a fast path, as long as it
+falls back to GPUArrays' method for the cases it doesn't handle and accepts AK's
+algorithm objects (`alg=AK.MergeSort()` etc.). There is no extension point for vendor
+reductions; `GPUArrays._ak_mapreducedim!` is internal.
+
+*Known regressions*:
+
+- Small reductions along `dims` can be slower than the back-ends' own kernels were (up to
+  about 4× for a 1000×10 array on CUDA), while large ones are faster. This is tracked in
+  [AcceleratedKernels.jl#135](https://github.com/JuliaGPU/AcceleratedKernels.jl/issues/135).
+
+JLArrays 0.4.0 is the matching release of the reference back-end.

Project.toml

@@ -28,10 +28,10 @@ JLD2Ext = "JLD2"
 
 [compat]
 AcceleratedKernels = "0.5"
-Adapt = "4.6.1"
+Adapt = "4.7.3"
 GPUArraysCore = "= 0.2.1"
 JLD2 = "0.4, 0.5, 0.6"
-KernelAbstractions = "0.9.28, 0.10"
+KernelAbstractions = "0.9.43, 0.10"
 LLVM = "3.9, 4, 5, 6, 7, 8, 9, 10"
 LinearAlgebra = "1"
 Printf = "1"

README.md

@@ -39,6 +39,14 @@ Some (CPU based) examples can be see in the testing library `JLArrays` (located
 
 ### Dense array support
 
+Define an array type that subtypes `AbstractGPUArray`, a device-side counterpart with
+`Adapt.adapt_structure` rules to convert to it, and a KernelAbstractions back-end returned
+by `KernelAbstractions.get_backend`. Broadcasting, `mapreduce` and the other reductions,
+sorting, scans, `findall`, logical indexing, `reverse` and most of linear algebra then work
+without further methods. See the [interface
+documentation](https://juliagpu.github.io/GPUArrays.jl/dev/interface/) and `JLArrays` for a
+complete example.
+
 ### Sparse array support (optional)
 
 `GPUArrays.jl` provides **device-side** array types for `CSC`, `CSR`, `COO`, and `BSR` matrices, as well as sparse vectors.
@@ -54,8 +62,6 @@ For **host-side** types, your custom sparse types should implement:
 - `csr_type(::Type{T})` - the compressed sparse row type for your backend. A `CuSparseMatrixCSC` would have a `csr_type` of `CuSparseMatrixCSR`. 
 - `coo_type(::Type{T})` - the coordinate sparse matrix type for your backend. A `CuSparseMatrixCSC` would have a `coo_type` of `CuSparseMatrixCOO`.
 
-To use `SparseArrays.findnz`, your host-side type **must** implement `sortperm`. This can be done with scalar indexing, but will be very slow.
-
 Additionally, you need to teach `GPUArrays.jl` how to translate your backend's specific types onto the device. `GPUArrays.jl` provides the device-side types:
 
 - `GPUSparseDeviceVector`

docs/src/interface.md

@@ -29,7 +29,7 @@ struct CustomBackend <: KernelAbstractions.GPU
 KernelAbstractions.get_backend(a::CA) where CA <: CustomArray = CustomBackend()
 ```
 
-There are numerous examples of potential interfaces for GPUArrays, such as with [JLArrays](https://github.com/JuliaGPU/GPUArrays.jl/blob/main/lib/JLArrays/src/JLArrays.jl), [CuArrays](https://github.com/JuliaGPU/CUDA.jl/blob/main/src/gpuarrays.jl), and [ROCArrays](https://github.com/JuliaGPU/AMDGPU.jl/blob/main/src/gpuarrays.jl).
+There are numerous examples of potential interfaces for GPUArrays, such as with [JLArrays](https://github.com/JuliaGPU/GPUArrays.jl/blob/main/lib/JLArrays/src/JLArrays.jl), [CuArrays](https://github.com/JuliaGPU/CUDA.jl/blob/main/CUDACore/src/array.jl), and [ROCArrays](https://github.com/JuliaGPU/AMDGPU.jl/blob/main/src/array.jl).
 
 ## Caching Allocator
 

lib/JLArrays/Project.toml

@@ -1,6 +1,6 @@
 name = "JLArrays"
 uuid = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
-version = "0.4.1"
+version = "0.4.0"
 authors = ["Tim Besard <tim.besard@gmail.com>"]
 
 [deps]

src/GPUArrays.jl

@@ -43,6 +43,5 @@ include("host/statistics.jl")
 include("host/sparse.jl")
 include("host/alloc_cache.jl")
 
-include("deprecated.jl")
 
 end # module

src/host/mapreduce.jl

@@ -19,6 +19,18 @@ Base.mapreducedim!(f, op, R::AnyGPUArray, A::AbstractArray) =
     _ak_mapreducedim!(f, op, R, A; backend=get_backend(R))
 Base.mapreducedim!(f, op, R::AnyGPUArray, A::Broadcast.Broadcasted) =
     _ak_mapreducedim!(f, op, R, A; backend=get_backend(R))
+# (LinearAlgebra's methods for adjoint and transposed matrices are ambiguous with these)
+for W in (:Adjoint, :Transpose)
+    @eval begin
+        Base.mapreducedim!(f, op::Union{typeof(&), typeof(|), typeof(+), typeof(Base.add_sum),
+                                        typeof(max), typeof(min), typeof(Base._extrema_rf)},
+                           R::AnyGPUArray, A::$W{<:Any, <:AbstractMatrix}) =
+            _ak_mapreducedim!(f, op, R, A; backend=get_backend(R))
+        Base.mapreducedim!(f::typeof(identity), op::Union{typeof(*), typeof(Base.mul_prod)},
+                           R::AnyGPUArray, A::$W{<:Union{Real, Complex}, <:AbstractMatrix}) =
+            _ak_mapreducedim!(f, op, R, A; backend=get_backend(R))
+    end
+end
 
 # `neutral_element` lives in GPUArraysCore, so that packages building on GPUArraysCore share it
 import GPUArraysCore: neutral_element

test/testsuite/accumulate.jl

@@ -61,4 +61,14 @@
     # Associative operators need not be commutative
     @test compare(A -> accumulate((a, b) -> a, A), AT, rand(1:10, 1000))
     @test compare(A -> accumulate((a, b) -> a, A; dims=2), AT, rand(1:10, 10, 100))
+    # composing affine maps x -> a*x + b is exact and sensitive to operand order, also in the
+    # combination of block aggregates
+    compose((a1, b1), (a2, b2)) = (a2 * a1, a2 * b1 + b2)
+    affine(dims...) = [(rand(-1:1), rand(-3:3)) for _ in CartesianIndices(dims)]
+    for n in (4, 1_000, 100_000)
+        @test compare(A -> accumulate(compose, A), AT, affine(n))
+        @test compare(A -> accumulate(compose, A; init=(1, 0)), AT, affine(n))
+    end
+    @test compare(A -> accumulate(compose, A; dims=2, init=(1, 0)), AT, affine(7, 10_000))
+    @test compare(A -> accumulate(compose, A; dims=1), AT, affine(10_000, 7))
 end

test/testsuite/indexing.jl

@@ -171,6 +171,12 @@ end
         # 1D
         @test compare(findfirst, AT, rand(Bool, 100))
         @test compare(x->findfirst(>(0.5f0), x), AT, rand(Float32, 100))
+        # predicates with fields, which the kernel's closure captures
+        let c = 0.5f0
+            @test compare(x->findfirst(y -> y > c, x), AT, rand(Float32, 100))
+        end
+        @test compare(x->findfirst(!isnan, x), AT, Float32[NaN, NaN, 1, NaN])
+        @test compare(x->findfirst(Base.Fix2(>, 0.5f0) ∘ abs, x), AT, rand(Float32, 100))
         @test compare(findfirst, AT, fill(false, 10))
         let x = fill(false, 10)
             @test findfirst(x) == findfirst(AT(x))
@@ -192,6 +198,12 @@ end
         # 1D
         @test compare(findlast, AT, rand(Bool, 100))
         @test compare(x->findlast(>(0.5f0), x), AT, rand(Float32, 100))
+        # predicates with fields, which the kernel's closure captures
+        let c = 0.5f0
+            @test compare(x->findlast(y -> y > c, x), AT, rand(Float32, 100))
+        end
+        @test compare(x->findlast(!isnan, x), AT, Float32[NaN, NaN, 1, NaN])
+        @test compare(x->findlast(Base.Fix2(>, 0.5f0) ∘ abs, x), AT, rand(Float32, 100))
         let x = fill(false, 10)
             @test findlast(x) == findlast(AT(x))
         end

test/testsuite/reductions.jl

@@ -39,6 +39,10 @@ end
         for t in [transpose, adjoint]
             @test compare((A,R)->Base.mapreducedim!(identity, +, R, A), AT, rand(range, (2,2)), t(zeros(ET, (2,))))
             @test compare((A,R)->Base.mapreducedim!(abs2, +, R, A), AT, rand(range, (3,2,10)), t(zeros(ET, (2,3))))
+            # ... and from them (LinearAlgebra has methods for these sources)
+            @test compare((A,R)->Base.mapreducedim!(identity, +, R, t(A)), AT, rand(range, (3,4)), zeros(ET, (4,)))
+            @test compare((A,R)->Base.mapreducedim!(abs2, max, R, t(A)), AT, rand(range, (3,4)), zeros(real(ET), (1,3)))
+            @test compare((A,R)->Base.mapreducedim!(identity, *, R, t(A)), AT, rand(range, (3,4)), ones(ET, (4,)))
 
             # TODO: reenable once https://github.com/JuliaGPU/Metal.jl/issues/907 is fixed
             # @test compare((A,R)->sum!(abs2, R, A), AT, rand(range, (3,2,10)), t(zeros(ET, (2,3))))
@@ -346,7 +350,7 @@ end
             @test same(f, m8)
         end
         # (Base's pairwise path reduces `Int8`s in `Int8`, which the values here do not overflow)
-        @test same(A -> reduce(+, A; dims, init=0.5), Int16.(m8))
+        @test same(A -> reduce(+, A; dims, init=0.5f0), Int16.(m8))
     end
     # ... and empty reduced dimensions: `init`, else Base's initial value or error
     e8 = zeros(Int8, 0, 3)

test/testsuite/sorting.jl

@@ -63,6 +63,10 @@ end
     @test_throws UndefKeywordError sort(AT(rand(Float32, 4, 4)))
     @test_throws UndefKeywordError sort!(AT(rand(Float32, 4, 4)))
     @test_throws UndefKeywordError sortperm(AT(rand(Float32, 4, 4)))
+    # a `dims` out of range (Base's error for `sort(A; dims=0)` differs between versions)
+    for (f, dims) in ((sort!, 3), (sort, 3), (sortperm, 0), (sortperm, 3))
+        @test_throws ArgumentError f(AT(rand(Float32, 4, 4)); dims)
+    end
 end
 
 @testsuite "sorting/sortperm" (AT, eltypes)->begin