From 3e0cca2a151afa79a4617dad13a509a0e6d3c0d5 Mon Sep 17 00:00:00 2001
From: Benjamin Brock <benjamin.brock@intel.com>
Date: Mon, 5 Jun 2023 11:49:18 -0700
Subject: [PATCH 1/3] Begin implementation of spec wrapper for GBTL

---
 src/examples/Makefile                         |  16 ++
 src/examples/test.cpp                         |  25 ++
 src/interfaces/spec/grb/detail/concepts.hpp   | 105 ++++++++
 src/interfaces/spec/grb/detail/cpos.hpp       | 107 ++++++++
 src/interfaces/spec/grb/detail/detail.hpp     |   5 +
 src/interfaces/spec/grb/detail/get.hpp        |  39 +++
 .../spec/grb/detail/matrix_traits.hpp         | 106 ++++++++
 .../spec/grb/detail/namespace_macros.hpp      |   5 +
 src/interfaces/spec/grb/detail/tag_invoke.hpp |  86 +++++++
 src/interfaces/spec/grb/grb.hpp               |   4 +
 src/interfaces/spec/grb/matrix.hpp            | 118 +++++++++
 src/interfaces/spec/grb/scalar_ref.hpp        |  45 ++++
 src/interfaces/spec/grb/util/hints.hpp        |  15 ++
 src/interfaces/spec/grb/util/index.hpp        |  83 +++++++
 src/interfaces/spec/grb/util/matrix_entry.hpp | 233 ++++++++++++++++++
 src/interfaces/spec/grb/util/util.hpp         |   6 +
 16 files changed, 998 insertions(+)
 create mode 100644 src/examples/Makefile
 create mode 100644 src/examples/test.cpp
 create mode 100644 src/interfaces/spec/grb/detail/concepts.hpp
 create mode 100644 src/interfaces/spec/grb/detail/cpos.hpp
 create mode 100644 src/interfaces/spec/grb/detail/detail.hpp
 create mode 100644 src/interfaces/spec/grb/detail/get.hpp
 create mode 100644 src/interfaces/spec/grb/detail/matrix_traits.hpp
 create mode 100644 src/interfaces/spec/grb/detail/namespace_macros.hpp
 create mode 100644 src/interfaces/spec/grb/detail/tag_invoke.hpp
 create mode 100644 src/interfaces/spec/grb/grb.hpp
 create mode 100644 src/interfaces/spec/grb/matrix.hpp
 create mode 100644 src/interfaces/spec/grb/scalar_ref.hpp
 create mode 100644 src/interfaces/spec/grb/util/hints.hpp
 create mode 100644 src/interfaces/spec/grb/util/index.hpp
 create mode 100644 src/interfaces/spec/grb/util/matrix_entry.hpp
 create mode 100644 src/interfaces/spec/grb/util/util.hpp

diff --git a/src/examples/Makefile b/src/examples/Makefile
new file mode 100644
index 0000000..5250c12
--- /dev/null
+++ b/src/examples/Makefile
@@ -0,0 +1,16 @@
+CXX=g++-13
+
+SOURCES += $(wildcard *.cpp)
+TARGETS := $(patsubst %.cpp, %, $(SOURCES))
+
+INCLUDE_FLAGS = -I../ -I../graphblas/platforms/sequential -I../interfaces/spec
+
+CXXFLAGS = -std=c++2b -O3 $(INCLUDE_FLAGS)
+
+all: $(TARGETS)
+
+%: %.cpp
+	$(CXX) $(CXXFLAGS) -o $@ $^ $(CXXFLAGS)
+
+clean:
+	rm -fv $(TARGETS)
diff --git a/src/examples/test.cpp b/src/examples/test.cpp
new file mode 100644
index 0000000..0c9f3c9
--- /dev/null
+++ b/src/examples/test.cpp
@@ -0,0 +1,25 @@
+#include <grb/grb.hpp>
+
+int main(int argc, char** argv) {
+  spec::matrix<int> m({100, 100});
+
+  // Write to missing element.
+  m[{4, 4}] = 12;
+
+  // Access present element.
+  int v = m[{4, 4}];
+  std::cout << v << std::endl;
+
+  // Access missing element.
+  int g = m[{4, 3}];
+  std::cout << g << std::endl;
+
+  // Write to present element.
+  m[{4, 3}] = 12;
+
+  g = m[{4, 3}];
+  std::cout << g << std::endl;
+
+
+  return 0;
+}
diff --git a/src/interfaces/spec/grb/detail/concepts.hpp b/src/interfaces/spec/grb/detail/concepts.hpp
new file mode 100644
index 0000000..ce50e42
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/concepts.hpp
@@ -0,0 +1,105 @@
+
+#pragma once
+
+#include <ranges>
+
+#include "matrix_traits.hpp"
+#include "cpos.hpp"
+#include "get.hpp"
+#include "namespace_macros.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename T, std::size_t I, typename U = GRB_SPEC_NAMESPACE::any>
+concept TupleElementGettable = requires(T tuple) {
+                                 {GRB_SPEC_NAMESPACE::get<I>(tuple)} -> std::convertible_to<U>;
+                               };
+
+template <typename T, typename... Args>
+concept TupleLike =
+  requires {
+    typename std::tuple_size<std::remove_cvref_t<T>>::type;
+    requires std::same_as<std::remove_cvref_t<decltype(std::tuple_size_v<std::remove_cvref_t<T>>)>, std::size_t>;
+  } &&
+  sizeof...(Args) == std::tuple_size_v<std::remove_cvref_t<T>> &&
+  []<std::size_t... I>(std::index_sequence<I...>) {
+    return (TupleElementGettable<T, I, Args> && ...);
+  }(std::make_index_sequence<std::tuple_size_v<std::remove_cvref_t<T>>>());
+
+template <typename Entry, typename T, typename I>
+concept MatrixEntry = TupleLike<Entry, GRB_SPEC_NAMESPACE::any, GRB_SPEC_NAMESPACE::any> &&
+                      requires(Entry entry) { {GRB_SPEC_NAMESPACE::get<0>(entry)} -> TupleLike<I, I>; } &&
+                      requires(Entry entry) { {GRB_SPEC_NAMESPACE::get<1>(entry)} -> std::convertible_to<T>; };
+
+template <typename Entry, typename T, typename I, typename U>
+concept MutableMatrixEntry = MatrixEntry<Entry, T, I> &&
+                             std::is_assignable_v<decltype(GRB_SPEC_NAMESPACE::get<1>(std::declval<Entry>())), U>;
+
+template <typename M>
+concept MatrixRange = std::ranges::sized_range<M> &&
+  requires(M matrix) {
+    typename container_traits<std::remove_cvref_t<M>>;
+    // typename GRB_SPEC_NAMESPACE::matrix_scalar_t<M>;
+    // typename GRB_SPEC_NAMESPACE::matrix_index_t<M>;
+    {std::declval<std::ranges::range_value_t<M>>()}
+      -> MatrixEntry<GRB_SPEC_NAMESPACE::matrix_scalar_t<M>,
+                     GRB_SPEC_NAMESPACE::matrix_index_t<M>>;
+    {GRB_SPEC_NAMESPACE::shape(matrix)} -> TupleLike<GRB_SPEC_NAMESPACE::matrix_index_t<M>,
+                                  GRB_SPEC_NAMESPACE::matrix_index_t<M>>;
+    {GRB_SPEC_NAMESPACE::find(matrix, {GRB_SPEC_NAMESPACE::matrix_index_t<M>{}, GRB_SPEC_NAMESPACE::matrix_index_t<M>{}})} -> std::convertible_to<std::ranges::iterator_t<M>>;
+  };
+
+template <typename M, typename T>
+concept MutableMatrixRange = GRB_SPEC_NAMESPACE::MatrixRange<M> &&
+                             GRB_SPEC_NAMESPACE::MutableMatrixEntry<std::ranges::range_reference_t<M>,
+                                                  GRB_SPEC_NAMESPACE::matrix_scalar_t<M>,
+                                                  GRB_SPEC_NAMESPACE::matrix_index_t<M>,
+                                                  T> &&
+  requires(M matrix, T value) {
+    {GRB_SPEC_NAMESPACE::insert(matrix, {{GRB_SPEC_NAMESPACE::matrix_index_t<M>{}, GRB_SPEC_NAMESPACE::matrix_index_t<M>{}}, value})}
+      -> std::same_as<std::pair<std::ranges::iterator_t<M>, bool>>;
+  } &&
+  std::is_constructible_v<GRB_SPEC_NAMESPACE::matrix_scalar_t<M>, T>;
+  
+template <typename M>
+concept MaskMatrixRange = MatrixRange<M> &&
+                          std::is_convertible_v<GRB_SPEC_NAMESPACE::matrix_scalar_t<M>, bool>;
+
+template <typename Entry, typename T, typename I>
+concept VectorEntry = TupleLike<Entry, GRB_SPEC_NAMESPACE::any, GRB_SPEC_NAMESPACE::any> &&
+                      requires(Entry entry) { {GRB_SPEC_NAMESPACE::get<0>(entry)} -> std::integral; } &&
+                      requires(Entry entry) { {GRB_SPEC_NAMESPACE::get<1>(entry)} -> std::convertible_to<T>; };
+
+template <typename Entry, typename T, typename I, typename U>
+concept MutableVectorEntry = VectorEntry<Entry, T, I> &&
+                             std::is_assignable_v<decltype(GRB_SPEC_NAMESPACE::get<1>(std::declval<Entry>())), U>;
+
+template <typename V>
+concept VectorRange = std::ranges::sized_range<V> &&
+  requires(V vector) {
+    typename GRB_SPEC_NAMESPACE::vector_scalar_t<V>;
+    typename GRB_SPEC_NAMESPACE::vector_index_t<V>;
+    {std::declval<std::ranges::range_value_t<V>>()}
+      -> VectorEntry<GRB_SPEC_NAMESPACE::vector_scalar_t<V>,
+                     GRB_SPEC_NAMESPACE::vector_index_t<V>>;
+    {GRB_SPEC_NAMESPACE::shape(vector)} -> std::same_as<GRB_SPEC_NAMESPACE::vector_index_t<V>>;
+    {GRB_SPEC_NAMESPACE::find(vector, GRB_SPEC_NAMESPACE::vector_index_t<V>{})} -> std::convertible_to<std::ranges::iterator_t<V>>;
+};
+
+template <typename V, typename T>
+concept MutableVectorRange = VectorRange<V> &&
+                             MutableVectorEntry<std::ranges::range_reference_t<V>,
+                                                GRB_SPEC_NAMESPACE::vector_scalar_t<V>,
+                                                GRB_SPEC_NAMESPACE::vector_index_t<V>,
+                                                T> &&
+  requires(V vector, T value) {
+    {GRB_SPEC_NAMESPACE::insert(vector, {GRB_SPEC_NAMESPACE::vector_index_t<V>{}, GRB_SPEC_NAMESPACE::vector_scalar_t<V>{}})}
+      -> std::same_as<std::pair<std::ranges::iterator_t<V>, bool>>;
+  } &&
+  std::is_constructible_v<GRB_SPEC_NAMESPACE::vector_scalar_t<V>, T>;
+
+template <typename M>
+concept MaskVectorRange = VectorRange<M> &&
+                          std::is_convertible_v<GRB_SPEC_NAMESPACE::vector_scalar_t<M>, bool>;
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/detail/cpos.hpp b/src/interfaces/spec/grb/detail/cpos.hpp
new file mode 100644
index 0000000..e54a08e
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/cpos.hpp
@@ -0,0 +1,107 @@
+#pragma once
+
+#include <type_traits>
+#include <concepts>
+#include <utility>
+#include <ranges>
+#include <any>
+#include "tag_invoke.hpp"
+#include "matrix_traits.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+// Helper concepts for CPOs.
+
+namespace {
+
+template <typename T>
+concept has_matrix_shape = requires(T t) { {t.shape()} -> std::same_as<typename std::tuple_element<0, std::ranges::range_value_t<T>>::type>; };
+
+template <typename T>
+concept has_vector_shape = requires(T t) { {t.shape()} -> std::same_as<std::ranges::range_value_t<T>>; };
+
+template <typename T>
+concept has_find_method = requires(T t) { {t.find(std::declval<typename container_traits<T>::key_type>())} -> std::same_as<typename container_traits<T>::iterator>; };
+
+template <typename T>
+concept has_insert_method = requires(T t) { {t.insert({std::declval<typename container_traits<T>::key_type>(), std::declval<container_scalar_t<T>>()})}; };
+
+template <typename T, typename M>
+concept has_insert_or_assign_method = requires(T t, M obj) { {t.insert_or_assign(std::declval<typename container_traits<T>::key_type>(), std::forward<M>(obj))}; };
+
+} // end anonymous
+
+inline constexpr struct shape_fn_ {
+  template <typename T>
+  auto operator()(T&& x) const
+    requires(GRB_SPEC_NAMESPACE::is_tag_invocable_v<shape_fn_, T> ||
+             has_matrix_shape<T> ||
+             has_vector_shape<T>)
+    {
+    if constexpr(GRB_SPEC_NAMESPACE::is_tag_invocable_v<shape_fn_, T>) {
+      return GRB_SPEC_NAMESPACE::tag_invoke(*this, std::forward<T>(x));
+    } else if constexpr(has_matrix_shape<T>) {
+      return std::forward<T>(x).shape();
+    } else if constexpr(has_vector_shape<T>) {
+      return std::forward<T>(x).shape();
+    }
+  }
+} shape{};
+
+inline constexpr struct size_fn_ {
+  template <typename T>
+  auto operator()(T&& x) const
+  requires(GRB_SPEC_NAMESPACE::is_tag_invocable_v<size_fn_, T> ||
+           std::ranges::sized_range<T>)
+  {
+    if constexpr(GRB_SPEC_NAMESPACE::is_tag_invocable_v<size_fn_, T>) {
+      return GRB_SPEC_NAMESPACE::tag_invoke(*this, std::forward<T>(x));
+    } else if constexpr(std::ranges::sized_range<T>) {
+      return std::ranges::size(std::forward<T>(x));
+    }
+  }
+} size{};
+
+inline constexpr struct find_fn_ {
+  template <typename T>
+  auto operator()(T&& x, typename GRB_SPEC_NAMESPACE::container_traits<T>::key_type key) const
+  requires(GRB_SPEC_NAMESPACE::is_tag_invocable_v<find_fn_, T, typename GRB_SPEC_NAMESPACE::container_traits<T>::key_type> ||
+           has_find_method<T>)
+  {
+    if constexpr(GRB_SPEC_NAMESPACE::is_tag_invocable_v<find_fn_, T, typename GRB_SPEC_NAMESPACE::container_traits<T>::key_type>) {
+      return GRB_SPEC_NAMESPACE::tag_invoke(*this, std::forward<T>(x), key);
+    } else if constexpr(has_find_method<T>) {
+      return std::forward<T>(x).find(key);
+    }
+  }
+} find{};
+
+inline constexpr struct insert_fn_ {
+  template <typename T>
+  auto operator()(T&& x, const container_value_t<T>& entry) const
+  requires(GRB_SPEC_NAMESPACE::is_tag_invocable_v<insert_fn_, T, const container_value_t<T>&> ||
+           has_insert_method<T>)
+  {
+    if constexpr(GRB_SPEC_NAMESPACE::is_tag_invocable_v<insert_fn_, T, const container_value_t<T>&>) {
+      return GRB_SPEC_NAMESPACE::tag_invoke(*this, std::forward<T>(x), entry);
+    } else if constexpr(has_insert_method<T>) {
+      return std::forward<T>(x).insert(entry);
+    }
+  }
+} insert{};
+
+inline constexpr struct insert_or_assign_fn_ {
+  template <typename T, typename M>
+  auto operator()(T&& x, typename GRB_SPEC_NAMESPACE::container_traits<T>::key_type key, M&& obj) const
+  requires(GRB_SPEC_NAMESPACE::is_tag_invocable_v<insert_or_assign_fn_, T, typename GRB_SPEC_NAMESPACE::container_traits<T>::key_type, M> ||
+           has_insert_or_assign_method<T, M>)
+  {
+    if constexpr(GRB_SPEC_NAMESPACE::is_tag_invocable_v<insert_or_assign_fn_, T, typename GRB_SPEC_NAMESPACE::container_traits<T>::key_type, M>) {
+      return GRB_SPEC_NAMESPACE::tag_invoke(*this, std::forward<T>(x), key, std::forward<M>(obj));
+    } else if constexpr(has_insert_or_assign_method<T, M>) {
+      return std::forward<T>(x).insert_or_assign(key, std::forward<M>(obj));
+    }
+  }
+} insert_or_assign{};
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/detail/detail.hpp b/src/interfaces/spec/grb/detail/detail.hpp
new file mode 100644
index 0000000..f17bf44
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/detail.hpp
@@ -0,0 +1,5 @@
+
+#pragma once
+
+#include "namespace_macros.hpp"
+#include "concepts.hpp"
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/detail/get.hpp b/src/interfaces/spec/grb/detail/get.hpp
new file mode 100644
index 0000000..64027a7
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/get.hpp
@@ -0,0 +1,39 @@
+#pragma once
+
+#include <tuple>
+#include <concepts>
+
+#include "namespace_macros.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename T, std::size_t I>
+concept method_gettable = requires(T tuple) {
+  {tuple. template get<I>()} -> std::convertible_to<GRB_SPEC_NAMESPACE::any>;
+};
+
+template <typename T, std::size_t I>
+concept adl_gettable = requires(T tuple) {
+  {get<I>(tuple)} -> std::convertible_to<GRB_SPEC_NAMESPACE::any>;
+};
+
+template <typename T, std::size_t I>
+concept std_gettable = requires(T tuple) {
+  {std::get<I>(tuple)} -> std::convertible_to<GRB_SPEC_NAMESPACE::any>;
+};
+
+template <std::size_t I, typename T>
+inline constexpr decltype(auto) get(T&& tuple)
+requires(method_gettable<T, I>)
+{
+  return std::forward<T>(tuple). template get<I>();
+}
+
+template <std::size_t I, typename T>
+inline constexpr decltype(auto) get(T&& tuple)
+requires(!method_gettable<T, I> && std_gettable<T, I>)
+{
+  return std::get<I>(std::forward<T>(tuple));
+}
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/detail/matrix_traits.hpp b/src/interfaces/spec/grb/detail/matrix_traits.hpp
new file mode 100644
index 0000000..547d62d
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/matrix_traits.hpp
@@ -0,0 +1,106 @@
+#pragma once
+
+#include <type_traits>
+#include <concepts>
+#include <utility>
+#include <ranges>
+#include <tuple>
+#include <any>
+
+#include "tag_invoke.hpp"
+#include "namespace_macros.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+using any = std::any;
+
+namespace {
+
+template <typename Container>
+struct get_index_type {
+private:
+  using value_type = std::ranges::range_value_t<Container>;
+  using key_type = typename std::tuple_element<0, value_type>::type;
+public:
+  using type = key_type;
+};
+
+template <typename Container>
+requires requires { typename std::tuple_element<0, typename std::tuple_element<0, std::ranges::range_value_t<Container>>::type>::type; }
+struct get_index_type<Container>
+{
+private:
+  using value_type = std::ranges::range_value_t<Container>;
+  using key_type = typename std::tuple_element<0, value_type>::type;
+public:
+  using type = typename std::tuple_element<0, key_type>::type;
+};
+}
+
+template <typename Container>
+struct container_traits {
+  using size_type = std::ranges::range_size_t<Container>;
+  using difference_type = std::ranges::range_difference_t<Container>;
+  using iterator = std::ranges::iterator_t<Container>;
+  using value_type = std::ranges::range_value_t<Container>;
+  using key_type = std::remove_cvref_t<typename std::tuple_element<0, value_type>::type>;
+  using map_type = std::remove_cvref_t<typename std::tuple_element<1, value_type>::type>;
+  using scalar_type = std::remove_cvref_t<typename std::tuple_element<1, value_type>::type>;
+  using index_type = std::remove_cvref_t<typename get_index_type<Container>::type>;
+  using reference = std::ranges::range_reference_t<Container>;
+};
+
+template <typename Container>
+struct matrix_traits {
+  using size_type = std::ranges::range_size_t<Container>;
+  using difference_type = std::ranges::range_difference_t<Container>;
+  using iterator = std::ranges::iterator_t<Container>;
+  using value_type = std::ranges::range_value_t<Container>;
+  using key_type = std::remove_cvref_t<typename std::tuple_element<0, value_type>::type>;
+  using map_type = std::remove_cvref_t<typename std::tuple_element<1, value_type>::type>;
+  using scalar_type = std::remove_cvref_t<typename std::tuple_element<1, value_type>::type>;
+  using index_type = std::remove_cvref_t<typename std::tuple_element<0, key_type>::type>;
+  using reference = std::ranges::range_reference_t<Container>;
+};
+
+template <typename Container>
+struct vector_traits {
+  using size_type = std::ranges::range_size_t<Container>;
+  using difference_type = std::ranges::range_difference_t<Container>;
+  using iterator = std::ranges::iterator_t<Container>;
+  using value_type = std::ranges::range_value_t<Container>;
+  using key_type = std::remove_cvref_t<typename std::tuple_element<0, value_type>::type>;
+  using map_type = std::remove_cvref_t<typename std::tuple_element<1, value_type>::type>;
+  using scalar_type = std::remove_cvref_t<typename std::tuple_element<1, value_type>::type>;
+  using index_type = key_type;
+  using reference = std::ranges::range_reference_t<Container>;
+};
+
+template <typename T>
+using container_scalar_t = typename container_traits<T>::scalar_type;
+
+template <typename T>
+using container_index_t = typename container_traits<T>::index_type;
+
+template <typename T>
+using container_value_t = typename container_traits<T>::value_type;
+
+template <typename T>
+using container_key_t = typename container_traits<T>::key_type;
+
+template <typename T>
+using matrix_scalar_t = typename container_traits<std::remove_cvref_t<T>>::scalar_type;
+
+template <typename T>
+using matrix_index_t = typename container_traits<std::remove_cvref_t<T>>::index_type;
+
+template <typename T>
+using vector_scalar_t = typename container_traits<std::remove_cvref_t<T>>::scalar_type;
+
+template <typename T>
+using vector_index_t = typename container_traits<std::remove_cvref_t<T>>::index_type;
+
+template <typename A, typename B, typename Fn>
+using elementwise_return_type_t = std::invoke_result_t<Fn, container_scalar_t<A>, container_scalar_t<B>>;
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/detail/namespace_macros.hpp b/src/interfaces/spec/grb/detail/namespace_macros.hpp
new file mode 100644
index 0000000..71053a9
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/namespace_macros.hpp
@@ -0,0 +1,5 @@
+
+#pragma once
+
+#define GRB_SPEC_NAMESPACE spec
+#define GBTL_NAMESPACE grb
diff --git a/src/interfaces/spec/grb/detail/tag_invoke.hpp b/src/interfaces/spec/grb/detail/tag_invoke.hpp
new file mode 100644
index 0000000..57fdadb
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/tag_invoke.hpp
@@ -0,0 +1,86 @@
+#pragma once
+
+#include <utility>
+#include <type_traits>
+
+#include "namespace_macros.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+namespace _tag_invoke {
+  void tag_invoke();
+
+  struct _fn {
+    template <typename CPO, typename... Args>
+    constexpr auto operator()(CPO cpo, Args&&... args) const
+        noexcept(noexcept(tag_invoke((CPO &&) cpo, (Args &&) args...)))
+        -> decltype(tag_invoke((CPO &&) cpo, (Args &&) args...)) {
+      return tag_invoke((CPO &&) cpo, (Args &&) args...);
+    }
+  };
+
+  template <typename CPO, typename... Args>
+  using tag_invoke_result_t = decltype(
+      tag_invoke(std::declval<CPO>(), std::declval<Args>()...));
+
+  using yes_type = char;
+  using no_type = char(&)[2];
+
+  template <typename CPO, typename... Args>
+  auto try_tag_invoke(int) //
+      noexcept(noexcept(tag_invoke(
+          std::declval<CPO>(), std::declval<Args>()...)))
+      -> decltype(static_cast<void>(tag_invoke(
+          std::declval<CPO>(), std::declval<Args>()...)), yes_type{});
+
+  template <typename CPO, typename... Args>
+  no_type try_tag_invoke(...) noexcept(false);
+
+  template <template <typename...> class T, typename... Args>
+  struct defer {
+    using type = T<Args...>;
+  };
+
+  struct empty {};
+}  // namespace _tag_invoke
+
+namespace _tag_invoke_cpo {
+  inline constexpr _tag_invoke::_fn tag_invoke{};
+}
+using namespace _tag_invoke_cpo;
+
+template <auto& CPO>
+using tag_t = std::remove_cvref_t<decltype(CPO)>;
+
+using _tag_invoke::tag_invoke_result_t;
+
+template <typename CPO, typename... Args>
+inline constexpr bool is_tag_invocable_v =
+    (sizeof(_tag_invoke::try_tag_invoke<CPO, Args...>(0)) ==
+     sizeof(_tag_invoke::yes_type));
+
+template <typename CPO, typename... Args>
+struct tag_invoke_result
+  : std::conditional_t<
+        is_tag_invocable_v<CPO, Args...>,
+        _tag_invoke::defer<tag_invoke_result_t, CPO, Args...>,
+        _tag_invoke::empty> 
+{};
+
+template <typename CPO, typename... Args>
+using is_tag_invocable = std::bool_constant<is_tag_invocable_v<CPO, Args...>>;
+
+template <typename CPO, typename... Args>
+inline constexpr bool is_nothrow_tag_invocable_v =
+    noexcept(_tag_invoke::try_tag_invoke<CPO, Args...>(0));
+
+template <typename CPO, typename... Args>
+using is_nothrow_tag_invocable =
+    std::bool_constant<is_nothrow_tag_invocable_v<CPO, Args...>>;
+
+template <typename CPO, typename... Args>
+concept tag_invocable =
+    (sizeof(_tag_invoke::try_tag_invoke<CPO, Args...>(0)) ==
+     sizeof(_tag_invoke::yes_type));
+
+} // end GRB_SPEC_NAMESPACE
diff --git a/src/interfaces/spec/grb/grb.hpp b/src/interfaces/spec/grb/grb.hpp
new file mode 100644
index 0000000..7e4a17f
--- /dev/null
+++ b/src/interfaces/spec/grb/grb.hpp
@@ -0,0 +1,4 @@
+
+#pragma once
+
+#include "matrix.hpp"
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/matrix.hpp b/src/interfaces/spec/grb/matrix.hpp
new file mode 100644
index 0000000..74aba03
--- /dev/null
+++ b/src/interfaces/spec/grb/matrix.hpp
@@ -0,0 +1,118 @@
+
+#pragma once
+
+#include <graphblas/graphblas.hpp>
+
+#include "detail/detail.hpp"
+#include "util/util.hpp"
+#include "scalar_ref.hpp"
+
+#include <concepts>
+#include <memory>
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename T,
+          std::integral I = std::size_t,
+          typename Hint = GRB_SPEC_NAMESPACE::sparse,
+          typename Allocator = std::allocator<T>>
+class matrix {
+public:
+  /// Type of scalar elements stored in the matrix.
+  using scalar_type = T;
+
+  /// Type used to reference the indices of elements stored in the matrix.
+  using index_type = I;
+
+  /// A tuple-like type containing:
+  /// 1. A tuple-like type with two `index_type` elements storing the row and column
+  ///    index of the element.
+  /// 2. An element of type `scalar_type`.
+  using value_type = matrix_entry<T, I>;
+
+  using key_type = index<I>;
+  using map_type = T;
+
+  /// Allocator type
+  using allocator_type = Allocator;
+
+  /// A large unsigned integral type
+  using size_type = std::size_t;
+
+  /// A large signed integral type
+  using difference_type = std::ptrdiff_t;
+
+  using hint_type = Hint;
+
+  using backend_type = GBTL_NAMESPACE::Matrix<T>;
+
+  using scalar_reference = scalar_ref<T, I, backend_type>;
+  using const_scalar_reference = scalar_ref<const T, I, const backend_type>;
+
+  /// Construct an empty matrix of dimension `shape[0]` x `shape[1]`
+  matrix(index<I> shape) : backend_(shape[0], shape[1]) {}
+
+  matrix(std::initializer_list<I> shape)
+    : backend_(*shape.begin(), *(shape.begin()+1)) {}
+
+  index<I> shape() const {
+    return index<I>(backend_.nrows(), backend_.ncols());
+  }
+
+  size_type size() const {
+    return backend_.nvals();
+  }
+
+  template <class M>
+  void insert_or_assign(key_type k, M&& obj) {
+    backend_.setElement(k[0], k[1], std::forward<M>(obj));
+  }
+
+  void insert(const value_type& value) {
+    auto&& [k, v] = value;
+
+    if (!backend_.hasElement(k[0], k[1])) {
+      backend_.setElement(k[0], k[1], v);
+    }
+  }
+
+  void insert(value_type&& value) {
+    auto&& [k, v] = value;
+
+    if (!backend_.hasElement(k[0], k[1])) {
+      backend_.setElement(k[0], k[1], std::move(v));
+    }
+  }
+
+  template <typename InputIt>
+  void insert(InputIt first, InputIt last) {
+    for (auto it = first; it != last; ++it) {
+      auto&& [k, v] = *it;
+      insert_or_assign(k, v);
+    }
+  }
+
+  scalar_reference operator[](key_type index) {
+    if (!backend_.hasElement(index[0], index[1])) {
+      backend_.setElement(index[0], index[1], T());
+    }
+    return scalar_reference(index, &backend_);
+  }
+
+  void clear() {
+    backend_.clear();
+  }
+
+  matrix() = default;
+  matrix(const Allocator& allocator) {}
+
+  matrix(const matrix&) = default;
+  matrix(matrix&&) = default;
+  matrix& operator=(const matrix&) = default;
+  matrix& operator=(matrix&&) = default;
+
+private:
+  backend_type backend_;
+};
+
+};
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/scalar_ref.hpp b/src/interfaces/spec/grb/scalar_ref.hpp
new file mode 100644
index 0000000..b21aa16
--- /dev/null
+++ b/src/interfaces/spec/grb/scalar_ref.hpp
@@ -0,0 +1,45 @@
+
+#pragma once
+
+#include <graphblas/graphblas.hpp>
+
+#include "detail/detail.hpp"
+#include "util/util.hpp"
+
+#include <concepts>
+#include <memory>
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename T, std::integral I, typename Matrix>
+class scalar_ref {
+public:
+  scalar_ref() = delete;
+  ~scalar_ref() = default;
+  scalar_ref(const scalar_ref &) = default;
+
+  scalar_ref(index<I> index, Matrix* matrix) : index_(index), matrix_(matrix) {}
+
+  operator T() const {
+    return matrix_->extractElement(index_[0], index_[1]);
+  }
+
+  scalar_ref operator=(const T &value) const
+    requires(!std::is_const_v<T>)
+  {
+    matrix_->setElement(index_[0], index_[1], value);
+    return *this;
+  }
+
+  scalar_ref operator=(const scalar_ref &other) const {
+    T value = other;
+    *this = value;
+    return *this;
+  }
+
+private:
+  index<I> index_;
+  Matrix* matrix_;
+};
+
+} // end GRB_SPEC_NAMESPACE
diff --git a/src/interfaces/spec/grb/util/hints.hpp b/src/interfaces/spec/grb/util/hints.hpp
new file mode 100644
index 0000000..a3c91e4
--- /dev/null
+++ b/src/interfaces/spec/grb/util/hints.hpp
@@ -0,0 +1,15 @@
+#pragma once
+
+#include "../detail/detail.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+struct sparse {};
+struct dense {};
+struct row {};
+struct column {};
+
+template <typename... Hints>
+struct compose {};
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/util/index.hpp b/src/interfaces/spec/grb/util/index.hpp
new file mode 100644
index 0000000..c56e96b
--- /dev/null
+++ b/src/interfaces/spec/grb/util/index.hpp
@@ -0,0 +1,83 @@
+
+#pragma once
+
+#include <tuple>
+#include <concepts>
+#include <limits>
+
+#include "../detail/namespace_macros.hpp"
+#include "../detail/concepts.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <std::integral T = std::size_t>
+class index {
+public:
+  using index_type = T;
+
+  using first_type = T;
+  using second_type = T;
+
+  constexpr index_type operator[](index_type dim) const noexcept {
+    if (dim == 0) {
+      return first;
+    } else {
+      return second;
+    }
+  }
+
+  template <std::integral U>
+  requires(std::numeric_limits<U>::max() >= std::numeric_limits<T>::max())
+  constexpr operator index<U>() const noexcept {
+    return index<U>(first, second);
+  }
+
+  constexpr index(index_type first, index_type second) : first(first), second(second) {}
+
+  template <GRB_SPEC_NAMESPACE::TupleLike<T, T> Tuple>
+  constexpr index(Tuple tuple) : first(GRB_SPEC_NAMESPACE::get<0>(tuple)), second(GRB_SPEC_NAMESPACE::get<1>(tuple)) {}
+
+  constexpr bool operator==(const index&) const noexcept = default;
+
+  template <std::size_t Index>
+  constexpr T get() const noexcept
+  requires(Index <= 1)
+  {
+    if constexpr(Index == 0) { return first; }
+    if constexpr(Index == 1) { return second; }
+  }
+
+  index() = default;
+  ~index() = default;
+  index(const index&) = default;
+  index& operator=(const index&) = default;
+  index(index&&) = default;
+  index& operator=(index&&) = default;
+
+  index_type first;
+  index_type second;
+};
+
+} // end GRB_SPEC_NAMESPACE
+
+namespace std {
+
+template <std::size_t Index, std::integral I>
+struct tuple_element<Index, GRB_SPEC_NAMESPACE::index<I>>
+  : tuple_element<Index, std::tuple<I, I>> 
+{};
+
+template <std::integral I>
+struct tuple_size<GRB_SPEC_NAMESPACE::index<I>>
+    : integral_constant<std::size_t, 2> {};
+
+template <std::size_t Index, std::integral I>
+inline constexpr I get(GRB_SPEC_NAMESPACE::index<I> index)
+requires(Index <= 1)
+{
+  if constexpr(Index == 0) { return index.first; }
+  if constexpr(Index == 1) { return index.second; }
+}
+
+} // end std
+
diff --git a/src/interfaces/spec/grb/util/matrix_entry.hpp b/src/interfaces/spec/grb/util/matrix_entry.hpp
new file mode 100644
index 0000000..8a6692a
--- /dev/null
+++ b/src/interfaces/spec/grb/util/matrix_entry.hpp
@@ -0,0 +1,233 @@
+#pragma once
+
+#include <type_traits>
+#include <concepts>
+#include <limits>
+
+#include "index.hpp"
+#include "../detail/concepts.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename T,
+          typename I = std::size_t>
+class matrix_entry {
+public:
+  using index_type = I;
+  using map_type = T;
+
+  matrix_entry(GRB_SPEC_NAMESPACE::index<I> index, const map_type& value) : value_(value), index_(index) {}
+  matrix_entry(GRB_SPEC_NAMESPACE::index<I> index, map_type&& value) : value_(std::move(value)), index_(index) {}
+
+  template <typename U>
+  requires(std::is_constructible_v<T, U>)
+  matrix_entry(GRB_SPEC_NAMESPACE::index<I> index, U&& value) : value_(std::forward<U>(value)), index_(index) {}
+
+  template <typename Entry>
+  matrix_entry(Entry&& entry) : index_(GRB_SPEC_NAMESPACE::get<0>(entry)),
+                                value_(GRB_SPEC_NAMESPACE::get<1>(entry))
+                                {}
+
+  template <std::size_t Index>
+  requires(Index <= 1)
+  auto get() const noexcept {
+    if constexpr(Index == 0) { return index(); }
+    if constexpr(Index == 1) { return value(); }
+  }
+
+  operator std::pair<std::pair<I, I>, T>() const noexcept {
+    return {{index_[0], index_[1]}, value_};
+  }
+
+  GRB_SPEC_NAMESPACE::index<I> index() const noexcept {
+    return index_;
+  }
+
+  map_type value() const noexcept {
+    return value_;
+  }
+
+  template <std::integral U>
+  requires(!std::is_same_v<I, U> &&
+           std::numeric_limits<U>::max() >= std::numeric_limits<I>::max())
+  operator matrix_entry<T, U>() const noexcept {
+    return matrix_entry<T, U>(index_, value_);
+  }
+
+  template <std::integral U>
+  requires(!std::is_const_v<T>   &&
+           !std::is_same_v<I, U> &&
+           std::numeric_limits<U>::max() >= std::numeric_limits<I>::max())
+  operator matrix_entry<std::add_const_t<T>, U>() const noexcept {
+    return matrix_entry<std::add_const_t<T>, U>(index_, value_);
+  }
+
+  bool operator<(const matrix_entry& other) const noexcept {
+    if (index()[0] < other.index()[0]) {
+      return true;
+    } else if (index()[0] == other.index()[0] &&
+               index()[1] < other.index()[1]) {
+      return true;
+    }
+    return false;
+  }
+
+  matrix_entry() = default;
+  ~matrix_entry() = default;
+
+  matrix_entry(const matrix_entry&) = default;
+  matrix_entry(matrix_entry&&) = default;
+  matrix_entry& operator=(const matrix_entry&) = default;
+  matrix_entry& operator=(matrix_entry&&) = default;
+
+private:
+  GRB_SPEC_NAMESPACE::index<I> index_;
+  map_type value_;
+};
+
+  
+} // end GRB_SPEC_NAMESPACE
+
+namespace std {
+
+template <typename T, typename I>
+requires(!std::is_const_v<T>)
+void swap(GRB_SPEC_NAMESPACE::matrix_entry<T, I> a, GRB_SPEC_NAMESPACE::matrix_entry<T, I> b) {
+  GRB_SPEC_NAMESPACE::matrix_entry<T, I> other = a;
+  a = b;
+  b = other;
+}
+
+template <std::size_t Index, typename T, typename I>
+struct tuple_element<Index, GRB_SPEC_NAMESPACE::matrix_entry<T, I>>
+  : tuple_element<Index, std::tuple<GRB_SPEC_NAMESPACE::index<I>, T>>
+{
+};
+
+template <typename T, typename I>
+struct tuple_size<GRB_SPEC_NAMESPACE::matrix_entry<T, I>>
+    : integral_constant<size_t, 2> {};
+    
+
+} // end std
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename T,
+          typename I = std::size_t,
+          typename TRef = T&>
+class matrix_ref {
+public:
+  using scalar_type = T;
+  using index_type = I;
+  
+  using key_type = GRB_SPEC_NAMESPACE::index<I>;
+  using map_type = T;
+
+  using scalar_reference = TRef;
+
+  using value_type = GRB_SPEC_NAMESPACE::matrix_entry<T, I>;
+
+  matrix_ref(GRB_SPEC_NAMESPACE::index<I> index, scalar_reference value) : index_(index), value_(value) {}
+
+  operator value_type() const noexcept {
+    return value_type(index_, value_);
+  }
+
+  operator std::pair<std::pair<I, I>, T>() const noexcept {
+    return {{index_[0], index_[1]}, value_};
+  }
+
+  template <std::size_t Index>
+  requires(Index <= 1)
+  decltype(auto) get() const noexcept
+  {
+    if constexpr(Index == 0) { return index(); }
+    if constexpr(Index == 1) { return value(); }
+  }
+
+  GRB_SPEC_NAMESPACE::index<I> index() const noexcept {
+    return index_;
+  }
+
+  scalar_reference value() const noexcept {
+    return value_;
+  }
+
+  template <std::integral U>
+  requires(!std::is_same_v<I, U> &&
+           std::numeric_limits<U>::max() >= std::numeric_limits<I>::max())
+  operator matrix_ref<T, U, TRef>() const noexcept {
+    return matrix_ref<T, U, TRef>(index_, value_);
+  }
+
+  template <std::integral U>
+  requires(!std::is_const_v<T>   &&
+           !std::is_same_v<I, U> &&
+           std::numeric_limits<U>::max() >= std::numeric_limits<I>::max())
+  operator matrix_ref<std::add_const_t<T>, U, TRef>() const noexcept {
+    return matrix_ref<std::add_const_t<T>, U, TRef>(index_, value_);
+  }
+
+  bool operator<(matrix_entry<T, I> other) const noexcept {
+    if (index()[0] < other.index()[0]) {
+      return true;
+    } else if (index()[0] == other.index()[0] &&
+               index()[1] < other.index()[1]) {
+      return true;
+    }
+    return false;
+  }
+
+  matrix_ref() = delete;
+  ~matrix_ref() = default;
+
+  matrix_ref(const matrix_ref&) = default;
+  matrix_ref& operator=(const matrix_ref&) = delete;
+  matrix_ref(matrix_ref&&) = default;
+  matrix_ref& operator=(matrix_ref&&) = default;
+
+private:
+  GRB_SPEC_NAMESPACE::index<I> index_;
+  scalar_reference value_;
+};
+
+  
+} // end GRB_SPEC_NAMESPACE
+
+namespace std {
+
+template <typename T, typename I, typename TRef>
+requires(!std::is_const_v<T>)
+void swap(GRB_SPEC_NAMESPACE::matrix_ref<T, I, TRef> a, GRB_SPEC_NAMESPACE::matrix_ref<T, I, TRef> b) {
+  GRB_SPEC_NAMESPACE::matrix_entry<T, I> other = a;
+  a = b;
+  b = other;
+}
+
+template <std::size_t Index, typename T, typename I, typename TRef>
+struct tuple_element<Index, GRB_SPEC_NAMESPACE::matrix_ref<T, I, TRef>>
+  : tuple_element<Index, std::tuple<GRB_SPEC_NAMESPACE::index<I>, TRef>> {};
+
+template <typename T, typename I, typename TRef>
+struct tuple_size<GRB_SPEC_NAMESPACE::matrix_ref<T, I, TRef>>
+    : integral_constant<std::size_t, 2> {};
+
+template <std::size_t Index, typename T, typename I, typename TRef>
+inline decltype(auto) get(GRB_SPEC_NAMESPACE::matrix_ref<T, I, TRef> ref)
+requires(Index <= 1)
+{
+  if constexpr(Index == 0) { return ref.index(); }
+  if constexpr(Index == 1) { return ref.value(); }
+}
+
+template <std::size_t Index, typename T, typename I>
+inline decltype(auto) get(GRB_SPEC_NAMESPACE::matrix_entry<T, I> entry)
+requires(Index <= 1)
+{
+  if constexpr(Index == 0) { return entry.index(); }
+  if constexpr(Index == 1) { return entry.value(); }
+}
+
+} // end std
+
diff --git a/src/interfaces/spec/grb/util/util.hpp b/src/interfaces/spec/grb/util/util.hpp
new file mode 100644
index 0000000..53ad3a1
--- /dev/null
+++ b/src/interfaces/spec/grb/util/util.hpp
@@ -0,0 +1,6 @@
+
+#pragma once
+
+#include "hints.hpp"
+#include "index.hpp"
+#include "matrix_entry.hpp"
\ No newline at end of file

From db5f10be4609b9a0c249baa60158e0190c2617f1 Mon Sep 17 00:00:00 2001
From: Benjamin Brock <benjamin.brock@intel.com>
Date: Mon, 5 Jun 2023 12:54:11 -0700
Subject: [PATCH 2/3] Implement `erase`

---
 src/interfaces/spec/grb/matrix.hpp | 32 ++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/src/interfaces/spec/grb/matrix.hpp b/src/interfaces/spec/grb/matrix.hpp
index 74aba03..c928bce 100644
--- a/src/interfaces/spec/grb/matrix.hpp
+++ b/src/interfaces/spec/grb/matrix.hpp
@@ -9,6 +9,7 @@
 
 #include <concepts>
 #include <memory>
+#include <limits>
 
 namespace GRB_SPEC_NAMESPACE {
 
@@ -63,6 +64,14 @@ class matrix {
     return backend_.nvals();
   }
 
+  bool empty() const {
+    return size() == 0;
+  }
+
+  std::size_t max_size() const {
+    return std::numeric_limits<I>::max() - 1;
+  }
+
   template <class M>
   void insert_or_assign(key_type k, M&& obj) {
     backend_.setElement(k[0], k[1], std::forward<M>(obj));
@@ -99,10 +108,33 @@ class matrix {
     return scalar_reference(index, &backend_);
   }
 
+  std::size_t erase(key_type index) {
+    std::size_t num_deleted = backend_.hasElement(index[0], index[1]);
+
+    if (num_deleted) {
+      backend_.removeElement(index[0], index[1]);
+    }
+    return num_deleted;
+  }
+
   void clear() {
     backend_.clear();
   }
 
+  // Missing
+
+  // at()
+
+  // Iteration
+  // begin()
+  // end()
+  // cbegin()
+  // cend()
+
+  // find()
+
+  // insert variants need to return iterators
+
   matrix() = default;
   matrix(const Allocator& allocator) {}
 

From 970e8526f906eec599254b3e5d3ad0c2d38aa741 Mon Sep 17 00:00:00 2001
From: Benjamin Brock <benjamin.brock@intel.com>
Date: Tue, 6 Jun 2023 13:15:44 -0700
Subject: [PATCH 3/3] Implement `ewise_intersection`

---
 src/examples/ewise.cpp                        |  20 +
 .../spec/grb/algorithms/algorithms.hpp        |   4 +
 .../grb/algorithms/ewise_intersection.hpp     |  30 ++
 .../spec/grb/algorithms/multiply.hpp          |  27 ++
 src/interfaces/spec/grb/detail/detail.hpp     |   3 +-
 .../spec/grb/detail/iterator_adaptor.hpp      | 183 +++++++++
 .../spec/grb/functional/functional.hpp        |   3 +
 .../spec/grb/functional/monoid_traits.hpp     |  75 ++++
 .../spec/grb/functional/op_definitions.hpp    | 357 ++++++++++++++++++
 src/interfaces/spec/grb/grb.hpp               |   4 +-
 src/interfaces/spec/grb/matrix.hpp            |   9 +
 src/interfaces/spec/grb/util/util.hpp         |   2 +-
 .../spec/grb/views/full_matrix_view.hpp       | 155 ++++++++
 src/interfaces/spec/grb/views/views.hpp       |   3 +
 14 files changed, 872 insertions(+), 3 deletions(-)
 create mode 100644 src/examples/ewise.cpp
 create mode 100644 src/interfaces/spec/grb/algorithms/algorithms.hpp
 create mode 100644 src/interfaces/spec/grb/algorithms/ewise_intersection.hpp
 create mode 100644 src/interfaces/spec/grb/algorithms/multiply.hpp
 create mode 100644 src/interfaces/spec/grb/detail/iterator_adaptor.hpp
 create mode 100644 src/interfaces/spec/grb/functional/functional.hpp
 create mode 100644 src/interfaces/spec/grb/functional/monoid_traits.hpp
 create mode 100644 src/interfaces/spec/grb/functional/op_definitions.hpp
 create mode 100644 src/interfaces/spec/grb/views/full_matrix_view.hpp
 create mode 100644 src/interfaces/spec/grb/views/views.hpp

diff --git a/src/examples/ewise.cpp b/src/examples/ewise.cpp
new file mode 100644
index 0000000..40e30f1
--- /dev/null
+++ b/src/examples/ewise.cpp
@@ -0,0 +1,20 @@
+#include <grb/grb.hpp>
+
+int main(int argc, char** argv) {
+  spec::matrix<int> a({10, 10});
+  spec::matrix<int> b({10, 10});
+  spec::matrix<int> c({10, 10});
+
+  a[{2, 3}] = 12;
+  b[{2, 3}] = 12;
+
+  a[{1, 8}] = 7;
+  b[{1, 8}] = 4;
+
+  a[{7, 3}] = 2;
+  b[{4, 3}] = 2;
+
+  spec::ewise_intersection(c, a, b, spec::plus{});
+
+  return 0;
+}
diff --git a/src/interfaces/spec/grb/algorithms/algorithms.hpp b/src/interfaces/spec/grb/algorithms/algorithms.hpp
new file mode 100644
index 0000000..fa83c95
--- /dev/null
+++ b/src/interfaces/spec/grb/algorithms/algorithms.hpp
@@ -0,0 +1,4 @@
+#pragma once
+
+// #include "multiply.hpp"
+#include "ewise_intersection.hpp"
diff --git a/src/interfaces/spec/grb/algorithms/ewise_intersection.hpp b/src/interfaces/spec/grb/algorithms/ewise_intersection.hpp
new file mode 100644
index 0000000..d26130f
--- /dev/null
+++ b/src/interfaces/spec/grb/algorithms/ewise_intersection.hpp
@@ -0,0 +1,30 @@
+#pragma once
+
+#include "../detail/detail.hpp"
+#include "../util/util.hpp"
+#include "../functional/functional.hpp"
+#include "../matrix.hpp"
+#include "../views/views.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+// NOTE: concepts are missing because `GRB_SPEC_NAMESPACE::matrix` does not
+//       satisfy iteration yet.
+template <typename A,
+          typename B,
+          typename Combine,
+          typename C,
+          typename M = GRB_SPEC_NAMESPACE::full_matrix_mask<>,
+          typename Accumulate = GRB_SPEC_NAMESPACE::take_right<>
+          >
+void ewise_intersection(C&& c, A&& a, B&& b, Combine&& combine, M&& mask = M{}, Accumulate&& acc = Accumulate{}, bool merge = false)
+{
+  auto merge_enum = (merge) ? GBTL_NAMESPACE::OutputControlEnum::MERGE : GBTL_NAMESPACE::OutputControlEnum::REPLACE;
+  if constexpr(std::is_same_v<M, GRB_SPEC_NAMESPACE::full_matrix_mask<>>) {
+    GBTL_NAMESPACE::eWiseAdd(c.backend_, GBTL_NAMESPACE::NoMask(), acc, combine, a.backend_, b.backend_, merge_enum);
+  } else {
+    GBTL_NAMESPACE::eWiseAdd(c.backend_, mask.backend_, acc, combine, a.backend_, b.backend_, merge_enum);
+  }
+}
+
+} // GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/algorithms/multiply.hpp b/src/interfaces/spec/grb/algorithms/multiply.hpp
new file mode 100644
index 0000000..6550b60
--- /dev/null
+++ b/src/interfaces/spec/grb/algorithms/multiply.hpp
@@ -0,0 +1,27 @@
+#pragma once
+
+#include "../detail/detail.hpp"
+#include "../util/util.hpp"
+#include "../matrix.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <MatrixRange A,
+          MatrixRange B,
+          BinaryOperator<GRB_SPEC_NAMESPACE::matrix_scalar_t<A>, GRB_SPEC_NAMESPACE::matrix_scalar_t<B>> Combine = GRB_SPEC_NAMESPACE::multiplies<>,
+          BinaryOperator<GRB_SPEC_NAMESPACE::elementwise_return_type_t<A, B, Combine>,
+                         GRB_SPEC_NAMESPACE::elementwise_return_type_t<A, B, Combine>,
+                         GRB_SPEC_NAMESPACE::elementwise_return_type_t<A, B, Combine>> Reduce = GRB_SPEC_NAMESPACE::plus<>,
+          MaskMatrixRange M = GRB_SPEC_NAMESPACE::full_matrix_mask<>>
+auto multiply(A&& a,
+              B&& b,
+              Reduce&& reduce = Reduce{},
+              Combine&& combine = Combine{},
+              M&& mask = GRB_SPEC_NAMESPACE::full_matrix_mask())
+{
+  using T = GRB_SPEC_NAMESPACE::elementwise_return_type_t<A, B, Combine>;
+  matrix<T> c(a.shape()[0], b.shape()[1]);
+  multiply(c, a, b, reduce, combine, mask);
+}
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/detail/detail.hpp b/src/interfaces/spec/grb/detail/detail.hpp
index f17bf44..587cad5 100644
--- a/src/interfaces/spec/grb/detail/detail.hpp
+++ b/src/interfaces/spec/grb/detail/detail.hpp
@@ -2,4 +2,5 @@
 #pragma once
 
 #include "namespace_macros.hpp"
-#include "concepts.hpp"
\ No newline at end of file
+#include "concepts.hpp"
+#include "iterator_adaptor.hpp"
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/detail/iterator_adaptor.hpp b/src/interfaces/spec/grb/detail/iterator_adaptor.hpp
new file mode 100644
index 0000000..4a3a50b
--- /dev/null
+++ b/src/interfaces/spec/grb/detail/iterator_adaptor.hpp
@@ -0,0 +1,183 @@
+#pragma once
+
+#include <iterator>
+#include <type_traits>
+
+namespace GRB_SPEC_NAMESPACE {
+
+namespace detail {
+
+template <typename Accessor> class iterator_adaptor {
+public:
+  using accessor_type = Accessor;
+  using const_accessor_type = typename Accessor::const_iterator_accessor;
+  using nonconst_accessor_type = typename Accessor::nonconst_iterator_accessor;
+
+  using difference_type = typename Accessor::difference_type;
+  using value_type = typename Accessor::value_type;
+  using iterator = iterator_adaptor<accessor_type>;
+  using const_iterator = iterator_adaptor<const_accessor_type>;
+  using reference = typename Accessor::reference;
+  using iterator_category = typename Accessor::iterator_category;
+
+  using nonconst_iterator = iterator_adaptor<nonconst_accessor_type>;
+
+  static_assert(std::is_same_v<iterator, iterator_adaptor<Accessor>>);
+
+  iterator_adaptor() = default;
+  ~iterator_adaptor() = default;
+  iterator_adaptor(const iterator_adaptor &) = default;
+  iterator_adaptor &operator=(const iterator_adaptor &) = default;
+
+  template <typename... Args>
+    requires(
+        sizeof...(Args) >= 1 &&
+        !((sizeof...(Args) == 1 &&
+           (std::is_same_v<nonconst_iterator, std::decay_t<Args>> || ...)) ||
+          (std::is_same_v<const_iterator, std::decay_t<Args>> || ...) ||
+          (std::is_same_v<nonconst_accessor_type, std::decay_t<Args>> || ...) ||
+          (std::is_same_v<const_accessor_type, std::decay_t<Args>> || ...)) &&
+        std::is_constructible_v<accessor_type, Args...>)
+  iterator_adaptor(Args &&...args) : accessor_(std::forward<Args>(args)...) {}
+
+  iterator_adaptor(const accessor_type &accessor) : accessor_(accessor) {}
+  iterator_adaptor(const const_accessor_type &accessor)
+    requires(!std::is_same_v<accessor_type, const_accessor_type>)
+      : accessor_(accessor) {}
+
+  operator const_iterator() const
+    requires(!std::is_same_v<iterator, const_iterator>)
+  {
+    return const_iterator(accessor_);
+  }
+
+  bool operator==(const_iterator other) const {
+    return accessor_ == other.accessor_;
+  }
+
+  bool operator!=(const_iterator other) const { return !(*this == other); }
+
+  bool operator<(const_iterator other) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    return accessor_ < other.accessor_;
+  }
+
+  bool operator<=(const_iterator other) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    return *this < other || *this == other;
+  }
+
+  bool operator>(const_iterator other) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    return !(*this <= other);
+  }
+
+  bool operator>=(const_iterator other) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    return !(*this < other);
+  }
+
+  reference operator*() const { return *accessor_; }
+
+  reference operator[](difference_type offset) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    return *(*this + offset);
+  }
+
+  iterator &operator+=(difference_type offset) noexcept
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    accessor_ += offset;
+    return *this;
+  }
+
+  iterator &operator-=(difference_type offset) noexcept
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    accessor_ += -offset;
+    return *this;
+  }
+
+  iterator operator+(difference_type offset) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    iterator other = *this;
+    other += offset;
+    return other;
+  }
+
+  iterator operator-(difference_type offset) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    iterator other = *this;
+    other += -offset;
+    return other;
+  }
+
+  difference_type operator-(const_iterator other) const
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    return accessor_ - other.accessor_;
+  }
+
+  iterator &operator++() noexcept
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    *this += 1;
+    return *this;
+  }
+
+  iterator &operator++() noexcept
+    requires(
+        !std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    ++accessor_;
+    return *this;
+  }
+
+  iterator operator++(int) noexcept {
+    iterator other = *this;
+    ++(*this);
+    return other;
+  }
+
+  iterator &operator--() noexcept
+    requires(
+        std::is_same_v<iterator_category, std::random_access_iterator_tag> ||
+        std::is_same_v<iterator_category, std::bidirectional_iterator_tag>)
+  {
+    *this += -1;
+    return *this;
+  }
+
+  iterator operator--(int) noexcept
+    requires(
+        std::is_same_v<iterator_category, std::random_access_iterator_tag> ||
+        std::is_same_v<iterator_category, std::bidirectional_iterator_tag>)
+  {
+    iterator other = *this;
+    --(*this);
+    return other;
+  }
+
+  friend iterator operator+(difference_type n, iterator iter)
+    requires(std::is_same_v<iterator_category, std::random_access_iterator_tag>)
+  {
+    return iter + n;
+  }
+
+private:
+  friend const_iterator;
+  friend nonconst_iterator;
+
+  accessor_type accessor_;
+};
+
+} // end detail
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/functional/functional.hpp b/src/interfaces/spec/grb/functional/functional.hpp
new file mode 100644
index 0000000..ebc8a47
--- /dev/null
+++ b/src/interfaces/spec/grb/functional/functional.hpp
@@ -0,0 +1,3 @@
+#pragma once
+
+#include "op_definitions.hpp"
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/functional/monoid_traits.hpp b/src/interfaces/spec/grb/functional/monoid_traits.hpp
new file mode 100644
index 0000000..5972092
--- /dev/null
+++ b/src/interfaces/spec/grb/functional/monoid_traits.hpp
@@ -0,0 +1,75 @@
+
+#pragma once
+
+#include <type_traits>
+#include <any>
+
+namespace GRB_SPEC_NAMESPACE {
+
+using any = std::any;
+
+template <typename Fn, typename T, typename U = T, typename V = T>
+inline constexpr bool is_binary_op_v = requires(Fn fn, T t, U u) {
+                                         {fn(t, u)} -> std::convertible_to<V>;
+                                       };
+
+template <typename Fn, typename T>
+inline constexpr bool has_identity_template_v = requires { {Fn:: template identity<T>()} -> std::same_as<T>; };
+
+template <typename Fn, typename T>
+inline constexpr bool has_identity_method_v = requires { {Fn::identity()} -> std::same_as<T>; };
+
+
+
+template <typename Fn, typename T>
+requires(
+         is_binary_op_v<Fn, T, T, T> &&
+         (has_identity_method_v<Fn, T> || has_identity_template_v<Fn, T>)
+         )
+class monoid_traits
+{
+public:
+  static constexpr T identity() noexcept {
+    if constexpr(has_identity_method_v<Fn, T>) {
+      return Fn::identity();
+    } else if constexpr(has_identity_template_v<Fn, T>) {
+      return Fn:: template identity<T>();
+    }
+  }
+};
+
+template <typename T>
+class monoid_traits<std::plus<T>, T>
+{
+public:
+  static constexpr T identity() noexcept {
+    return T(0);
+  }
+};
+
+template <typename T>
+class monoid_traits<std::plus<void>, T>
+{
+public:
+  static constexpr T identity() noexcept {
+    return T(0);
+  }
+};
+
+template <typename Fn, typename T>
+inline constexpr bool has_identity_v = requires { {GRB_SPEC_NAMESPACE::monoid_traits<Fn, T>::identity()} -> std::same_as<T>; };
+
+template <typename Fn, typename T>
+inline constexpr bool is_monoid_v = is_binary_op_v<Fn, T, T, T> &&
+                                    has_identity_v<Fn, T>;
+
+template <typename Fn, typename T, typename U = T, typename V = GRB_SPEC_NAMESPACE::any>
+concept BinaryOperator = requires(Fn fn, T t, U u) {
+                           {fn(t, u)} -> std::convertible_to<V>;
+                         };
+
+template <typename Fn, typename T>
+concept Monoid = BinaryOperator<Fn, T, T, T> &&
+                 requires { {GRB_SPEC_NAMESPACE::monoid_traits<Fn, T>::identity()} -> std::same_as<T>; };
+
+} // end GRB_SPEC_NAMESPACE
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/functional/op_definitions.hpp b/src/interfaces/spec/grb/functional/op_definitions.hpp
new file mode 100644
index 0000000..590efd3
--- /dev/null
+++ b/src/interfaces/spec/grb/functional/op_definitions.hpp
@@ -0,0 +1,357 @@
+#pragma once
+
+#include <utility>
+#include <limits>
+#include <type_traits>
+#include "monoid_traits.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename Fn, typename T = void>
+class unary_op_impl_ {
+public:
+  constexpr T operator()(const T& value) const {
+    return Fn{}(value);
+  }
+};
+
+template <typename Fn>
+class unary_op_impl_<Fn, void> {
+public:
+  template <typename T>
+  constexpr auto operator()(T&& value) const {
+    return Fn{}(std::forward<T>(value));
+  }
+};
+
+template <typename Fn, typename T = void, typename U = T, typename V = void>
+class binary_op_impl_ {
+public:
+  constexpr V operator()(const T& lhs, const U& rhs) const {
+    return Fn{}(lhs, rhs);
+  }
+
+  template <typename X>
+  static constexpr X identity()
+  requires(has_identity_template_v<Fn, X>)
+  {
+    return Fn:: template identity<X>();
+  }
+};
+
+template <typename Fn, typename T, typename U>
+class binary_op_impl_<Fn, T, U, void> {
+public:
+  constexpr auto operator()(const T& lhs, const U& rhs) const {
+    return Fn{}(lhs, rhs);
+  }
+
+  static constexpr T identity()
+  requires(std::is_same_v<T, U> &&
+           has_identity_template_v<Fn, T>)
+  {
+    return Fn:: template identity<T>();
+  }
+};
+
+template <typename Fn>
+class binary_op_impl_<Fn, void, void, void>
+{
+public:
+  template <typename T, typename U>
+  constexpr auto operator()(T&& lhs, U&& rhs) const
+  {
+    return Fn{}(lhs, rhs);
+  }
+
+  template <typename T>
+  static constexpr T identity()
+  requires(has_identity_template_v<Fn, T>)
+  {
+    return Fn:: template identity<T>();
+  }
+};
+
+class plus_impl_ {
+public:
+  template <typename T, typename U>
+  constexpr auto operator()(T&& lhs, U&& rhs) const {
+    return std::forward<T>(lhs) + std::forward<U>(rhs);
+  }
+
+  template <typename T>
+  static constexpr T identity() {
+    return T(0);
+  }
+};
+
+/// Binary Operator to perform subtraction. Uses the `-` operator.
+struct minus_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return a - b;
+  }
+};
+
+/// Binary Operator to perform multiplication. Uses the `*` operator.
+struct multiplies_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return a * b;
+  }
+
+  template <typename T>
+  static constexpr T identity() {
+    return T(1);
+  }
+};
+
+/// Binary Operator to perform division. Uses the `/` operator.
+struct divides_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return a / b;
+  }
+};
+
+// omit std::functional comparisons
+
+/// Binary Operator to perform max, returning the greater of the two values,
+/// or the first element if they are equal.  Uses the `<` operator.
+
+template <typename T, typename U>
+using larger_max_integral_t =
+  std::conditional_t<std::cmp_less(std::numeric_limits<T>::max(),
+                                   std::numeric_limits<U>::max()),
+                     U, T>;
+
+struct max_impl_ {
+  template <std::integral T, std::integral U>
+  constexpr auto operator()(const T& a, const U& b) const
+    -> std::conditional_t<
+        std::cmp_less(std::numeric_limits<T>::max(),
+                      std::numeric_limits<U>::max()),
+        U, T
+       >
+  {
+    if (std::cmp_less(a, b)) {
+      return b;
+    } else {
+      return a;
+    }
+  }
+
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const
+    -> std::conditional_t<
+        std::numeric_limits<T>::max() < std::numeric_limits<U>::max(),
+        U, T
+       >
+  requires(!(std::is_integral_v<T> && std::is_integral_v<U>))
+  {
+    if (a < b) {
+      return b;
+    } else {
+      return a;
+    }
+  }
+
+  template <typename T>
+  static constexpr T identity()
+  requires(std::numeric_limits<T>::is_specialized())
+  {
+    return std::min(std::numeric_limits<T>::lowest(), -std::numeric_limits<T>::infinity());
+  }
+};
+
+
+/// Binary Operator to perform min, returning the lesser of the two values,
+/// or the first element if they are equal.  Uses the `<` operator.
+struct min_impl_ {
+  template <std::integral T, std::integral U>
+  constexpr auto operator()(const T& a, const U& b) const
+    -> std::conditional_t<
+         std::cmp_less(std::numeric_limits<U>::lowest(),
+                       std::numeric_limits<T>::lowest()),
+         U, T
+       >
+  {
+    if (std::cmp_less(b, a)) {
+      return b;
+    } else {
+      return a;
+    }
+  }
+
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const
+    -> std::conditional_t<
+         std::numeric_limits<U>::lowest() < std::numeric_limits<T>::lowest(),
+         U, T
+       >
+  requires(!(std::is_integral_v<T> && std::is_integral_v<U>))
+  {
+    if (b < a) {
+      return b;
+    } else {
+      return a;
+    }
+  }
+
+  template <typename T>
+  static constexpr T identity()
+  requires(std::numeric_limits<T>::is_specialized())
+  {
+    return std::max(std::numeric_limits<T>::max(), std::numeric_limits<T>::infinity());
+  }
+};
+
+/// Binary Operator to perform modulus, uses the `%` operator.
+struct modulus_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return a % b;
+  }
+};
+
+/// The binary operator `grb::plus`, which forms a monoid
+/// on integral types.
+/*
+template <typename T = void, typename U = T, typename V = void>
+using plus = binary_op_impl_<plus_impl_, T, U, V>;
+*/
+
+template <typename T = void, typename U = T, typename V = void>
+struct plus : public binary_op_impl_<plus_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct minus : public binary_op_impl_<minus_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct multiplies : public binary_op_impl_<multiplies_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct times : public multiplies<T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct divides : public binary_op_impl_<divides_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct max : public binary_op_impl_<max_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct min : public binary_op_impl_<min_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct modulus : public binary_op_impl_<modulus_impl_, T, U, V> {};
+
+// Unary operators
+
+struct negate_impl_ {
+  template <typename T>
+  constexpr auto operator()(const T& a) const {
+    return -a;
+  }
+};
+
+struct logical_not_impl_ {
+  template <typename T>
+  constexpr auto operator()(const T& a) const {
+    return !a;
+  }
+};
+
+template <typename T = void>
+struct negate : public unary_op_impl_<negate_impl_, T> {};
+
+template <typename T = void>
+struct logical_not : public unary_op_impl_<logical_not_impl_, T> {};
+
+// Logical operators
+
+struct logical_and_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return a && b;
+  }
+
+  template <typename T>
+  static constexpr T identity() {
+    return T(true);
+  }
+};
+
+struct logical_or_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return a || b;
+  }
+
+  template <typename T>
+  static constexpr T identity() {
+    return T(false);
+  }
+};
+
+struct logical_xor_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return (a || b) && !(a && b);
+  }
+
+  template <typename T>
+  static constexpr T identity() {
+    return T(false);
+  }
+};
+
+struct logical_xnor_impl_ {
+  template <typename T, typename U>
+  constexpr auto operator()(const T& a, const U& b) const {
+    return !((a || b) && !(a && b));
+  }
+};
+
+template <typename T = void, typename U = T, typename V = void>
+struct logical_and : binary_op_impl_<logical_and_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct logical_or : binary_op_impl_<logical_or_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct logical_xor : binary_op_impl_<logical_xor_impl_, T, U, V> {};
+
+template <typename T = void, typename U = T, typename V = void>
+struct logical_xnor : binary_op_impl_<logical_xnor_impl_, T, U, V> {};
+
+template <typename T = void>
+struct take_left {
+  T operator()(const T& left, const T& right) {
+    return left;
+  }
+};
+
+template <>
+struct take_left<void> {
+  template <typename T, typename U>
+  T operator()(const T& left, const U& right) const {
+    return left;
+  }
+};
+
+template <typename T = void>
+struct take_right {
+  T operator()(const T& left, const T& right) const {
+    return right;
+  }
+};
+
+template <>
+struct take_right<void> {
+  template <typename T, typename U>
+  U operator()(const T& left, const U& right) const {
+    return right;
+  }
+};
+
+} // end GRB_SPEC_NAMESPACE
diff --git a/src/interfaces/spec/grb/grb.hpp b/src/interfaces/spec/grb/grb.hpp
index 7e4a17f..080f4ec 100644
--- a/src/interfaces/spec/grb/grb.hpp
+++ b/src/interfaces/spec/grb/grb.hpp
@@ -1,4 +1,6 @@
 
 #pragma once
 
-#include "matrix.hpp"
\ No newline at end of file
+#include "matrix.hpp"
+#include "functional/functional.hpp"
+#include "algorithms/algorithms.hpp"
\ No newline at end of file
diff --git a/src/interfaces/spec/grb/matrix.hpp b/src/interfaces/spec/grb/matrix.hpp
index c928bce..dc8f31e 100644
--- a/src/interfaces/spec/grb/matrix.hpp
+++ b/src/interfaces/spec/grb/matrix.hpp
@@ -143,6 +143,15 @@ class matrix {
   matrix& operator=(const matrix&) = default;
   matrix& operator=(matrix&&) = default;
 
+template <typename A,
+          typename B,
+          typename Combine,
+          typename C,
+          typename M,
+          typename Accumulate
+          >
+friend void ewise_intersection(C&& c, A&& a, B&& b, Combine&& combine, M&& mask, Accumulate&& acc, bool merge);
+
 private:
   backend_type backend_;
 };
diff --git a/src/interfaces/spec/grb/util/util.hpp b/src/interfaces/spec/grb/util/util.hpp
index 53ad3a1..46c4141 100644
--- a/src/interfaces/spec/grb/util/util.hpp
+++ b/src/interfaces/spec/grb/util/util.hpp
@@ -3,4 +3,4 @@
 
 #include "hints.hpp"
 #include "index.hpp"
-#include "matrix_entry.hpp"
\ No newline at end of file
+#include "matrix_entry.hpp"
diff --git a/src/interfaces/spec/grb/views/full_matrix_view.hpp b/src/interfaces/spec/grb/views/full_matrix_view.hpp
new file mode 100644
index 0000000..a948617
--- /dev/null
+++ b/src/interfaces/spec/grb/views/full_matrix_view.hpp
@@ -0,0 +1,155 @@
+#pragma once
+
+#include <limits>
+#include <type_traits>
+#include <concepts>
+
+#include "../detail/detail.hpp"
+#include "../util/util.hpp"
+
+namespace GRB_SPEC_NAMESPACE {
+
+template <typename T, std::integral I>
+class full_matrix_accessor {
+public:
+  using value_type = GRB_SPEC_NAMESPACE::matrix_entry<T, I>;
+  using difference_type = std::ptrdiff_t;
+  using iterator_accessor = full_matrix_accessor;
+  using const_iterator_accessor = full_matrix_accessor;
+  using nonconst_iterator_accessor = full_matrix_accessor;
+  using reference = GRB_SPEC_NAMESPACE::matrix_entry<T, I>;
+  using iterator_category = std::random_access_iterator_tag;
+
+  constexpr full_matrix_accessor() noexcept = default;
+  constexpr ~full_matrix_accessor() noexcept = default;
+  constexpr full_matrix_accessor(const full_matrix_accessor&) noexcept = default;
+  constexpr full_matrix_accessor& operator=(const full_matrix_accessor&) noexcept = default;
+
+  constexpr full_matrix_accessor(I i, I j, I n, T value) noexcept : i_(i), j_(j), n_(n), value_(value) {}
+
+  constexpr full_matrix_accessor& operator+=(difference_type offset) noexcept {
+    auto new_offset = flat_offset() + offset;
+    i_ = new_offset / n_;
+    j_ = new_offset % n_;
+    return *this;
+  }
+
+  constexpr difference_type operator-(const const_iterator_accessor& other) const noexcept {
+    return difference_type(flat_offset()) - difference_type(other.flat_offset());
+  }
+
+  constexpr bool operator==(const const_iterator_accessor& other) const noexcept {
+    return i_ == other.i_ && j_ == other.j_;
+  }
+
+  constexpr bool operator<(const const_iterator_accessor& other) const noexcept {
+    return i_ < other.i_ || (i_ == other.i_ && j_ < other.j_);
+  }
+
+  constexpr reference operator*() const noexcept {
+    return GRB_SPEC_NAMESPACE::matrix_entry<T, I>({i_, j_}, value_);
+  }
+
+private:
+
+  std::size_t flat_offset() const noexcept {
+    return std::size_t(i_)*n_ + j_;
+  }
+
+  I i_, j_;
+  I n_;
+  T value_;
+};
+
+template <typename T, std::integral I = std::size_t>
+using full_matrix_iterator = GRB_SPEC_NAMESPACE::detail::iterator_adaptor<full_matrix_accessor<T, I>>;
+
+template <typename T, typename I = std::size_t>
+class full_matrix
+{
+public: 
+  using scalar_type = T;
+
+  using index_type = I;
+
+  using value_type = GRB_SPEC_NAMESPACE::matrix_entry<T, I>;
+
+  using key_type = GRB_SPEC_NAMESPACE::index<I>;
+  using map_type = T;
+
+  using size_type = std::size_t;
+
+  using difference_type = std::ptrdiff_t;
+
+  using scalar_reference = scalar_type;
+
+  using iterator = full_matrix_iterator<T, I>;
+  using const_iterator = iterator;
+
+
+  full_matrix(GRB_SPEC_NAMESPACE::index<I> shape = {std::numeric_limits<I>::max(), std::numeric_limits<I>::max()},
+              T value = T())
+    : shape_(shape), value_(value)
+  {}
+
+  iterator begin() const noexcept {
+    return iterator(0, 0, shape()[1], value_);
+  }
+
+  iterator end() const noexcept {
+    return iterator(shape()[0], 0, shape()[1], value_);
+  }
+
+  GRB_SPEC_NAMESPACE::index<I> shape() const noexcept {
+    return shape_;
+  }
+
+  size_type size() const noexcept {
+    return shape()[0]*shape()[1];
+  }
+
+  scalar_reference operator[](GRB_SPEC_NAMESPACE::index<I> index) const noexcept {
+    return value_;
+  }
+
+  iterator find(key_type key) const noexcept {
+    if constexpr(std::is_signed_v<I>) {
+      if (key[0] < 0 || key[1] < 0) {
+        return end();
+      }
+    }
+
+    if (key[0] < shape()[0] && key[1] < shape()[1]) {
+      return iterator(key[0], key[1], shape()[1], value_);
+    } else {
+      return end();
+    }
+  }
+
+private:
+  T value_;
+  GRB_SPEC_NAMESPACE::index<I> shape_;
+};
+
+template <typename I = std::size_t>
+class full_matrix_mask : public full_matrix<bool, I>
+{
+public:
+
+  full_matrix_mask(GRB_SPEC_NAMESPACE::index<I> shape = {std::numeric_limits<I>::max(), std::numeric_limits<I>::max()})
+    : full_matrix<bool, I>(shape, true)
+  {}
+};
+
+
+template <typename I = std::size_t>
+class empty_matrix_mask : public full_matrix<bool, I>
+{
+public:
+
+  empty_matrix_mask(GRB_SPEC_NAMESPACE::index<I> shape = {std::numeric_limits<I>::max(), std::numeric_limits<I>::max()})
+    : full_matrix<bool, I>(shape, false)
+  {}
+};
+
+} // end GRB_SPEC_NAMESPACE
diff --git a/src/interfaces/spec/grb/views/views.hpp b/src/interfaces/spec/grb/views/views.hpp
new file mode 100644
index 0000000..f97ff4d
--- /dev/null
+++ b/src/interfaces/spec/grb/views/views.hpp
@@ -0,0 +1,3 @@
+#pragma once
+
+#include "full_matrix_view.hpp"
\ No newline at end of file