diff --git a/.gitignore b/.gitignore
index 29db1c4..bfa7ba2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 *~
 /build/
+/build_03/
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 597f6fe..5bb9719 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,40 +1,55 @@
 cmake_minimum_required(VERSION 3.7.0)
 project(µSpectre_prototype)
 
 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 add_compile_options(-Wall -Wextra -Werror)
 
 
 
 enable_testing()
 
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake)
-
+if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+  # using Clang
+elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+  # using GCC
+  string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type )
+  if ("release" STREQUAL "${build_type}" )
+    add_compile_options(-march=native)
+    add_compile_options(-mavx)
+    add_compile_options(-mfma)
+    add_compile_options(-msse)
+  endif()
+elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+  # using Intel C++
+elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC")
+  # using Visual Studio C++
+endif()
 
 find_package(Boost COMPONENTS unit_test_framework REQUIRED)
 find_package(Eigen3 3.3 REQUIRED NO_MODULE)
 
 include_directories(
   tests
   src
   ${EIGEN3_INCLUDE_DIRS}
   )
 
 #find_package(MPI REQUIRED)
 find_package(FFTW REQUIRED)
 #find_package(FFTWMPI REQUIRED)
 
 add_subdirectory(
   "${CMAKE_SOURCE_DIR}/src/"
   )
 
 #build tests
 file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/test_*.cc")
 
 add_executable(main_test_suite tests/main_test_suite.cc ${TEST_SRCS})
 target_link_libraries(main_test_suite ${Boost_LIBRARIES} Eigen3::Eigen materials)
 add_test(main_test_suite main_test_suite)
 
 
diff --git a/src/common/units.hh b/src/common/units.hh
index 96910f4..81dc67d 100644
--- a/src/common/units.hh
+++ b/src/common/units.hh
@@ -1,229 +1,312 @@
 /**
  * file   units.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   02 May 2017
  *
  * @brief  Implements physical units without runtime overhead. This is inspi-
  *         red by section 28.7 of Stroustrup's "The C++ Programming Language",
  *         but extended to rational powers of basic units (fracture mechanics!)
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include <ratio>
 #include <limits>
 #include <iostream>
 #include <boost/preprocessor/seq/for_each.hpp>
 #include "common.hh"
 #include <Eigen/src/Core/functors/BinaryFunctors.h>
+#include <Eigen/Dense>
 
 
 #ifndef UNITS_H
 #define UNITS_H
 
 namespace muSpectre {
 
   namespace units {
 
     //----------------------------------------------------------------------------//
     template <typename m__ratio, typename kg_ratio, typename s__ratio>
     struct unit {
       using m  = m__ratio;
       using kg = kg_ratio;
       using s  = s__ratio;
     };
 
     //----------------------------------------------------------------------------//
     template <typename unit1, typename unit2>
     struct uplus {
       using type = unit<std::ratio_add<typename unit1::m,  typename unit2::m >,
                         std::ratio_add<typename unit1::kg, typename unit2::kg>,
                         std::ratio_add<typename unit1::s,  typename unit2::s >>;
     };
 
     //----------------------------------------------------------------------------//
     template<typename unit1, typename unit2>
     using unit_plus = typename uplus<unit1, unit2>::type;
 
     //----------------------------------------------------------------------------//
     template <typename unit1, typename unit2>
     struct uminus {
       using type = unit<std::ratio_subtract<typename unit1::m,  typename unit2::m >,
                         std::ratio_subtract<typename unit1::kg, typename unit2::kg>,
                         std::ratio_subtract<typename unit1::s,  typename unit2::s >>;
     };
 
     //----------------------------------------------------------------------------//
     template<typename unit1, typename unit2>
     using unit_minus = typename uminus<unit1, unit2>::type;
 
-    using ratio1 = std::ratio<1,1>;
-    using ratio0 = std::ratio<0,1>;
+    using ratio1 = std::ratio<1>;
+    using ratio0 = std::ratio<0>;
 
     using nodim = unit<ratio0, ratio0, ratio0>;
     using m  = unit<ratio1, ratio0, ratio0>;
     using kg = unit<ratio0, ratio1, ratio0>;
     using s  = unit<ratio0, ratio0, ratio1>;
 
     using si_speed = unit_minus<m, s>;
     using si_acceleration = unit_minus<si_speed, s>;
     using si_area = unit_plus<m, m>;
 
     using N  = unit_plus<kg, si_acceleration>;
     using Pa = unit_minus<N, si_area>;
     using J  = unit_plus<N,m>;
 
   }  // units
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
   struct Quantity {
     T val;
     explicit constexpr Quantity(T r): val{r} {}
-    explicit constexpr Quantity(): val{T()} {}
+    explicit Quantity(size_t rows, size_t cols): val(rows, cols){}
+    explicit Quantity(size_t i, size_t j, size_t k, size_t l): val(i, j, k, l){}
+    explicit constexpr Quantity(): val(){}
+  };
+
+  //----------------------------------------------------------------------------//
+  template<typename T, typename U, size_t Rows>
+  struct Quantity<Eigen::Matrix<T, Rows, Eigen::Dynamic>, U> {
+    Eigen::Matrix<T, Rows, Eigen::Dynamic> val;
+    Quantity(int dynrows, int dyncols): val(dynrows, dyncols) {}
+    constexpr Quantity() = default;
+  };
+
+
+  //----------------------------------------------------------------------------//
+  template<typename U>
+  template<typename T, size_t Rows, size_t Cols>
+  struct Quantity<Eigen::Matrix<T, Rows, Cols>, U> {
+    Eigen::Matrix<T, Rows, Cols> val;
+    Quantity(int dynrows, int dyncols): val(dynrows, dyncols) {}
+    constexpr Quantity() = default;
   };
 
   //----------------------------------------------------------------------------//
   template<typename T, typename U>
   std::ostream & operator<<(std::ostream & os, const Quantity<T, U> Q) {
     os << Q.val
        << "m^"  << U::m ::num/U::m ::den
        << "kg^" << U::kg::num/U::kg::den
        << "s^"  << U::s ::num/U::s ::den;
     return os;
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
   bool operator==(Quantity<T,U> x, Quantity<T,U> y) {
     return x.val == y.val;
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
-  Quantity<T,U> operator+(Quantity<T,U> x, Quantity<T,U> y) {
+  Quantity<T,U> operator+(const Quantity<T,U>& x, const Quantity<T,U> & y) {
     return Quantity<T,U>(x.val + y.val);
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
-  Quantity<T,U> operator-(Quantity<T,U> x, Quantity<T,U> y) {
+  Quantity<T,U>& operator+=(Quantity<T,U>& x, const Quantity<T,U>& y) {
+    x.val += y.val;
+    return x;
+  }
+
+  //----------------------------------------------------------------------------//
+  template <typename T, typename U>
+  Quantity<T,U> operator-(const Quantity<T,U> & x, const Quantity<T,U> & y) {
     return Quantity<T,U>(x.val - y.val);
   }
 
+  //----------------------------------------------------------------------------//
+  template <typename T, typename U>
+  Quantity<T,U>& operator-=(Quantity<T,U> & x, const Quantity<T,U> & y) {
+    x.val -= y.val;
+    return x;
+  }
+
   //----------------------------------------------------------------------------//
   template <typename T, typename U1, typename U2>
-  Quantity<T, units::unit_plus<U1, U2>> operator*(Quantity<T, U1> x, Quantity<T, U2> y) {
+  Quantity<T, units::unit_plus<U1, U2>> operator*(const Quantity<T, U1> & x,
+                                                  const Quantity<T, U2> & y) {
     return Quantity<T, units::unit_plus<U1, U2>>{x.val*y.val};
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U1, typename U2>
-  Quantity<T, units::unit_minus<U1, U2>> operator/(Quantity<T, U1> x, Quantity<T, U2> y) {
+  Quantity<T, units::unit_minus<U1, U2>> operator/(const Quantity<T, U1> & x,
+                                                   const Quantity<T, U2> & y) {
     return Quantity<T, units::unit_minus<U1, U2>>{x.val/y.val};
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
   Quantity<T,U> operator*(Quantity<T,U> x, Real y) {
     return Quantity<T,U>{x.val*y};
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
   Quantity<T,U> operator*(Real x, Quantity<T,U> y) {
     return Quantity<T,U>{y.val*x};
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
   Quantity<T,U> operator/(Quantity<T,U> x, Real y) {
     return Quantity<T,U>{x.val/y};
   }
 
   //----------------------------------------------------------------------------//
   template <typename T, typename U>
   Quantity<T,U> operator/(Real x, Quantity<T, U> y) {
     return Quantity<T, units::nodim>{x}/y;
   }
 
   //----------------------------------------------------------------------------//
   //! string literals for convenient use of units
 #define LITERAL(literal_string) \
   auto operator""_##literal_string(long double r) { \
     return Quantity<Real, units::literal_string> {Real(r)};  \
   } \
   auto operator""_##literal_string(unsigned long long r) { \
     return Quantity<Real, units::literal_string> {Real(r)}; \
   }
   // Warning: no literals seem possible for complex scalars
 
 #define LITERALS (m)(kg)(s)(N)(Pa)(J)
 #define MACRO(r, dummy, literal_string) LITERAL(literal_string)
 
   BOOST_PP_SEQ_FOR_EACH(MACRO, _, LITERALS)
 
 
 }  // muSpectre
 
   /** Add num_traits for use in eigen arrays
    */
   namespace Eigen {
 
     /* ---------------------------------------------------------------------- */
     template<>
     template<typename T, typename U>
     struct NumTraits<muSpectre::Quantity<T, U>>
       : NumTraits<T>
     {};
 
     /* ---------------------------------------------------------------------- */
     template<>
     template<typename T, typename U>
     struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U>, muSpectre::Real, internal::scalar_product_op<muSpectre::Quantity<T, U>, muSpectre::Real>>
     {
       typedef muSpectre::Quantity<T,U> ReturnType;
     };
 
     /* ---------------------------------------------------------------------- */
     template<>
     template<typename T, typename U>
-    struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U>, muSpectre::Real, internal::scalar_quotient_op<muSpectre::Quantity<T, U>, muSpectre::Real>>
+    struct ScalarBinaryOpTraits<muSpectre::Real, muSpectre::Quantity<T, U>, internal::scalar_product_op<muSpectre::Real, muSpectre::Quantity<T, U>>>
+    {
+      typedef muSpectre::Quantity<T,U> ReturnType;
+    };
+
+    //----------------------------------------------------------------------------//
+    //! required for matrix multiplications with same units
+    template<>
+    template<typename U1>
+    struct ScalarBinaryOpTraits<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>, internal::scalar_product_op<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>>>
+    {
+      typedef muSpectre::Quantity<muSpectre::Real,muSpectre::units::unit_plus<U1, U1>> ReturnType;
+    };
+
+    //----------------------------------------------------------------------------//
+    //! required for matrix multiplications with differing units
+    template<>
+    template<typename T, typename U1, typename U2>
+    struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>, internal::scalar_product_op<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>>>
+    {
+      typedef muSpectre::Quantity<T,muSpectre::units::unit_plus<U1, U2>> ReturnType;
+    };
+
+    /* ---------------------------------------------------------------------- */
+    template<>
+    template<typename T, typename U>
+    struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U>, muSpectre::Real,
+                                internal::scalar_quotient_op<muSpectre::Quantity<T, U>, muSpectre::Real>>
     {
       typedef muSpectre::Quantity<T,U> ReturnType;
     };
 
     /* ---------------------------------------------------------------------- */
     template<>
     template<typename T, typename U>
     struct ScalarBinaryOpTraits<muSpectre::Real, muSpectre::Quantity<T, U>,
-                                internal::scalar_product_op<muSpectre::Real,
-                                  muSpectre::Quantity<T, U>>>
+                                internal::scalar_quotient_op<muSpectre::Real,  muSpectre::Quantity<T, U>>>
     {
       typedef muSpectre::Quantity<T,U> ReturnType;
     };
 
+    //----------------------------------------------------------------------------//
+    //! required for matrix multiplications with same units
+    template<>
+    template<typename U1>
+    struct ScalarBinaryOpTraits<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>,
+                                internal::scalar_quotient_op<muSpectre::Quantity<muSpectre::Real, U1>, muSpectre::Quantity<muSpectre::Real, U1>>>
+    {
+      typedef muSpectre::Quantity<muSpectre::Real,muSpectre::units::nodim> ReturnType;
+    };
+
+    //----------------------------------------------------------------------------//
+    //! required for matrix multiplications with differing units
+    template<>
+    template<typename T, typename U1, typename U2>
+    struct ScalarBinaryOpTraits<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>,
+                                internal::scalar_quotient_op<muSpectre::Quantity<T, U1>, muSpectre::Quantity<T, U2>>>
+    {
+      typedef muSpectre::Quantity<T,muSpectre::units::unit_minus<U1, U2>> ReturnType;
+    };
+
 
   }  // Eigen
 
 
 #endif /* UNITS_H */
diff --git a/src/common/voigt_conversion.hh b/src/common/voigt_conversion.hh
index e1f4ad2..34e076b 100644
--- a/src/common/voigt_conversion.hh
+++ b/src/common/voigt_conversion.hh
@@ -1,38 +1,64 @@
 /**
  * file   voigt_conversion.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   02 May 2017
  *
  * @brief  utilities to transform vector notation arrays into voigt notation
  *         arrays and vice-versa
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
+#include "common/common.hh"
 #include <Eigen/Dense>
 #include <unsupported/Eigen/CXX11/Tensor>
 
 namespace muSpectre {
 
-  
+  template<Dim_t dim>
+  class VoigtConversion
+  {
+  public:
+    VoigtConversion();
+    virtual ~VoigtConversion();
+
+    template<class Tens, class Voigt>
+    inline static void sym_fourth_to_voigt(const Tens & t, Voigt & v);
+    template<class Tens, class Voigt>
+    inline static void sym_fourth_to_vect(const Voigt & v, Tens & t);
+    
+
+  public:
+  const static Dim_t size;
+  // matrix of vector index I as function of tensor indices i,j
+  const static Dim_t mat[dim][dim];
+  // array of matrix indices ij as function of vector index I
+  const static Dim_t sym_mat[dim][dim];
+  // array of matrix indices ij as function of vector index I
+  const static Dim_t vec[dim * dim][2];
+  // factors to multiply the strain by for voigt notation
+  const static Real factors[(dim * (dim - 1) / 2 + dim)];
+
+  };
+
 
 }  // muSpectre
diff --git a/src/materials/CMakeLists.txt b/src/materials/CMakeLists.txt
index caef772..aab3453 100644
--- a/src/materials/CMakeLists.txt
+++ b/src/materials/CMakeLists.txt
@@ -1,9 +1,8 @@
 set (materials_SRC
   material.cc
   material_hyper_elastic.cc
   )
 
 add_library(materials ${materials_SRC})
 target_link_libraries(materials Eigen3::Eigen)
-message("CMAKE_SOURCE_DIR = ${CMAKE_SOURCE_DIR}")
 
diff --git a/tests/test_units.cc b/tests/test_units.cc
index 7860693..de2bb5e 100644
--- a/tests/test_units.cc
+++ b/tests/test_units.cc
@@ -1,94 +1,166 @@
 /**
  * file   test_units.cc
  *
  * @author Till Junge <till.junge@epfl.ch>
  *
  * @date   02 May 2017
  *
  * @brief  Test the compile-time unit enforcement
  *
  * @section LICENCE
  *
  * Copyright (C) 2017 Till Junge
  *
  * µSpectre is free software; you can redistribute it and/or
  * modify it under the terms of the GNU General Public License as
  * published by the Free Software Foundation, either version 3, or (at
  * your option) any later version.
  *
  * µSpectre is distributed in the hope that it will be useful, but
  * WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  * General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with GNU Emacs; see the file COPYING. If not, write to the
  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  * Boston, MA 02111-1307, USA.
  */
 
 #include "boost/test/unit_test.hpp"
 #include <Eigen/Dense>
 #include <chrono>
 
 #include "common/units.hh"
 
 namespace muSpectre{
 
   BOOST_AUTO_TEST_SUITE(units)
 
   BOOST_AUTO_TEST_CASE (unit_test) {
     auto length = 20._m;
     auto mass = 15._kg;
     BOOST_CHECK_EQUAL(length*mass,mass*length);
     BOOST_CHECK_NE(length, 2*length);
     BOOST_CHECK_EQUAL(1._kg*1._m/1._s/1._s, 1._N);
     BOOST_CHECK_EQUAL(1._N*1._m, 1_J);
     BOOST_CHECK_EQUAL(1_N/1._m/1._m, 1._Pa);
   }
 
-  BOOST_AUTO_TEST_CASE(eigen_compatibility) {
+  BOOST_AUTO_TEST_CASE(eigen_basic_compatibility) {
     Eigen::Matrix<Quantity<Real,units::nodim>, 2, 2> mat;
     mat = mat.Random();
     Eigen::Matrix<Quantity<Real,units::nodim>, 2, 2> mat2;
     mat2 = mat;
     BOOST_CHECK_EQUAL(mat2, mat);
     BOOST_CHECK_NE(mat2, mat*mat);
     BOOST_CHECK_EQUAL(mat+mat, mat*2);
     BOOST_CHECK_EQUAL(mat+mat, 2*mat);
     BOOST_CHECK_EQUAL((mat+mat)/2, mat);
   }
 
+  BOOST_AUTO_TEST_CASE(eigen_compatibility) {
+    Eigen::Matrix<Quantity<Real,units::m>, 2, 2> mat;
+    mat = mat.Random();
+    Eigen::Matrix<Quantity<Real,units::m>, 2, 2> mat2;
+    mat2 = mat;
+    Eigen::Matrix<Quantity<Real,units::si_area>, 2, 2> mat3;
+    mat3 = mat3.Random();
+    Eigen::Matrix<Quantity<Real,units::N>, 2, 2> mat4;
+    mat4 = mat4.Random();
+    Eigen::Matrix<Quantity<Real,units::kg>, 2, 2> mat5;
+    mat5 = mat5.Random();
+    Eigen::Matrix<Quantity<Real,units::si_acceleration>, 2, 2> mat6;
+    mat6 = mat6.Random();
+    BOOST_CHECK_EQUAL(mat2, mat);
+    BOOST_CHECK_NE(mat3, mat*mat);
+    BOOST_CHECK_EQUAL(mat+mat, mat*2);
+    BOOST_CHECK_EQUAL(mat+mat, 2*mat);
+    BOOST_CHECK_EQUAL((mat+mat)/2, mat);
+    BOOST_CHECK_NE(mat4, mat5*mat6);
+  }
 
-  template <typename T>
-  Real timed_multiplication(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& mat, Uint nb_rep = 100) {
-    using mtype =Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
-    Eigen::Matrix<decltype(T()*T()), Eigen::Dynamic, Eigen::Dynamic> result(mat.rows(), mat.cols());
-    Real tot_duration{0};
-    for (Uint rep = 0; rep < nb_rep; ++rep) {
-      mat = mtype::Random(mat.rows(), mat.cols());
-      auto start = std::chrono::high_resolution_clock::now();
-      result = mat * mat;
-      std::chrono::duration<Real> duration = std::chrono::high_resolution_clock::now() - start;
-      tot_duration += duration.count();
-    }
-    
-    return tot_duration;
+  BOOST_AUTO_TEST_CASE(eigen_dynamic) {
+    const auto dyn = Eigen::Dynamic;
+    Eigen::Matrix<Quantity<Real,units::nodim>, dyn, dyn> mat(2, 2);
+    mat = mat.Random(2, 2);
+    Eigen::Matrix<Quantity<Real,units::nodim>, dyn, dyn> mat2(2, 2);
+    mat2 = mat;
+    Eigen::Matrix<Quantity<Real,units::nodim>, dyn, dyn> mat3(2, 2);
+    mat3 = mat3.Random(2, 2);
+    Eigen::Matrix<Quantity<Real,units::nodim>, dyn, dyn> mat4(2, 2);
+    mat4 = mat4.Random(2, 2);
+    Eigen::Matrix<Quantity<Real,units::nodim>, dyn, dyn> mat5(2, 2);
+    mat5 = mat5.Random(2, 2);
+    Eigen::Matrix<Quantity<Real,units::nodim>, dyn, dyn> mat6(2, 2);
+    mat6 = mat6.Random(2, 2);
+    BOOST_CHECK_EQUAL(mat2, mat);
+    BOOST_CHECK_EQUAL(mat+mat, mat*2);
+    BOOST_CHECK_EQUAL(mat+mat, 2*mat);
+    BOOST_CHECK_EQUAL((mat+mat)/2, mat);
+    BOOST_CHECK_NE(mat3, mat*mat);
+    BOOST_CHECK_NE(mat4, mat5*mat6);
+
+  }
+
+  BOOST_AUTO_TEST_CASE(eigen_compatibility_dynamic) {
+    Quantity<Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>, units::m> mat(2, 2);
+    mat.val = mat.val.Random(2, 2);
+    Quantity<Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>, units::m> mat2(2, 2);
+    mat2.val = mat.val;
+    Quantity<Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>, units::si_area> mat3(2, 2);
+    mat3.val = mat3.val.Random(2, 2);
+    Quantity<Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>, units::N> mat4(2, 2);
+    mat4.val = mat4.val.Random(2, 2);
+    Quantity<Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>, units::kg> mat5(2, 2);
+    mat5.val = mat5.val.Random(2, 2);
+    Quantity<Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>, units::si_acceleration> mat6(2, 2);
+    mat6.val = mat6.val.Random(2, 2);
+    BOOST_CHECK_EQUAL(mat2, mat);
+    BOOST_CHECK_EQUAL(mat+mat, mat*2);
+    BOOST_CHECK_EQUAL(mat+mat, 2*mat);
+    BOOST_CHECK_EQUAL((mat+mat)/2, mat);
+    BOOST_CHECK_NE(mat3, mat*mat);
+    BOOST_CHECK_NE(mat4, mat5*mat6);
   }
-  
+
   BOOST_AUTO_TEST_CASE(eigen_speed) {
-    const size_t size{100};
+    const size_t size{20};
     using rtype = Eigen::MatrixXd;
-    using utype = Eigen::Matrix<Quantity<Real, units::m>, Eigen::Dynamic, Eigen::Dynamic>;
+    using utype = Quantity<Eigen::MatrixXd, units::m>;
     rtype rmat(size, size);
     utype umat(size, size);
 
+    Real r_duration = 0, u_duration = 0;
 
-    auto t_r = timed_multiplication(rmat);
-    auto t_u = timed_multiplication(umat);
+    rtype r_result(size, size);
+    Quantity<Eigen::MatrixXd, units::si_area> u_result(size, size);
+
+    const size_t nb_rep = 50;
+
+    std::chrono::time_point<std::chrono::high_resolution_clock> start;
+    std::chrono::duration<Real> duration;
+    for (Uint rep = 0; rep < nb_rep; ++rep) {
+      umat.val = umat.val.Random(size, size);
+      rmat = rmat.Random(size, size);
+      start = std::chrono::high_resolution_clock::now();
+      u_result = umat * umat;
+      duration = std::chrono::high_resolution_clock::now() - start;
+      u_duration += duration.count();
+
+      start = std::chrono::high_resolution_clock::now();
+      r_result = rmat * rmat;
+      duration = std::chrono::high_resolution_clock::now() - start;
+      r_duration += duration.count();
+    }
+
+    for (Uint rep = 0; rep < nb_rep; ++rep) {
+    }
 
-    BOOST_CHECK_EQUAL(t_r, t_u);
+    // just to sound the alarms in case the units make the execution time explode
+    // (should not have *any* runtime cost)
+    BOOST_CHECK_LT(u_duration, r_duration*1.1);
   }
 
   BOOST_AUTO_TEST_SUITE_END()
-}  // muSp
+}  // muSpectre