Implement addition between AbstractDeviceSparseMatrix and DenseMatrix (#21)

* Initial plan

* Implement addition between AbstractDeviceSparseMatrix and DenseMatrix

Co-authored-by: albertomercurio <61953577+albertomercurio@users.noreply.github.com>

* Apply code formatting

* Add benchmarks for sparse + dense matrix addition

Co-authored-by: albertomercurio <61953577+albertomercurio@users.noreply.github.com>

---------

Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>
Co-authored-by: albertomercurio <61953577+albertomercurio@users.noreply.github.com>

Author: Copilot

benchmarks/matrix_benchmarks.jl

@@ -153,3 +153,53 @@ function benchmark_three_arg_dot!(
 
     return nothing
 end
+
+"""
+    benchmark_sparse_dense_add!(SUITE, array_constructor, array_type_name; N=10000, T=Float64)
+
+Benchmark sparse + dense matrix addition for CSC, CSR, and COO formats.
+
+# Arguments
+- `SUITE`: The BenchmarkGroup to add benchmarks to
+- `array_constructor`: Function to construct arrays (e.g., `Array`, `JLArray`)
+- `array_type_name`: String name for the array type (for display)
+
+# Keyword Arguments
+- `N`: Size of the matrix (default: 10000)
+- `T`: Element type (default: Float64)
+"""
+function benchmark_sparse_dense_add!(
+    SUITE,
+    array_constructor,
+    array_type_name;
+    N = 10000,
+    T = Float64,
+)
+    # Create sparse matrix with 1% density
+    sm_csc_std = sprand(T, N, N, 0.01)
+
+    # Convert to different formats
+    sm_csc = DeviceSparseMatrixCSC(sm_csc_std)
+    sm_csr = DeviceSparseMatrixCSR(sm_csc_std)
+    sm_coo = DeviceSparseMatrixCOO(sm_csc_std)
+
+    # Adapt to device
+    dsm_csc = adapt(array_constructor, sm_csc)
+    dsm_csr = adapt(array_constructor, sm_csr)
+    dsm_coo = adapt(array_constructor, sm_coo)
+
+    # Create dense matrix
+    dense_mat = adapt(array_constructor, randn(T, N, N))
+
+    # Level 3: Format (CSC, CSR, COO - will be plotted together)
+    SUITE["Sparse + Dense Addition"][array_type_name]["CSC"] =
+        @benchmarkable $dsm_csc + $dense_mat
+
+    SUITE["Sparse + Dense Addition"][array_type_name]["CSR"] =
+        @benchmarkable $dsm_csr + $dense_mat
+
+    SUITE["Sparse + Dense Addition"][array_type_name]["COO"] =
+        @benchmarkable $dsm_coo + $dense_mat
+
+    return nothing
+end

benchmarks/runbenchmarks.jl

@@ -20,6 +20,7 @@ benchmark_vector_sparse_dense_dot!(SUITE, Array, "Array")
 benchmark_matrix_vector_mul!(SUITE, Array, "Array")
 benchmark_matrix_matrix_mul!(SUITE, Array, "Array")
 benchmark_three_arg_dot!(SUITE, Array, "Array")
+benchmark_sparse_dense_add!(SUITE, Array, "Array")
 
 # Run benchmarks for JLArrays
 println("Running benchmarks for JLArrays...")
@@ -28,6 +29,7 @@ benchmark_vector_sparse_dense_dot!(SUITE, jl, "JLArray")
 benchmark_matrix_vector_mul!(SUITE, jl, "JLArray")
 benchmark_matrix_matrix_mul!(SUITE, jl, "JLArray")
 benchmark_three_arg_dot!(SUITE, jl, "JLArray")
+benchmark_sparse_dense_add!(SUITE, jl, "JLArray")
 
 # Tune and run benchmarks
 println("\nTuning benchmarks...")

src/core.jl

@@ -55,3 +55,56 @@ trans_adj_wrappers(fmt) = (
     ),
     (T -> :(Adjoint{$T,<:$fmt{$T}}), true, true, A -> :(parent($A)), T -> :($T)),
 )
+
+# Generic addition between AbstractDeviceSparseMatrix and DenseMatrix
+"""
+    +(A::AbstractDeviceSparseMatrix, B::DenseMatrix)
+
+Add a sparse matrix `A` to a dense matrix `B`, returning a dense matrix.
+All backends must be compatible.
+
+# Examples
+```jldoctest
+julia> using DeviceSparseArrays, SparseArrays
+
+julia> A = DeviceSparseMatrixCSC(sparse([1, 2], [1, 2], [1.0, 2.0], 3, 3));
+
+julia> B = ones(3, 3);
+
+julia> C = A + B;
+
+julia> collect(C)
+3×3 Matrix{Float64}:
+ 2.0  1.0  1.0
+ 1.0  3.0  1.0
+ 1.0  1.0  1.0
+```
+"""
+function Base.:+(A::AbstractDeviceSparseMatrix, B::DenseMatrix)
+    size(A) == size(B) || throw(
+        DimensionMismatch(
+            "dimensions must match: A has dims $(size(A)), B has dims $(size(B))",
+        ),
+    )
+
+    backend_A = get_backend(A)
+    backend_B = get_backend(B)
+
+    backend_A == backend_B || throw(ArgumentError("Both arrays must have the same backend"))
+
+    # Create a copy of B to avoid modifying the input
+    C = copy(B)
+
+    # Add the sparse values to C
+    _add_sparse_to_dense!(C, A)
+
+    return C
+end
+
+"""
+    +(B::DenseMatrix, A::AbstractDeviceSparseMatrix)
+
+Add a dense matrix `B` to a sparse matrix `A`, returning a dense matrix.
+This is the commutative version of `A + B`.
+"""
+Base.:+(B::DenseMatrix, A::AbstractDeviceSparseMatrix) = A + B

src/matrix_coo/matrix_coo.jl

@@ -320,3 +320,14 @@ for (wrapa, transa, conja, unwrapa, whereT1) in trans_adj_wrappers(:DeviceSparse
         return sum(block_results)
     end
 end
+
+# Helper function for adding DeviceSparseMatrixCOO to dense matrix
+function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCOO)
+    backend = get_backend(A)
+    nnz_val = nnz(A)
+
+    kernel! = kernel_add_sparse_to_dense_coo!(backend)
+    kernel!(C, getrowind(A), getcolind(A), getnzval(A); ndrange = (nnz_val,))
+
+    return C
+end

src/matrix_coo/matrix_coo_kernels.jl

@@ -127,3 +127,15 @@ end
         block_results[group_id] = sum
     end
 end
+
+# Kernel for adding sparse matrix to dense matrix (COO format)
+@kernel inbounds=true function kernel_add_sparse_to_dense_coo!(
+    C,
+    @Const(rowind),
+    @Const(colind),
+    @Const(nzval),
+)
+    i = @index(Global)
+
+    @inbounds C[rowind[i], colind[i]] += nzval[i]
+end

src/matrix_csc/matrix_csc.jl

@@ -307,3 +307,14 @@ for (wrapa, transa, conja, unwrapa, whereT1) in trans_adj_wrappers(:DeviceSparse
         return sum(block_results)
     end
 end
+
+# Helper function for adding DeviceSparseMatrixCSC to dense matrix
+function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCSC)
+    backend = get_backend(A)
+    n = size(A, 2)
+
+    kernel! = kernel_add_sparse_to_dense_csc!(backend)
+    kernel!(C, getcolptr(A), getrowval(A), getnzval(A); ndrange = (n,))
+
+    return C
+end

src/matrix_csc/matrix_csc_kernels.jl

@@ -129,3 +129,17 @@ end
         block_results[group_id] = sum
     end
 end
+
+# Kernel for adding sparse matrix to dense matrix (CSC format)
+@kernel inbounds=true function kernel_add_sparse_to_dense_csc!(
+    C,
+    @Const(colptr),
+    @Const(rowval),
+    @Const(nzval),
+)
+    col = @index(Global)
+
+    @inbounds for j = colptr[col]:(colptr[col+1]-1)
+        C[rowval[j], col] += nzval[j]
+    end
+end

src/matrix_csr/matrix_csr.jl

@@ -320,3 +320,14 @@ for (wrapa, transa, conja, unwrapa, whereT1) in trans_adj_wrappers(:DeviceSparse
         return sum(block_results)
     end
 end
+
+# Helper function for adding DeviceSparseMatrixCSR to dense matrix
+function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCSR)
+    backend = get_backend(A)
+    m = size(A, 1)
+
+    kernel! = kernel_add_sparse_to_dense_csr!(backend)
+    kernel!(C, getrowptr(A), getcolval(A), getnzval(A); ndrange = (m,))
+
+    return C
+end

src/matrix_csr/matrix_csr_kernels.jl

@@ -129,3 +129,17 @@ end
         block_results[group_id] = sum
     end
 end
+
+# Kernel for adding sparse matrix to dense matrix (CSR format)
+@kernel inbounds=true function kernel_add_sparse_to_dense_csr!(
+    C,
+    @Const(rowptr),
+    @Const(colval),
+    @Const(nzval),
+)
+    row = @index(Global)
+
+    @inbounds for j = rowptr[row]:(rowptr[row+1]-1)
+        C[row, colval[j]] += nzval[j]
+    end
+end

test/shared/matrix_coo.jl

@@ -252,4 +252,29 @@ function shared_test_linearalgebra_matrix_coo(
             end
         end
     end
+
+    @testset "Sparse + Dense Matrix Addition" begin
+        for T in (int_types..., float_types..., complex_types...)
+            m, n = 50, 40
+            A = sprand(T, m, n, 0.1)
+            B = rand(T, m, n)
+
+            dA = adapt(op, DeviceSparseMatrixCOO(A))
+            dB = op(B)
+
+            # Test sparse + dense
+            result = dA + dB
+            expected = Matrix(A) + B
+            @test collect(result) ≈ expected
+
+            # Test dense + sparse (commutative)
+            result2 = dB + dA
+            @test collect(result2) ≈ expected
+
+            # Test dimension mismatch
+            B_wrong = rand(T, m + 1, n)
+            dB_wrong = op(B_wrong)
+            @test_throws DimensionMismatch dA + dB_wrong
+        end
+    end
 end