Modernisation of Vector implementation

This modernises the implementation of Vector, principally by writing the operations using STL algorithms for the most part. Because Vector is now essentially a thin wrapper around Array combined with the appropriate STL algorithms, it can now be a header-only class.
gambitproject · Dec 18, 2024 · 4502065 · 4502065
1 parent c76d074
commit 4502065
Show file tree

Hide file tree

Showing 8 changed files with 129 additions and 264 deletions.
diff --git a/Makefile.am b/Makefile.am
@@ -248,7 +248,6 @@ core_SOURCES = \
 	src/core/array.h \
 	src/core/list.h \
 	src/core/vector.h \
-	src/core/vector.imp \
 	src/core/recarray.h \
 	src/core/matrix.h \
 	src/core/matrix.imp \
@@ -258,7 +257,6 @@ core_SOURCES = \
 	src/core/integer.h \
 	src/core/rational.cc \
 	src/core/rational.h \
-	src/core/vector.cc \
 	src/core/matrix.cc \
 	src/core/sqmatrix.cc \
 	src/core/function.cc \

diff --git a/src/core/array.h b/src/core/array.h
@@ -81,6 +81,12 @@ template <class T> class Array {
       m_index++;
       return *this;
     }
+    iterator operator++(int)
+    {
+      auto ret = *this;
+      m_index++;
+      return ret;
+    }
     bool operator==(const iterator &it) const
     {
       return (m_array == it.m_array) && (m_index == it.m_index);
@@ -108,6 +114,12 @@ template <class T> class Array {
       m_index++;
       return *this;
     }
+    const_iterator operator++(int)
+    {
+      auto ret = *this;
+      m_index++;
+      return ret;
+    }
     bool operator==(const const_iterator &it) const
     {
       return (m_array == it.m_array) && (m_index == it.m_index);

diff --git a/src/core/matrix.imp b/src/core/matrix.imp
@@ -169,7 +169,7 @@ template <class T> void Matrix<T>::CMultiply(const Vector<T> &in, Vector<T> &out
     T sum = (T)0;
 
     T *src1 = this->data[i] + this->mincol;
-    T *src2 = in.data + this->mincol;
+    auto src2 = in.begin();
     int j = this->maxcol - this->mincol + 1;
     while (j--) {
       sum += *(src1++) * *(src2++);
@@ -217,7 +217,7 @@ template <class T> void Matrix<T>::RMultiply(const Vector<T> &in, Vector<T> &out
     T k = in[i];
 
     T *src = this->data[i] + this->mincol;
-    T *dst = out.data + this->mincol;
+    auto dst = out.begin();
     int j = this->maxcol - this->mincol + 1;
     while (j--) {
       *(dst++) += *(src++) * k;

diff --git a/src/core/vector.cc b/src/core/vector.cc
diff --git a/src/core/vector.h b/src/core/vector.h
@@ -23,13 +23,20 @@
 #ifndef LIBGAMBIT_VECTOR_H
 #define LIBGAMBIT_VECTOR_H
 
+#include <numeric>
+
 namespace Gambit {
 
 template <class T> class Matrix;
 
 /// A mathematical vector: a list of numbers with the standard math operators
 template <class T> class Vector : public Array<T> {
-  friend class Matrix<T>;
+private:
+  // check vector for identical boundaries
+  bool Check(const Vector<T> &v) const
+  {
+    return (v.mindex == this->mindex && v.maxdex == this->maxdex);
+  }
 
 public:
   /** Create a vector of length len, starting at 1 */
@@ -42,35 +49,119 @@ template <class T> class Vector : public Array<T> {
   ~Vector() override = default;
 
   /** Assignment operator: requires vectors to be of same length */
-  Vector<T> &operator=(const Vector<T> &V);
+  Vector<T> &operator=(const Vector<T> &V)
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    Array<T>::operator=(V);
+    return *this;
+  }
   /** Assigns the value c to all components of the vector */
-  Vector<T> &operator=(T c);
-
-  Vector<T> operator+(const Vector<T> &V) const;
-  Vector<T> &operator+=(const Vector<T> &V);
-
-  Vector<T> operator-();
-  Vector<T> operator-(const Vector<T> &V) const;
-  Vector<T> &operator-=(const Vector<T> &V);
-
-  Vector<T> operator*(T c) const;
-  Vector<T> &operator*=(T c);
-  T operator*(const Vector<T> &V) const;
-
-  Vector<T> operator/(T c) const;
-
-  bool operator==(const Vector<T> &V) const;
+  Vector<T> &operator=(const T &c)
+  {
+    std::fill(this->begin(), this->end(), c);
+    return *this;
+  }
+
+  Vector<T> operator+(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), tmp.begin(), std::plus<>());
+    return tmp;
+  }
+
+  Vector<T> &operator+=(const Vector<T> &V)
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), this->begin(), std::plus<>());
+    return *this;
+  }
+
+  Vector<T> operator-() const
+  {
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(tmp.cbegin(), tmp.cend(), tmp.begin(), std::negate<>());
+    return tmp;
+  }
+
+  Vector<T> operator-(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), tmp.begin(), std::minus<>());
+    return tmp;
+  }
+
+  Vector<T> &operator-=(const Vector<T> &V)
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    std::transform(this->cbegin(), this->cend(), V.cbegin(), this->begin(), std::minus<>());
+    return *this;
+  }
+
+  Vector<T> operator*(const T &c) const
+  {
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), tmp.begin(), [&](const T &v) { return v * c; });
+    return tmp;
+  }
+
+  Vector<T> &operator*=(const T &c)
+  {
+    std::transform(this->cbegin(), this->cend(), this->begin(), [&](const T &v) { return v * c; });
+    return *this;
+  }
+
+  T operator*(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    return std::inner_product(this->begin(), this->end(), V.begin(), static_cast<T>(0));
+  }
+
+  Vector<T> operator/(const T &c) const
+  {
+    Vector<T> tmp(this->mindex, this->maxdex);
+    std::transform(this->cbegin(), this->cend(), tmp.begin(), [&](const T &v) { return v / c; });
+    return tmp;
+  }
+
+  bool operator==(const Vector<T> &V) const
+  {
+    if (!Check(V)) {
+      throw DimensionException();
+    }
+    return Array<T>::operator==(V);
+  }
   bool operator!=(const Vector<T> &V) const { return !(*this == V); }
 
   /** Tests if all components of the vector are equal to a constant c */
-  bool operator==(T c) const;
-  bool operator!=(T c) const;
+  bool operator==(const T &c) const
+  {
+    return std::all_of(this->begin(), this->end(), [&](const T &v) { return v == c; });
+  }
+  bool operator!=(const T &c) const
+  {
+    return std::any_of(this->begin(), this->end(), [&](const T &v) { return v != c; });
+  }
 
   // square of length
-  T NormSquared() const;
-
-  // check vector for identical boundaries
-  bool Check(const Vector<T> &v) const;
+  T NormSquared() const
+  {
+    return std::accumulate(this->begin(), this->end(), static_cast<T>(0),
+                           [](const T &t, const T &v) { return t + v * v; });
+  }
 };
 
 } // end namespace Gambit