SuiteSparseGraphBLAS.jl and Krylov.jl

[amontoison]

SuiteSparseGraphBLAS.jl could provide efficient sparse matrix - dense vector products for Krylov.jl but we can't easily use it here for the following reasons:

• The type of vector that we use to store variables in a KrylovSolver must be a subtype of DenseVector and GBVector is an AbstractSparseVector (even if GBVector is dense vector internally);
• BLAS1 routines (dot, axpy!, axpby!, norm, etc...) are not all supported by the GBVector type.

However, I don't think that it's relevant to store variables inside the Krylov workspaces with GBVector.
We want to use it because mul!(y, GBMatrix, x) is only working if x and y are GBVector.
Is it possible to make it working with x::Vector and y::Vector or not?

If it's not possible, we could just define an operator that take as input a SparseMatrixCSC and create a structure with a GBMatrix and two auxiliary GBVector such that we can do:

function mul!(y::Vector{T}, A::OpGraphBLAS{T}, x::Vector{T})
  A.gbx .= x
  mul!(A.gby, A.gbA, A.gpx)
  y .= A.gpy
end

In the documentation, I did some benchmarks with MKLSparse, the default mul! of SparseArrays and a multhreaded mul! implemented julia.
It will be great to add some benchmarks with SuiteSparseGraphBLAS.jl and how to use it in Krylov.jl if we solve the previous problems.

@wimmerer

[amontoison]

We could reuse the following script for the benchmarks. It uses real life matrices from the SuiteSparseMatrixCollection.
I generated the graphic in the documentation with it.

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

[rayegun]

I'll work on benchmarking. I definitely have y = A::GBMatrix * b::Vector implemented, we can then unpack y into a dense vector in O(1) time. I do not currently have mul!(y::Vector, A::GBMatrix, b::Vector). That may not be possible at this juncture since I believe "shallow" arrays cannot be output.

Let me run some stability tests and check with Tim Davis though. I may be able to "cheat" here a bit and tell SuiteSparse:GraphBLAS that the output vector is not "shallow". This would be doable, but slightly dangerous. If we pass an accumulator (effectively y::Vector .+= A::GBMatrix * b::Vector then I can do this reasonably safely), but I'm not sure what the runtime looks like in that case.

(For clarity, GraphBLAS is very particular about the patterns of vectors and matrices. If A has an empty row then y will have an empty index. It would no longer be a dense vector. This poses a problem from a memory perspective, GraphBLAS might try to free y and rebuild it as bitmap, and from a performance perspective for your usecase).