array: deal with alignment

src/array.hpp

@@ -181,7 +181,7 @@ namespace detail
 
 /*! Simple data-by-reference strided array class a la Fortran
  */
-template <typename Scalar, typename Shape>
+template <typename Scalar, typename Shape, size_t Alignment=0>
 class Array : public std::conditional_t<Shape::fixed, detail::fixed_base<Shape>, detail::dynamic_base<Shape> >
 {
 private:
@@ -196,12 +196,28 @@ private:
             return (n == Dynamic) ? Eigen::Dynamic : n;
         }
 
+    static constexpr int eigen_alignment()
+        {
+            switch (alignment()) {
+            case 128: return Eigen::Aligned128;
+            case 64: return Eigen::Aligned64;
+            case 32: return Eigen::Aligned32;
+            case 16: return Eigen::Aligned16;
+            case 8: return Eigen::Aligned8;
+            default: return Eigen::Unaligned;
+            }
+        };
+
 public:
     Array(Scalar *data, size_t size, const std::array<size_t, Shape::ndim> shape)
         : Base(shape), data_(data)
         {
             static_assert(!Shape::fixed, "array shape must not be compile-time fixed");
             Base::data_bounds_check(size);
+#ifndef NO_BOUNDS_CHECK
+            if (reinterpret_cast<intptr_t>(data) % alignment() != 0)
+                throw std::runtime_error("data is not aligned");
+#endif
         }
 
     Array(Scalar *data, size_t size)
@@ -209,6 +225,10 @@ public:
         {
             static_assert(Shape::fixed, "array shape must be compile-time fixed");
             Base::data_bounds_check(size);
+#ifndef NO_BOUNDS_CHECK
+            if (reinterpret_cast<intptr_t>(data) % alignment() != 0)
+                throw std::runtime_error("data is not aligned");
+#endif
         }
 
     template <typename... Idx>
@@ -229,6 +249,18 @@ public:
 
     Scalar* data() { return data_; }
 
+    static constexpr int alignment()
+        {
+            if (Alignment > 0)
+                return Alignment;
+            else if (sizeof(Scalar) % 16 == 0)
+                return 16;
+            else if (sizeof(Scalar) % 8 == 0)
+                return 8;
+            else
+                return 1;
+        };
+
     template <typename... Idx>
     size_t index(Idx... idxs) const
         {
@@ -281,7 +313,7 @@ public:
         }
 
     typedef Eigen::Matrix<Scalar, eigen_shape<0>(), eigen_shape<1>(), Eigen::RowMajor> EigenMatrix;
-    typedef Eigen::Map<EigenMatrix> EigenMap;
+    typedef Eigen::Map<EigenMatrix, eigen_alignment()> EigenMap;
 
     EigenMap matrix() const
         {

src/main.cpp

@@ -38,7 +38,7 @@ int main()
 
     mat = mat * mat * mat;
 
-    static_assert(std::is_same<decltype(mat), Eigen::Map<Eigen::Matrix<double, 2, 2, Eigen::RowMajor> > >::value, "wrong type");
+    static_assert(std::is_same<decltype(mat), Eigen::Map<Eigen::Matrix<double, 2, 2, Eigen::RowMajor>, Eigen::Aligned8 > >::value, "wrong type");
 
     dump(Q0);
     return 0;