Use hash_accumulator for outer product, use sparse_intersection for

first phase of inner product.

Author: BenBrock

include/spblas/algorithms/detail/sparse_dot_product.hpp

@@ -33,6 +33,23 @@ std::optional<T> sparse_dot_product(__backend::spa_accumulator<T, I>& acc,
   }
 }
 
+template <typename Set, typename A, typename B>
+bool sparse_intersection(Set&& set, A&& a, B&& b) {
+  set.clear();
+
+  for (auto&& [i, v] : a) {
+    set.insert(i);
+  }
+
+  for (auto&& [i, v] : b) {
+    if (set.contains(i)) {
+      return true;
+    }
+  }
+
+  return false;
+}
+
 } // namespace __detail
 
 } // namespace spblas

include/spblas/algorithms/detail/spgemm/spgemm_innerproduct.hpp

@@ -83,16 +83,16 @@ operation_info_t multiply_compute(A&& a, B&& b, C&& c) {
 
   O nnz = 0;
 
-  __backend::spa_accumulator<T, I> dot_product_acc(__backend::shape(a)[1]);
+  __backend::spa_set<I> dot_product_acc(__backend::shape(a)[1]);
 
   for (auto&& [i, a_row] : __backend::rows(a)) {
     if (!__ranges::empty(a_row)) {
       for (auto&& [j, b_column] : __backend::columns(b)) {
         if (!__ranges::empty(b_column)) {
           auto v =
-              __detail::sparse_dot_product<T>(dot_product_acc, a_row, b_column);
+              __detail::sparse_intersection(dot_product_acc, a_row, b_column);
 
-          if (v.has_value()) {
+          if (v) {
             nnz++;
           }
         }

include/spblas/algorithms/detail/spgemm/spgemm_outerproduct.hpp

@@ -7,7 +7,7 @@
 #include <spblas/algorithms/detail/sparse_dot_product.hpp>
 #include <spblas/algorithms/transposed.hpp>
 #include <spblas/backend/csr_builder.hpp>
-#include <spblas/backend/spa_accumulator.hpp>
+#include <spblas/backend/hash_accumulator.hpp>
 #include <spblas/detail/operation_info_t.hpp>
 
 namespace spblas {
@@ -30,7 +30,7 @@ void multiply(A&& a, B&& b, C&& c) {
   using T = tensor_scalar_t<C>;
   using I = tensor_index_t<C>;
 
-  std::vector<__backend::spa_accumulator<T, I>> row_accumulators;
+  std::vector<__backend::hash_accumulator<T, I>> row_accumulators;
 
   for (std::size_t i = 0; i < __backend::shape(c)[0]; i++) {
     row_accumulators.emplace_back(__backend::shape(c)[1]);
@@ -80,7 +80,7 @@ operation_info_t multiply_compute(A&& a, B&& b, C&& c) {
   using I = tensor_index_t<C>;
   using O = tensor_offset_t<C>;
 
-  std::vector<__backend::spa_accumulator<T, I>> row_accumulators;
+  std::vector<__backend::hash_accumulator<T, I>> row_accumulators;
 
   for (std::size_t i = 0; i < __backend::shape(c)[0]; i++) {
     row_accumulators.emplace_back(__backend::shape(c)[1]);

include/spblas/backend/hash_accumulator.hpp

@@ -0,0 +1,92 @@
+#pragma once
+
+#include <functional>
+#include <span>
+#include <tuple>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+#include <spblas/detail/ranges.hpp>
+
+namespace spblas {
+
+namespace __backend {
+
+template <typename T, std::integral I>
+class hash_accumulator {
+public:
+  hash_accumulator(I count) {}
+
+  T& operator[](I pos) {
+    return hash_[pos];
+  }
+
+  bool contains(I pos) {
+    return hash_.contains(pos);
+  }
+
+  void clear() {
+    hash_.clear();
+  }
+
+  I size() const {
+    return hash_.size();
+  }
+
+  bool empty() {
+    return hash_.empty();
+  }
+
+  void sort() {}
+
+  auto get() {
+    std::vector<std::pair<I, T>> values(hash_.begin(), hash_.end());
+
+    std::sort(values.begin(), values.end(), [](auto&& a, auto&& b) {
+      return std::get<0>(a) < std::get<0>(b);
+    });
+
+    return values;
+  }
+
+private:
+  std::unordered_map<I, T> hash_;
+};
+
+template <std::integral T>
+class hash_set {
+public:
+  hash_set(T count) {}
+
+  void insert(T key) {
+    set_.insert(key);
+  }
+
+  bool contains(T key) {
+    return set_.contains(key);
+  }
+
+  void clear() {
+    set_.clear();
+  }
+
+  T size() const {
+    return set_.size();
+  }
+
+  bool empty() {
+    return set_.empty();
+  }
+
+  auto get() const {
+    return __ranges::views::all(set_);
+  }
+
+private:
+  std::unordered_set<T> set_;
+};
+
+} // namespace __backend
+
+} // namespace spblas

include/spblas/backend/spa_accumulator.hpp

@@ -2,7 +2,6 @@
 
 #include <functional>
 #include <span>
-#include <spblas/detail/ranges.hpp>
 #include <tuple>
 #include <vector>
 

include/spblas/concepts.hpp

@@ -11,13 +11,14 @@ namespace spblas {
 /*
   The following types fulfill the matrix concept:
   - Instantiations of csr_view<...>
+  - Instantiations of csc_view<...>
   - Instantiations of mdspan<...> with rank 2
   - Instantiations of scaled_view<T> where M is a matrix
 */
 
 template <typename M>
 concept matrix =
-    __detail::is_csr_view_v<M> ||
+    __detail::is_csr_view_v<M> || __detail::is_csc_view_v<M> ||
     __detail::is_matrix_instantiation_of_mdspan_v<M> || __detail::matrix<M>;
 
 /*