From 53ed2f1de35f29b6177b4a2811b6ab2bf8bea72c Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Wed, 3 Aug 2022 16:39:59 +0300
Subject: [PATCH] array: use pointers, not offsets

---
 src/array.hpp | 55 ++++++++++++++++++++-------------------------------
 src/main.cpp  |  4 ++--
 2 files changed, 23 insertions(+), 36 deletions(-)

diff --git a/src/array.hpp b/src/array.hpp
index 575be9b..e223b14 100644
--- a/src/array.hpp
+++ b/src/array.hpp
@@ -90,19 +90,16 @@ namespace detail
     class fixed_base
     {
     protected:
-        size_t offset_;
-
-        template <typename Storage>
-        void data_bounds_check(Storage& data) const
+        void data_bounds_check(size_t size) const
             {
 #ifndef NO_BOUNDS_CHECK
-                if (Shape::size > data.size())
+                if (Shape::size > size)
                     throw std::out_of_range("data array too small");
 #endif
-        }
+            }
 
     public:
-        fixed_base(size_t offset=0) : offset_(offset)
+        fixed_base()
             {
                 static_assert(Shape::fixed, "array shape must be compile-time fixed");
             }
@@ -112,7 +109,6 @@ namespace detail
 
         constexpr std::array<size_t, Shape::ndim> shape() const { return Shape::shape(); }
         constexpr std::array<size_t, Shape::ndim> strides() const { return Shape::strides(); }
-        constexpr size_t offset() const { return offset_; }
         constexpr size_t size() const { return Shape::size; }
     };
 
@@ -124,10 +120,8 @@ namespace detail
     protected:
         std::array<size_t, Shape::ndim> shape_;
         std::array<size_t, Shape::ndim> strides_;
-        size_t offset_;
 
-        template <typename Storage>
-        void data_bounds_check(Storage& data) const
+        void data_bounds_check(size_t size) const
             {
 #ifndef NO_BOUNDS_CHECK
                 size_t total = 1;
@@ -137,15 +131,14 @@ namespace detail
                     if (shape_[i] != dim(i))
                         throw std::out_of_range("mismatch with fixed shape");
                 }
-                total += offset_;
-                if (total > data.size())
+                if (total > size)
                     throw std::out_of_range("data array too small");
 #endif
         }
 
     public:
-        dynamic_base(const std::array<size_t, Shape::ndim> shape, size_t offset)
-            : shape_(shape), offset_(offset)
+        dynamic_base(const std::array<size_t, Shape::ndim> shape)
+            : shape_(shape)
             {
                 static_assert(!Shape::fixed, "array shape must not be compile-time fixed");
 
@@ -181,7 +174,6 @@ namespace detail
 
         const std::array<size_t, Shape::ndim>& shape() const { return shape_; }
         const std::array<size_t, Shape::ndim>& strides() const { return strides_; }
-        size_t offset() const { return offset_; }
         size_t size() const
             {
                 size_t size = 1;
@@ -194,13 +186,13 @@ namespace detail
 
 /*! Simple data-by-reference strided array class a la Fortran
  */
-template <typename Scalar, typename Shape, typename Storage = std::vector<Scalar> >
+template <typename Scalar, typename Shape>
 class Array : public std::conditional_t<Shape::fixed, detail::fixed_base<Shape>, detail::dynamic_base<Shape> >
 {
 private:
     typedef std::conditional_t<Shape::fixed, detail::fixed_base<Shape>, detail::dynamic_base<Shape> > Base;
 
-    Storage& data_;
+    Scalar *data_;
 
     template <size_t axis>
     static constexpr Eigen::Index eigen_shape()
@@ -210,16 +202,16 @@ private:
         }
 
 public:
-    Array(Storage& data, const std::array<size_t, Shape::ndim> shape, const size_t offset=0)
-        : Base(shape, offset), data_(data)
+    Array(Scalar *data, size_t size, const std::array<size_t, Shape::ndim> shape)
+        : Base(shape), data_(data)
         {
-            Base::data_bounds_check(data);
+            Base::data_bounds_check(size);
         }
 
-    Array(Storage& data, const size_t offset=0)
-        : Base(offset), data_(data)
+    Array(Scalar *data, size_t size)
+        : Base(), data_(data)
         {
-            Base::data_bounds_check(data);
+            Base::data_bounds_check(size);
         }
 
     template <typename... Idx>
@@ -238,7 +230,7 @@ public:
             return data_[index(idxs...)];
         }
 
-    Storage& data() { return data_; }
+    Scalar* data() { return data_; }
 
     template <typename... Idx>
     size_t index(Idx... idxs) const
@@ -247,7 +239,7 @@ public:
                           "number of indices must be smaller than the number of dimensions");
 
             const std::array<size_t, Shape::ndim> m{idxs...};
-            size_t idx = Base::offset();
+            size_t idx = 0;
 
             for (size_t i = 0; i < sizeof...(idxs); ++i) {
                 idx += Base::stride(i) * m[i];
@@ -257,10 +249,6 @@ public:
 #endif
             }
 
-#ifndef NO_BOUNDS_CHECK
-            if (idx >= data_.size())
-                throw std::out_of_range("total index out of bounds");
-#endif
             return idx;
         }
 
@@ -275,7 +263,7 @@ public:
                           "number of indices must be small than the number of dimensions");
 
             size_t offset = index(idxs...);
-            return {data_, offset};
+            return {data_ + offset, Base::size() - offset};
         }
 
     /*! Extract partial matrix (compile-time dynamic size) */
@@ -292,7 +280,7 @@ public:
             std::array<size_t, Shape::ndim - nidxs> shape;
             for (size_t i = nidxs; i < Shape::ndim; ++i)
                 shape[i - nidxs] = Base::dim(i);
-            return {data_, shape, offset};
+            return {data_ + offset, Base::size() - offset, shape};
         }
 
     typedef Eigen::Matrix<Scalar, eigen_shape<0>(), eigen_shape<1>(), Eigen::RowMajor> EigenMatrix;
@@ -301,8 +289,7 @@ public:
     EigenMap matrix() const
         {
             static_assert(Shape::ndim == 2, "matrix must be two-dimensional");
-            return {(Scalar *)data_.data() + Base::offset(),
-                static_cast<Eigen::Index>(Base::dim(0)), static_cast<Eigen::Index>(Base::dim(1))};
+            return {data_, static_cast<Eigen::Index>(Base::dim(0)), static_cast<Eigen::Index>(Base::dim(1))};
         }
 
     operator EigenMap() const { return matrix(); }
diff --git a/src/main.cpp b/src/main.cpp
index 1a1444e..5ffac0b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -17,7 +17,7 @@ void dump(Array<Scalar,Shape> Q)
         }
     }
 
-    for (size_t i = 0; i < Q.data().size(); ++i) {
+    for (size_t i = 0; i < Q.size(); ++i) {
         std::cout << Q.data()[i] << std::endl;
     }
 }
@@ -25,7 +25,7 @@ void dump(Array<Scalar,Shape> Q)
 int main()
 {
     std::vector<double> v(7*2*2*2);
-    Array<double,Shape<Dynamic,Dynamic,2,2> > Q0(v,{7,2,2,2});
+    Array<double,Shape<Dynamic,Dynamic,2,2> > Q0(v.data(), v.size(), {7,2,2,2});
     auto Q = Q0.part(0u,1u);
 
     S_2(Q);
-- 
GitLab