Add some sparse ops and rename others to have more explicit names

doc/extending/other_ops.rst

@@ -96,7 +96,7 @@ You can extract the 4 fields with
 :func:`PyTensor.sparse.basic.csm_shape` to extract the individual
 fields.
 
-You can look at the `AddSD` sparse `Op` for an example with C code. It implements
+You can look at the `AddSparseDense` sparse `Op` for an example with C code. It implements
 the addition of a sparse matrix with a dense matrix.
 
 Sparse Tests

pytensor/link/numba/dispatch/sparse/basic.py

@@ -18,6 +18,7 @@
     CSMProperties,
     DenseFromSparse,
     Diag,
+    EnsureSortedIndices,
     GetItem2d,
     GetItem2Lists,
     GetItem2ListsGrad,
@@ -27,6 +28,7 @@
     HStack,
     RowScaleCSC,
     SparseFromDense,
+    SumDuplicates,
     Transpose,
     VStack,
 )
@@ -69,6 +71,57 @@ def csm_constructor(data, indices, indptr, shape):
     return csm_constructor
 
 
+@register_funcify_default_op_cache_key(EnsureSortedIndices)
+def numba_funcify_EnsureSortedIndices(op, node, **kwargs):
+    inplace = op.inplace
+
+    @numba_basic.numba_njit
+    def ensure_sorted_indices(x):
+        if inplace:
+            x.sort_indices()
+            return x
+        return x.sorted_indices()
+
+    return ensure_sorted_indices
+
+
+@register_funcify_default_op_cache_key(SumDuplicates)
+def numba_funcify_SumDuplicates(op, node, **kwargs):
+    format = node.outputs[0].type.format
+
+    @numba_basic.numba_njit
+    def sum_duplicates(x):
+        result = x.sorted_indices()
+        data = result.data
+        indices = result.indices.view(np.uint32)
+        indptr = result.indptr.view(np.uint32)
+        new_indptr = np.empty_like(indptr)
+        new_indptr[0] = 0
+        nnz = 0
+        for major in range(len(indptr) - 1):
+            for ptr in range(indptr[major], indptr[major + 1]):
+                index = indices[ptr]
+                if nnz > new_indptr[major] and indices[nnz - 1] == index:
+                    data[nnz - 1] += data[ptr]
+                else:
+                    data[nnz] = data[ptr]
+                    indices[nnz] = index
+                    nnz += 1
+            new_indptr[major + 1] = nnz
+
+        components = (
+            data[:nnz],
+            indices[:nnz].view(np.int32),
+            new_indptr.view(np.int32),
+        )
+
+        if format == "csr":
+            return sp.sparse.csr_matrix(components, shape=x.shape)
+        return sp.sparse.csc_matrix(components, shape=x.shape)
+
+    return sum_duplicates
+
+
 @register_funcify_default_op_cache_key(Cast)
 def numba_funcify_Cast(op, node, **kwargs):
     inp_dtype = node.inputs[0].type.dtype