Benchmarks MatVec

mklsparse.jl

using BenchmarkTools, SuiteSparseMatrixCollection, MatrixMarket
using SparseArrays, LinearAlgebra, Printf, JLD, Base.Threads

ssmc = ssmc_db(verbose=false)
dataset = ssmc[(ssmc.binary .== false) .& (10000 .≤ ssmc.nrows .≤ 20000) .& (ssmc.numerical_symmetry .== 1.0), :]
# dataset = ssmc[(ssmc.binary .== false) .& (1000 .≤ ssmc.nrows .≤ 10000), :]
paths = fetch_ssmc(dataset, format="MM")
names = dataset[!,:name]

function threaded_mul!(y::Vector{T}, A::SparseMatrixCSC{T}, x::Vector{T}) where T <: Number
  A.m == A.n || error("A is not a square matrix!")
  @threads for i = 1 : A.n
    tmp = zero(T)
    @inbounds for j = A.colptr[i] : (A.colptr[i+1] - 1)
      tmp += A.nzval[j] * x[A.rowval[j]]
    end
    @inbounds y[i] = tmp
  end
  return y
end

# nb_pbs = length(paths)  # 796 matrices...
nb_pbs = 15
benchmark_sparse = true
benchmark_mkl = true
benchmark_julia = true
time_sparse = zeros(nb_pbs)
time_mkl = zeros(nb_pbs)
time_julia = zeros(nb_pbs)

if benchmark_sparse
  for i = 1 : nb_pbs
    name = dataset[!,:name][i]
    path = paths[i]
    @printf("SparseArrays -- %3d / %3d -- %s\n", i, nb_pbs, name)
    A = Float64.(MatrixMarket.mmread(path * "/$name.mtx"))
    n, m = size(A)
    x = rand(m)
    y = rand(n)
    time_sparse[i] = @belapsed for k in 1:1000 mul!($y, $A, $x) end
  end
  save("time_sparse.jld", "time_sparse", time_sparse)
else
  time_sparse = load("time_sparse.jld", "time_sparse")
end

using MKLSparse

if benchmark_mkl
  for i = 1 : nb_pbs
    name = names[i]
    path = paths[i]
    @printf("MKLSparse -- %3d / %3d -- %s\n", i, nb_pbs, name)
    A = MatrixMarket.mmread(path * "/$name.mtx")
    n, m = size(A)
    x = rand(m)
    y = rand(n)
    time_mkl[i] = @belapsed for k in 1:1000 mul!($y, $A, $x) end
  end
  save("time_mkl.jld", "time_mkl", time_mkl)
else
  time_mkl = load("time_mkl.jld", "time_mkl")
end

if benchmark_julia
  for i = 1 : nb_pbs
    name = dataset[!,:name][i]
    path = paths[i]
    @printf("Julia -- %3d / %3d -- %s\n", i, nb_pbs, name)
    A = Float64.(MatrixMarket.mmread(path * "/$name.mtx"))
    n, m = size(A)
    x = rand(m)
    y = rand(n)
    time_julia[i] = @belapsed for k in 1:1000 threaded_mul!($y, $A, $x) end
  end
  save("time_julia.jld", "time_julia", time_julia)
else
  time_julia = load("time_julia.jld", "time_julia")
end

using PlotlyJS

p = plot([
  bar(name="mul! SparseArrays", x=names[1:nb_pbs], y=time_sparse),
  bar(name="mul! MKLSparse", x=names[1:nb_pbs], y=time_mkl),
  bar(name="threaded_mul!", x=names[1:nb_pbs], y=time_julia),
  ],
)
relayout!(p, barmode="group")
savefig(p, "julia_vs_mkl_eps", format="eps")
savefig(p, "julia_vs_mkl_svg", format="svg")