From fcbfc9f859bde84b48fe41685f01771ad7c322c7 Mon Sep 17 00:00:00 2001
From: Prayush Kumar <prayush.kumar@ligo.org>
Date: Mon, 21 May 2018 11:37:50 +0530
Subject: [PATCH 1/4] Add class ModalVector for spectral coefficients,
 attempting to resolve #378

---
 src/DataStructures/CMakeLists.txt             |   1 +
 src/DataStructures/ModalVector.cpp            | 132 +++++
 src/DataStructures/ModalVector.hpp            | 485 ++++++++++++++++++
 tests/Unit/DataStructures/CMakeLists.txt      |   1 +
 .../Unit/DataStructures/Test_ModalVector.cpp  | 469 +++++++++++++++++
 5 files changed, 1088 insertions(+)
 create mode 100644 src/DataStructures/ModalVector.cpp
 create mode 100644 src/DataStructures/ModalVector.hpp
 create mode 100644 tests/Unit/DataStructures/Test_ModalVector.cpp

diff --git a/src/DataStructures/CMakeLists.txt b/src/DataStructures/CMakeLists.txt
index 744a928e2e2e..9874fbff7e4a 100644
--- a/src/DataStructures/CMakeLists.txt
+++ b/src/DataStructures/CMakeLists.txt
@@ -8,6 +8,7 @@ set(LIBRARY_SOURCES
     Index.cpp
     IndexIterator.cpp
     LeviCivitaIterator.cpp
+    ModalVector.cpp
     SliceIterator.cpp
     StripeIterator.cpp
     VariablesHelpers.cpp
diff --git a/src/DataStructures/ModalVector.cpp b/src/DataStructures/ModalVector.cpp
new file mode 100644
index 000000000000..bf1e1d2be2d7
--- /dev/null
+++ b/src/DataStructures/ModalVector.cpp
@@ -0,0 +1,132 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "DataStructures/ModalVector.hpp"
+
+#include <algorithm>
+#include <pup.h>
+#include <pup_stl.h>
+
+#include "Utilities/StdHelpers.hpp"
+
+ModalVector::ModalVector(const size_t size, const double value)
+    : size_(size),
+      owned_data_(size_, value),
+      data_(owned_data_.data(), size_) {}
+
+template <class T, Requires<cpp17::is_same_v<T, double>>>
+ModalVector::ModalVector(std::initializer_list<T> list)
+    : size_(list.size()),
+      owned_data_(std::move(list)),
+      data_(owned_data_.data(), size_) {}
+
+ModalVector::ModalVector(double* start, size_t size)
+    : size_(size), owned_data_(0), data_(start, size_), owning_(false) {}
+
+/// \cond HIDDEN_SYMBOLS
+ModalVector::ModalVector(const ModalVector& rhs) {
+  size_ = rhs.size();
+  if (rhs.is_owning()) {
+    owned_data_ = rhs.owned_data_;
+  } else {
+    owned_data_ = InternalStorage_t(rhs.begin(), rhs.end());
+  }
+  data_ = decltype(data_){owned_data_.data(), size_};
+}
+
+ModalVector& ModalVector::operator=(const ModalVector& rhs) {
+  if (this == &rhs) {
+    return *this;
+  }
+  if (owning_) {
+    size_ = rhs.size();
+    if (rhs.is_owning()) {
+      owned_data_ = rhs.owned_data_;
+    } else {
+      owned_data_ = InternalStorage_t(rhs.begin(), rhs.end());
+    }
+    data_ = decltype(data_){owned_data_.data(), size_};
+  } else {
+    ASSERT(rhs.size() == size(), "Must copy into same size, not "
+                                     << rhs.size() << " into " << size());
+    std::copy(rhs.begin(), rhs.end(), begin());
+  }
+  return *this;
+}
+
+ModalVector::ModalVector(ModalVector&& rhs) noexcept {
+  size_ = rhs.size_;
+  owned_data_ = std::move(rhs.owned_data_);
+  // clang-tidy: move trivially copyable type, future proof in case impl
+  // changes
+  data_ = std::move(rhs.data_);  // NOLINT
+  owning_ = rhs.owning_;
+
+  rhs.owning_ = true;
+  rhs.size_ = 0;
+  rhs.data_ = decltype(rhs.data_){};
+}
+
+ModalVector&
+ModalVector::operator=(ModalVector&& rhs) noexcept {
+  if (this == &rhs) {
+    return *this;
+  }
+  if (owning_) {
+    size_ = rhs.size_;
+    owned_data_ = std::move(rhs.owned_data_);
+    // clang-tidy: move trivially copyable type, future proof in case impl
+    // changes
+    data_ = std::move(rhs.data_);  // NOLINT
+    owning_ = rhs.owning_;
+  } else {
+    ASSERT(rhs.size() == size(), "Must copy into same size, not "
+                                     << rhs.size() << " into " << size());
+    std::copy(rhs.begin(), rhs.end(), begin());
+  }
+  rhs.owning_ = true;
+  rhs.size_ = 0;
+  rhs.data_ = decltype(rhs.data_){};
+  return *this;
+}
+/// \endcond
+
+void ModalVector::pup(PUP::er& p) noexcept {  // NOLINT
+  p | size_;
+  if (p.isUnpacking()) {
+    owning_ = true;
+    p | owned_data_;
+    data_ = decltype(data_){owned_data_.data(), size_};
+  } else {
+    if (not owning_) {
+      owned_data_ =
+          InternalStorage_t(data_.data(), data_.data() + size_);  // NOLINT
+      p | owned_data_;
+      owned_data_.clear();
+    } else {
+      p | owned_data_;
+    }
+  }
+}
+
+std::ostream& operator<<(std::ostream& os, const ModalVector& d) {
+  // This function is inside the detail namespace StdHelpers.hpp
+  StdHelpers_detail::print_helper(os, d.begin(), d.end());
+  return os;
+}
+
+/// Equivalence operator for ModalVector
+bool operator==(const ModalVector& lhs, const ModalVector& rhs) {
+  return lhs.size() == rhs.size() and
+         std::equal(lhs.begin(), lhs.end(), rhs.begin());
+}
+
+/// Inequivalence operator for ModalVector
+bool operator!=(const ModalVector& lhs, const ModalVector& rhs) {
+  return not(lhs == rhs);
+}
+
+/// \cond
+template
+ModalVector::ModalVector(std::initializer_list<double> list);
+/// \endcond
diff --git a/src/DataStructures/ModalVector.hpp b/src/DataStructures/ModalVector.hpp
new file mode 100644
index 000000000000..f0c3c54eecac
--- /dev/null
+++ b/src/DataStructures/ModalVector.hpp
@@ -0,0 +1,485 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+/// \file
+/// Defines class Data.
+
+#pragma once
+
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <initializer_list>
+#include <limits>
+#include <ostream>
+#include <type_traits>
+#include <vector>
+
+#include "DataStructures/DataVector.hpp"
+#include "ErrorHandling/Assert.hpp"
+#include "Utilities/ForceInline.hpp"
+#include "Utilities/Gsl.hpp"
+//~ #include "Utilities/PointerVector.hpp" // IWYU complains on double includes
+#include "Utilities/Requires.hpp"
+
+/// \cond HIDDEN_SYMBOLS
+namespace PUP {
+class er;
+}  // namespace PUP
+
+// clang-tidy: no using declarations in header files
+//             We want the std::abs to be used
+using std::abs;  // NOLINT
+/// \endcond
+
+/*!
+ * \ingroup DataStructuresGroup
+ * \brief A class for storing spectral coefficients on a mesh.
+ *
+ * A ModalVector holds an array of spectral coefficients, and can be
+ * either owning (the array is deleted when the Data goes out of scope) or
+ * non-owning, meaning it just has a pointer to an array.
+ *
+ * Only basic mathematical operations are supported with ModalVectors. In
+ * addition to the addition, subtraction, multiplication, division, etc. there
+ * are the following element-wise operations:
+ *
+ * - abs
+ * - max
+ * - min
+ */
+class ModalVector {
+  /// \cond HIDDEN_SYMBOLS
+  static constexpr int private_asserts() {
+    static_assert(std::is_nothrow_move_constructible<ModalVector>::value,
+                  "Missing move semantics");
+    return 0;
+  }
+  /// \endcond
+ public:
+  using value_type = double;
+  using allocator_type = std::allocator<value_type>;
+  using size_type = size_t;
+  using difference_type = std::ptrdiff_t;
+
+  // The type alias ElementType is needed because blaze::IsInvertible<T> is not
+  // SFINAE friendly to an ElementType type alias
+  using ElementType = double;
+
+ private:
+  /// The type of the "pointer" used internally
+  using InternalModalVector_t = PointerVector<double>;
+  /// The type used to store the data in
+  using InternalStorage_t = std::vector<double, allocator_type>;
+
+ public:
+  /// Create with the given size and value.
+  ///
+  /// \param size number of values
+  /// \param value the value to initialize each element.
+  explicit ModalVector(
+  size_t size, double value = std::numeric_limits<double>::signaling_NaN());
+
+  /// Create a non-owning ModalVector that points to `start`
+  ModalVector(double* start, size_t size);
+
+  /// Create from an initializer list of doubles. All elements in the
+  /// `std::initializer_list` must have decimal points
+  // cppcheck-suppress syntaxError
+  template <class T, Requires<cpp17::is_same_v<T, double>> = nullptr>
+  ModalVector(std::initializer_list<T> list);
+
+  /// Empty ModalVector
+  ModalVector() = default;
+  /// \cond HIDDEN_SYMBOLS
+  ~ModalVector() = default;
+
+  ModalVector(const ModalVector& rhs);
+  ModalVector(ModalVector&& rhs) noexcept;
+  ModalVector& operator=(const ModalVector& rhs);
+  ModalVector& operator=(ModalVector&& rhs) noexcept;
+
+  // This is a converting constructor. clang-tidy complains that it's not
+  // explicit, but we want it to allow conversion.
+  // clang-tidy: mark as explicit (we want conversion to ModalVector)
+  template <typename VT, bool VF>
+  ModalVector(const blaze::Vector<VT, VF>& expression);  // NOLINT
+
+  template <typename VT, bool VF>
+  ModalVector& operator=(const blaze::Vector<VT, VF>& expression);
+  /// \endcond
+
+  /// Number of values stored
+  size_t size() const noexcept { return size_; }
+
+  // @{
+  /// Set the ModalVector to be a reference to another ModalVector
+  /// object
+  void set_data_ref(gsl::not_null<ModalVector*> rhs) noexcept {
+    set_data_ref(rhs->data(), rhs->size_);
+  }
+  void set_data_ref(double* start, size_t size) noexcept {
+    size_ = size;
+    owned_data_ = decltype(owned_data_){};
+    data_ = decltype(data_){start, size_};
+    owning_ = false;
+  }
+  // @}
+
+  /// Returns true if the class owns the data
+  bool is_owning() const noexcept { return owning_; }
+
+  // @{
+  /// Access ith element
+  double& operator[](const size_type i) noexcept {
+    ASSERT(i < size_, "i = " << i << ", size = " << size_);
+    // clang-tidy: do not use pointer arithmetic
+    return data_[i];  // NOLINT
+  }
+  const double& operator[](const size_type i) const noexcept {
+    ASSERT(i < size_, "i = " << i << ", size = " << size_);
+    // clang-tidy: do not use pointer arithmetic
+    return data_[i];  // NOLINT
+  }
+  // @}
+
+  // @{
+  /// Access to the pointer
+  double* data() noexcept { return data_.data(); }
+  const double* data() const noexcept { return data_.data(); }
+  // @}
+
+  // @{
+  /// Returns iterator to beginning of data
+  decltype(auto) begin() noexcept { return data_.begin(); }
+  decltype(auto) begin() const noexcept { return data_.begin(); }
+  // @}
+  // @{
+  /// Returns iterator to end of data
+  decltype(auto) end() noexcept { return data_.end(); }
+  decltype(auto) end() const noexcept { return data_.end(); }
+  // @}
+
+  /// Serialization for Charm++
+  // clang-tidy: google-runtime-references
+  void pup(PUP::er& p) noexcept;  // NOLINT
+
+  // @{
+  /// See the Blaze library documentation for details on these functions since
+  /// they merely forward to Blaze.
+  /// https://bitbucket.org/blaze-lib/blaze/overview
+  ModalVector& operator=(const double& rhs) noexcept {
+    data_ = rhs;
+    return *this;
+  }
+
+  ModalVector& operator+=(const ModalVector& rhs) noexcept {
+    data_ += rhs.data_;
+    return *this;
+  }
+  template <typename VT, bool VF>
+  ModalVector& operator+=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ += rhs;
+    return *this;
+  }
+  ModalVector& operator+=(const double& rhs) noexcept {
+    data_ += rhs;
+    return *this;
+  }
+
+  ModalVector& operator-=(const ModalVector& rhs) noexcept {
+    data_ -= rhs.data_;
+    return *this;
+  }
+  template <typename VT, bool VF>
+  ModalVector& operator-=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ -= rhs;
+    return *this;
+  }
+  ModalVector& operator-=(const double& rhs) noexcept {
+    data_ -= rhs;
+    return *this;
+  }
+
+  template <typename VT, bool VF>
+  ModalVector& operator*=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ *= rhs;
+    return *this;
+  }
+  ModalVector& operator*=(const double& rhs) noexcept {
+    data_ *= rhs;
+    return *this;
+  }
+  ModalVector& operator*=(const ModalVector& rhs) noexcept {
+    data_ *= rhs.data_;
+    return *this;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV object in this function
+  ModalVector& operator*=(const DataVector& rhs) noexcept {
+    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    data_ *= _rhs_local;
+    return *this;
+  }
+
+  template <typename VT, bool VF>
+  ModalVector& operator/=(const blaze::Vector<VT, VF>& rhs) noexcept {
+    data_ /= rhs;
+    return *this;
+  }
+  ModalVector& operator/=(const double& rhs) noexcept {
+    data_ /= rhs;
+    return *this;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV object in this function
+  ModalVector& operator/=(const DataVector& rhs) noexcept {
+    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    data_ /= _rhs_local;
+    return *this;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const ModalVector& lhs, const ModalVector& rhs) noexcept {
+    return lhs.data_ + rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
+    return ~lhs + rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ + ~rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ + rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
+      const double& lhs, const ModalVector& rhs) noexcept {
+    return lhs + rhs.data_;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& lhs, const ModalVector& rhs) noexcept {
+    return lhs.data_ - rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
+    return ~lhs - rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ - ~rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ - rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const double& lhs, const ModalVector& rhs) noexcept {
+    return lhs - rhs.data_;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
+      const ModalVector& rhs) noexcept {
+    return -rhs.data_;
+  };
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ * rhs;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const double& lhs, const ModalVector& rhs) noexcept {
+    return lhs * rhs.data_;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const ModalVector& rhs) noexcept {
+    return lhs.data_ * rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
+    return ~lhs * rhs.data_;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ * ~rhs;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV objects in the
+  // next 2 functions?
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const ModalVector& lhs, const DataVector& rhs) noexcept {
+    blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    return lhs * _rhs_local;
+  }
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
+      const DataVector& lhs, const ModalVector& rhs) noexcept {
+    blaze::DynamicVector<double> _lhs_local(lhs.size(), lhs.data());
+    return rhs * _lhs_local;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
+      const ModalVector& lhs, const double& rhs) noexcept {
+    return lhs.data_ / rhs;
+  }
+  template <typename VT, bool VF>
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
+      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
+    return lhs.data_ / ~rhs;
+  }
+  // FIXME PK: Can we avoid the instantiation of new DV object in this function
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
+      const ModalVector& lhs, const DataVector& rhs) noexcept {
+    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
+    return lhs / _rhs_local;
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) min(
+      const ModalVector& t) noexcept {
+    return min(t.data_);
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) max(
+      const ModalVector& t) noexcept {
+    return max(t.data_);
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) abs(
+      const ModalVector& t) noexcept {
+    return abs(t.data_);
+  }
+
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) fabs(
+      const ModalVector& t) noexcept {
+    return abs(t.data_);
+  }
+  // @}
+
+  /// If less than zero returns zero, otherwise returns one
+  SPECTRE_ALWAYS_INLINE friend decltype(auto) step_function(
+      const ModalVector& t) noexcept {
+    return step_function(t.data_);
+  }
+
+ private:
+  /// \cond HIDDEN_SYMBOLS
+  size_t size_ = 0;
+  InternalStorage_t owned_data_;
+  InternalModalVector_t data_;
+  bool owning_{true};
+  /// \endcond
+};
+
+/// Output operator for ModalVector
+std::ostream& operator<<(std::ostream& os, const ModalVector& d);
+
+/// Equivalence operator for ModalVector
+bool operator==(const ModalVector& lhs, const ModalVector& rhs);
+
+/// Inequivalence operator for ModalVector
+bool operator!=(const ModalVector& lhs, const ModalVector& rhs);
+
+template <typename T, size_t Dim>
+std::array<ModalVector, Dim> operator+(
+    const std::array<T, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) + gsl::at(rhs, i);
+  }
+  return result;
+}
+template <typename U, size_t Dim>
+std::array<ModalVector, Dim> operator+(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<U, Dim>& rhs) noexcept {
+  return rhs + lhs;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim> operator+(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) + gsl::at(rhs, i);
+  }
+  return result;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim>& operator+=(
+    std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(lhs, i) += gsl::at(rhs, i);
+  }
+  return lhs;
+}
+
+template <typename T, size_t Dim>
+std::array<ModalVector, Dim> operator-(
+    const std::array<T, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) - gsl::at(rhs, i);
+  }
+  return result;
+}
+template <typename U, size_t Dim>
+std::array<ModalVector, Dim> operator-(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<U, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) - gsl::at(rhs, i);
+  }
+  return result;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim> operator-(
+    const std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  std::array<ModalVector, Dim> result;
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(result, i) = gsl::at(lhs, i) - gsl::at(rhs, i);
+  }
+  return result;
+}
+template <size_t Dim>
+std::array<ModalVector, Dim>& operator-=(
+    std::array<ModalVector, Dim>& lhs,
+    const std::array<ModalVector, Dim>& rhs) noexcept {
+  for (size_t i = 0; i < Dim; i++) {
+    gsl::at(lhs, i) -= gsl::at(rhs, i);
+  }
+  return lhs;
+}
+
+// FIXME PK: What is the purpose of next 2 functions copied from DataVector?
+/// \cond HIDDEN_SYMBOLS
+template <typename VT, bool VF>
+ModalVector::ModalVector(const blaze::Vector<VT, VF>& expression)
+    : size_((~expression).size()),
+      owned_data_((~expression).size()),
+      data_(owned_data_.data(), (~expression).size()) {
+  data_ = expression;
+}
+
+template <typename VT, bool VF>
+ModalVector& ModalVector::operator=(
+    const blaze::Vector<VT, VF>& expression) {
+  if (owning_ and (~expression).size() != size()) {
+    size_ = (~expression).size();
+    owned_data_ = InternalStorage_t(size_);
+    data_ = decltype(data_){owned_data_.data(), size_};
+  } else if (not owning_) {
+    ASSERT((~expression).size() == size(), "Must copy into same size");
+  }
+  data_ = expression;
+  return *this;
+}
+/// \endcond
diff --git a/tests/Unit/DataStructures/CMakeLists.txt b/tests/Unit/DataStructures/CMakeLists.txt
index 85c05e6cbf84..020f1bbbd6f3 100644
--- a/tests/Unit/DataStructures/CMakeLists.txt
+++ b/tests/Unit/DataStructures/CMakeLists.txt
@@ -12,6 +12,7 @@ set(LIBRARY_SOURCES
   Test_Index.cpp
   Test_IndexIterator.cpp
   Test_LeviCivitaIterator.cpp
+  Test_ModalVector.cpp
   Test_OrientVariablesOnSlice.cpp
   Test_SliceIterator.cpp
   Test_StripeIterator.cpp
diff --git a/tests/Unit/DataStructures/Test_ModalVector.cpp b/tests/Unit/DataStructures/Test_ModalVector.cpp
new file mode 100644
index 000000000000..ef09225c9628
--- /dev/null
+++ b/tests/Unit/DataStructures/Test_ModalVector.cpp
@@ -0,0 +1,469 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "tests/Unit/TestingFramework.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <numeric>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/ModalVector.hpp"
+#include "ErrorHandling/Error.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/StdArrayHelpers.hpp"
+#include "Utilities/StdHelpers.hpp"  // IWYU pragma: keep
+#include "tests/Unit/TestHelpers.hpp"
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector",
+                  "[DataStructures][Unit]") {
+  ModalVector a{2};
+  CHECK(a.size() == 2);
+  ModalVector b{2, 10.0};
+  CHECK(b.size() == 2);
+  for (size_t i = 0; i < b.size(); ++i) {
+    INFO(i);
+    CHECK(b[i] == 10.0);
+  }
+
+  ModalVector t(5, 10.0);
+  CHECK(t.size() == 5);
+  for (size_t i = 0; i < t.size(); ++i) {
+    INFO(i);
+    CHECK(t[i] == 10.0);
+  }
+  for (const auto& p : t) {
+    CHECK(p == 10.0);
+  }
+  for (auto& p : t) {
+    CHECK(p == 10.0);
+  }
+  ModalVector t2{1.43, 2.83, 3.94, 7.85};
+  CHECK(t2.size() == 4);
+  CHECK(t2.is_owning());
+  CHECK(t2[0] == 1.43);
+  CHECK(t2[1] == 2.83);
+  CHECK(t2[2] == 3.94);
+  CHECK(t2[3] == 7.85);
+  test_copy_semantics(t);
+  auto t_copy = t;
+  CHECK(t_copy.is_owning());
+  test_move_semantics(std::move(t), t_copy);
+  ModalVector t_move_assignment = std::move(t_copy);
+  CHECK(t_move_assignment.is_owning());
+  ModalVector t_move_constructor = std::move(t_move_assignment);
+  CHECK(t_move_constructor.is_owning());
+}
+
+SPECTRE_TEST_CASE("Unit.Serialization.ModalVector",
+                  "[DataStructures][Unit][Serialization]") {
+  const size_t npts = 10;
+  ModalVector t(npts), tgood(npts);
+  std::iota(t.begin(), t.end(), 1.2);
+  std::iota(tgood.begin(), tgood.end(), 1.2);
+  CHECK(tgood == t);
+  CHECK(t.is_owning());
+  CHECK(tgood.is_owning());
+  const ModalVector serialized = serialize_and_deserialize(t);
+  CHECK(tgood == t);
+  CHECK(serialized == tgood);
+  CHECK(serialized.is_owning());
+  CHECK(serialized.data() != t.data());
+  CHECK(t.is_owning());
+}
+
+SPECTRE_TEST_CASE("Unit.Serialization.ModalVector_Ref",
+                  "[DataStructures][Unit][Serialization]") {
+  const size_t npts = 10;
+  ModalVector t(npts);
+  std::iota(t.begin(), t.end(), 4.3);
+  ModalVector t2;
+  t2.set_data_ref(&t);
+  CHECK(t.is_owning());
+  CHECK_FALSE(t2.is_owning());
+  CHECK(t2 == t);
+  const ModalVector serialized = serialize_and_deserialize(t);
+  CHECK(t2 == t);
+  CHECK(serialized == t2);
+  CHECK(serialized.is_owning());
+  CHECK(serialized.data() != t.data());
+  CHECK(t.is_owning());
+  const ModalVector serialized2 = serialize_and_deserialize(t2);
+  CHECK(t2 == t);
+  CHECK(serialized2 == t2);
+  CHECK(serialized2.is_owning());
+  CHECK(serialized2.data() != t2.data());
+  CHECK_FALSE(t2.is_owning());
+}
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector_Ref",
+                  "[DataStructures][Unit]") {
+  ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector t;
+  t.set_data_ref(&data);
+  CHECK(not t.is_owning());
+  CHECK(data.is_owning());
+  CHECK(t.data() == data.data());
+  CHECK(t.size() == 4);
+  CHECK(t[0] == 1.43);
+  CHECK(t[1] == 2.83);
+  CHECK(t[2] == 3.94);
+  CHECK(t[3] == 7.85);
+  test_copy_semantics(t);
+  ModalVector t_copy;
+  t_copy.set_data_ref(&t);
+  test_move_semantics(std::move(t), t_copy);
+  ModalVector t_move_assignment = std::move(t_copy);
+  CHECK(not t_move_assignment.is_owning());
+  ModalVector t_move_constructor = std::move(t_move_assignment);
+  CHECK(not t_move_constructor.is_owning());
+  {
+    ModalVector data_2{1.43, 2.83, 3.94, 7.85};
+    ModalVector data_2_ref;
+    data_2_ref.set_data_ref(&data_2);
+    ModalVector data_3{2.43, 3.83, 4.94, 8.85};
+    data_2_ref = std::move(data_3);
+    CHECK(data_2[0] == 2.43);
+    CHECK(data_2[1] == 3.83);
+    CHECK(data_2[2] == 4.94);
+    CHECK(data_2[3] == 8.85);
+// Intentionally testing self-move
+#ifdef __clang__
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wself-move"
+#endif  // defined(__clang__)
+    data_2_ref = std::move(data_2_ref);
+#ifdef __clang__
+#pragma GCC diagnostic pop
+#endif  // defined(__clang__)
+    CHECK(data_2[0] == 2.43);
+    CHECK(data_2[1] == 3.83);
+    CHECK(data_2[2] == 4.94);
+    CHECK(data_2[3] == 8.85);
+    ModalVector owned_data;
+    // clang-tidy: false positive, used after it was moved
+    owned_data = data_2_ref;  // NOLINT
+    CHECK(owned_data[0] == 2.43);
+    CHECK(owned_data[1] == 3.83);
+    CHECK(owned_data[2] == 4.94);
+    CHECK(owned_data[3] == 8.85);
+    // Test operator!=
+    CHECK_FALSE(owned_data != data_2_ref);
+  }
+}
+
+// [[OutputRegex, Must copy into same size]]
+[[noreturn]] SPECTRE_TEST_CASE(
+                "Unit.DataStructures.ModalVector.ref_diff_size",
+                "[DataStructures][Unit]") {
+  ASSERTION_TEST();
+#ifdef SPECTRE_DEBUG
+  ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector data_ref;
+  data_ref.set_data_ref(&data);
+  ModalVector data2{1.43, 2.83, 3.94};
+  data_ref = data2;
+  ERROR("Failed to trigger ASSERT in an assertion test");
+#endif
+}
+
+// [[OutputRegex, Must copy into same size]]
+[[noreturn]] SPECTRE_TEST_CASE(
+    "Unit.DataStructures.ModalVector.move_ref_diff_size",
+    "[DataStructures][Unit]") {
+  ASSERTION_TEST();
+#ifdef SPECTRE_DEBUG
+  ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector data_ref;
+  data_ref.set_data_ref(&data);
+  ModalVector data2{1.43, 2.83, 3.94};
+  data_ref = std::move(data2);
+  ERROR("Failed to trigger ASSERT in an assertion test");
+#endif
+}
+
+namespace {
+template <typename T1, typename T2>
+void check_vectors(const T1& t1, const T2& t2) {
+  CHECK_ITERABLE_APPROX(t1, t2);
+}
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathAfterMove",
+                  "[Unit][DataStructures]") {
+  ModalVector m0(10, 3.0), m1(10, 9.0);
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{};
+    b = std::move(a);
+    b = m0 + m1;
+    check_vectors(b, ModalVector(10, 12.0));
+    // clang-tidy: use after move (intentional here)
+    CHECK(a.size() == 0);  // NOLINT
+    CHECK(a.is_owning());
+    a = m0 * m1;
+    check_vectors(a, ModalVector(10, 27.0));
+    check_vectors(b, ModalVector(10, 12.0));
+  }
+
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{};
+    b = std::move(a);
+    a = m0 + m1;
+    check_vectors(b, ModalVector(10, 2.0));
+    check_vectors(a, ModalVector(10, 12.0));
+  }
+
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{std::move(a)};
+    b = m0 + m1;
+    CHECK(b.size() == 10);
+    check_vectors(b, ModalVector(10, 12.0));
+    // clang-tidy: use after move (intentional here)
+    CHECK(a.size() == 0);  // NOLINT
+    CHECK(a.is_owning());
+    a = m0 * m1;
+    check_vectors(a, ModalVector(10, 27.0));
+    check_vectors(b, ModalVector(10, 12.0));
+  }
+
+  {
+    ModalVector a(10, 2.0);
+    ModalVector b{std::move(a)};
+    a = m0 + m1;
+    check_vectors(b, ModalVector(10, 2.0));
+    check_vectors(a, ModalVector(10, 12.0));
+  }
+}
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
+                  "[Unit][DataStructures]") {
+  const size_t num_pts = 19;
+  ModalVector val{1., -2., 3., -4., 8., 12., -14.};
+  ModalVector nine(num_pts, 9.0);
+  ModalVector one(num_pts, 1.0);
+
+  // Test unary minus
+  check_vectors(-nine, ModalVector(num_pts, -9.0));
+
+  check_vectors(nine + 2.0, ModalVector(num_pts, 11.0));
+  check_vectors(2.0 + nine, ModalVector(num_pts, 11.0));
+  check_vectors(nine - 2.0, ModalVector(num_pts, 7.0));
+  check_vectors(2.0 - nine, ModalVector(num_pts, -7.0));
+  check_vectors(nine + nine, ModalVector(num_pts, 18.0));
+  check_vectors(nine + (one * nine), ModalVector(num_pts, 18.0));
+  check_vectors((one * nine) + nine, ModalVector(num_pts, 18.0));
+  check_vectors(nine - ModalVector(num_pts, 8.0), one);
+  check_vectors(nine - (one * nine), ModalVector(num_pts, 0.0));
+  check_vectors((one * nine) - nine, ModalVector(num_pts, 0.0));
+
+  //~ check_vectors(ModalVector(num_pts, -1.0 / 9.0), -one / nine);
+  //~ check_vectors(ModalVector(num_pts,-8.0 / 9.0),-(nine - one) / nine);
+  check_vectors(ModalVector(num_pts, 18.0), (one / 0.5) * nine);
+  //~ check_vectors(ModalVector(num_pts, 1.0), 9.0 / nine);
+  //~ check_vectors(ModalVector(num_pts, 1.0), (one * 9.0) / nine);
+
+  CHECK(-14 == min(val));
+  CHECK(12 == max(val));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, abs(val));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, fabs(val));
+
+  check_vectors(step_function(
+                ModalVector{-12.3, 2.0, -4.0, 0.0, 7.0, -8.0}),
+                ModalVector{0.0, 1.0, 0.0, 1.0, 1.0, 0.0});
+
+  check_vectors(ModalVector(num_pts, 81.0), nine * nine);
+  check_vectors(ModalVector(num_pts, 81.0), nine * (nine * one));
+  check_vectors(ModalVector(num_pts, 81.0), (nine * nine) * one);
+  check_vectors(ModalVector(num_pts, 81.0), 9.0 * nine);
+  check_vectors(ModalVector(num_pts, 81.0), nine * 9.0);
+  check_vectors(ModalVector(num_pts, 1.0), nine / 9.0);
+  //~ check_vectors(ModalVector(num_pts, 1.0), nine / (one * 9.0));
+  //~ check_vectors(ModalVector(num_pts, 9.0), (one * nine) / one);
+
+  ModalVector dummy(nine * nine * 1.0);
+  check_vectors(ModalVector(num_pts, 81.0), dummy);
+
+  // Test assignment
+  ModalVector test_81(num_pts, -1.0);
+  test_81 = nine * nine;
+  check_vectors(ModalVector(num_pts, 81.0), test_81);
+  CHECK(test_81.is_owning());
+  ModalVector test_81_ref(test_81.data(), test_81.size());
+  test_81_ref = 0.0;
+  test_81 = 0.0;
+  test_81_ref = nine * nine;
+  check_vectors(ModalVector(num_pts, 81.0), test_81);
+  CHECK(test_81.is_owning());
+  check_vectors(ModalVector(num_pts, 81.0), test_81_ref);
+  CHECK_FALSE(test_81_ref.is_owning());
+  ModalVector second_81(num_pts);
+  second_81 = test_81;
+  check_vectors(ModalVector(num_pts, 81.0), second_81);
+  CHECK(second_81.is_owning());
+  test_81_ref = 0.0;
+  check_vectors(ModalVector(num_pts, 0.0), test_81_ref);
+  test_81_ref = second_81;
+  check_vectors(ModalVector(num_pts, 81.0), test_81_ref);
+  second_81 = 0.0;
+  test_81_ref.set_data_ref(&test_81);
+  second_81 = std::move(test_81_ref);
+  check_vectors(ModalVector(num_pts, 81.0), second_81);
+  CHECK_FALSE(second_81.is_owning());
+  second_81 = 0.0;
+  check_vectors(ModalVector(num_pts, 0.0), second_81);
+  test_81 = 81.0;
+  check_vectors(ModalVector(num_pts, 81.0), second_81);
+  CHECK_FALSE(second_81.is_owning());
+
+  test_81 = 81.0;
+  ModalVector test_081;
+  test_081.set_data_ref(&test_81);
+  check_vectors(ModalVector(num_pts, 81.0), test_081);
+  CHECK_FALSE(test_081.is_owning());
+
+  ModalVector test_assignment(num_pts, 7.0);
+  test_assignment += nine;
+  check_vectors(ModalVector(num_pts, 16.0), test_assignment);
+  test_assignment += 3.0;
+  check_vectors(ModalVector(num_pts, 19.0), test_assignment);
+  test_assignment += (nine * nine);
+  check_vectors(ModalVector(num_pts, 100.0), test_assignment);
+
+  test_assignment = 7.0;
+  test_assignment -= nine;
+  check_vectors(ModalVector(num_pts, -2.0), test_assignment);
+  test_assignment -= 3.0;
+  check_vectors(ModalVector(num_pts, -5.0), test_assignment);
+  test_assignment -= nine * nine;
+  check_vectors(ModalVector(num_pts, -86.0), test_assignment);
+
+  test_assignment = 2.0;
+  test_assignment *= 3.0;
+  check_vectors(ModalVector(num_pts, 6.0), test_assignment);
+  test_assignment *= nine;
+  check_vectors(ModalVector(num_pts, 54.0), test_assignment);
+  test_assignment = 1.0;
+  test_assignment *= nine * nine;
+  check_vectors(ModalVector(num_pts, 81.0), test_assignment);
+
+  test_assignment = 2.0;
+  test_assignment /= 2.0;
+  check_vectors(ModalVector(num_pts, 1.0), test_assignment);
+  //~ test_assignment /= ModalVector(num_pts, 0.5);
+  //~ check_vectors(ModalVector(num_pts, 2.0), test_assignment);
+  //~ test_assignment /=
+  //~       (ModalVector(num_pts, 2.0) * ModalVector(num_pts, 3.0));
+  //~ check_vectors(ModalVector(num_pts, 1.0 / 3.0), test_assignment);
+
+  // Test assignment where the RHS is an expression that contains the LHS
+  ModalVector x(num_pts, 4.);
+  x *= x;
+  check_vectors(ModalVector(num_pts, 16.0), x);
+  //~ x /= x;
+  //~ check_vectors(ModalVector(num_pts, 1.0), x);
+
+  // Test addition of arrays of ModalVectors to arrays of doubles.
+  const ModalVector t1{0.0, 1.0, 2.0, 3.0};
+  const ModalVector t2{-0.1, -0.2, -0.3, -0.4};
+  const ModalVector t3{5.0, 4.0, 3.0, 2.0};
+  const ModalVector e1{10.0, 11.0, 12.0, 13.0};
+  const ModalVector e2{19.9, 19.8, 19.7, 19.6};
+  const ModalVector e3{35.0, 34.0, 33.0, 32.0};
+  const ModalVector e4{-10.0, -9.0, -8.0, -7.0};
+  const ModalVector e5{-20.1, -20.2, -20.3, -20.4};
+  const ModalVector e6{-25.0, -26.0, -27.0, -28.0};
+  const ModalVector e7{10.0, 12.0, 14.0, 16.0};
+  const ModalVector e8{19.8, 19.6, 19.4, 19.2};
+  const ModalVector e9{40.0, 38.0, 36.0, 34.0};
+  const ModalVector e10{-10.0, -10.0, -10.0, -10.0};
+  const ModalVector e11{-20.0, -20.0, -20.0, -20.0};
+  const ModalVector e12{-30.0, -30.0, -30.0, -30.0};
+
+  const std::array<double, 3> point{{10.0, 20.0, 30.0}};
+  std::array<ModalVector, 3> points{{t1, t2, t3}};
+  const std::array<ModalVector, 3> expected{{e1, e2, e3}};
+  const std::array<ModalVector, 3> expected2{{e4, e5, e6}};
+  std::array<ModalVector, 3> dvectors1{{t1, t2, t3}};
+  const std::array<ModalVector, 3> dvectors2{{e1, e2, e3}};
+  const std::array<ModalVector, 3> expected3{{e7, e8, e9}};
+  const std::array<ModalVector, 3> expected4{{e10, e11, e12}};
+  CHECK(points + point == expected);
+  CHECK(point + points == expected);
+  CHECK(points - point == expected2);
+  CHECK(point - points == -expected2);
+
+  points += point;
+  CHECK(points == expected);
+  points -= point;
+  points -= point;
+  CHECK(points == expected2);
+
+  CHECK_ITERABLE_APPROX(dvectors1 + dvectors2, expected3);
+  CHECK_ITERABLE_APPROX(dvectors1 - dvectors2, expected4);
+
+  dvectors1 += dvectors2;
+  CHECK_ITERABLE_APPROX(dvectors1, expected3);
+  dvectors1 -= dvectors2;
+  dvectors1 -= dvectors2;
+  CHECK_ITERABLE_APPROX(dvectors1, expected4);
+
+  // Test calculation of magnitude of ModalVector
+  const std::array<ModalVector, 1> d1{{ModalVector{-2.5, 3.4}}};
+  const ModalVector expected_d1{2.5, 3.4};
+  const auto magnitude_d1 = magnitude(d1);
+  CHECK_ITERABLE_APPROX(expected_d1, magnitude_d1);
+  const std::array<ModalVector, 2> d2{{ModalVector(2, 3.),
+                                             ModalVector(2, 4.)}};
+  const ModalVector expected_d2(2, 5.);
+  const auto magnitude_d2 = magnitude(d2);
+  CHECK_ITERABLE_APPROX(expected_d2, magnitude_d2);
+  const std::array<ModalVector, 3> d3{
+      {ModalVector(2, 3.), ModalVector(2, -4.),
+       ModalVector(2, 12.)}};
+  const ModalVector expected_d3(2, 13.);
+  const auto magnitude_d3 = magnitude(d3);
+  CHECK_ITERABLE_APPROX(expected_d3, magnitude_d3);
+}
+
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathWithDataVector",
+                  "[Unit][DataStructures]") {
+  const size_t num_pts = 23;
+  ModalVector nine(num_pts, 9.0);
+  CHECK(nine.is_owning());
+  DataVector eight(num_pts, 8.0);
+
+  // Test math with DataVector
+  ModalVector test_72(num_pts, -1.0);
+  test_72 = nine * eight;
+  check_vectors(ModalVector(num_pts, 72.0), test_72);
+  CHECK(test_72.is_owning());
+  test_72 = eight * nine;
+  check_vectors(ModalVector(num_pts, 72.0), test_72);
+  CHECK(test_72.is_owning());
+  test_72 = nine;
+  test_72 *= eight;
+  check_vectors(ModalVector(num_pts, 72.0), test_72);
+  CHECK(test_72.is_owning());
+  test_72 /= eight;
+  check_vectors(ModalVector(num_pts, 9.0), test_72);
+  CHECK(test_72.is_owning());
+}
+
+// [[OutputRegex, Must copy into same size]]
+[[noreturn]] SPECTRE_TEST_CASE(
+    "Unit.DataStructures.ModalVector.ExpressionAssignError",
+    "[Unit][DataStructures]") {
+  ASSERTION_TEST();
+#ifdef SPECTRE_DEBUG
+  ModalVector one(10, 1.0);
+  ModalVector one_ref(one.data(), one.size());
+  ModalVector one_b(2, 1.0);
+  one_ref = (one_b * one_b);
+  ERROR("Failed to trigger ASSERT in an assertion test");
+#endif
+}

From f14f6f06763635ece6f061299a8f742ed24fd7df Mon Sep 17 00:00:00 2001
From: Prayush Kumar <prayush.kumar@ligo.org>
Date: Mon, 4 Jun 2018 23:33:02 +0530
Subject: [PATCH 2/4] addressing review comments re: documentation, formatting;
 removed unused tests and step_function

---
 src/DataStructures/ModalVector.hpp            | 31 +++++++++++--------
 .../Unit/DataStructures/Test_ModalVector.cpp  | 19 ++----------
 2 files changed, 20 insertions(+), 30 deletions(-)

diff --git a/src/DataStructures/ModalVector.hpp b/src/DataStructures/ModalVector.hpp
index f0c3c54eecac..443d069fa091 100644
--- a/src/DataStructures/ModalVector.hpp
+++ b/src/DataStructures/ModalVector.hpp
@@ -2,7 +2,7 @@
 // See LICENSE.txt for details.
 
 /// \file
-/// Defines class Data.
+/// Defines class ModalVector.
 
 #pragma once
 
@@ -19,7 +19,7 @@
 #include "ErrorHandling/Assert.hpp"
 #include "Utilities/ForceInline.hpp"
 #include "Utilities/Gsl.hpp"
-//~ #include "Utilities/PointerVector.hpp" // IWYU complains on double includes
+#include "Utilities/PointerVector.hpp" // IWYU pragma: keep
 #include "Utilities/Requires.hpp"
 
 /// \cond HIDDEN_SYMBOLS
@@ -37,16 +37,26 @@ using std::abs;  // NOLINT
  * \brief A class for storing spectral coefficients on a mesh.
  *
  * A ModalVector holds an array of spectral coefficients, and can be
- * either owning (the array is deleted when the Data goes out of scope) or
- * non-owning, meaning it just has a pointer to an array.
+ * either owning (the array is deleted when the ModalVector goes out of scope)
+ * or non-owning, meaning it just has a pointer to an array.
  *
  * Only basic mathematical operations are supported with ModalVectors. In
- * addition to the addition, subtraction, multiplication, division, etc. there
+ * addition to addition, subtraction, multiplication, division, there
  * are the following element-wise operations:
  *
  * - abs
  * - max
  * - min
+ *
+ * In order to allow filtering, multiplication (*, *=) and division (/, /=)
+ * operations with a DataVectors (holding filters) is supported.
+ *
+ * Note : Because of no tagging feature yet, expressions containing purely
+ * DataVectors will evaluate and get copy-constructed to a ModalVector.
+ * Disallowing this would also mean disallowing expressions such as
+ * ModalVector coeffs_filtered(coeffs * filter) - where:
+ *  - coeffs is a ModalVector
+ *  - filter is a DataVector
  */
 class ModalVector {
   /// \cond HIDDEN_SYMBOLS
@@ -77,8 +87,9 @@ class ModalVector {
   ///
   /// \param size number of values
   /// \param value the value to initialize each element.
-  explicit ModalVector(
-  size_t size, double value = std::numeric_limits<double>::signaling_NaN());
+  explicit
+  ModalVector(size_t size,
+              double value = std::numeric_limits<double>::signaling_NaN());
 
   /// Create a non-owning ModalVector that points to `start`
   ModalVector(double* start, size_t size);
@@ -359,11 +370,6 @@ class ModalVector {
   }
   // @}
 
-  /// If less than zero returns zero, otherwise returns one
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) step_function(
-      const ModalVector& t) noexcept {
-    return step_function(t.data_);
-  }
 
  private:
   /// \cond HIDDEN_SYMBOLS
@@ -459,7 +465,6 @@ std::array<ModalVector, Dim>& operator-=(
   return lhs;
 }
 
-// FIXME PK: What is the purpose of next 2 functions copied from DataVector?
 /// \cond HIDDEN_SYMBOLS
 template <typename VT, bool VF>
 ModalVector::ModalVector(const blaze::Vector<VT, VF>& expression)
diff --git a/tests/Unit/DataStructures/Test_ModalVector.cpp b/tests/Unit/DataStructures/Test_ModalVector.cpp
index ef09225c9628..6ef417964720 100644
--- a/tests/Unit/DataStructures/Test_ModalVector.cpp
+++ b/tests/Unit/DataStructures/Test_ModalVector.cpp
@@ -261,29 +261,19 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
   check_vectors(nine - (one * nine), ModalVector(num_pts, 0.0));
   check_vectors((one * nine) - nine, ModalVector(num_pts, 0.0));
 
-  //~ check_vectors(ModalVector(num_pts, -1.0 / 9.0), -one / nine);
-  //~ check_vectors(ModalVector(num_pts,-8.0 / 9.0),-(nine - one) / nine);
   check_vectors(ModalVector(num_pts, 18.0), (one / 0.5) * nine);
-  //~ check_vectors(ModalVector(num_pts, 1.0), 9.0 / nine);
-  //~ check_vectors(ModalVector(num_pts, 1.0), (one * 9.0) / nine);
 
   CHECK(-14 == min(val));
   CHECK(12 == max(val));
   check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, abs(val));
   check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, fabs(val));
 
-  check_vectors(step_function(
-                ModalVector{-12.3, 2.0, -4.0, 0.0, 7.0, -8.0}),
-                ModalVector{0.0, 1.0, 0.0, 1.0, 1.0, 0.0});
-
   check_vectors(ModalVector(num_pts, 81.0), nine * nine);
   check_vectors(ModalVector(num_pts, 81.0), nine * (nine * one));
   check_vectors(ModalVector(num_pts, 81.0), (nine * nine) * one);
   check_vectors(ModalVector(num_pts, 81.0), 9.0 * nine);
   check_vectors(ModalVector(num_pts, 81.0), nine * 9.0);
   check_vectors(ModalVector(num_pts, 1.0), nine / 9.0);
-  //~ check_vectors(ModalVector(num_pts, 1.0), nine / (one * 9.0));
-  //~ check_vectors(ModalVector(num_pts, 9.0), (one * nine) / one);
 
   ModalVector dummy(nine * nine * 1.0);
   check_vectors(ModalVector(num_pts, 81.0), dummy);
@@ -354,18 +344,11 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
   test_assignment = 2.0;
   test_assignment /= 2.0;
   check_vectors(ModalVector(num_pts, 1.0), test_assignment);
-  //~ test_assignment /= ModalVector(num_pts, 0.5);
-  //~ check_vectors(ModalVector(num_pts, 2.0), test_assignment);
-  //~ test_assignment /=
-  //~       (ModalVector(num_pts, 2.0) * ModalVector(num_pts, 3.0));
-  //~ check_vectors(ModalVector(num_pts, 1.0 / 3.0), test_assignment);
 
   // Test assignment where the RHS is an expression that contains the LHS
   ModalVector x(num_pts, 4.);
   x *= x;
   check_vectors(ModalVector(num_pts, 16.0), x);
-  //~ x /= x;
-  //~ check_vectors(ModalVector(num_pts, 1.0), x);
 
   // Test addition of arrays of ModalVectors to arrays of doubles.
   const ModalVector t1{0.0, 1.0, 2.0, 3.0};
@@ -445,6 +428,8 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathWithDataVector",
   test_72 = eight * nine;
   check_vectors(ModalVector(num_pts, 72.0), test_72);
   CHECK(test_72.is_owning());
+  test_72 = test_72 / eight;
+  check_vectors(ModalVector(num_pts, 9.0), test_72);
   test_72 = nine;
   test_72 *= eight;
   check_vectors(ModalVector(num_pts, 72.0), test_72);

From 16451c6517c0a7df4e198f511d476d784cc0c54f Mon Sep 17 00:00:00 2001
From: Prayush Kumar <prayush.kumar@ligo.org>
Date: Fri, 27 Jul 2018 00:11:10 -0400
Subject: [PATCH 3/4] updating the PR to include the ResultType tagging feature
 in PointerVector

---
 src/DataStructures/ModalVector.cpp            | 121 ++--
 src/DataStructures/ModalVector.hpp            | 520 +++++++++---------
 .../Unit/DataStructures/Test_ModalVector.cpp  | 249 ++++++---
 3 files changed, 465 insertions(+), 425 deletions(-)

diff --git a/src/DataStructures/ModalVector.cpp b/src/DataStructures/ModalVector.cpp
index bf1e1d2be2d7..cd99f18f76ff 100644
--- a/src/DataStructures/ModalVector.cpp
+++ b/src/DataStructures/ModalVector.cpp
@@ -5,107 +5,91 @@
 
 #include <algorithm>
 #include <pup.h>
-#include <pup_stl.h>
 
 #include "Utilities/StdHelpers.hpp"
 
 ModalVector::ModalVector(const size_t size, const double value)
-    : size_(size),
-      owned_data_(size_, value),
-      data_(owned_data_.data(), size_) {}
+    : owned_data_(size, value) {
+  reset_pointer_vector();
+}
+
+ModalVector::ModalVector(double* start, size_t size) noexcept
+    : BaseType(start, size), owned_data_(0), owning_(false) {}
 
 template <class T, Requires<cpp17::is_same_v<T, double>>>
 ModalVector::ModalVector(std::initializer_list<T> list)
-    : size_(list.size()),
-      owned_data_(std::move(list)),
-      data_(owned_data_.data(), size_) {}
-
-ModalVector::ModalVector(double* start, size_t size)
-    : size_(size), owned_data_(0), data_(start, size_), owning_(false) {}
+    : owned_data_(std::move(list)) {
+  reset_pointer_vector();
+}
 
 /// \cond HIDDEN_SYMBOLS
-ModalVector::ModalVector(const ModalVector& rhs) {
-  size_ = rhs.size();
+// clang-tidy: calling a base constructor other than the copy constructor.
+//             We reset the base class in reset_pointer_vector after calling its
+//             default constructor
+ModalVector::ModalVector(const ModalVector& rhs) : BaseType{} {  // NOLINT
   if (rhs.is_owning()) {
     owned_data_ = rhs.owned_data_;
   } else {
-    owned_data_ = InternalStorage_t(rhs.begin(), rhs.end());
+    owned_data_.assign(rhs.begin(), rhs.end());
   }
-  data_ = decltype(data_){owned_data_.data(), size_};
+  reset_pointer_vector();
 }
 
 ModalVector& ModalVector::operator=(const ModalVector& rhs) {
-  if (this == &rhs) {
-    return *this;
-  }
-  if (owning_) {
-    size_ = rhs.size();
-    if (rhs.is_owning()) {
-      owned_data_ = rhs.owned_data_;
+  if (this != &rhs) {
+    if (owning_) {
+      if (rhs.is_owning()) {
+        owned_data_ = rhs.owned_data_;
+      } else {
+        owned_data_.assign(rhs.begin(), rhs.end());
+      }
+      reset_pointer_vector();
     } else {
-      owned_data_ = InternalStorage_t(rhs.begin(), rhs.end());
+      ASSERT(rhs.size() == size(), "Must copy into same size, not "
+                                       << rhs.size() << " into " << size());
+      std::copy(rhs.begin(), rhs.end(), begin());
     }
-    data_ = decltype(data_){owned_data_.data(), size_};
-  } else {
-    ASSERT(rhs.size() == size(), "Must copy into same size, not "
-                                     << rhs.size() << " into " << size());
-    std::copy(rhs.begin(), rhs.end(), begin());
   }
   return *this;
 }
 
 ModalVector::ModalVector(ModalVector&& rhs) noexcept {
-  size_ = rhs.size_;
   owned_data_ = std::move(rhs.owned_data_);
-  // clang-tidy: move trivially copyable type, future proof in case impl
-  // changes
-  data_ = std::move(rhs.data_);  // NOLINT
+  ~*this = ~rhs;  // PointerVector is trivially copyable
   owning_ = rhs.owning_;
 
   rhs.owning_ = true;
-  rhs.size_ = 0;
-  rhs.data_ = decltype(rhs.data_){};
+  rhs.reset();
 }
 
-ModalVector&
-ModalVector::operator=(ModalVector&& rhs) noexcept {
-  if (this == &rhs) {
-    return *this;
-  }
-  if (owning_) {
-    size_ = rhs.size_;
-    owned_data_ = std::move(rhs.owned_data_);
-    // clang-tidy: move trivially copyable type, future proof in case impl
-    // changes
-    data_ = std::move(rhs.data_);  // NOLINT
-    owning_ = rhs.owning_;
-  } else {
-    ASSERT(rhs.size() == size(), "Must copy into same size, not "
-                                     << rhs.size() << " into " << size());
-    std::copy(rhs.begin(), rhs.end(), begin());
+ModalVector& ModalVector::operator=(ModalVector&& rhs) noexcept {
+  if (this != &rhs) {
+    if (owning_) {
+      owned_data_ = std::move(rhs.owned_data_);
+      ~*this = ~rhs;  // PointerVector is trivially copyable
+      owning_ = rhs.owning_;
+    } else {
+      ASSERT(rhs.size() == size(), "Must copy into same size, not "
+                                       << rhs.size() << " into " << size());
+      std::copy(rhs.begin(), rhs.end(), begin());
+    }
+    rhs.owning_ = true;
+    rhs.reset();
   }
-  rhs.owning_ = true;
-  rhs.size_ = 0;
-  rhs.data_ = decltype(rhs.data_){};
   return *this;
 }
 /// \endcond
 
 void ModalVector::pup(PUP::er& p) noexcept {  // NOLINT
-  p | size_;
-  if (p.isUnpacking()) {
-    owning_ = true;
-    p | owned_data_;
-    data_ = decltype(data_){owned_data_.data(), size_};
-  } else {
-    if (not owning_) {
-      owned_data_ =
-          InternalStorage_t(data_.data(), data_.data() + size_);  // NOLINT
-      p | owned_data_;
-      owned_data_.clear();
-    } else {
-      p | owned_data_;
+  auto my_size = size();
+  p | my_size;
+  if (my_size > 0) {
+    if (p.isUnpacking()) {
+      owning_ = true;
+      owned_data_.resize(my_size);
+      reset_pointer_vector();
     }
+    PUParray(p, data(), size());
   }
 }
 
@@ -116,17 +100,16 @@ std::ostream& operator<<(std::ostream& os, const ModalVector& d) {
 }
 
 /// Equivalence operator for ModalVector
-bool operator==(const ModalVector& lhs, const ModalVector& rhs) {
+bool operator==(const ModalVector& lhs, const ModalVector& rhs) noexcept {
   return lhs.size() == rhs.size() and
          std::equal(lhs.begin(), lhs.end(), rhs.begin());
 }
 
 /// Inequivalence operator for ModalVector
-bool operator!=(const ModalVector& lhs, const ModalVector& rhs) {
+bool operator!=(const ModalVector& lhs, const ModalVector& rhs) noexcept {
   return not(lhs == rhs);
 }
 
 /// \cond
-template
-ModalVector::ModalVector(std::initializer_list<double> list);
+template ModalVector::ModalVector(std::initializer_list<double> list);
 /// \endcond
diff --git a/src/DataStructures/ModalVector.hpp b/src/DataStructures/ModalVector.hpp
index 443d069fa091..2099b878f2c6 100644
--- a/src/DataStructures/ModalVector.hpp
+++ b/src/DataStructures/ModalVector.hpp
@@ -15,12 +15,13 @@
 #include <type_traits>
 #include <vector>
 
-#include "DataStructures/DataVector.hpp"
 #include "ErrorHandling/Assert.hpp"
 #include "Utilities/ForceInline.hpp"
 #include "Utilities/Gsl.hpp"
+#include "Utilities/MakeWithValue.hpp"
 #include "Utilities/PointerVector.hpp" // IWYU pragma: keep
 #include "Utilities/Requires.hpp"
+#include "Utilities/TMPL.hpp"          // for list
 
 /// \cond HIDDEN_SYMBOLS
 namespace PUP {
@@ -32,6 +33,39 @@ class er;
 using std::abs;  // NOLINT
 /// \endcond
 
+// IWYU doesn't like that we want PointerVector.hpp to expose Blaze and also
+// have ModalVector.hpp to expose PointerVector.hpp without including Blaze
+// directly in ModalVector.hpp
+//
+// IWYU pragma: no_include <blaze/math/dense/DenseVector.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecAddExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecDivExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecMultExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecDVecSubExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecMapExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecScalarDivExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DVecScalarMultExpr.h>
+// IWYU pragma: no_include <blaze/math/expressions/DenseVector.h>
+// IWYU pragma: no_include <blaze/math/expressions/Vector.h>
+// IWYU pragma: no_include <blaze/math/typetraits/IsVector.h>
+// IWYU pragma: no_include <blaze/math/expressions/Forward.h>
+// IWYU pragma: no_include <blaze/math/AlignmentFlag.h>
+// IWYU pragma: no_include <blaze/math/PaddingFlag.h>
+// IWYU pragma: no_include <blaze/math/traits/AddTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/DivTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/MultTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/SubTrait.h>
+// IWYU pragma: no_include <blaze/system/TransposeFlag.h>
+// IWYU pragma: no_include <blaze/math/traits/BinaryMapTrait.h>
+// IWYU pragma: no_include <blaze/math/traits/UnaryMapTrait.h>
+// IWYU pragma: no_include <blaze/math/typetraits/TransposeFlag.h>
+
+// IWYU pragma: no_forward_declare blaze::DenseVector
+// IWYU pragma: no_forward_declare blaze::UnaryMapTrait
+// IWYU pragma: no_forward_declare blaze::BinaryMapTrait
+// IWYU pragma: no_forward_declare blaze::IsVector
+// IWYU pragma: no_forward_declare blaze::TransposeFlag
+
 /*!
  * \ingroup DataStructuresGroup
  * \brief A class for storing spectral coefficients on a mesh.
@@ -44,21 +78,17 @@ using std::abs;  // NOLINT
  * addition to addition, subtraction, multiplication, division, there
  * are the following element-wise operations:
  *
- * - abs
+ * - abs/fabs/magnitude
  * - max
  * - min
  *
  * In order to allow filtering, multiplication (*, *=) and division (/, /=)
- * operations with a DataVectors (holding filters) is supported.
+ * operations with a DenseVectors (holding filters) is supported.
  *
- * Note : Because of no tagging feature yet, expressions containing purely
- * DataVectors will evaluate and get copy-constructed to a ModalVector.
- * Disallowing this would also mean disallowing expressions such as
- * ModalVector coeffs_filtered(coeffs * filter) - where:
- *  - coeffs is a ModalVector
- *  - filter is a DataVector
  */
-class ModalVector {
+class ModalVector
+    : public PointerVector<double, blaze::unaligned, blaze::unpadded,
+                           blaze::defaultTransposeFlag, ModalVector> {
   /// \cond HIDDEN_SYMBOLS
   static constexpr int private_asserts() {
     static_assert(std::is_nothrow_move_constructible<ModalVector>::value,
@@ -71,18 +101,32 @@ class ModalVector {
   using allocator_type = std::allocator<value_type>;
   using size_type = size_t;
   using difference_type = std::ptrdiff_t;
+  using BaseType = PointerVector<value_type, blaze::unaligned, blaze::unpadded,
+                                 blaze::defaultTransposeFlag, ModalVector>;
+  static constexpr bool transpose_flag = blaze::defaultTransposeFlag;
+
+  using BaseType::ElementType;
+  using TransposeType = ModalVector;
+  using CompositeType = const ModalVector&;
+
+  /// Iterators etc obtained from base class PointerVector
+  using BaseType::operator[];
+  using BaseType::begin;
+  using BaseType::cbegin;
+  using BaseType::cend;
+  using BaseType::data;
+  using BaseType::end;
+  using BaseType::size;
+
+  /// Do we need an upcast from ModalVector to PointerVector?
+  // @{
+  // Upcast to `BaseType`
+  const BaseType& operator~() const noexcept {
+    return static_cast<const BaseType&>(*this);
+  }
+  BaseType& operator~() noexcept { return static_cast<BaseType&>(*this); }
+  // @}
 
-  // The type alias ElementType is needed because blaze::IsInvertible<T> is not
-  // SFINAE friendly to an ElementType type alias
-  using ElementType = double;
-
- private:
-  /// The type of the "pointer" used internally
-  using InternalModalVector_t = PointerVector<double>;
-  /// The type used to store the data in
-  using InternalStorage_t = std::vector<double, allocator_type>;
-
- public:
   /// Create with the given size and value.
   ///
   /// \param size number of values
@@ -92,47 +136,58 @@ class ModalVector {
               double value = std::numeric_limits<double>::signaling_NaN());
 
   /// Create a non-owning ModalVector that points to `start`
-  ModalVector(double* start, size_t size);
+  ModalVector(double* start, size_t size) noexcept;
 
-  /// Create from an initializer list of doubles. All elements in the
-  /// `std::initializer_list` must have decimal points
-  // cppcheck-suppress syntaxError
+  /// Create from an initializer list of doubles (must have decimal points)
   template <class T, Requires<cpp17::is_same_v<T, double>> = nullptr>
   ModalVector(std::initializer_list<T> list);
 
   /// Empty ModalVector
-  ModalVector() = default;
+  ModalVector() noexcept = default;
   /// \cond HIDDEN_SYMBOLS
   ~ModalVector() = default;
 
+  // Initialize ModalVector with expressions involving itself
   ModalVector(const ModalVector& rhs);
   ModalVector(ModalVector&& rhs) noexcept;
   ModalVector& operator=(const ModalVector& rhs);
   ModalVector& operator=(ModalVector&& rhs) noexcept;
 
+  /// Catch DenseVector expressions into ModalVectors, only if their return
+  /// types are defined as ModalVectors
+  //
   // This is a converting constructor. clang-tidy complains that it's not
   // explicit, but we want it to allow conversion.
   // clang-tidy: mark as explicit (we want conversion to ModalVector)
   template <typename VT, bool VF>
-  ModalVector(const blaze::Vector<VT, VF>& expression);  // NOLINT
+  ModalVector(const blaze::DenseVector<VT, VF>& expression) noexcept;  // NOLINT
 
   template <typename VT, bool VF>
-  ModalVector& operator=(const blaze::Vector<VT, VF>& expression);
+  ModalVector& operator=(const blaze::DenseVector<VT, VF>& expression) noexcept;
   /// \endcond
 
-  /// Number of values stored
-  size_t size() const noexcept { return size_; }
+  ModalVector& operator=(const double& rhs) noexcept {
+    ~*this = rhs;
+    return *this;
+  }
+
+  // Using expression macro, defined in PointerVector.hpp, these lines overload
+  // assignment operators { +=, -=, *=, /= } for LHS = ModalVector and RHS =
+  // (i) ModalVector, (ii) blaze::DenseVector (to hold spectral filters), and
+  // (iii) double (i.e. the ElementType)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(+=, ModalVector)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(-=, ModalVector)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(*=, ModalVector)
+  MAKE_EXPRESSION_MATH_ASSIGN_PV(/=, ModalVector)
 
   // @{
-  /// Set the ModalVector to be a reference to another ModalVector
-  /// object
+  /// Set ModalVector to be a reference to another ModalVector object
   void set_data_ref(gsl::not_null<ModalVector*> rhs) noexcept {
-    set_data_ref(rhs->data(), rhs->size_);
+    set_data_ref(rhs->data(), rhs->size());
   }
   void set_data_ref(double* start, size_t size) noexcept {
-    size_ = size;
     owned_data_ = decltype(owned_data_){};
-    data_ = decltype(data_){start, size_};
+    (~*this).reset(start, size);
     owning_ = false;
   }
   // @}
@@ -140,254 +195,162 @@ class ModalVector {
   /// Returns true if the class owns the data
   bool is_owning() const noexcept { return owning_; }
 
-  // @{
-  /// Access ith element
-  double& operator[](const size_type i) noexcept {
-    ASSERT(i < size_, "i = " << i << ", size = " << size_);
-    // clang-tidy: do not use pointer arithmetic
-    return data_[i];  // NOLINT
-  }
-  const double& operator[](const size_type i) const noexcept {
-    ASSERT(i < size_, "i = " << i << ", size = " << size_);
-    // clang-tidy: do not use pointer arithmetic
-    return data_[i];  // NOLINT
-  }
-  // @}
-
-  // @{
-  /// Access to the pointer
-  double* data() noexcept { return data_.data(); }
-  const double* data() const noexcept { return data_.data(); }
-  // @}
-
-  // @{
-  /// Returns iterator to beginning of data
-  decltype(auto) begin() noexcept { return data_.begin(); }
-  decltype(auto) begin() const noexcept { return data_.begin(); }
-  // @}
-  // @{
-  /// Returns iterator to end of data
-  decltype(auto) end() noexcept { return data_.end(); }
-  decltype(auto) end() const noexcept { return data_.end(); }
-  // @}
-
   /// Serialization for Charm++
   // clang-tidy: google-runtime-references
   void pup(PUP::er& p) noexcept;  // NOLINT
 
-  // @{
-  /// See the Blaze library documentation for details on these functions since
-  /// they merely forward to Blaze.
-  /// https://bitbucket.org/blaze-lib/blaze/overview
-  ModalVector& operator=(const double& rhs) noexcept {
-    data_ = rhs;
-    return *this;
+ private:
+  SPECTRE_ALWAYS_INLINE void reset_pointer_vector() noexcept {
+    reset(owned_data_.data(), owned_data_.size());
   }
 
-  ModalVector& operator+=(const ModalVector& rhs) noexcept {
-    data_ += rhs.data_;
-    return *this;
-  }
-  template <typename VT, bool VF>
-  ModalVector& operator+=(const blaze::Vector<VT, VF>& rhs) noexcept {
-    data_ += rhs;
-    return *this;
-  }
-  ModalVector& operator+=(const double& rhs) noexcept {
-    data_ += rhs;
-    return *this;
-  }
+  /// \cond HIDDEN_SYMBOLS
+  std::vector<value_type, allocator_type> owned_data_;
+  bool owning_{true};
+  /// \endcond
+};
 
-  ModalVector& operator-=(const ModalVector& rhs) noexcept {
-    data_ -= rhs.data_;
-    return *this;
-  }
-  template <typename VT, bool VF>
-  ModalVector& operator-=(const blaze::Vector<VT, VF>& rhs) noexcept {
-    data_ -= rhs;
-    return *this;
-  }
-  ModalVector& operator-=(const double& rhs) noexcept {
-    data_ -= rhs;
-    return *this;
-  }
+/// Output operator for ModalVector
+std::ostream& operator<<(std::ostream& os, const ModalVector& d);
 
-  template <typename VT, bool VF>
-  ModalVector& operator*=(const blaze::Vector<VT, VF>& rhs) noexcept {
-    data_ *= rhs;
-    return *this;
-  }
-  ModalVector& operator*=(const double& rhs) noexcept {
-    data_ *= rhs;
-    return *this;
-  }
-  ModalVector& operator*=(const ModalVector& rhs) noexcept {
-    data_ *= rhs.data_;
-    return *this;
-  }
-  // FIXME PK: Can we avoid the instantiation of new DV object in this function
-  ModalVector& operator*=(const DataVector& rhs) noexcept {
-    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
-    data_ *= _rhs_local;
-    return *this;
-  }
+/// Equivalence operator for ModalVector
+bool operator==(const ModalVector& lhs, const ModalVector& rhs) noexcept;
 
-  template <typename VT, bool VF>
-  ModalVector& operator/=(const blaze::Vector<VT, VF>& rhs) noexcept {
-    data_ /= rhs;
-    return *this;
-  }
-  ModalVector& operator/=(const double& rhs) noexcept {
-    data_ /= rhs;
-    return *this;
-  }
-  // FIXME PK: Can we avoid the instantiation of new DV object in this function
-  ModalVector& operator/=(const DataVector& rhs) noexcept {
-    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
-    data_ /= _rhs_local;
-    return *this;
-  }
+/// Inequivalence operator for ModalVector
+bool operator!=(const ModalVector& lhs, const ModalVector& rhs) noexcept;
 
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
-      const ModalVector& lhs, const ModalVector& rhs) noexcept {
-    return lhs.data_ + rhs.data_;
-  }
-  template <typename VT, bool VF>
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
-      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
-    return ~lhs + rhs.data_;
-  }
-  template <typename VT, bool VF>
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
-      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
-    return lhs.data_ + ~rhs;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
-      const ModalVector& lhs, const double& rhs) noexcept {
-    return lhs.data_ + rhs;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator+(
-      const double& lhs, const ModalVector& rhs) noexcept {
-    return lhs + rhs.data_;
-  }
+/// \cond
+// Used for comparing ModalVector to an expression
+template <typename VT, bool VF>
+bool operator==(const ModalVector& lhs,
+                const blaze::DenseVector<VT, VF>& rhs) noexcept {
+  return lhs == ModalVector(rhs);
+}
 
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
-      const ModalVector& lhs, const ModalVector& rhs) noexcept {
-    return lhs.data_ - rhs.data_;
-  }
-  template <typename VT, bool VF>
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
-      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
-    return ~lhs - rhs.data_;
-  }
-  template <typename VT, bool VF>
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
-      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
-    return lhs.data_ - ~rhs;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
-      const ModalVector& lhs, const double& rhs) noexcept {
-    return lhs.data_ - rhs;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
-      const double& lhs, const ModalVector& rhs) noexcept {
-    return lhs - rhs.data_;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator-(
-      const ModalVector& rhs) noexcept {
-    return -rhs.data_;
-  };
-
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
-      const ModalVector& lhs, const double& rhs) noexcept {
-    return lhs.data_ * rhs;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
-      const double& lhs, const ModalVector& rhs) noexcept {
-    return lhs * rhs.data_;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
-      const ModalVector& lhs, const ModalVector& rhs) noexcept {
-    return lhs.data_ * rhs.data_;
-  }
-  template <typename VT, bool VF>
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
-      const blaze::Vector<VT, VF>& lhs, const ModalVector& rhs) noexcept {
-    return ~lhs * rhs.data_;
-  }
-  template <typename VT, bool VF>
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
-      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
-    return lhs.data_ * ~rhs;
-  }
-  // FIXME PK: Can we avoid the instantiation of new DV objects in the
-  // next 2 functions?
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
-      const ModalVector& lhs, const DataVector& rhs) noexcept {
-    blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
-    return lhs * _rhs_local;
-  }
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator*(
-      const DataVector& lhs, const ModalVector& rhs) noexcept {
-    blaze::DynamicVector<double> _lhs_local(lhs.size(), lhs.data());
-    return rhs * _lhs_local;
-  }
+template <typename VT, bool VF>
+bool operator!=(const ModalVector& lhs,
+                const blaze::DenseVector<VT, VF>& rhs) noexcept {
+  return not(lhs == rhs);
+}
 
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
-      const ModalVector& lhs, const double& rhs) noexcept {
-    return lhs.data_ / rhs;
-  }
-  template <typename VT, bool VF>
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
-      const ModalVector& lhs, const blaze::Vector<VT, VF>& rhs) noexcept {
-    return lhs.data_ / ~rhs;
-  }
-  // FIXME PK: Can we avoid the instantiation of new DV object in this function
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) operator/(
-      const ModalVector& lhs, const DataVector& rhs) noexcept {
-    const blaze::DynamicVector<double> _rhs_local(rhs.size(), rhs.data());
-    return lhs / _rhs_local;
-  }
+template <typename VT, bool VF>
+bool operator==(const blaze::DenseVector<VT, VF>& lhs,
+                const ModalVector& rhs) noexcept {
+  return ModalVector(lhs) == rhs;
+}
 
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) min(
-      const ModalVector& t) noexcept {
-    return min(t.data_);
-  }
+template <typename VT, bool VF>
+bool operator!=(const blaze::DenseVector<VT, VF>& lhs,
+                const ModalVector& rhs) noexcept {
+  return not(lhs == rhs);
+}
+/// \endcond
 
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) max(
-      const ModalVector& t) noexcept {
-    return max(t.data_);
-  }
+// Specialize the Blaze type traits to correctly handle ModalVector
+namespace blaze {
+template <>
+struct IsVector<ModalVector> : std::true_type {};
 
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) abs(
-      const ModalVector& t) noexcept {
-    return abs(t.data_);
-  }
+template <>
+struct TransposeFlag<ModalVector> : BoolConstant<
+                ModalVector::transpose_flag> {};
 
-  SPECTRE_ALWAYS_INLINE friend decltype(auto) fabs(
-      const ModalVector& t) noexcept {
-    return abs(t.data_);
-  }
-  // @}
+template <>
+struct AddTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
 
+template <>
+struct AddTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
 
- private:
-  /// \cond HIDDEN_SYMBOLS
-  size_t size_ = 0;
-  InternalStorage_t owned_data_;
-  InternalModalVector_t data_;
-  bool owning_{true};
-  /// \endcond
+template <>
+struct AddTrait<double, ModalVector> {
+  using Type = ModalVector;
 };
 
-/// Output operator for ModalVector
-std::ostream& operator<<(std::ostream& os, const ModalVector& d);
+template <>
+struct SubTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
 
-/// Equivalence operator for ModalVector
-bool operator==(const ModalVector& lhs, const ModalVector& rhs);
+template <>
+struct SubTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
 
-/// Inequivalence operator for ModalVector
-bool operator!=(const ModalVector& lhs, const ModalVector& rhs);
+template <>
+struct SubTrait<double, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct MultTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct MultTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
+
+template <>
+struct MultTrait<double, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct DivTrait<ModalVector, ModalVector> {
+  using Type = ModalVector;
+};
+
+template <>
+struct DivTrait<ModalVector, double> {
+  using Type = ModalVector;
+};
+
+// Forbid math operations in this specialization of UnaryMap traits for
+// ModalVector that are unlikely to be used on spectral coefficients
+template <typename Operator>
+struct UnaryMapTrait<ModalVector, Operator> {
+  static_assert(not tmpl::list_contains_v<tmpl::list<//blaze::Max, blaze::Min,
+                //blaze::Sqrt,
+                blaze::Cbrt,
+                blaze::InvSqrt, blaze::InvCbrt,
+                blaze::Acos, blaze::Acosh, blaze::Cos, blaze::Cosh,
+                blaze::Asin, blaze::Asinh, blaze::Sin, blaze::Sinh,
+                blaze::Atan, blaze::Atan2, blaze::Atanh,
+                blaze::Tan, blaze::Tanh, blaze::Hypot,
+                blaze::Exp, blaze::Exp2, blaze::Exp10,
+                blaze::Log, blaze::Log2, blaze::Log10,
+                blaze::Erf, blaze::Erfc, blaze::StepFunction
+                >, Operator>,
+                "This operation is not permitted on a ModalVector."
+                "Only unary operation permitted are: abs, fabs, sqrt.");
+  using Type = ModalVector;
+};
+
+// Forbid math operations in this specialization of BinaryMap traits for
+// ModalVector that are unlikely to be used on spectral coefficients
+template <typename Operator>
+struct BinaryMapTrait<ModalVector, ModalVector, Operator> {
+  static_assert(not tmpl::list_contains_v<tmpl::list<blaze::Max, blaze::Min
+                >, Operator>,
+                "This operation is not permitted on a ModalVector."
+                "Only unary operation are permitted: abs, fabs, sqrt.");
+  using Type = ModalVector;
+};
+
+}  // namespace blaze
+
+SPECTRE_ALWAYS_INLINE decltype(auto) fabs(const ModalVector& t) noexcept {
+  return abs(~t);
+}
+
+SPECTRE_ALWAYS_INLINE decltype(auto) abs(const ModalVector& t) noexcept {
+  return abs(~t);
+}
 
 template <typename T, size_t Dim>
 std::array<ModalVector, Dim> operator+(
@@ -467,24 +430,39 @@ std::array<ModalVector, Dim>& operator-=(
 
 /// \cond HIDDEN_SYMBOLS
 template <typename VT, bool VF>
-ModalVector::ModalVector(const blaze::Vector<VT, VF>& expression)
-    : size_((~expression).size()),
-      owned_data_((~expression).size()),
-      data_(owned_data_.data(), (~expression).size()) {
-  data_ = expression;
+ModalVector::ModalVector(const blaze::DenseVector<VT, VF>& expression) noexcept
+    : owned_data_((~expression).size()) {
+  static_assert(cpp17::is_same_v<typename VT::ResultType, ModalVector>,
+                "You are attempting to assign the result of an expression that "
+                "is not a ModalVector to a ModalVector.");
+  reset_pointer_vector();
+  ~*this = expression;
 }
 
 template <typename VT, bool VF>
 ModalVector& ModalVector::operator=(
-    const blaze::Vector<VT, VF>& expression) {
+    const blaze::DenseVector<VT, VF>& expression) noexcept {
+  static_assert(cpp17::is_same_v<typename VT::ResultType, ModalVector>,
+                "You are attempting to assign the result of an expression that "
+                "is not a ModalVector to a ModalVector.");
   if (owning_ and (~expression).size() != size()) {
-    size_ = (~expression).size();
-    owned_data_ = InternalStorage_t(size_);
-    data_ = decltype(data_){owned_data_.data(), size_};
+    owned_data_.resize((~expression).size());
+    reset_pointer_vector();
   } else if (not owning_) {
     ASSERT((~expression).size() == size(), "Must copy into same size");
   }
-  data_ = expression;
+  ~*this = expression;
   return *this;
 }
 /// \endcond
+
+namespace MakeWithValueImpls {
+/// \brief Returns a ModalVector the same size as `input`, with each element
+/// equal to `value`.
+template <>
+SPECTRE_ALWAYS_INLINE ModalVector
+MakeWithValueImpl<ModalVector, ModalVector>::apply(const ModalVector& input,
+                                                 const double value) {
+  return ModalVector(input.size(), value);
+}
+}  // namespace MakeWithValueImpls
diff --git a/tests/Unit/DataStructures/Test_ModalVector.cpp b/tests/Unit/DataStructures/Test_ModalVector.cpp
index 6ef417964720..288e43b0a8be 100644
--- a/tests/Unit/DataStructures/Test_ModalVector.cpp
+++ b/tests/Unit/DataStructures/Test_ModalVector.cpp
@@ -3,20 +3,26 @@
 
 #include "tests/Unit/TestingFramework.hpp"
 
-#include <algorithm>
-#include <array>
-#include <cmath>
-#include <cstddef>
-#include <numeric>
+#include <algorithm>                         // for move
+#include <array>                             // for array, operator==
+#include <cmath>                             // for abs
+#include <cstddef>                           // for size_t
+#include <functional>                        // for std::reference_wrapper
+#include <numeric>                           // for iota
 
-#include "DataStructures/DataVector.hpp"
 #include "DataStructures/ModalVector.hpp"
 #include "ErrorHandling/Error.hpp"
+#include "Utilities/DereferenceWrapper.hpp"
 #include "Utilities/Gsl.hpp"
+#include "Utilities/Requires.hpp"
 #include "Utilities/StdArrayHelpers.hpp"
 #include "Utilities/StdHelpers.hpp"  // IWYU pragma: keep
 #include "tests/Unit/TestHelpers.hpp"
 
+// IWYU wants to include DataVector.hpp to access the function blaze::smpAssign
+//
+// IWYU pragma: no_include "DataStructures/DataVector.hpp"
+
 SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector",
                   "[DataStructures][Unit]") {
   ModalVector a{2};
@@ -59,7 +65,7 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector",
 
 SPECTRE_TEST_CASE("Unit.Serialization.ModalVector",
                   "[DataStructures][Unit][Serialization]") {
-  const size_t npts = 10;
+  const size_t npts = 100;
   ModalVector t(npts), tgood(npts);
   std::iota(t.begin(), t.end(), 1.2);
   std::iota(tgood.begin(), tgood.end(), 1.2);
@@ -76,7 +82,7 @@ SPECTRE_TEST_CASE("Unit.Serialization.ModalVector",
 
 SPECTRE_TEST_CASE("Unit.Serialization.ModalVector_Ref",
                   "[DataStructures][Unit][Serialization]") {
-  const size_t npts = 10;
+  const size_t npts = 11;
   ModalVector t(npts);
   std::iota(t.begin(), t.end(), 4.3);
   ModalVector t2;
@@ -101,10 +107,10 @@ SPECTRE_TEST_CASE("Unit.Serialization.ModalVector_Ref",
 SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector_Ref",
                   "[DataStructures][Unit]") {
   ModalVector data{1.43, 2.83, 3.94, 7.85};
+  CHECK(data.is_owning());
   ModalVector t;
   t.set_data_ref(&data);
   CHECK(not t.is_owning());
-  CHECK(data.is_owning());
   CHECK(t.data() == data.data());
   CHECK(t.size() == 4);
   CHECK(t[0] == 1.43);
@@ -125,6 +131,7 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector_Ref",
     data_2_ref.set_data_ref(&data_2);
     ModalVector data_3{2.43, 3.83, 4.94, 8.85};
     data_2_ref = std::move(data_3);
+    // check that data_2 now contain data_3's values
     CHECK(data_2[0] == 2.43);
     CHECK(data_2[1] == 3.83);
     CHECK(data_2[2] == 4.94);
@@ -145,6 +152,7 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector_Ref",
     ModalVector owned_data;
     // clang-tidy: false positive, used after it was moved
     owned_data = data_2_ref;  // NOLINT
+    CHECK(owned_data.is_owning());
     CHECK(owned_data[0] == 2.43);
     CHECK(owned_data[1] == 3.83);
     CHECK(owned_data[2] == 4.94);
@@ -156,15 +164,17 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector_Ref",
 
 // [[OutputRegex, Must copy into same size]]
 [[noreturn]] SPECTRE_TEST_CASE(
-                "Unit.DataStructures.ModalVector.ref_diff_size",
-                "[DataStructures][Unit]") {
+    "Unit.DataStructures.ModalVector.ref_diff_size",
+    "[DataStructures][Unit]") {
   ASSERTION_TEST();
 #ifdef SPECTRE_DEBUG
   ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector data2{1.43, 2.83, 5.94};
   ModalVector data_ref;
   data_ref.set_data_ref(&data);
-  ModalVector data2{1.43, 2.83, 3.94};
+  CHECK(data_ref[2] == 3.94);
   data_ref = data2;
+  CHECK(data_ref[2] == 5.94);
   ERROR("Failed to trigger ASSERT in an assertion test");
 #endif
 }
@@ -176,9 +186,9 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector_Ref",
   ASSERTION_TEST();
 #ifdef SPECTRE_DEBUG
   ModalVector data{1.43, 2.83, 3.94, 7.85};
+  ModalVector data2{1.43, 2.83, 3.94};
   ModalVector data_ref;
   data_ref.set_data_ref(&data);
-  ModalVector data2{1.43, 2.83, 3.94};
   data_ref = std::move(data2);
   ERROR("Failed to trigger ASSERT in an assertion test");
 #endif
@@ -191,6 +201,7 @@ void check_vectors(const T1& t1, const T2& t2) {
 }
 }  // namespace
 
+/// Tests of ModalVector math
 SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathAfterMove",
                   "[Unit][DataStructures]") {
   ModalVector m0(10, 3.0), m1(10, 9.0);
@@ -240,64 +251,121 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathAfterMove",
   }
 }
 
-SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
-                  "[Unit][DataStructures]") {
-  const size_t num_pts = 19;
-  ModalVector val{1., -2., 3., -4., 8., 12., -14.};
-  ModalVector nine(num_pts, 9.0);
+namespace {
+enum class UseRefWrap { None, Cref, Ref };
+
+template <UseRefWrap Wrap, class T,
+          Requires<Wrap == UseRefWrap::Cref> = nullptr>
+decltype(auto) wrap(const T& t) noexcept {
+  return std::cref(t);
+}
+template <UseRefWrap Wrap, class T,
+          Requires<Wrap == UseRefWrap::Ref> = nullptr>
+decltype(auto) wrap(T& t) noexcept {
+  return std::ref(t);
+}
+template <UseRefWrap Wrap, class T,
+          Requires<Wrap == UseRefWrap::None> = nullptr>
+decltype(auto) wrap(const T& t) noexcept {
+  return t;
+}
+
+// Wrap is used to wrap values in a std::reference_wrapper using std::cref and
+// std::ref, or to not wrap at all. This is done to verify that all math
+// operations work transparently with a `std::reference_wrapper` too.
+template <UseRefWrap WrapLeftOp, UseRefWrap WrapRightOp>
+void test_ModalVector_math() noexcept {
+  constexpr size_t num_pts = 19;
+  ModalVector val{1., 2., 3., -4., 8., 12., -14.};
   ModalVector one(num_pts, 1.0);
+  ModalVector eight(num_pts, 8.0);
+  ModalVector nine(num_pts, 9.0);
 
   // Test unary minus
-  check_vectors(-nine, ModalVector(num_pts, -9.0));
-
-  check_vectors(nine + 2.0, ModalVector(num_pts, 11.0));
-  check_vectors(2.0 + nine, ModalVector(num_pts, 11.0));
-  check_vectors(nine - 2.0, ModalVector(num_pts, 7.0));
-  check_vectors(2.0 - nine, ModalVector(num_pts, -7.0));
-  check_vectors(nine + nine, ModalVector(num_pts, 18.0));
-  check_vectors(nine + (one * nine), ModalVector(num_pts, 18.0));
-  check_vectors((one * nine) + nine, ModalVector(num_pts, 18.0));
-  check_vectors(nine - ModalVector(num_pts, 8.0), one);
-  check_vectors(nine - (one * nine), ModalVector(num_pts, 0.0));
-  check_vectors((one * nine) - nine, ModalVector(num_pts, 0.0));
-
-  check_vectors(ModalVector(num_pts, 18.0), (one / 0.5) * nine);
-
-  CHECK(-14 == min(val));
-  CHECK(12 == max(val));
-  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, abs(val));
-  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.}, fabs(val));
-
-  check_vectors(ModalVector(num_pts, 81.0), nine * nine);
-  check_vectors(ModalVector(num_pts, 81.0), nine * (nine * one));
-  check_vectors(ModalVector(num_pts, 81.0), (nine * nine) * one);
-  check_vectors(ModalVector(num_pts, 81.0), 9.0 * nine);
-  check_vectors(ModalVector(num_pts, 81.0), nine * 9.0);
-  check_vectors(ModalVector(num_pts, 1.0), nine / 9.0);
-
-  ModalVector dummy(nine * nine * 1.0);
+  check_vectors(-wrap<WrapLeftOp>(nine), ModalVector(num_pts, -9.0));
+
+  check_vectors(wrap<WrapLeftOp>(nine) + 2.0, ModalVector(num_pts, 11.0));
+  check_vectors(2.0 + wrap<WrapLeftOp>(nine), ModalVector(num_pts, 11.0));
+  check_vectors(wrap<WrapLeftOp>(nine) - 2.0, ModalVector(num_pts, 7.0));
+  check_vectors(2.0 - wrap<WrapLeftOp>(nine), ModalVector(num_pts, -7.0));
+  check_vectors(wrap<WrapLeftOp>(nine) + wrap<WrapRightOp>(nine),
+                ModalVector(num_pts, 18.0));
+  check_vectors(wrap<WrapLeftOp>(nine) +
+                    (wrap<WrapLeftOp>(one) * wrap<WrapRightOp>(nine)),
+                ModalVector(num_pts, 18.0));
+  check_vectors((wrap<WrapLeftOp>(one) * wrap<WrapLeftOp>(nine)) +
+                    wrap<WrapRightOp>(nine),
+                ModalVector(num_pts, 18.0));
+  check_vectors(wrap<WrapLeftOp>(nine) - ModalVector(num_pts, 8.0),
+                ModalVector(num_pts, 1.0));
+  check_vectors(wrap<WrapLeftOp>(nine) -
+                    (wrap<WrapRightOp>(one) * wrap<WrapLeftOp>(nine)),
+                ModalVector(num_pts, 0.0));
+  check_vectors((wrap<WrapLeftOp>(one) * wrap<WrapLeftOp>(nine)) -
+                    wrap<WrapRightOp>(nine),
+                ModalVector(num_pts, 0.0));
+
+  check_vectors(ModalVector(num_pts, -1.0 / 9.0),
+                -one / wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, -8.0 / 9.0),
+                -(wrap<WrapLeftOp>(nine) - wrap<WrapRightOp>(one)) /
+                    wrap<WrapLeftOp>(nine));
+  check_vectors(ModalVector(num_pts, 18.0),
+                (wrap<WrapLeftOp>(one) / 0.5) * wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, 1.0), 9.0 / wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, 1.0),
+                (wrap<WrapLeftOp>(one) * 9.0) / wrap<WrapRightOp>(nine));
+
+  check_vectors(ModalVector(num_pts, 81.0),
+                wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(nine));
+  check_vectors(ModalVector(num_pts, 81.0),
+                wrap<WrapLeftOp>(nine) *
+                    (wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(one)));
+  check_vectors(ModalVector(num_pts, 81.0),
+                (wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(nine)) *
+                    wrap<WrapLeftOp>(one));
+  check_vectors(ModalVector(num_pts, 81.0), 9.0 * wrap<WrapLeftOp>(nine));
+  check_vectors(ModalVector(num_pts, 81.0), wrap<WrapLeftOp>(nine) * 9.0);
+  check_vectors(ModalVector(num_pts, 1.0), wrap<WrapLeftOp>(nine) / 9.0);
+  check_vectors(ModalVector(num_pts, 1.0),
+                wrap<WrapLeftOp>(nine) / (wrap<WrapRightOp>(one) * 9.0));
+  check_vectors(ModalVector(num_pts, 9.0),
+                (wrap<WrapLeftOp>(one) * wrap<WrapLeftOp>(nine)) /
+                    wrap<WrapRightOp>(one));
+
+  CHECK(-14 == min(wrap<WrapLeftOp>(val)));
+  CHECK(12 == max(wrap<WrapLeftOp>(val)));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.},
+                abs(wrap<WrapLeftOp>(val)));
+  check_vectors(ModalVector{1., 2., 3., 4., 8., 12., 14.},
+                fabs(wrap<WrapLeftOp>(val)));
+
+  check_vectors(sqrt(wrap<WrapLeftOp>(nine)), ModalVector(num_pts, 3.0));
+
+  ModalVector dummy(wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(nine) * 1.0);
   check_vectors(ModalVector(num_pts, 81.0), dummy);
 
   // Test assignment
   ModalVector test_81(num_pts, -1.0);
-  test_81 = nine * nine;
+  // cannot assign to a reference_wrapper using an expr :(
+  test_81 = wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(nine);
   check_vectors(ModalVector(num_pts, 81.0), test_81);
   CHECK(test_81.is_owning());
   ModalVector test_81_ref(test_81.data(), test_81.size());
   test_81_ref = 0.0;
   test_81 = 0.0;
-  test_81_ref = nine * nine;
+  test_81_ref = wrap<WrapLeftOp>(nine) * wrap<WrapRightOp>(nine);
   check_vectors(ModalVector(num_pts, 81.0), test_81);
   CHECK(test_81.is_owning());
   check_vectors(ModalVector(num_pts, 81.0), test_81_ref);
   CHECK_FALSE(test_81_ref.is_owning());
   ModalVector second_81(num_pts);
-  second_81 = test_81;
+  second_81 = wrap<WrapLeftOp>(test_81);
   check_vectors(ModalVector(num_pts, 81.0), second_81);
   CHECK(second_81.is_owning());
   test_81_ref = 0.0;
   check_vectors(ModalVector(num_pts, 0.0), test_81_ref);
-  test_81_ref = second_81;
+  test_81_ref = wrap<WrapLeftOp>(second_81);
   check_vectors(ModalVector(num_pts, 81.0), test_81_ref);
   second_81 = 0.0;
   test_81_ref.set_data_ref(&test_81);
@@ -317,39 +385,61 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
   CHECK_FALSE(test_081.is_owning());
 
   ModalVector test_assignment(num_pts, 7.0);
-  test_assignment += nine;
+  test_assignment += wrap<WrapLeftOp>(nine);
   check_vectors(ModalVector(num_pts, 16.0), test_assignment);
   test_assignment += 3.0;
   check_vectors(ModalVector(num_pts, 19.0), test_assignment);
-  test_assignment += (nine * nine);
+  test_assignment += (wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(nine));
   check_vectors(ModalVector(num_pts, 100.0), test_assignment);
 
   test_assignment = 7.0;
-  test_assignment -= nine;
+  test_assignment -= wrap<WrapLeftOp>(nine);
   check_vectors(ModalVector(num_pts, -2.0), test_assignment);
   test_assignment -= 3.0;
   check_vectors(ModalVector(num_pts, -5.0), test_assignment);
-  test_assignment -= nine * nine;
+  test_assignment -= wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(nine);
   check_vectors(ModalVector(num_pts, -86.0), test_assignment);
 
   test_assignment = 2.0;
   test_assignment *= 3.0;
   check_vectors(ModalVector(num_pts, 6.0), test_assignment);
-  test_assignment *= nine;
+  test_assignment *= wrap<WrapLeftOp>(nine);
   check_vectors(ModalVector(num_pts, 54.0), test_assignment);
   test_assignment = 1.0;
-  test_assignment *= nine * nine;
+  test_assignment *= wrap<WrapLeftOp>(nine) * wrap<WrapLeftOp>(nine);
   check_vectors(ModalVector(num_pts, 81.0), test_assignment);
 
   test_assignment = 2.0;
   test_assignment /= 2.0;
   check_vectors(ModalVector(num_pts, 1.0), test_assignment);
+  test_assignment /= ModalVector(num_pts, 0.5);
+  check_vectors(ModalVector(num_pts, 2.0), test_assignment);
+  test_assignment /= (ModalVector(num_pts, 2.0) * ModalVector(num_pts, 3.0));
+  check_vectors(ModalVector(num_pts, 1.0 / 3.0), test_assignment);
 
   // Test assignment where the RHS is an expression that contains the LHS
   ModalVector x(num_pts, 4.);
-  x *= x;
-  check_vectors(ModalVector(num_pts, 16.0), x);
+  x += sqrt(wrap<WrapLeftOp>(x));
+  check_vectors(ModalVector(num_pts, 6.0), x);
+  x -= sqrt(wrap<WrapLeftOp>(x) - 2.0);
+  check_vectors(ModalVector(num_pts, 4.0), x);
+  x = sqrt(wrap<WrapLeftOp>(x));
+  check_vectors(ModalVector(num_pts, 2.0), x);
+  x *= wrap<WrapLeftOp>(x);
+  check_vectors(ModalVector(num_pts, 4.0), x);
+  x /= wrap<WrapLeftOp>(x);
+  check_vectors(ModalVector(num_pts, 1.0), x);
+
+  // Test composition of constant expressions with ModalVector math member
+  // functions
+  x = ModalVector(num_pts, 2.);
+  check_vectors(ModalVector(num_pts, 1.4142135623730951),
+                sqrt(abs(wrap<WrapLeftOp>(x))));
+  check_vectors(ModalVector(num_pts, 2.),
+                sqrt(wrap<WrapLeftOp>(x) * wrap<WrapRightOp>(x)));
+}
 
+void test_ModalVector_array_math() noexcept {
   // Test addition of arrays of ModalVectors to arrays of doubles.
   const ModalVector t1{0.0, 1.0, 2.0, 3.0};
   const ModalVector t2{-0.1, -0.2, -0.3, -0.4};
@@ -400,43 +490,31 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
   const ModalVector expected_d1{2.5, 3.4};
   const auto magnitude_d1 = magnitude(d1);
   CHECK_ITERABLE_APPROX(expected_d1, magnitude_d1);
-  const std::array<ModalVector, 2> d2{{ModalVector(2, 3.),
-                                             ModalVector(2, 4.)}};
+  const std::array<ModalVector, 2> d2{{ModalVector(2, 3.), ModalVector(2, 4.)}};
   const ModalVector expected_d2(2, 5.);
   const auto magnitude_d2 = magnitude(d2);
   CHECK_ITERABLE_APPROX(expected_d2, magnitude_d2);
   const std::array<ModalVector, 3> d3{
-      {ModalVector(2, 3.), ModalVector(2, -4.),
-       ModalVector(2, 12.)}};
+      {ModalVector(2, 3.), ModalVector(2, -4.), ModalVector(2, 12.)}};
   const ModalVector expected_d3(2, 13.);
   const auto magnitude_d3 = magnitude(d3);
   CHECK_ITERABLE_APPROX(expected_d3, magnitude_d3);
 }
+}  // namespace
 
-SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathWithDataVector",
+SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.Math",
                   "[Unit][DataStructures]") {
-  const size_t num_pts = 23;
-  ModalVector nine(num_pts, 9.0);
-  CHECK(nine.is_owning());
-  DataVector eight(num_pts, 8.0);
-
-  // Test math with DataVector
-  ModalVector test_72(num_pts, -1.0);
-  test_72 = nine * eight;
-  check_vectors(ModalVector(num_pts, 72.0), test_72);
-  CHECK(test_72.is_owning());
-  test_72 = eight * nine;
-  check_vectors(ModalVector(num_pts, 72.0), test_72);
-  CHECK(test_72.is_owning());
-  test_72 = test_72 / eight;
-  check_vectors(ModalVector(num_pts, 9.0), test_72);
-  test_72 = nine;
-  test_72 *= eight;
-  check_vectors(ModalVector(num_pts, 72.0), test_72);
-  CHECK(test_72.is_owning());
-  test_72 /= eight;
-  check_vectors(ModalVector(num_pts, 9.0), test_72);
-  CHECK(test_72.is_owning());
+  test_ModalVector_math<UseRefWrap::Cref, UseRefWrap::Cref>();
+  test_ModalVector_math<UseRefWrap::Cref, UseRefWrap::Ref>();
+  test_ModalVector_math<UseRefWrap::Cref, UseRefWrap::None>();
+  test_ModalVector_math<UseRefWrap::Ref, UseRefWrap::Cref>();
+  test_ModalVector_math<UseRefWrap::Ref, UseRefWrap::Ref>();
+  test_ModalVector_math<UseRefWrap::Ref, UseRefWrap::None>();
+  test_ModalVector_math<UseRefWrap::None, UseRefWrap::Cref>();
+  test_ModalVector_math<UseRefWrap::None, UseRefWrap::Ref>();
+  test_ModalVector_math<UseRefWrap::None, UseRefWrap::None>();
+
+  test_ModalVector_array_math();
 }
 
 // [[OutputRegex, Must copy into same size]]
@@ -452,3 +530,4 @@ SPECTRE_TEST_CASE("Unit.DataStructures.ModalVector.MathWithDataVector",
   ERROR("Failed to trigger ASSERT in an assertion test");
 #endif
 }
+

From dc43352d68fedbb51260bded72693fc78843a862 Mon Sep 17 00:00:00 2001
From: Prayush Kumar <prayush.kumar@ligo.org>
Date: Fri, 27 Jul 2018 02:09:09 -0400
Subject: [PATCH 4/4] fix IWYU's complaint

---
 src/DataStructures/ModalVector.hpp | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/DataStructures/ModalVector.hpp b/src/DataStructures/ModalVector.hpp
index 2099b878f2c6..5e3d80fc77e9 100644
--- a/src/DataStructures/ModalVector.hpp
+++ b/src/DataStructures/ModalVector.hpp
@@ -59,6 +59,7 @@ using std::abs;  // NOLINT
 // IWYU pragma: no_include <blaze/math/traits/BinaryMapTrait.h>
 // IWYU pragma: no_include <blaze/math/traits/UnaryMapTrait.h>
 // IWYU pragma: no_include <blaze/math/typetraits/TransposeFlag.h>
+// IWYU pragma: no_include "DataStructures/DataVector.hpp"
 
 // IWYU pragma: no_forward_declare blaze::DenseVector
 // IWYU pragma: no_forward_declare blaze::UnaryMapTrait