Roadmap

[gdalle]

Tip

Active development is happening over in https://github.com/AntoineBut/GPUGraphs.jl for the time being.

24/02/2025

Practical

• Implement sparse matrix storage with KernelAbstractions.jl
• Implement matrix-vector multiplication (LinearAlgebra.mul!)
• Implement BFS using abstractions from GPUArrays.jl whenever possible
• Add tests on simple graphs
• (Possibly) Include conversion utilities to native graph formats (e.g. CuSparse)

Theoretical

• Start reviewing graph algorithms on GPU and GraphBLAS specification
• Start thinking about a more complicated algorithm than BFS as an end goal (e.g. "community detection", "stochastic block model")
• Read up on Graph Neural Networks in PyTorch and JAX: which libraries? what do they use to encode graphs? where are they fast/slow?
• Same for GraphNeuralNetworks.jl

[AntoineBut]

[gdalle]

04/03/2025

Objectives of the project

1. Implement a subset of the GraphBLAS API with GPU backend (inspired by SuiteSparseGraphBLAS.jl)
2. Implement a few algorithms to test it, find cases where GPUs are useful

Next few weeks

1. Create data structures
2. Add basic API functions

For next week: have a basic package with a CI test suite on AntoineBut/GPUGraphs.jl

People to meet

• Raye Kimmerer
• Willow Ahrens

[AntoineBut]

Basic functionality priority order:

Priority 1:

• Sparse Matrix storage structure (GPUGraphMatrix)
• Basic constructors for GPUGraphMatrix
• Sparse Matrix x Sparse Matrix multiplication
• Sparse Matrix x Dense vector multiplication
• Element-wise Sparse Matrix and vector map (. brodacasting)
• Element-wise Sparse Vector add and Mul

Priority 2:

• (Lazy) Transpose operator
• Array view (readonly shallow copy)
• Sparse Matrix x Sparse vector multiplication (not sure if this is relevant)
• Some more elaborate constructors / conversion from Graphs.jl and SuiteSparse:GraphBLAS.jl objects

Priority 3:

• Masking for GraphBLAS operations
• All remaining operators on https://graphblas.juliasparse.org/dev/#GraphBLAS-Operations

@gdalle I'll get started on the elements with "priority 1" for the next weeks (alongside the basic test / CI pipeline), any feedback on this first draft on priorities is of course welcome :))

[gdalle]

Looks good overall, although I would probably switch 2 and 3.

[gdalle]

11/03/2025

Work done

• Basic sparse types based on GPU arrays
• Basic tests

Todo

• Transfer back from GPU to CPU
• Same behavior for getindex to allow/disallow scalar access with a macro
• First operators

[gdalle]

18/03/2025

Work done

• Make KA work on CPU backend
• Struggle with first kernel and scalar indexing (see below)
• KernelAbstractions doesn't accept anything but primitive types

Kernel explanation

using KernelAbstractions

@kernel function spmv!(y::AbstractVector, M::AbstractSparseMatrix, x::AbstractVector)
    # option 1
    #=
    blockIdx = @index(Group)
    blockDimx = @groupsize()
    threadIdx = @index(Local)
    row = blockIdx * blockDimx + threadIdx # modulo 1-indexing
    =#
    # option 2
    row = @index(Global)
    for i in M.rowptr[row]:M.rowptr[row+1]
        col = M.colval[i]
        value = M.nzval[i]
        sum += x[col] * value
    end
    y[row] += sum
end

# alloc y, x, M
nb_threads_by_block = 32 # (can be automatically set)
grid_size = (size(M, 1),)  # not necessarily x32
spmv_configured! = spmv!(backend, nb_threads_by_block)
spmv_configured!(y, M, x; ndrange=grid_size)

Custom operation example

using KernelAbstractions

struct Operation{O}
    op::O
end

@kernel function mul2_kernel(A, op::Operation)
    I = @index(Global)
    A[I] = op.op(2, A[I])
end

dev = CPU()
A = ones(1024, 1024)
ev = mul2_kernel(dev, 64)(A, Operation(*), ndrange=size(A))
synchronize(dev)
all(A .== 2.0)

[gdalle]

25/03/2025

Work done

• Matrix-vector multiplication works
• KernelAbstractions on CPU fails
• Benchmark against CPU version
• In some cases the GPU version is faster on bigger matrices than on smaller matrices, when Int32 indexing fails

Todo

• Implement padded version of CSR format
• Benchmark with $p \propto 1/n$ inside sprand instead of $p$ constant
• Initialize at zero
• Reuse SuiteSparseGraphBLAS API for accumulation (optional kwarg to mul!)
• Instead of SemiRing, define something which adapts to element types
  • monoid_neutral(::Type{T}, ::typeof(+)) = zero(T)
  • monoid_neutral(::Type{T}, ::typeof(*)) = one(T)
  • monoid_neutral(::Type{T}, ::typeof(max)) = typemin(T)

[AntoineBut]

Looking back at

struct Operation{O}
  op::O
end
[...]
ev = mul2_kernel(dev, 64)(A, Operation(*), ndrange=size(A))

from last week, is there any particular reason to wrap the *::Function in a custom Operation Object? I currently just don't see why we would need this.

[gdalle]

No, I was just following your initial draft, but as discussed yesterday, directly passing things like * or max is completely fine. As long as we have specialized methods for the neutral elements.

[gdalle]

01/04/2025

Work done

• Implemented CSR format with padding (ELL)
• Compared performance on various matrices from the SuiteSparse collection
• ELL not efficient when there is only one non-zero in a row
• Explore profiling of Apple GPU (contribute docs to Metal.jl?)
• NLP-KKT case might be favorable for naive implem because fixed degree means strided memory access, even if not coalesced they may be quite good
• Behavior is different on Bool vs Float32
• Theory: memory accesses in b are limiting

Todo

• Start taking notes about all of this
• Test theory about b: replace it with constant values.
• Try tiling or find some literature about less basic matrix-vector mult

[gdalle]

08/04/2025

Work done

• Wrote little to no code
• Benchmarked to confirm that memory accesses inside b are a big problem
• Read some papers, everyones tries to optimize memory access in A but no one cares about b

Ideas

• Process segments of b one after the other by cutting the matrix into vertical bands to prevent it from getting too large
• Will work better for dense matrices
• Problem is symmetric for A*b = c between b and c: switching A to CSC doesn't save us
• Slicing A into k vertical bands multiplies the amount of padding in ELL (hence the matrix density) by k
• Sorting rows by density (JDS format) can be useful
• Even in BFS, the vector b doesn't remain sparse for small-world graphs
• If b is sparse and fits in shared memory, bingo

Todo

• Implement CSC for the matrix
• Implement matrix-vector product where the vector is sparse

[gdalle]

15/04/2025

Work done

• Slowed down by exams
• Coded CSC storage
• Struggling with numerical tests up to floating point precision
• Look up atomic operations

Ideas

• Check that norm of the difference is small (that's what isapprox does?)
• How to accumulate into c for CSC? Check out low-level interface of Atomix.jl
• Beware of functions passed as arguments, see https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
• Type promotion didn't work with SuiteSparseGraphBLAS, will it work here?
• Will masking work here?

Goals

• Give up on further optimizations with bands, not very promising
• Add elementwise stuff
• Aim for a working BFS in two weeks with a dense vector
  • at least distances
  • ideally parents
• Adapt index-related operators like secondi by extending add(x, y) into add(x, y, some_relevant_indices)

[AntoineBut]

29/04/2025

Work done

• Made custom functions to use as addition, multiplication and accumulator in SpMV kernel. These functions all take the same arguments : (x, y, row_x, col_x, row_y, col_y) which allows for functions like second_i to be used as well.
• Implemented a first BFS returning distances from the source node.
• Started working on masking feature, encountered some difficulties.
• Looked briefly into the question of type promotion : looks like it is working nicely. I just need to rewrite some stuff to enable it. Currently we have a single type parameter for input/output vectors and matrix elements:

# current version : 
function gpu_spmv!(
    C::AV,
    A::SparseGPUMatrixCSR{Tv,Ti},
    B::AV;
    mul::Function = GPUGraphs_mul,
    add::Function = GPUGraphs_add,
    accum::Function = GPUGraphs_second,
    mask::Union{SparseGPUVector{Bool,Ti}, Nothing} = nothing,
) where {Tv,Ti,AV<:AbstractVector{Tv}}
# new version should look like : 
function gpu_spmv!(
    C::ResVec,
    A::SparseGPUMatrixCSR{Tv,Ti},
    B::InputVec;
    mul::Function = GPUGraphs_mul,
    add::Function = GPUGraphs_add,
    accum::Function = GPUGraphs_second,
    mask::Union{SparseGPUVector{Bool,Ti}, Nothing} = nothing,
) where {Tv,Ti<:Integer, ResType<:Number, InputType<:Number, ResVec<:AbstractVector{ResType}, InputVec<:AbstractVector{InputType}}

Ideas

• Rencently stumbled upon this performance tip. This is probably something I can change in the BFS code as I did not use the @. broadcasting syntax.

Todo

• Implement BFS_parents now that the index-based functions are available.
• Make a clear benchmark of BFS performances, identify bottleneck(s).
• Read some literature about masking, even though there appears to be no simple solution to our problem (according to this paper)
• Rewrite function header to allow type promotion

[gdalle]

06/05/2025

Work done

• Implemented masking with dense mask vector
• GPU BFS 3x faster without masking, 6x with masking, compared with CPU Graphs.jl
• Sofar only distances, but functions adapted to work with indices
• Type promotion works
• Profiled BFS: most of the time is spent in the sparse matrix-vector multiplication

Remarks

• No need to prefix every function with GPUGraphs
• Can we optimize the reduction inside the kernel?
• Stop once we have found a minimum, greedy and/or
• Hard to have a sparse vector as mask, not worth parallelizing anyway
• In the regime where parallelization is useful, the vector of frontier is dense
• Reduction transfers to the CPU (https://github.com/JuliaGPU/GPUArrays.jl/blob/e8e9b031613f31818e75a6c7f8745788fb80b71f/src/host/mapreduce.jl#L81)

Todo

• Try ELL version of storage for BFS
• Try out larger matrices, same ones as last year with parametric size / density
• Update benchmark scripts to fetch matrices automatically and run on NVIDIA GPUs
• Compare performance with SuiteSparseGraphBLAS
• Bellman-Ford algorithm as perspective
• Give up on sparse masks and sparse vectors