diff --git a/benchmarks/matrix_benchmarks.jl b/benchmarks/matrix_benchmarks.jl
index ef9f64d..3d241e6 100644
--- a/benchmarks/matrix_benchmarks.jl
+++ b/benchmarks/matrix_benchmarks.jl
@@ -228,6 +228,66 @@ function benchmark_sparse_dense_add!(
     return nothing
 end
 
+"""
+    benchmark_sparse_sparse_add!(SUITE, array_constructor, array_type_name; N=10000, T=Float64)
+
+Benchmark sparse + sparse 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_sparse_add!(
+    SUITE,
+    array_constructor,
+    array_type_name;
+    N = 10000,
+    T = Float64,
+)
+    # Create two sparse matrices with 1% density
+    sm_a_csc_std = sprand(T, N, N, 0.01)
+    sm_b_csc_std = sprand(T, N, N, 0.01)
+
+    # Convert to different formats
+    sm_a_csc = DeviceSparseMatrixCSC(sm_a_csc_std)
+    sm_b_csc = DeviceSparseMatrixCSC(sm_b_csc_std)
+    sm_a_csr = DeviceSparseMatrixCSR(sm_a_csc_std)
+    sm_b_csr = DeviceSparseMatrixCSR(sm_b_csc_std)
+    sm_a_coo = DeviceSparseMatrixCOO(sm_a_csc_std)
+    sm_b_coo = DeviceSparseMatrixCOO(sm_b_csc_std)
+
+    # Adapt to device
+    dsm_a_csc = adapt(array_constructor, sm_a_csc)
+    dsm_b_csc = adapt(array_constructor, sm_b_csc)
+    dsm_a_csr = adapt(array_constructor, sm_a_csr)
+    dsm_b_csr = adapt(array_constructor, sm_b_csr)
+    dsm_a_coo = adapt(array_constructor, sm_a_coo)
+    dsm_b_coo = adapt(array_constructor, sm_b_coo)
+
+    # Level 3: Format (CSC, CSR, COO - will be plotted together)
+    SUITE["Sparse + Sparse Addition"][array_type_name]["CSC"] = @benchmarkable begin
+        $dsm_a_csc + $dsm_b_csc
+        _synchronize_backend($dsm_a_csc)
+    end
+
+    SUITE["Sparse + Sparse Addition"][array_type_name]["CSR"] = @benchmarkable begin
+        $dsm_a_csr + $dsm_b_csr
+        _synchronize_backend($dsm_a_csr)
+    end
+
+    SUITE["Sparse + Sparse Addition"][array_type_name]["COO"] = @benchmarkable begin
+        $dsm_a_coo + $dsm_b_coo
+        _synchronize_backend($dsm_a_coo)
+    end
+
+    return nothing
+end
+
 """
     benchmark_kron!(SUITE, array_constructor, array_type_name; N=100, T=Float64)
 
diff --git a/benchmarks/runbenchmarks.jl b/benchmarks/runbenchmarks.jl
index ca78d25..763e879 100644
--- a/benchmarks/runbenchmarks.jl
+++ b/benchmarks/runbenchmarks.jl
@@ -26,6 +26,7 @@ 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")
+benchmark_sparse_sparse_add!(SUITE, Array, "Array")
 benchmark_kron!(SUITE, Array, "Array")
 benchmark_conversions!(SUITE, Array, "Array")
 
@@ -37,6 +38,7 @@ benchmark_matrix_vector_mul!(SUITE, JLArray, "JLArray")
 benchmark_matrix_matrix_mul!(SUITE, JLArray, "JLArray")
 benchmark_three_arg_dot!(SUITE, JLArray, "JLArray")
 benchmark_sparse_dense_add!(SUITE, JLArray, "JLArray")
+benchmark_sparse_sparse_add!(SUITE, JLArray, "JLArray")
 benchmark_kron!(SUITE, JLArray, "JLArray")
 benchmark_conversions!(SUITE, JLArray, "JLArray")
 
diff --git a/src/matrix_coo/matrix_coo.jl b/src/matrix_coo/matrix_coo.jl
index 82f22a8..e526040 100644
--- a/src/matrix_coo/matrix_coo.jl
+++ b/src/matrix_coo/matrix_coo.jl
@@ -328,6 +328,110 @@ function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCOO)
     return C
 end
 
+"""
+    +(A::DeviceSparseMatrixCOO, B::DeviceSparseMatrixCOO)
+
+Add two sparse matrices in COO format. Both matrices must have the same dimensions
+and be on the same backend (device).
+
+The result is a COO matrix with entries from both A and B properly merged,
+with duplicate entries (same row and column) combined by summing their values.
+
+# Examples
+```jldoctest
+julia> using DeviceSparseArrays, SparseArrays
+
+julia> A = DeviceSparseMatrixCOO(sparse([1, 2], [1, 2], [1.0, 2.0], 2, 2));
+
+julia> B = DeviceSparseMatrixCOO(sparse([1, 2], [2, 1], [3.0, 4.0], 2, 2));
+
+julia> C = A + B;
+
+julia> collect(C)
+2×2 Matrix{Float64}:
+ 1.0  3.0
+ 4.0  2.0
+```
+"""
+function Base.:+(A::DeviceSparseMatrixCOO, B::DeviceSparseMatrixCOO)
+    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 matrices must have the same backend"))
+
+    m, n = size(A)
+    Ti = eltype(getrowind(A))
+    Tv = promote_type(eltype(nonzeros(A)), eltype(nonzeros(B)))
+
+    # Concatenate the coordinate arrays
+    nnz_A = nnz(A)
+    nnz_B = nnz(B)
+    nnz_concat = nnz_A + nnz_B
+
+    # Allocate concatenated arrays
+    rowind_concat = similar(getrowind(A), nnz_concat)
+    colind_concat = similar(getcolind(A), nnz_concat)
+    nzval_concat = similar(nonzeros(A), Tv, nnz_concat)
+
+    # Copy entries from A and B
+    rowind_concat[1:nnz_A] .= getrowind(A)
+    colind_concat[1:nnz_A] .= getcolind(A)
+    nzval_concat[1:nnz_A] .= nonzeros(A)
+    rowind_concat[(nnz_A+1):end] .= getrowind(B)
+    colind_concat[(nnz_A+1):end] .= getcolind(B)
+    nzval_concat[(nnz_A+1):end] .= nonzeros(B)
+
+    # Sort by (row, col) using keys similar to COO->CSC conversion
+    backend = backend_A
+    keys = similar(rowind_concat, Ti, nnz_concat)
+    kernel_make_keys! = kernel_make_csc_keys!(backend)
+    kernel_make_keys!(keys, rowind_concat, colind_concat, m; ndrange = (nnz_concat,))
+
+    # Sort using AcceleratedKernels
+    perm = _sortperm_AK(keys)
+
+    # Apply permutation to get sorted arrays
+    rowind_sorted = rowind_concat[perm]
+    colind_sorted = colind_concat[perm]
+    nzval_sorted = nzval_concat[perm]
+
+    # Mark unique entries (first occurrence of each (row, col) pair)
+    keep_mask = similar(rowind_sorted, Bool, nnz_concat)
+    kernel_mark! = kernel_mark_unique_coo!(backend)
+    kernel_mark!(keep_mask, rowind_sorted, colind_sorted, nnz_concat; ndrange = (nnz_concat,))
+
+    # Compute write indices using cumsum
+    write_indices = _cumsum_AK(keep_mask)
+    nnz_final = allowed_getindex(write_indices, nnz_concat)
+
+    # Allocate final arrays
+    rowind_C = similar(getrowind(A), nnz_final)
+    colind_C = similar(getcolind(A), nnz_final)
+    nzval_C = similar(nonzeros(A), Tv, nnz_final)
+
+    # Compact: merge duplicates by summing
+    kernel_compact! = kernel_compact_coo!(backend)
+    kernel_compact!(
+        rowind_C,
+        colind_C,
+        nzval_C,
+        rowind_sorted,
+        colind_sorted,
+        nzval_sorted,
+        write_indices,
+        nnz_concat;
+        ndrange = (nnz_concat,),
+    )
+
+    return DeviceSparseMatrixCOO(m, n, rowind_C, colind_C, nzval_C)
+end
+
 """
     kron(A::DeviceSparseMatrixCOO, B::DeviceSparseMatrixCOO)
 
diff --git a/src/matrix_coo/matrix_coo_kernels.jl b/src/matrix_coo/matrix_coo_kernels.jl
index 2184223..e412fe2 100644
--- a/src/matrix_coo/matrix_coo_kernels.jl
+++ b/src/matrix_coo/matrix_coo_kernels.jl
@@ -181,3 +181,55 @@ end
         nzval_C[idx] = val_A * val_B
     end
 end
+
+# Kernel for marking duplicate entries in sorted COO format
+# Returns a mask where mask[i] = true if entry i should be kept (first occurrence or sum)
+@kernel inbounds=true function kernel_mark_unique_coo!(
+    keep_mask,
+    @Const(rowind),
+    @Const(colind),
+    @Const(nnz_total),
+)
+    i = @index(Global)
+
+    if i == 1
+        # Always keep the first entry
+        keep_mask[i] = true
+    elseif i <= nnz_total
+        # Keep if different from previous entry
+        keep_mask[i] = (rowind[i] != rowind[i-1] || colind[i] != colind[i-1])
+    end
+end
+
+# Kernel for compacting COO by summing duplicate entries
+@kernel inbounds=true function kernel_compact_coo!(
+    rowind_out,
+    colind_out,
+    nzval_out,
+    @Const(rowind_in),
+    @Const(colind_in),
+    @Const(nzval_in),
+    @Const(write_indices),
+    @Const(nnz_in),
+)
+    i = @index(Global)
+
+    if i <= nnz_in
+        out_idx = write_indices[i]
+        
+        # If this is a new entry (or first of duplicates), write it
+        if i == 1 || (rowind_in[i] != rowind_in[i-1] || colind_in[i] != colind_in[i-1])
+            rowind_out[out_idx] = rowind_in[i]
+            colind_out[out_idx] = colind_in[i]
+            
+            # Sum all duplicates
+            val_sum = nzval_in[i]
+            j = i + 1
+            while j <= nnz_in && rowind_in[j] == rowind_in[i] && colind_in[j] == colind_in[i]
+                val_sum += nzval_in[j]
+                j += 1
+            end
+            nzval_out[out_idx] = val_sum
+        end
+    end
+end
diff --git a/src/matrix_csc/matrix_csc.jl b/src/matrix_csc/matrix_csc.jl
index 68091a0..8a67f3d 100644
--- a/src/matrix_csc/matrix_csc.jl
+++ b/src/matrix_csc/matrix_csc.jl
@@ -315,6 +315,92 @@ function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCSC)
     return C
 end
 
+"""
+    +(A::DeviceSparseMatrixCSC, B::DeviceSparseMatrixCSC)
+
+Add two sparse matrices in CSC format. Both matrices must have the same dimensions
+and be on the same backend (device).
+
+# Examples
+```jldoctest
+julia> using DeviceSparseArrays, SparseArrays
+
+julia> A = DeviceSparseMatrixCSC(sparse([1, 2], [1, 2], [1.0, 2.0], 2, 2));
+
+julia> B = DeviceSparseMatrixCSC(sparse([1, 2], [2, 1], [3.0, 4.0], 2, 2));
+
+julia> C = A + B;
+
+julia> collect(C)
+2×2 Matrix{Float64}:
+ 1.0  3.0
+ 4.0  2.0
+```
+"""
+function Base.:+(A::DeviceSparseMatrixCSC, B::DeviceSparseMatrixCSC)
+    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 matrices must have the same backend"))
+
+    m, n = size(A)
+    Ti = eltype(getcolptr(A))
+    Tv = promote_type(eltype(nonzeros(A)), eltype(nonzeros(B)))
+
+    # Count non-zeros per column
+    nnz_per_col = similar(getcolptr(A), n)
+    fill!(nnz_per_col, zero(Ti))
+
+    backend = backend_A
+    kernel_count! = kernel_count_nnz_per_col_csc!(backend)
+    kernel_count!(
+        nnz_per_col,
+        getcolptr(A),
+        getrowval(A),
+        getcolptr(B),
+        getrowval(B);
+        ndrange = (n,),
+    )
+
+    # Build colptr for result matrix using cumsum
+    # colptr_C[i+1] = 1 + sum(nnz_per_col[1:i])
+    cumsum_nnz = _cumsum_AK(nnz_per_col)
+    colptr_C = similar(getcolptr(A), n + 1)
+    # Set colptr_C[2:end] to cumsum + 1
+    colptr_C[2:end] .= cumsum_nnz
+    colptr_C[2:end] .+= one(Ti)
+    # Set colptr_C[1] to 1 using broadcasting
+    colptr_C[1:1] .= one(Ti)
+
+    # Allocate result arrays
+    nnz_total = allowed_getindex(colptr_C, n + 1) - one(Ti)
+    rowval_C = similar(getrowval(A), nnz_total)
+    nzval_C = similar(nonzeros(A), Tv, nnz_total)
+
+    # Merge the two matrices
+    kernel_merge! = kernel_merge_csc!(backend)
+    kernel_merge!(
+        rowval_C,
+        nzval_C,
+        colptr_C,
+        getcolptr(A),
+        getrowval(A),
+        nonzeros(A),
+        getcolptr(B),
+        getrowval(B),
+        nonzeros(B);
+        ndrange = (n,),
+    )
+
+    return DeviceSparseMatrixCSC(m, n, colptr_C, rowval_C, nzval_C)
+end
+
 """
     kron(A::DeviceSparseMatrixCSC, B::DeviceSparseMatrixCSC)
 
diff --git a/src/matrix_csc/matrix_csc_kernels.jl b/src/matrix_csc/matrix_csc_kernels.jl
index e1e29bf..12fdb93 100644
--- a/src/matrix_csc/matrix_csc_kernels.jl
+++ b/src/matrix_csc/matrix_csc_kernels.jl
@@ -143,3 +143,102 @@ end
         C[rowval[j], col] += nzval[j]
     end
 end
+
+# Kernel for counting non-zeros per column when adding two CSC matrices
+@kernel inbounds=true function kernel_count_nnz_per_col_csc!(
+    nnz_per_col,
+    @Const(colptr_A),
+    @Const(rowval_A),
+    @Const(colptr_B),
+    @Const(rowval_B),
+)
+    col = @index(Global)
+
+    i_A = colptr_A[col]
+    i_B = colptr_B[col]
+    end_A = colptr_A[col+1]
+    end_B = colptr_B[col+1]
+
+    count = 0
+    while i_A < end_A && i_B < end_B
+        row_A = rowval_A[i_A]
+        row_B = rowval_B[i_B]
+        if row_A < row_B
+            count += 1
+            i_A += 1
+        elseif row_A > row_B
+            count += 1
+            i_B += 1
+        else  # row_A == row_B
+            count += 1
+            i_A += 1
+            i_B += 1
+        end
+    end
+
+    # Add remaining entries from A
+    count += (end_A - i_A)
+    # Add remaining entries from B
+    count += (end_B - i_B)
+
+    nnz_per_col[col] = count
+end
+
+# Kernel for merging two CSC matrices (addition)
+@kernel inbounds=true function kernel_merge_csc!(
+    rowval_C,
+    nzval_C,
+    @Const(colptr_C),
+    @Const(colptr_A),
+    @Const(rowval_A),
+    @Const(nzval_A),
+    @Const(colptr_B),
+    @Const(rowval_B),
+    @Const(nzval_B),
+)
+    col = @index(Global)
+
+    i_A = colptr_A[col]
+    i_B = colptr_B[col]
+    i_C = colptr_C[col]
+    end_A = colptr_A[col+1]
+    end_B = colptr_B[col+1]
+
+    while i_A < end_A && i_B < end_B
+        row_A = rowval_A[i_A]
+        row_B = rowval_B[i_B]
+        if row_A < row_B
+            rowval_C[i_C] = row_A
+            nzval_C[i_C] = nzval_A[i_A]
+            i_A += 1
+            i_C += 1
+        elseif row_A > row_B
+            rowval_C[i_C] = row_B
+            nzval_C[i_C] = nzval_B[i_B]
+            i_B += 1
+            i_C += 1
+        else  # row_A == row_B
+            rowval_C[i_C] = row_A
+            nzval_C[i_C] = nzval_A[i_A] + nzval_B[i_B]
+            i_A += 1
+            i_B += 1
+            i_C += 1
+        end
+    end
+
+    # Copy remaining entries from A
+    while i_A < end_A
+        rowval_C[i_C] = rowval_A[i_A]
+        nzval_C[i_C] = nzval_A[i_A]
+        i_A += 1
+        i_C += 1
+    end
+
+    # Copy remaining entries from B
+    while i_B < end_B
+        rowval_C[i_C] = rowval_B[i_B]
+        nzval_C[i_C] = nzval_B[i_B]
+        i_B += 1
+        i_C += 1
+    end
+end
diff --git a/src/matrix_csr/matrix_csr.jl b/src/matrix_csr/matrix_csr.jl
index fb97aec..9dd8035 100644
--- a/src/matrix_csr/matrix_csr.jl
+++ b/src/matrix_csr/matrix_csr.jl
@@ -312,6 +312,92 @@ function _add_sparse_to_dense!(C::DenseMatrix, A::DeviceSparseMatrixCSR)
     return C
 end
 
+"""
+    +(A::DeviceSparseMatrixCSR, B::DeviceSparseMatrixCSR)
+
+Add two sparse matrices in CSR format. Both matrices must have the same dimensions
+and be on the same backend (device).
+
+# Examples
+```jldoctest
+julia> using DeviceSparseArrays, SparseArrays
+
+julia> A = DeviceSparseMatrixCSR(sparse([1, 2], [1, 2], [1.0, 2.0], 2, 2));
+
+julia> B = DeviceSparseMatrixCSR(sparse([1, 2], [2, 1], [3.0, 4.0], 2, 2));
+
+julia> C = A + B;
+
+julia> collect(C)
+2×2 Matrix{Float64}:
+ 1.0  3.0
+ 4.0  2.0
+```
+"""
+function Base.:+(A::DeviceSparseMatrixCSR, B::DeviceSparseMatrixCSR)
+    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 matrices must have the same backend"))
+
+    m, n = size(A)
+    Ti = eltype(getrowptr(A))
+    Tv = promote_type(eltype(nonzeros(A)), eltype(nonzeros(B)))
+
+    # Count non-zeros per row
+    nnz_per_row = similar(getrowptr(A), m)
+    fill!(nnz_per_row, zero(Ti))
+
+    backend = backend_A
+    kernel_count! = kernel_count_nnz_per_row_csr!(backend)
+    kernel_count!(
+        nnz_per_row,
+        getrowptr(A),
+        getcolval(A),
+        getrowptr(B),
+        getcolval(B);
+        ndrange = (m,),
+    )
+
+    # Build rowptr for result matrix using cumsum
+    # rowptr_C[i+1] = 1 + sum(nnz_per_row[1:i])
+    cumsum_nnz = _cumsum_AK(nnz_per_row)
+    rowptr_C = similar(getrowptr(A), m + 1)
+    # Set rowptr_C[2:end] to cumsum + 1
+    rowptr_C[2:end] .= cumsum_nnz
+    rowptr_C[2:end] .+= one(Ti)
+    # Set rowptr_C[1] to 1 using broadcasting
+    rowptr_C[1:1] .= one(Ti)
+
+    # Allocate result arrays
+    nnz_total = allowed_getindex(rowptr_C, m + 1) - one(Ti)
+    colval_C = similar(getcolval(A), nnz_total)
+    nzval_C = similar(nonzeros(A), Tv, nnz_total)
+
+    # Merge the two matrices
+    kernel_merge! = kernel_merge_csr!(backend)
+    kernel_merge!(
+        colval_C,
+        nzval_C,
+        rowptr_C,
+        getrowptr(A),
+        getcolval(A),
+        nonzeros(A),
+        getrowptr(B),
+        getcolval(B),
+        nonzeros(B);
+        ndrange = (m,),
+    )
+
+    return DeviceSparseMatrixCSR(m, n, rowptr_C, colval_C, nzval_C)
+end
+
 """
     kron(A::DeviceSparseMatrixCSR, B::DeviceSparseMatrixCSR)
 
diff --git a/src/matrix_csr/matrix_csr_kernels.jl b/src/matrix_csr/matrix_csr_kernels.jl
index ac6ae5a..1aa7cea 100644
--- a/src/matrix_csr/matrix_csr_kernels.jl
+++ b/src/matrix_csr/matrix_csr_kernels.jl
@@ -143,3 +143,102 @@ end
         C[row, colval[j]] += nzval[j]
     end
 end
+
+# Kernel for counting non-zeros per row when adding two CSR matrices
+@kernel inbounds=true function kernel_count_nnz_per_row_csr!(
+    nnz_per_row,
+    @Const(rowptr_A),
+    @Const(colval_A),
+    @Const(rowptr_B),
+    @Const(colval_B),
+)
+    row = @index(Global)
+
+    i_A = rowptr_A[row]
+    i_B = rowptr_B[row]
+    end_A = rowptr_A[row+1]
+    end_B = rowptr_B[row+1]
+
+    count = 0
+    while i_A < end_A && i_B < end_B
+        col_A = colval_A[i_A]
+        col_B = colval_B[i_B]
+        if col_A < col_B
+            count += 1
+            i_A += 1
+        elseif col_A > col_B
+            count += 1
+            i_B += 1
+        else  # col_A == col_B
+            count += 1
+            i_A += 1
+            i_B += 1
+        end
+    end
+
+    # Add remaining entries from A
+    count += (end_A - i_A)
+    # Add remaining entries from B
+    count += (end_B - i_B)
+
+    nnz_per_row[row] = count
+end
+
+# Kernel for merging two CSR matrices (addition)
+@kernel inbounds=true function kernel_merge_csr!(
+    colval_C,
+    nzval_C,
+    @Const(rowptr_C),
+    @Const(rowptr_A),
+    @Const(colval_A),
+    @Const(nzval_A),
+    @Const(rowptr_B),
+    @Const(colval_B),
+    @Const(nzval_B),
+)
+    row = @index(Global)
+
+    i_A = rowptr_A[row]
+    i_B = rowptr_B[row]
+    i_C = rowptr_C[row]
+    end_A = rowptr_A[row+1]
+    end_B = rowptr_B[row+1]
+
+    while i_A < end_A && i_B < end_B
+        col_A = colval_A[i_A]
+        col_B = colval_B[i_B]
+        if col_A < col_B
+            colval_C[i_C] = col_A
+            nzval_C[i_C] = nzval_A[i_A]
+            i_A += 1
+            i_C += 1
+        elseif col_A > col_B
+            colval_C[i_C] = col_B
+            nzval_C[i_C] = nzval_B[i_B]
+            i_B += 1
+            i_C += 1
+        else  # col_A == col_B
+            colval_C[i_C] = col_A
+            nzval_C[i_C] = nzval_A[i_A] + nzval_B[i_B]
+            i_A += 1
+            i_B += 1
+            i_C += 1
+        end
+    end
+
+    # Copy remaining entries from A
+    while i_A < end_A
+        colval_C[i_C] = colval_A[i_A]
+        nzval_C[i_C] = nzval_A[i_A]
+        i_A += 1
+        i_C += 1
+    end
+
+    # Copy remaining entries from B
+    while i_B < end_B
+        colval_C[i_C] = colval_B[i_B]
+        nzval_C[i_C] = nzval_B[i_B]
+        i_B += 1
+        i_C += 1
+    end
+end
diff --git a/test/shared/matrix_coo.jl b/test/shared/matrix_coo.jl
index 861a2fc..0e5e4d3 100644
--- a/test/shared/matrix_coo.jl
+++ b/test/shared/matrix_coo.jl
@@ -278,6 +278,37 @@ function shared_test_linearalgebra_matrix_coo(
         end
     end
 
+    @testset "Sparse + Sparse Matrix Addition" begin
+        for T in (int_types..., float_types..., complex_types...)
+            m, n = 50, 40
+            A = sprand(T, m, n, 0.1)
+            B = sprand(T, m, n, 0.15)
+
+            dA = adapt(op, DeviceSparseMatrixCOO(A))
+            dB = adapt(op, DeviceSparseMatrixCOO(B))
+
+            # Test sparse + sparse
+            result = dA + dB
+            expected = A + B
+            @test collect(result) ≈ Matrix(expected)
+            @test result isa DeviceSparseMatrixCOO
+
+            # Test with overlapping entries
+            A_overlap = sparse([1, 2, 3], [1, 2, 3], T[1, 2, 3], m, n)
+            B_overlap = sparse([1, 2, 4], [1, 2, 4], T[4, 5, 6], m, n)
+            dA_overlap = adapt(op, DeviceSparseMatrixCOO(A_overlap))
+            dB_overlap = adapt(op, DeviceSparseMatrixCOO(B_overlap))
+            result_overlap = dA_overlap + dB_overlap
+            expected_overlap = A_overlap + B_overlap
+            @test collect(result_overlap) ≈ Matrix(expected_overlap)
+
+            # Test dimension mismatch
+            B_wrong = sprand(T, m + 1, n, 0.1)
+            dB_wrong = adapt(op, DeviceSparseMatrixCOO(B_wrong))
+            @test_throws DimensionMismatch dA + dB_wrong
+        end
+    end
+
     @testset "Kronecker Product" begin
         for T in (int_types..., float_types..., complex_types...)
             # Test with rectangular matrices
diff --git a/test/shared/matrix_csc.jl b/test/shared/matrix_csc.jl
index ee4768b..cc5a297 100644
--- a/test/shared/matrix_csc.jl
+++ b/test/shared/matrix_csc.jl
@@ -276,6 +276,37 @@ function shared_test_linearalgebra_matrix_csc(
         end
     end
 
+    @testset "Sparse + Sparse Matrix Addition" begin
+        for T in (int_types..., float_types..., complex_types...)
+            m, n = 50, 40
+            A = sprand(T, m, n, 0.1)
+            B = sprand(T, m, n, 0.15)
+
+            dA = adapt(op, DeviceSparseMatrixCSC(A))
+            dB = adapt(op, DeviceSparseMatrixCSC(B))
+
+            # Test sparse + sparse
+            result = dA + dB
+            expected = A + B
+            @test collect(result) ≈ Matrix(expected)
+            @test result isa DeviceSparseMatrixCSC
+
+            # Test with overlapping entries
+            A_overlap = sparse([1, 2, 3], [1, 2, 3], T[1, 2, 3], m, n)
+            B_overlap = sparse([1, 2, 4], [1, 2, 4], T[4, 5, 6], m, n)
+            dA_overlap = adapt(op, DeviceSparseMatrixCSC(A_overlap))
+            dB_overlap = adapt(op, DeviceSparseMatrixCSC(B_overlap))
+            result_overlap = dA_overlap + dB_overlap
+            expected_overlap = A_overlap + B_overlap
+            @test collect(result_overlap) ≈ Matrix(expected_overlap)
+
+            # Test dimension mismatch
+            B_wrong = sprand(T, m + 1, n, 0.1)
+            dB_wrong = adapt(op, DeviceSparseMatrixCSC(B_wrong))
+            @test_throws DimensionMismatch dA + dB_wrong
+        end
+    end
+
     @testset "Kronecker Product" begin
         if array_type != "JLArray"
             for T in (int_types..., float_types..., complex_types...)
diff --git a/test/shared/matrix_csr.jl b/test/shared/matrix_csr.jl
index 2be572d..4ea1a41 100644
--- a/test/shared/matrix_csr.jl
+++ b/test/shared/matrix_csr.jl
@@ -275,6 +275,37 @@ function shared_test_linearalgebra_matrix_csr(
         end
     end
 
+    @testset "Sparse + Sparse Matrix Addition" begin
+        for T in (int_types..., float_types..., complex_types...)
+            m, n = 50, 40
+            A = sprand(T, m, n, 0.1)
+            B = sprand(T, m, n, 0.15)
+
+            dA = adapt(op, DeviceSparseMatrixCSR(A))
+            dB = adapt(op, DeviceSparseMatrixCSR(B))
+
+            # Test sparse + sparse
+            result = dA + dB
+            expected = A + B
+            @test collect(result) ≈ Matrix(expected)
+            @test result isa DeviceSparseMatrixCSR
+
+            # Test with overlapping entries
+            A_overlap = sparse([1, 2, 3], [1, 2, 3], T[1, 2, 3], m, n)
+            B_overlap = sparse([1, 2, 4], [1, 2, 4], T[4, 5, 6], m, n)
+            dA_overlap = adapt(op, DeviceSparseMatrixCSR(A_overlap))
+            dB_overlap = adapt(op, DeviceSparseMatrixCSR(B_overlap))
+            result_overlap = dA_overlap + dB_overlap
+            expected_overlap = A_overlap + B_overlap
+            @test collect(result_overlap) ≈ Matrix(expected_overlap)
+
+            # Test dimension mismatch
+            B_wrong = sprand(T, m + 1, n, 0.1)
+            dB_wrong = adapt(op, DeviceSparseMatrixCSR(B_wrong))
+            @test_throws DimensionMismatch dA + dB_wrong
+        end
+    end
+
     @testset "Kronecker Product" begin
         if array_type != "JLArray"
             for T in (int_types..., float_types..., complex_types...)