array: more checks

src/action.hpp

@@ -8,8 +8,8 @@
 template <typename Scalar>
 inline Scalar S_2(Array<Scalar,2> Q)
 {
-    Q(0u,0u) = 123;
-    Q(1u,1u) = 123;
+    Q(0u,0u) = 1;
+    Q(0u,1u) = 2;
     return Scalar(0);
 }
 

src/array.hpp

@@ -4,12 +4,13 @@
 #include <vector>
 #include <array>
 #include <type_traits>
+#include <initializer_list>
 #include <stdexcept>
 
 
-/*! Simple data-by-reference strided array class
+/*! Simple data-by-reference strided array class a la Fortran
  */
-template <typename Scalar, int NDim, typename Storage = std::vector<Scalar> >
+template <typename Scalar, size_t NDim, typename Storage = std::vector<Scalar> >
 class Array
 {
 private:
@@ -19,24 +20,36 @@ private:
     size_t offset_ = 0;
 
 public:
-    Array(Array<Scalar,NDim>&& array)
-        : data_(array.data_), shape_(array.shape_), stride_(array.stride_), offset_(array.offset_)
-        {}
-
-    Array(Storage& data, const std::array<size_t, NDim> shape, const std::array<size_t, NDim> stride, const size_t offset=0)
+    Array(Storage& data, const std::array<size_t, NDim> shape, const std::array<size_t, NDim> stride, const size_t offset)
         : data_(data), shape_(shape), stride_(stride), offset_(offset)
         {}
 
-    Array(Storage& data, const std::array<size_t, NDim> shape)
+    Array(Storage& data, const std::array<size_t, NDim> shape, const bool row_major=true)
         : data_(data), shape_(shape), offset_(0)
         {
-            // row-major strides
-            for (size_t i = NDim; i > 0; --i) {
-                if (i == NDim)
-                    stride_[i-1] = 1;
-                else
-                    stride_[i-1] = shape[i] * stride_[i];
+            if (row_major) {
+                for (size_t i = NDim; i > 0; --i) {
+                    if (i == NDim)
+                        stride_[i-1] = 1;
+                    else
+                        stride_[i-1] = shape[i] * stride_[i];
+                }
+            } else {
+                for (size_t i = 0; i < NDim; ++i) {
+                    if (i == 0)
+                        stride_[i] = 1;
+                    else
+                        stride_[i] = shape[i-1] * stride_[i-1];
+                }
             }
+
+#ifndef NO_BOUNDS_CHECK
+            size_t total = 1;
+            for (size_t i = 0; i < NDim; ++i)
+                total *= shape[i];
+            if (total > data.size())
+                throw std::out_of_range("data array too small");
+#endif
         }
 
     template <typename... Idx>
@@ -48,7 +61,7 @@ public:
     template <size_t axis>
     Array<Scalar, NDim-1> slice(size_t pos=0) const
         {
-            static_assert(axis < NDim, "invalid slice index");
+            static_assert(axis < NDim, "invalid axis");
 
             std::array<size_t, NDim-1> shape;
             std::array<size_t, NDim-1> stride;
@@ -70,7 +83,7 @@ public:
     template <size_t axis>
     Array<Scalar, NDim> slice(size_t begin, size_t end) const
         {
-            static_assert(axis < NDim, "invalid slice index");
+            static_assert(axis < NDim, "invalid axis");
 
             std::array<size_t, NDim> shape;
             std::array<size_t, NDim> stride;
@@ -90,12 +103,20 @@ public:
         }
 
     Storage& data() { return data_; }
-    
+
     template <size_t axis>
-    const size_t shape() const { return shape_[axis]; }
+    const size_t shape() const
+        {
+            static_assert(axis < NDim, "invalid axis");
+            return shape_[axis];
+        }
 
     template <size_t axis>
-    const size_t stride() const { return stride_[axis]; }
+    const size_t stride() const
+        {
+            static_assert(axis < NDim, "invalid axis");
+            return stride_[axis];
+        }
 
     size_t offset() const { return offset_; }
 
@@ -116,6 +137,10 @@ public:
 #endif
             }
 
+#ifndef NO_BOUNDS_CHECK
+            if (idx >= data_.size())
+                throw std::out_of_range("total index out of bounds");
+#endif
             return idx;
         }
 };

src/main.cpp

 #include "common.hpp"
+#include "array.hpp"
 #include "action.hpp"
 
 
+template <typename Scalar=double>
+void dump(Array<Scalar,3> Q)
+{
+    for (size_t i = 0; i < Q.template shape<0>(); ++i) {
+        for (size_t j = 0; j < Q.template shape<1>(); ++j) {
+            for (size_t k = 0; k < Q.template shape<2>(); ++k) {
+                std::cout << i << "," << j << "," << k << " = " << Q(i,j,k) << std::endl;
+            }
+        }
+    }
+
+    for (size_t i = 0; i < Q.data().size(); ++i) {
+        std::cout << Q.data()[i] << std::endl;
+    }
+}
+
 int main()
 {
-    std::vector<ADComplex> v(2*2*2);
-    Array<ADComplex,3> Q(v,{2,2,2});
+    std::vector<double> v(2*2*2);
+    Array<double,3> Q(v,{2,2,2},false);
 
     S_2(Q.slice<2>());
 
-    for (size_t i = 0; i < Q.shape<0>(); ++i) {
-        for (size_t j = 0; j < Q.shape<1>(); ++j) {
-            for (size_t k = 0; k < Q.shape<2>(); ++k) {
-                std::cout << Q(i,j,k) << std::endl;
-            }
-        }
-    }
+    dump<>(Q);
 
     return 0;
 }