Slow 2D sum

[Cvikli]

Coming from this observation #issue 1323.

Having a code:

using BenchmarkTools
using CUDA
function mysum(X,Y,n) 
	I = (blockIdx().x - 1) * blockDim().x + threadIdx().x
	I > n && return
	cY = CUDA.Const(Y)
	@inbounds for j in 1:n
		X[I] += cY[I,j]
	end
	return
end

# MYSUM PERF. TEST.
nth=512
N = 500; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 1000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 2000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 4000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 6000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 8000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 10000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 12000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 

# BUILTIN PERF. TEST.
N = 500; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))
N = 1000; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))
N = 2000; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))
N = 4000; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))
N = 6000; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))
N = 8000; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))
N = 10000; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))
N = 12000; A = CUDA.randn(Float32, N, N); (@btime CUDA.@sync sum($A, dims=2))

Gives us this results:

  27.702 μs (44 allocations: 1.84 KiB)
  80.011 μs (34 allocations: 1.44 KiB)
  155.594 μs (168 allocations: 5.62 KiB)
  323.131 μs (40 allocations: 1.62 KiB)
  576.200 μs (314 allocations: 10.19 KiB)
  885.445 μs (630 allocations: 20.06 KiB)
  1.408 ms (988 allocations: 31.25 KiB)
  1.965 ms (1230 allocations: 38.81 KiB)

  61.396 μs (134 allocations: 4.84 KiB)
  144.424 μs (266 allocations: 8.95 KiB)
  333.410 μs (321 allocations: 10.67 KiB)
  968.642 μs (1165 allocations: 37.05 KiB)
  2.979 ms (1567 allocations: 49.61 KiB)
  5.656 ms (2234 allocations: 70.45 KiB)
  8.874 ms (10531 allocations: 329.73 KiB)
  13.106 ms (19971 allocations: 624.73 KiB)

Possibly: 200-1000% speedup.

Why are we faster with "mysum"? Is it the CUDA.const? Or the specialisation on dims=2? Or we don't see something?

[Cvikli]

Or... the preallocation! Damn, I am comparing different operation. :(

I will be back with the right measurement.

[Cvikli]

This shows a little bit better speed without preallocation:

N = 500; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))
N = 1000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))
N = 2000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))
N = 4000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))
N = 6000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))
N = 8000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))
N = 10000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))
N = 12000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync sum!($B, $A))

New results:

  41.969 μs (73 allocations: 2.47 KiB)
  95.652 μs (152 allocations: 4.94 KiB)
  246.136 μs (207 allocations: 6.66 KiB)
  921.868 μs (1572 allocations: 49.31 KiB)
  2.781 ms (2101 allocations: 65.84 KiB)
  5.428 ms (5666 allocations: 177.25 KiB)
  8.370 ms (10369 allocations: 324.22 KiB)
 12.587 ms (16901 allocations: 528.34 KiB)

But actually still notable difference.

[Cvikli]

Don't want to mix things, but running the bulitin sum(..., dims=2) and mean(..., dims=2) gives nearly exactly the same speed every testcases. Does this help for you to understand what can be the problem?

[maleadt]

running the bulitin sum(..., dims=2) and mean(..., dims=2) gives nearly exactly the same speed every testcases

That's expected, right? mean is just sum divided by length.

[Cvikli]

Yes, definitely. I just wanted to make sure that the sum is really important as there are really multiple aggregation algorithm that can depend on that.

[maleadt]

I wonder if it's just a launch configuration issue again, because generally CUDA.jl uses the max threads returned by the occupancy API, likely 1024 threads. Could you compare against that?

[Cvikli]

You mean test with 64/256/512/1024 threads and compare each testcases?

[maleadt]

Yeah. nth=1024 would compare with CUDA.jl's launch configuration.

[Cvikli]

I got the same results. The code I ran:

nth=1024
N = 500; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 1000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 2000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 4000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 6000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 8000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 10000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 
N = 12000; A = CUDA.randn(Float32, N, N); B = CUDA.randn(Float32, N, 1); (@btime CUDA.@sync @cuda threads=nth blocks=cld(N*N, nth)  mysum($B,$A,$N)) 

[Tuebel]

Hi,

my idea would be to turn the summation into a matrix vector multiplication, where the vector consists of ones. I have implemented and benchmarked a few different combinations:

1. Reducing only one dimension via multiplication and then using sum
2. Reducing both dimensions via multiplication
3. Testing 1 again with a preallocated vector

Here is the code:

using CUDA
M_rand = CUDA.rand(Float32, 1000, 1001)
M_cpu = Array(M_rand)
v_ones = CUDA.ones(Float32, 1, size(M_rand)[1])
M_ones = CUDA.ones(Float32, size(M_rand))

function dot_sum(M)
    v = CUDA.ones(Float32, 1, size(M)[1])
    sum(v * M)
end
function dot_sum_full(M)
    v = CUDA.ones(Float32, 1, size(M)[1])
    reduced = v * M
    v2 = CUDA.ones(Float32, length(reduced))
    # Avoid scalar indexing warning
    reduce(+, reduced * v2)
end
function dot_sum_prealloc(v, M)
    sum(v * M)
end

using LinearAlgebra
@cuassert sum(M_rand) ≈ dot_sum(M_rand) ≈ dot_sum_full(M_rand) ≈ dot_sum_prealloc(v_ones, M_rand) ≈ dot(M_ones, M_rand)

using BenchmarkTools
@benchmark CUDA.@sync dot_sum(M_rand)
@benchmark CUDA.@sync dot_sum_full(M_rand)
@benchmark CUDA.@sync dot_sum_prealloc(v_ones, M_rand)

@benchmark CUDA.@sync sum(M_rand)
@benchmark CUDA.@sync dot(M_ones, M_rand)
@benchmark sum(M_cpu)

And the results

# @benchmark CUDA.@sync dot_sum(M_rand)
BenchmarkTools.Trial: 9906 samples with 1 evaluation.
 Range (min … max):  175.332 μs … 262.126 ms  ┊ GC (min … max): 0.00% … 9.63%
 Time  (median):     203.573 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   496.464 μs ±   2.749 ms  ┊ GC (mean ± σ):  0.51% ± 0.10%

  █▅▄▂▁                      ▁                      ▂▁          ▁
  █████▇▆▄▅▄▃▄▆▆▆█▇▇█▇▆▆▆▆▇▇███▇▅▄▄▄▁▁▄▄▅▄▃▃▅▃▄▆▇█▆▆████▇▆▅▃▄▆▆ █
  175 μs        Histogram: log(frequency) by time       3.79 ms <

 Memory estimate: 5.22 KiB, allocs estimate: 112.

# @benchmark CUDA.@sync dot_sum_full(M_rand)
BenchmarkTools.Trial: 9739 samples with 1 evaluation.
 Range (min … max):  178.480 μs … 8.318 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     199.840 μs             ┊ GC (median):    0.00%
 Time  (mean ± σ):   506.466 μs ± 1.029 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▄▂▁                           ▁▁                           ▁
  ████▇▆▅▆▇█▇▇▇▆▇███▇▆▅▅▅▅▅▄▄▅▆▅▇███▇▆▆▅▅▅▅▄▅▅▅▄▅▄▆▅▆▆▆▅▄▁▅▆▅ █
  178 μs       Histogram: log(frequency) by time       6.1 ms <

 Memory estimate: 3.62 KiB, allocs estimate: 95.

# @benchmark CUDA.@sync dot_sum_prealloc(v_ones, M_rand)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  167.418 μs … 157.410 ms  ┊ GC (min … max): 0.00% … 11.92%
 Time  (median):     182.596 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   441.800 μs ±   1.753 ms  ┊ GC (mean ± σ):  0.42% ±  0.12%

  █▄▂                                          ▁▁               ▁
  ████▆▅▅▄▄▁▄▁▅▄▇▇▆▇▇▇▆▆▆▆▆▇▇██▇▇▅▅▄▅▅▅▆▅▄▆▇▇▇▆██▇█▇▆▄▄▄▃▁▁▃▄▃▃ █
  167 μs        Histogram: log(frequency) by time       4.21 ms <

 Memory estimate: 5.06 KiB, allocs estimate: 106.

# @benchmark CUDA.@sync sum(M_rand)
BenchmarkTools.Trial: 3969 samples with 1 evaluation.
 Range (min … max):  492.841 μs … 12.605 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     713.582 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):     1.245 ms ±  1.029 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇▇█▇▆▄▁      ▁▁▁▂▂▁        ▁      ▁▁  ▁▂▃▃▂▂▁                ▁
  ███████▇█▆▄▆████████▇▆▇██████▇▆▇▆▇███████████▇▃▃▃▅▇█▆▅▄▃▆█▅▇ █
  493 μs        Histogram: log(frequency) by time      4.24 ms <

 Memory estimate: 4.59 KiB, allocs estimate: 82.

# @benchmark CUDA.@sync dot(M_ones, M_rand)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  248.969 μs …   4.300 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     262.834 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   473.993 μs ± 704.821 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▄▁▂           ▁ ▁                              ▁          ▂▁ ▁
  ████▄█▇▃▄▁▄▃▅▃▅█▇█▇▅▅▅▅▅▅▆▅▄▆▆▁▃▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▅█▇▄▅▄▃▅▄▅▅▆██ █
  249 μs        Histogram: log(frequency) by time       3.24 ms <

 Memory estimate: 32 bytes, allocs estimate: 2.

# @benchmark sum(M_cpu)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  104.205 μs …  1.220 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     198.371 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   229.274 μs ± 94.172 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▂█▇▇█▇▆▄▃▃▁                                                
  ▁▁▂█████████████▇▇▇▆▆▆▆▅▆▅▄▅▅▄▄▄▄▄▃▃▃▃▃▂▃▃▂▂▂▂▂▂▂▂▂▂▁▂▁▁▁▁▁▁ ▄
  104 μs          Histogram: frequency by time          539 μs <

 Memory estimate: 16 bytes, allocs estimate: 1.

So I'm getting CPU-like speeds with the implementations, which in my eyes is pretty good for a reduce operation on the GPU.
And compared to the existing sum it is ~2.5 faster for the mean, but with a higher standard deviation.