Interesting timings

[Cvikli]

Dear CUDA developers,
I just noticed something strange. What can cause these timings? (Why does it decrease over time?)

This basic code timings:

using BenchmarkTools
using CUDA
N= 500; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
N= 1000; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
N= 2000; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
N= 4000; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
N= 6000; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
N= 8000; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
N= 10000; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
N= 12000; A = CUDA.randn(Float32, N, N); @btime sum($A, dims=2)
;

Results:

  7.228 μs (64 allocations: 2.66 KiB)
  8.189 μs (82 allocations: 3.20 KiB)
  8.270 μs (83 allocations: 3.23 KiB)
  8.307 μs (83 allocations: 3.23 KiB)
  6.517 μs (43 allocations: 1.89 KiB)
  6.505 μs (43 allocations: 1.89 KiB)
  6.379 μs (43 allocations: 1.89 KiB)
  6.635 μs (43 allocations: 1.89 KiB)

Sounds like an interesting anomaly.

versioninfo()
Julia Version 1.7.0
Commit 3bf9d17731 (2021-11-30 12:12 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-12.0.1 (ORCJIT, skylake)

CUDA.versioninfo()
CUDA toolkit 11.4.1, artifact installation
CUDA driver 11.5.0
NVIDIA driver 495.29.5

Libraries:

• CUBLAS: 11.5.4
• CURAND: 10.2.5
• CUFFT: 10.5.1
• CUSOLVER: 11.2.0
• CUSPARSE: 11.6.0
• CUPTI: 14.0.0
• NVML: 11.0.0+495.29.5
• CUDNN: 8.20.2 (for CUDA 11.4.0)
• CUTENSOR: 1.3.0 (for CUDA 11.2.0)

Toolchain:

• Julia: 1.7.0
• LLVM: 12.0.1
• PTX ISA support: 3.2, 4.0, 4.1, 4.2, 4.3, 5.0, 6.0, 6.1, 6.3, 6.4, 6.5, 7.0
• Device capability support: sm_35, sm_37, sm_50, sm_52, sm_53, sm_60, sm_61, sm_62, sm_70, sm_72, sm_75, sm_80

1 device:
0: NVIDIA GeForce GTX 1050 (sm_61, 310.438 MiB / 3.945 GiB available)

[Cvikli]

This sounds interesting, using @benchmark tool I get this:

using BenchmarkTools
using CUDA
N = 500; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2))
N = 1000; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2))
N = 2000; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2))
N = 4000; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2))
N = 6000; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2))
N = 8000; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2))
N = 10000; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2))
N = 12000; A = CUDA.randn(Float32, N, N); display(@benchmark sum($A, dims=2));

Results:

BenchmarkTools.Trial: 
  memory estimate:  2.66 KiB
  allocs estimate:  64
  --------------
  minimum time:     7.946 μs (0.00% GC)
  median time:      293.935 μs (0.00% GC)
  mean time:        294.296 μs (0.31% GC)
  maximum time:     13.970 ms (37.28% GC)
  --------------
  samples:          5647
  evals/sample:     3
BenchmarkTools.Trial: 
  memory estimate:  3.20 KiB
  allocs estimate:  82
  --------------
  minimum time:     8.984 μs (0.00% GC)
  median time:      766.087 μs (0.00% GC)
  mean time:        767.465 μs (0.00% GC)
  maximum time:     2.916 ms (0.00% GC)
  --------------
  samples:          6487
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  3.23 KiB
  allocs estimate:  83
  --------------
  minimum time:     9.058 μs (0.00% GC)
  median time:      2.143 ms (0.00% GC)
  mean time:        2.081 ms (0.00% GC)
  maximum time:     3.195 ms (0.00% GC)
  --------------
  samples:          801
  evals/sample:     3
BenchmarkTools.Trial: 
  memory estimate:  3.23 KiB
  allocs estimate:  83
  --------------
  minimum time:     10.125 μs (0.00% GC)
  median time:      8.150 ms (0.00% GC)
  mean time:        7.454 ms (0.00% GC)
  maximum time:     8.962 ms (0.00% GC)
  --------------
  samples:          224
  evals/sample:     3
BenchmarkTools.Trial: 
  memory estimate:  1.89 KiB
  allocs estimate:  43
  --------------
  minimum time:     6.985 μs (0.00% GC)
  median time:      2.979 ms (0.00% GC)
  mean time:        2.764 ms (0.00% GC)
  maximum time:     3.929 ms (0.00% GC)
  --------------
  samples:          362
  evals/sample:     5
BenchmarkTools.Trial: 
  memory estimate:  1.89 KiB
  allocs estimate:  43
  --------------
  minimum time:     11.775 μs (0.00% GC)
  median time:      7.150 ms (0.00% GC)
  mean time:        6.414 ms (0.00% GC)
  maximum time:     8.047 ms (0.00% GC)
  --------------
  samples:          195
  evals/sample:     4
BenchmarkTools.Trial: 
  memory estimate:  1.89 KiB
  allocs estimate:  43
  --------------
  minimum time:     7.156 μs (0.00% GC)
  median time:      9.378 ms (0.00% GC)
  mean time:        7.958 ms (0.00% GC)
  maximum time:     10.279 ms (0.00% GC)
  --------------
  samples:          158
  evals/sample:     4
BenchmarkTools.Trial: 
  memory estimate:  1.89 KiB
  allocs estimate:  43
  --------------
  minimum time:     7.674 μs (0.00% GC)
  median time:      11.546 ms (0.00% GC)
  mean time:        10.846 ms (0.00% GC)
  maximum time:     13.873 ms (0.00% GC)
  --------------
  samples:          116
  evals/sample:     4

So basically about 1000 times slower version. Not sure what does this mean?
Mean times of our measures respectively:

500:         294us
1000:        767us
2000:      2.081ms
4000:      7.454ms
6000:      2.764ms
8000:      6.414ms
10000:     7.958ms
12000:    10.846ms

Also that 6000 case was outstanding without reason. What could be the explain behind these?

[Cvikli]

Interesting note also:

If I write a custom hardcoded, preallocated sum function, it performs pretty well and hold us with my requirements:

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

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

Results:

BenchmarkTools.Trial: 
  memory estimate:  6.08 KiB
  allocs estimate:  179
  --------------
  minimum time:     249.119 μs (0.00% GC)
  median time:      253.769 μs (0.00% GC)
  mean time:        274.501 μs (0.63% GC)
  maximum time:     5.091 ms (97.14% GC)
  --------------
  samples:          10000
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  4.58 KiB
  allocs estimate:  133
  --------------
  minimum time:     501.477 μs (0.00% GC)
  median time:      505.177 μs (0.00% GC)
  mean time:        539.152 μs (0.53% GC)
  maximum time:     4.519 ms (90.12% GC)
  --------------
  samples:          9232
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  5.52 KiB
  allocs estimate:  163
  --------------
  minimum time:     1.028 ms (0.00% GC)
  median time:      1.034 ms (0.00% GC)
  mean time:        1.095 ms (0.51% GC)
  maximum time:     4.830 ms (78.14% GC)
  --------------
  samples:          4554
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  27.33 KiB
  allocs estimate:  861
  --------------
  minimum time:     2.092 ms (0.00% GC)
  median time:      2.110 ms (0.00% GC)
  mean time:        2.280 ms (0.46% GC)
  maximum time:     7.295 ms (67.63% GC)
  --------------
  samples:          2189
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  4.83 KiB
  allocs estimate:  141
  --------------
  minimum time:     3.380 ms (0.00% GC)
  median time:      3.417 ms (0.00% GC)
  mean time:        3.627 ms (0.45% GC)
  maximum time:     8.932 ms (0.00% GC)
  --------------
  samples:          1377
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  153.52 KiB
  allocs estimate:  4899
  --------------
  minimum time:     4.526 ms (0.00% GC)
  median time:      4.554 ms (0.00% GC)
  mean time:        4.796 ms (0.43% GC)
  maximum time:     8.382 ms (49.48% GC)
  --------------
  samples:          1042
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  189.89 KiB
  allocs estimate:  6063
  --------------
  minimum time:     6.448 ms (0.00% GC)
  median time:      6.492 ms (0.00% GC)
  mean time:        6.829 ms (0.41% GC)
  maximum time:     9.916 ms (34.97% GC)
  --------------
  samples:          732
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  620.14 KiB
  allocs estimate:  19831
  --------------
  minimum time:     13.025 ms (0.00% GC)
  median time:      13.521 ms (0.00% GC)
  mean time:        13.816 ms (0.41% GC)
  maximum time:     16.263 ms (21.11% GC)
  --------------
  samples:          362
  evals/sample:     1

Respectively compared to previous version:

500:     253us  - 294us
1000:    505us  - 767us
2000:  1.034ms - 2.081ms
4000:  2.110ms - 7.454ms
6000:  3.417ms - 2.764ms
8000:  4.554ms - 6.414ms
10000: 6.492ms - 7.958ms
12000:13.521ms - 10.846ms

Thes diffs are really strange... not sure what is going on under the hood that provide these differences. Any idea on this?

[Cvikli]

I just realised I forgot to put @inbounds and CUDA.Const(...) to the appropriate tags.
The latest one looked like this:

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

With this respectively:

Case | mysum with @inbounds + CUDA.Const(...)  | mysum with @inbounds |   mysum  | std sum
500:                     50us                  -        152us         -    274us -   294us
1000:                    93us                  -        284us         -    539us -   767us
2000:                   216us                  -        611us         -  1.095ms - 2.081ms
4000:                   697us                  -      1.309ms         -  2.280ms - 7.454ms
6000:                 1.548ms                  -      2.236ms         -  3.627ms - 2.764ms
8000:                 2.795ms                  -      3.298ms         -  4.796ms - 6.414ms
10000:                4.922ms                  -      5.031ms         -  6.829ms - 7.958ms
12000:                7.028ms                  -      9.312ms         - 13.816ms - 10.846ms

Each of them looks really interesting. Sounds like a 2-10x speedup with the custom mysum aritmetic.
Also the question with the @btime sounds like a bad measurement, as those values differ so much. (Maybe there were problem with syncronization?) If so the only question is that, can this speedup be real compared to the standard code?

[maleadt]

You're forgetting to synchronize the GPU. Always do @benchmark CUDA.@sync. In your first example, sum(; dims=...) returns an array, so can complete asynchronously (as opposed to a sum returning a scalar), necessitating synchronization.

[Cvikli]

Indeed, I missed the CUDA.@sync . Interestingly, the results didn't really change. I guess the @benchmark does it somehow indirectly... but the @btime definitely needs that I see.

Anyway, it is interesting to see that there is this 2-10x speedup.

[maleadt]

Interestingly, the results didn't really change.

Then you are doing something wrong:

julia> N= 500; A = CUDA.randn(Float32, N, N); @btime CUDA.@sync sum($A, dims=2);
  72.071 μs (40 allocations: 1.94 KiB)

julia> N= 1000; A = CUDA.randn(Float32, N, N); @btime CUDA.@sync sum($A, dims=2);
  147.723 μs (41 allocations: 1.95 KiB)

julia> N= 2000; A = CUDA.randn(Float32, N, N); @btime CUDA.@sync sum($A, dims=2);
  342.341 μs (41 allocations: 1.95 KiB)

There can't a speed-up when the input size changes, that doesn't make any sense. If it clocks in at only around 10us, you're measuring kernel launch time.

[Cvikli]

I just saw you comment. I think it worth to look at, my latest measures already isn't 10us and still there is a speedup. But please ignore me and continue your work if I am too much on this CUDA topic. With your guide I got in the other issue, I will be able to write my code in 2 kernel and won't bother with fusing the whole into one kernel.

So I run this:

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))

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)

Rel. Diff/speedup is about: 2-10x in certain cases.

I ran this on my workstation with 1080 TI. (That is why the timings changed.)

[maleadt]

If you think our sum implementation is slow, you can open a different issue.

[Cvikli]

Okay!