Commit missing files

Author: BenBrock

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

@@ -126,4 +126,92 @@ operation_info_t multiply_compute(A&& a, B&& b, C&& c) {
   return info;
 }
 
+// C = AB
+// CSR * CSR -> CSC
+// SpGEMM (Gustavson's Algorithm, scattered)
+template <matrix A, matrix B, matrix C>
+  requires(__backend::row_iterable<A> && __backend::row_iterable<B> &&
+           __detail::is_csc_view_v<C>)
+void multiply(A&& a, B&& b, C&& c) {
+  log_trace("");
+  if (__backend::shape(a)[0] != __backend::shape(c)[0] ||
+      __backend::shape(b)[1] != __backend::shape(c)[1] ||
+      __backend::shape(a)[1] != __backend::shape(b)[0]) {
+    throw std::invalid_argument(
+        "multiply: matrix dimensions are incompatible.");
+  }
+
+  using T = tensor_scalar_t<C>;
+  using I = tensor_index_t<C>;
+
+  __backend::spa_accumulator<T, I> c_row(__backend::shape(c)[1]);
+
+  std::vector<std::vector<std::pair<I, T>>> columns(__backend::shape(c)[1]);
+
+  for (auto&& [i, a_row] : __backend::rows(a)) {
+    c_row.clear();
+    for (auto&& [k, a_v] : a_row) {
+      for (auto&& [j, b_v] : __backend::lookup_row(b, k)) {
+        c_row[j] += a_v * b_v;
+      }
+    }
+    for (auto&& [j, v] : c_row.get()) {
+      columns[j].push_back({i, v});
+    }
+  }
+
+  __backend::csc_builder c_builder(c);
+
+  for (std::size_t j = 0; j < columns.size(); j++) {
+    auto&& column = columns[j];
+    std::sort(column.begin(), column.end(),
+              [](auto&& a, auto&& b) { return a.first < b.first; });
+
+    try {
+      c_builder.insert_column(j, column);
+    } catch (...) {
+      throw std::runtime_error("multiply: SpGEMM ran out of memory.");
+    }
+  }
+  c.update(c.values(), c.colptr(), c.rowind(), c.shape(),
+           c.colptr()[c.shape()[1]]);
+}
+
+// C = AB
+// CSR * CSR -> CSC
+// SpGEMM (Gustavson's Algorithm, scattered)
+template <matrix A, matrix B, matrix C>
+  requires(__backend::row_iterable<A> && __backend::row_iterable<B> &&
+           __detail::is_csc_view_v<C>)
+operation_info_t multiply_compute(A&& a, B&& b, C&& c) {
+  log_trace("");
+  if (__backend::shape(a)[0] != __backend::shape(c)[0] ||
+      __backend::shape(b)[1] != __backend::shape(c)[1] ||
+      __backend::shape(a)[1] != __backend::shape(b)[0]) {
+    throw std::invalid_argument(
+        "multiply: matrix dimensions are incompatible.");
+  }
+
+  using T = tensor_scalar_t<C>;
+  using I = tensor_index_t<C>;
+  using O = tensor_offset_t<C>;
+
+  O nnz = 0;
+  __backend::spa_set<I> c_row(__backend::shape(c)[1]);
+
+  for (auto&& [i, a_row] : __backend::rows(a)) {
+    c_row.clear();
+
+    for (auto&& [k, a_v] : a_row) {
+      for (auto&& [j, b_v] : __backend::lookup_row(b, k)) {
+        c_row.insert(j);
+      }
+    }
+
+    nnz += c_row.size();
+  }
+
+  return operation_info_t{__backend::shape(c), nnz};
+}
+
 } // namespace spblas

include/spblas/vendor/armpl/detail/armpl.hpp

@@ -124,6 +124,22 @@ template <>
 inline constexpr auto export_spmat_csr<std::complex<double>> =
     &armpl_spmat_export_csr_z;
 
+template <typename T>
+armpl_status_t (*export_spmat_csc)(armpl_const_spmat_t, armpl_int_t,
+                                   armpl_int_t*, armpl_int_t*,
+                                   const armpl_int_t**, const armpl_int_t**,
+                                   const T**);
+template <>
+inline constexpr auto export_spmat_csc<float> = &armpl_spmat_export_csc_s;
+template <>
+inline constexpr auto export_spmat_csc<double> = &armpl_spmat_export_csc_d;
+template <>
+inline constexpr auto export_spmat_csc<std::complex<float>> =
+    &armpl_spmat_export_csc_c;
+template <>
+inline constexpr auto export_spmat_csc<std::complex<double>> =
+    &armpl_spmat_export_csc_z;
+
 } // namespace __armpl
 
 } // namespace spblas

include/spblas/vendor/armpl/detail/detail.hpp

@@ -2,3 +2,4 @@
 
 #include "armpl.hpp"
 #include "create_matrix_handle.hpp"
+#include "export_matrix_handle.hpp"

include/spblas/vendor/armpl/detail/export_matrix_handle.hpp

@@ -1,31 +1,54 @@
 #pragma once
 
+#include <spblas/detail/operation_info_t.hpp>
+
 #include <armpl_sparse.h>
 #include <spblas/detail/view_inspectors.hpp>
 
 namespace spblas {
 
 namespace __armpl {
 
-template <matrix M, typename O>
+template <matrix M>
   requires __detail::is_csr_view_v<M>
-void export_matrix_handle(operation_info_t& info, M&& m, armpl_spmat_t m_handle) {
+void export_matrix_handle(operation_info_t& info, M&& matrix,
+                          armpl_spmat_t matrix_handle) {
   auto nnz = info.result_nnz();
   armpl_int_t m, n;
   armpl_int_t *rowptr, *colind;
   tensor_scalar_t<M>* values;
-  __armpl::export_spmat_csr<tensor_scalar_t<M>>(m_handle, 0, &m, &n, &rowptr,
-                                                &colind, &values);
+  __armpl::export_spmat_csr<tensor_scalar_t<M>>(matrix_handle, 0, &m, &n,
+                                                &rowptr, &colind, &values);
 
-  std::copy(values, values + nnz, m.values().begin());
-  std::copy(colind, colind + nnz, m.colind().begin());
-  std::copy(rowptr, rowptr + m + 1, m.rowptr().begin());
+  std::copy(values, values + nnz, matrix.values().begin());
+  std::copy(colind, colind + nnz, matrix.colind().begin());
+  std::copy(rowptr, rowptr + m + 1, matrix.rowptr().begin());
 
   free(values);
   free(rowptr);
   free(colind);
 }
 
+template <matrix M>
+  requires __detail::is_csc_view_v<M>
+void export_matrix_handle(operation_info_t& info, M&& matrix,
+                          armpl_spmat_t matrix_handle) {
+  auto nnz = info.result_nnz();
+  armpl_int_t m, n;
+  armpl_int_t *colptr, *rowind;
+  tensor_scalar_t<M>* values;
+  __armpl::export_spmat_csc<tensor_scalar_t<M>>(matrix_handle, 0, &m, &n,
+                                                &rowind, &colptr, &values);
+
+  std::copy(values, values + nnz, matrix.values().begin());
+  std::copy(rowind, rowind + nnz, matrix.rowind().begin());
+  std::copy(colptr, colptr + n + 1, matrix.colptr().begin());
+
+  free(values);
+  free(colptr);
+  free(rowind);
 }
 
-}
+} // namespace __armpl
+
+} // namespace spblas

include/spblas/vendor/armpl/multiply_impl.hpp

@@ -87,7 +87,7 @@ void multiply(A&& a, B&& b, C&& c) {
 template <matrix A, matrix B, matrix C>
   requires((__detail::has_csr_base<A> || __detail::has_csc_base<A>) &&
            (__detail::has_csr_base<B> || __detail::has_csc_base<B>) &&
-           __detail::is_csr_view_v<C>)
+           (__detail::is_csr_view_v<C> || __detail::is_csc_view_v<C>) )
 operation_info_t multiply_compute(A&& a, B&& b, C&& c) {
   log_trace("");
   auto a_base = __detail::get_ultimate_base(a);
@@ -117,42 +117,12 @@ operation_info_t multiply_compute(A&& a, B&& b, C&& c) {
 template <matrix A, matrix B, matrix C>
   requires((__detail::has_csr_base<A> || __detail::has_csc_base<A>) &&
            (__detail::has_csr_base<B> || __detail::has_csc_base<B>) &&
-           __detail::is_csr_view_v<C>)
+           (__detail::is_csr_view_v<C> || __detail::is_csc_view_v<C>) )
 void multiply_fill(operation_info_t& info, A&& a, B&& b, C&& c) {
   log_trace("");
-  auto a_handle = info.state_.a_handle;
-  auto b_handle = info.state_.b_handle;
   auto c_handle = info.state_.c_handle;
 
-  armpl_int_t m, n;
-  auto nnz = info.result_nnz();
-  armpl_int_t *rowptr, *colind;
-  tensor_scalar_t<C>* values;
-  __armpl::export_spmat_csr<tensor_scalar_t<C>>(c_handle, 0, &m, &n, &rowptr,
-                                                &colind, &values);
-
-  std::copy(values, values + nnz, c.values().begin());
-  std::copy(colind, colind + nnz, c.colind().begin());
-  std::copy(rowptr, rowptr + m + 1, c.rowptr().begin());
-
-  free(values);
-  free(rowptr);
-  free(colind);
-}
-
-template <matrix A, matrix B, matrix C>
-  requires __detail::has_csc_base<A> && __detail::has_csc_base<B> &&
-           __detail::is_csc_view_v<C>
-operation_info_t multiply_compute(A&& a, B&& b, C&& c) {
-  return multiply_compute(transposed(b), transposed(a), transposed(c));
-}
-
-template <matrix A, matrix B, matrix C>
-  requires((__detail::has_csr_base<A> || __detail::has_csc_base<A>) &&
-           (__detail::has_csr_base<B> || __detail::has_csc_base<B>) &&
-           __detail::is_csc_view_v<C>)
-void multiply_fill(operation_info_t& info, A&& a, B&& b, C&& c) {
-  multiply_fill(info, transposed(b), transposed(a), transposed(c));
+  __armpl::export_matrix_handle(info, c, c_handle);
 }
 
 } // namespace spblas