diff --git a/src/array.hpp b/src/array.hpp
index bbdec4f37f03210b2624b2d909e82847eb810823..575be9b4a09a62e612e8afc2bdb9098b2d929975 100644
--- a/src/array.hpp
+++ b/src/array.hpp
@@ -113,6 +113,7 @@ 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; }
     };
 
     /*! Base class for compile-time dynamic size arrays
@@ -178,9 +179,16 @@ namespace detail
                 return n == Dynamic ? strides_[axis] : n;
             }
 
-        constexpr std::array<size_t, Shape::ndim>& shape() const { return shape_; }
-        constexpr std::array<size_t, Shape::ndim>& strides() const { return strides_; }
-        constexpr size_t offset() const { return offset_; }
+        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;
+                for (size_t i = 0; i < Shape::ndim; ++i)
+                    size *= dim(i);
+                return size;
+            }
     };
 }
 
diff --git a/src/main.cpp b/src/main.cpp
index 8efa8a7dc404e5725eaf14191ca3ddfa75ba7270..1a1444e2a10cff562da76d6337763fc90e2fb8c0 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,12 +6,13 @@ using namespace std::complex_literals;
 
 template <typename Scalar, typename Shape>
 void dump(Array<Scalar,Shape> Q)
-#if 1
 {
     for (size_t i = 0; i < Q.dim(0); ++i) {
         for (size_t j = 0; j < Q.dim(1); ++j) {
             for (size_t k = 0; k < Q.dim(2); ++k) {
-                std::cout << i << "," << j << "," << k << " = " << Q(i,j,k) << std::endl;
+                for (size_t l = 0; l < Q.dim(3); ++l) {
+                    std::cout << i << "," << j << "," << k << "," << l << " = " << Q(i,j,k,l) << std::endl;
+                }
             }
         }
     }
@@ -20,28 +21,25 @@ void dump(Array<Scalar,Shape> Q)
         std::cout << Q.data()[i] << std::endl;
     }
 }
-#else
-;
-#endif
 
 int main()
 {
-    std::vector<double> v(2*2*2);
-    Array<double,Shape<Dynamic,2,2> > Q0(v,{2,2,2});
-    auto Q = Q0.part(0u);
+    std::vector<double> v(7*2*2*2);
+    Array<double,Shape<Dynamic,Dynamic,2,2> > Q0(v,{7,2,2,2});
+    auto Q = Q0.part(0u,1u);
 
     S_2(Q);
 
+    std::cout << "size = " << Q.size() << std::endl;
+
     dump(Q0);
 
     auto mat = Q.matrix();
 
     mat = mat * mat * mat;
 
-    if constexpr(true) {
-    }
+    static_assert(std::is_same<decltype(mat), Eigen::Map<Eigen::Matrix<double, 2, 2, Eigen::RowMajor> > >::value, "wrong type");
 
     dump(Q0);
-
     return 0;
 }