Merge remote-tracking branch 'upstream/master' into QR_views

Author: evelyne-ringoot

.buildkite/pipeline.yml

@@ -58,10 +58,6 @@ steps:
                 Pkg.develop(; path=pwd())
                 Pkg.develop(; name="Metal")
 
-                println("+++ :julia: Building support library")
-                include(joinpath(Pkg.devdir(), "Metal", "deps", "build_ci.jl"))
-                Pkg.activate()
-
                 println("+++ :julia: Running tests")
                 Pkg.test("Metal"; coverage=true)'
     agents:

Project.toml

@@ -1,6 +1,6 @@
 name = "GPUArrays"
 uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
-version = "8.6.3"
+version = "8.6.6"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -16,6 +16,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 [compat]
 Adapt = "2.0, 3.0"
 GPUArraysCore = "= 0.1.4"
-LLVM = "3.9, 4"
+LLVM = "3.9, 4, 5"
 Reexport = "1"
 julia = "1.6"

src/host/broadcast.jl

@@ -51,53 +51,22 @@ end
     bc′ = Broadcast.preprocess(dest, bc)
 
     # grid-stride kernel
-    function broadcast_kernel(ctx, dest, ::Val{Is}, bc′, nelem) where Is
-        j = 0
-        while j < nelem
-            j += 1
-
-            i = @linearidx(dest, j)
-
-            # cartesian indexing is slow, so avoid it if possible
-            if isa(IndexStyle(dest), IndexCartesian) || isa(IndexStyle(bc′), IndexCartesian)
-                # this performs an integer division, which is expensive. to make it possible
-                # for the compiler to optimize it away, we put the iterator in the type
-                # domain so that the indices are available at compile time. note that LLVM
-                # only seems to replace pow2 divisions (with bitshifts), but other back-ends
-                # may be smarter and replace arbitrary divisions by bit operations.
-                #
-                # also see maleadt/StaticCartesian.jl, which implements this in Julia,
-                # but does not result in an additional speed-up on tested back-ends.
-                #
-                # in addition, we use @inbounds to avoid bounds checks, but we also need to
-                # inform the compiler about the bounds that we are assuming. this is done
-                # using the assume intrinsic, and in case of Metal yields a 8x speed-up.
-                assume(1 <= i <= length(Is))
-                I = @inbounds Is[i]
-            end
-
-            val = if isa(IndexStyle(bc′), IndexCartesian)
-                @inbounds bc′[I]
-            else
-                @inbounds bc′[i]
-            end
-
-            if isa(IndexStyle(dest), IndexCartesian)
-                @inbounds dest[I] = val
-            else
-                @inbounds dest[i] = val
-            end
+    function broadcast_kernel(ctx, dest, bc′, nelem)
+        i = 0
+        while i < nelem
+            i += 1
+            I = @cartesianidx(dest, i)
+            @inbounds dest[I] = bc′[I]
         end
         return
     end
     elements = length(dest)
     elements_per_thread = typemax(Int)
-    Is = CartesianIndices(dest)
-    heuristic = launch_heuristic(backend(dest), broadcast_kernel, dest, Val(Is), bc′, 1;
+    heuristic = launch_heuristic(backend(dest), broadcast_kernel, dest, bc′, 1;
                                  elements, elements_per_thread)
     config = launch_configuration(backend(dest), heuristic;
                                   elements, elements_per_thread)
-    gpu_call(broadcast_kernel, dest, Val(Is), bc′, config.elements_per_thread;
+    gpu_call(broadcast_kernel, dest, bc′, config.elements_per_thread;
              threads=config.threads, blocks=config.blocks)
 
     return dest

src/host/math.jl

@@ -2,7 +2,7 @@
 
 function Base.clamp!(A::AnyGPUArray, low, high)
     gpu_call(A, low, high) do ctx, A, low, high
-        I = @linearidx A
+        I = @cartesianidx A
         A[I] = clamp(A[I], low, high)
         return
     end

test/testsuite/math.jl

@@ -17,7 +17,7 @@ end
 @testsuite "math/power" (AT, eltypes)->begin
     for ET in eltypes
         for p in 0:5
-            compare(x->x^p, AT, rand(ET, 2,2))
+            @test compare(x->x^p, AT, rand(ET, 2,2))
         end
     end
 end

test/testsuite/reductions.jl

@@ -15,6 +15,22 @@
     end
 end
 
+@testsuite "reductions/mapreducedim!_large" (AT, eltypes)->begin
+    @testset "$ET" for ET in eltypes
+        # Skip smaller floating types due to precision issues
+        if ET in (Float16, ComplexF16)
+            continue
+        end
+
+        range = ET <: Real ? (ET(1):ET(10)) : ET
+        # Reduce larger array sizes to test multiple-element reading in certain implementations
+        for (sz,red) in [(1000000,)=>(1,), (5000,500)=>(1,1), (500,5000)=>(1,1),
+                         (500,5000)=>(500,1), (5000,500)=>(1,500)]
+            @test compare((A,R)->Base.mapreducedim!(identity, +, R, A), AT, rand(range, sz), zeros(ET, red))
+        end
+    end
+end
+
 @testsuite "reductions/reducedim!" (AT, eltypes)->begin
     @testset "$ET" for ET in eltypes
         range = ET <: Real ? (ET(1):ET(10)) : ET