Skip to content

Instantly share code, notes, and snippets.

@carstenbauer

carstenbauer/bench.jl

Created August 9, 2021 18:22
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save carstenbauer/6cbc3d2c32422472a579bf61b5d62015 to your computer and use it in GitHub Desktop.
Save carstenbauer/6cbc3d2c32422472a579bf61b5d62015 to your computer and use it in GitHub Desktop.
Download ZIP
SAXPY benchmark CUDA.jl (broadcasting)
Raw
bench.jl
print("Loading modules...")
using BenchmarkTools
using CUDA
using DataFrames
using CSV
using Printf
println("done!")
const a = 3.1415f0
function saxpy_julia!(z,a,x,y)
z .= a .* x .+ y
return z
end
df = DataFrame(n=Int[], var"GFLOPS/s"=Float64[], var"GB/s"=Float64[], time=Float64[])
println("Running measurements...")
for i in 8:128
n = 1024*1024*i
x, y, z = CUDA.ones(n), CUDA.ones(n), CUDA.zeros(n)
t_saxpy = @belapsed CUDA.@sync saxpy_julia!($z,$a,$x,$y)
gflops = 2.0 * n * (1000)^(-3) / t_saxpy
bandwidth = 3.0 * sizeof(Float32) * n * (1000)^(-3) / t_saxpy
@printf("saxpy (julia): n= %12d %7.3f GFLOP/s %7.3f GB/s %7.3f s\n",
n, gflops, bandwidth, t_saxpy);
push!(df, (n, gflops, bandwidth, t_saxpy))
flush(stdout) # force immediate printing to .out logfile
# free memory (to be safe)
x, y, z = nothing, nothing, nothing
GC.gc(true)
end
println("done!")
print("Writing results to disk...")
CSV.write("bench_results.csv", df)
println("done!")
@carstenbauer
Copy link
Author

Some numbers from the Julia benchmark above:

saxpy (julia): n=      8388608 105.363 GFLOP/s 632.180 GB/s   0.000 s
saxpy (julia): n=      9437184 104.407 GFLOP/s 626.445 GB/s   0.000 s
saxpy (julia): n=     10485760 122.029 GFLOP/s 732.173 GB/s   0.000 s
saxpy (julia): n=     11534336 106.108 GFLOP/s 636.649 GB/s   0.000 s
saxpy (julia): n=     12582912 106.191 GFLOP/s 637.144 GB/s   0.000 s
saxpy (julia): n=     13631488 106.750 GFLOP/s 640.497 GB/s   0.000 s
saxpy (julia): n=     14680064 106.272 GFLOP/s 637.633 GB/s   0.000 s
saxpy (julia): n=     15728640 125.742 GFLOP/s 754.453 GB/s   0.000 s
saxpy (julia): n=     16777216 107.945 GFLOP/s 647.671 GB/s   0.000 s
saxpy (julia): n=     17825792 107.061 GFLOP/s 642.367 GB/s   0.000 s
saxpy (julia): n=     18874368 107.441 GFLOP/s 644.648 GB/s   0.000 s
saxpy (julia): n=     19922944 107.643 GFLOP/s 645.860 GB/s   0.000 s
saxpy (julia): n=     20971520 126.920 GFLOP/s 761.521 GB/s   0.000 s
saxpy (julia): n=     22020096 107.749 GFLOP/s 646.496 GB/s   0.000 s
saxpy (julia): n=     23068672 107.915 GFLOP/s 647.490 GB/s   0.000 s
saxpy (julia): n=     24117248 107.621 GFLOP/s 645.728 GB/s   0.000 s
saxpy (julia): n=     25165824 107.616 GFLOP/s 645.697 GB/s   0.000 s
saxpy (julia): n=     26214400 127.738 GFLOP/s 766.428 GB/s   0.000 s
saxpy (julia): n=     27262976 107.963 GFLOP/s 647.781 GB/s   0.001 s
saxpy (julia): n=     28311552 108.033 GFLOP/s 648.199 GB/s   0.001 s
saxpy (julia): n=     29360128 107.715 GFLOP/s 646.288 GB/s   0.001 s
saxpy (julia): n=     30408704 107.333 GFLOP/s 643.995 GB/s   0.001 s

The same from a C / CUBLAS version which seems to be around 20% faster in most cases:

saxpy (cuBLAS): n=      8388608 123.188 GFLOP/s 739.128 GB/s
saxpy (cuBLAS): n=      9437184 123.705 GFLOP/s 742.228 GB/s
saxpy (cuBLAS): n=     10485760 124.878 GFLOP/s 749.268 GB/s
saxpy (cuBLAS): n=     11534336 125.156 GFLOP/s 750.933 GB/s
saxpy (cuBLAS): n=     12582912 127.337 GFLOP/s 764.021 GB/s
saxpy (cuBLAS): n=     13631488 126.781 GFLOP/s 760.686 GB/s
saxpy (cuBLAS): n=     14680064 127.431 GFLOP/s 764.587 GB/s
saxpy (cuBLAS): n=     15728640 127.469 GFLOP/s 764.813 GB/s
saxpy (cuBLAS): n=     16777216 127.502 GFLOP/s 765.012 GB/s
saxpy (cuBLAS): n=     17825792 128.472 GFLOP/s 770.834 GB/s
saxpy (cuBLAS): n=     18874368 128.000 GFLOP/s 768.000 GB/s
saxpy (cuBLAS): n=     19922944 129.276 GFLOP/s 775.654 GB/s
saxpy (cuBLAS): n=     20971520 129.620 GFLOP/s 777.722 GB/s
saxpy (cuBLAS): n=     22020096 129.542 GFLOP/s 777.253 GB/s
saxpy (cuBLAS): n=     23068672 130.220 GFLOP/s 781.318 GB/s
saxpy (cuBLAS): n=     24117248 129.052 GFLOP/s 774.312 GB/s
saxpy (cuBLAS): n=     25165824 130.032 GFLOP/s 780.190 GB/s
saxpy (cuBLAS): n=     26214400 130.612 GFLOP/s 783.673 GB/s
saxpy (cuBLAS): n=     27262976 129.557 GFLOP/s 777.343 GB/s
saxpy (cuBLAS): n=     28311552 129.803 GFLOP/s 778.817 GB/s
saxpy (cuBLAS): n=     29360128 130.922 GFLOP/s 785.534 GB/s
saxpy (cuBLAS): n=     30408704 130.819 GFLOP/s 784.916 GB/s

@carstenbauer
Copy link
Author

What is particularly curious is that for some values of n, we seem to essentially match the C performance in a spike-like fashion…. For example, see n = 15728640 where we have 764.813 GB/s for C and 754.453 GB/s for Julia.

Interestingly, the positive spikes seem to appear with a constant frequency: every 5th run is fast?!?

@carstenbauer
Copy link
Author

carstenbauer commented Aug 9, 2021 •
edited

Alright, using an explicit kernel closes the gap to C.

print("Loading modules...")
using BenchmarkTools
using CUDA
using DataFrames
using CSV
using Printf
println("done!")

const a = 3.1415f0

function saxpy_gpu_kernel!(z, a, x, y)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(z)
        @inbounds z[i] = a * x[i] + y[i]
    end
    return nothing
end

function saxpy_gpu!(z, a, x, y; nthreads, nblocks)
    CUDA.@sync @cuda(
        threads = nthreads,
        blocks = nblocks,
        saxpy_gpu_kernel!(z, a, x, y)
    )
end

# query how many threads per block are available on the GPU
nthreads = CUDA.attribute(
    device(),
    CUDA.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK
)

df = DataFrame(n=Int[], var"GFLOPS/s"=Float64[], var"GB/s"=Float64[], time=Float64[])

println("Running measurements...")
for i in 8:20 # 128
    n = 1024*1024*i
    x, y, z = CUDA.ones(n), CUDA.ones(n), CUDA.zeros(n)
    # compute how many blocks we need to use for the given `dim`
    nblocks = cld(n, nthreads)

    t_saxpy = @belapsed saxpy_gpu!($z, $a, $x, $y; nblocks=$nblocks, nthreads=$nthreads)

    gflops = 2.0 * n * (1000)^(-3) / t_saxpy
    bandwidth = 3.0 * sizeof(Float32) * n * (1000)^(-3) / t_saxpy
    @printf("saxpy (julia): n= %12d %7.3f GFLOP/s %7.3f GB/s %7.3f s\n",
                 n, gflops, bandwidth, t_saxpy);
    push!(df, (n, gflops, bandwidth, t_saxpy))


    flush(stdout) # force immediate printing to .out logfile
    # free memory (to be safe)
    x, y, z = nothing, nothing, nothing
    GC.gc(true)
end
println("done!")
print("Writing results to disk...")
CSV.write("bench_results.csv", df)
println("done!")

Output:

saxpy (julia): n=      8388608 125.526 GFLOP/s 753.158 GB/s   0.000 s
saxpy (julia): n=      9437184 124.185 GFLOP/s 745.109 GB/s   0.000 s
saxpy (julia): n=     10485760 126.514 GFLOP/s 759.086 GB/s   0.000 s
saxpy (julia): n=     11534336 127.279 GFLOP/s 763.674 GB/s   0.000 s
saxpy (julia): n=     12582912 126.933 GFLOP/s 761.601 GB/s   0.000 s
saxpy (julia): n=     13631488 127.728 GFLOP/s 766.366 GB/s   0.000 s
saxpy (julia): n=     14680064 128.188 GFLOP/s 769.127 GB/s   0.000 s
saxpy (julia): n=     15728640 128.661 GFLOP/s 771.967 GB/s   0.000 s
saxpy (julia): n=     16777216 128.307 GFLOP/s 769.844 GB/s   0.000 s
saxpy (julia): n=     17825792 129.633 GFLOP/s 777.799 GB/s   0.000 s
saxpy (julia): n=     18874368 129.656 GFLOP/s 777.934 GB/s   0.000 s
saxpy (julia): n=     19922944 129.482 GFLOP/s 776.895 GB/s   0.000 s
saxpy (julia): n=     20971520 129.418 GFLOP/s 776.510 GB/s   0.000 s

@carstenbauer
Copy link
Author

carstenbauer commented Aug 10, 2021 •
edited

Benchmark results:

saxpy_julia

Benchmark code (without the plotting part):

print("Loading modules...");
flush(stdout);
using CUDA
using DataFrames
using CSV
using Printf
println("done!");
flush(stdout);

# kernel definition
function saxpy_gpu_kernel!(z, a, x, y)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(z)
        @inbounds z[i] = a * x[i] + y[i]
    end
    return nothing
end

# calling the kernel
function saxpy_kernel!(z, a, x, y; nthreads, nblocks)
    @cuda(threads = nthreads, blocks = nblocks, saxpy_gpu_kernel!(z, a, x, y))
    return nothing
end

# high-level broadcasting version
function saxpy_broadcasting!(z, a, x, y; nthreads, nblocks)
    z .= a .* x .+ y
    return z
end

function run_benchmarks(saxpy!)
    # for storing the results
    df = DataFrame(; n=Int[], var"GFLOP/s"=Float64[], var"GB/s"=Float64[], kind=String[])
    fname = string(saxpy!)

    # query how many threads per block are available on the GPU
    nthreads = CUDA.attribute(device(), CUDA.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)

    for i in 8:128
        # vector length:
        # always a multiple of nthreads per block (1024)
        # such that blocks are fully busy (no remainder).
        n = 1024 * 1024 * i
        # arbitrary constant
        a = 3.1415f0
        # allocate GPU memory
        x, y, z = CUDA.ones(n), CUDA.ones(n), CUDA.zeros(n)
        # compute how many blocks we need to use for the given `dim`
        nblocks = cld(n, nthreads)
        # benchmark: minimum time of 10 trials
        t_saxpy = 1e6
        for j in 1:10
            t = CUDA.@elapsed saxpy!(z, a, x, y; nthreads, nblocks)
            t_saxpy = min(t_saxpy, t)
        end
        # print and save results
        flops = 2.0 * n * (1000)^(-3) / t_saxpy
        bandwidth = 3.0 * sizeof(Float32) * n * (1000)^(-3) / t_saxpy
        @printf(
            "%s (julia): n= %12d %7.3f GFLOP/s %7.3f GB/s \n",
            fname,
            n,
            flops,
            bandwidth,
        )
        # force immediate printing to .out logfile
        flush(stdout)
        push!(df, (n, flops, bandwidth, fname))
        # explicitly free memory (to be safe)
        x, y, z = nothing, nothing, nothing
        GC.gc(true)
    end
    return df
end

println("Running measurements (kernel)...");
flush(stdout);
df_kernel = run_benchmarks(saxpy_kernel!)
df_broadcasting = run_benchmarks(saxpy_broadcasting!)
df = vcat(df_kernel, df_broadcasting)

println("done!")
print("Writing results to disk...")
CSV.write("bench_results.csv", df)
println("done!")

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment