diff --git a/src/common/aka_random_generator.hh b/src/common/aka_random_generator.hh
index 816173a6a..3fada453f 100644
--- a/src/common/aka_random_generator.hh
+++ b/src/common/aka_random_generator.hh
@@ -1,261 +1,265 @@
 /**
  * @file   aka_random_generator.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Feb 21 2013
  * @date last modification: Wed Nov 11 2015
  *
  * @brief  generic random generator
  *
  * @section LICENSE
  *
  * Copyright  (©)  2014,  2015 EPFL  (Ecole Polytechnique  Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <random>
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_AKA_RANDOM_GENERATOR_HH__
 #define __AKANTU_AKA_RANDOM_GENERATOR_HH__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /* List of available distributions                                            */
 /* -------------------------------------------------------------------------- */
 // clang-format off
 #define AKANTU_RANDOM_DISTRIBUTION_TYPES                \
   ((uniform      , std::uniform_real_distribution ))    \
   ((exponential  , std::exponential_distribution  ))    \
   ((gamma        , std::gamma_distribution        ))    \
   ((weibull      , std::weibull_distribution      ))    \
   ((extreme_value, std::extreme_value_distribution))    \
   ((normal       , std::normal_distribution       ))    \
   ((lognormal    , std::lognormal_distribution    ))    \
   ((chi_squared  , std::chi_squared_distribution  ))    \
   ((cauchy       , std::cauchy_distribution       ))    \
   ((fisher_f     , std::fisher_f_distribution     ))    \
   ((student_t    , std::student_t_distribution    ))
 // clang-format on
 
 #define AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(elem) BOOST_PP_CAT(_rdt_, elem)
 #define AKANTU_RANDOM_DISTRIBUTION_PREFIX(s, data, elem)                       \
   AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(BOOST_PP_TUPLE_ELEM(2, 0, elem))
 
 enum RandomDistributionType {
   BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_RANDOM_DISTRIBUTION_PREFIX, _,
                                            AKANTU_RANDOM_DISTRIBUTION_TYPES)),
   _rdt_not_defined
 };
 
 /* -------------------------------------------------------------------------- */
 /* Generator                                                                  */
 /* -------------------------------------------------------------------------- */
 template <typename T> class RandomGenerator {
   /* ------------------------------------------------------------------------ */
 public:
   inline T operator()() { return generator(); }
 
   /// function to print the contain of the class
   virtual void printself(std::ostream & stream, int) const {
     stream << "RandGenerator [seed=" << _seed << "]";
   }
 
   /* ------------------------------------------------------------------------ */
 public:
   static void seed(long int s) {
     _seed = s;
     generator.seed(_seed);
   }
   static long int seed() { return _seed; }
 
   static constexpr T min() { return generator.min(); }
   static constexpr T max() { return generator.max(); }
 
   /* ------------------------------------------------------------------------ */
 private:
   static long int _seed;
   static std::default_random_engine generator;
 };
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 #undef AKANTU_RANDOM_DISTRIBUTION_PREFIX
 
 #define AKANTU_RANDOM_DISTRIBUTION_TYPE_PRINT_CASE(r, data, elem)              \
   case AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(                                \
       BOOST_PP_TUPLE_ELEM(2, 0, elem)): {                                      \
     stream << BOOST_PP_STRINGIZE(AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(      \
         BOOST_PP_TUPLE_ELEM(2, 0, elem)));                                     \
     break;                                                                     \
   }
 
 inline std::ostream & operator<<(std::ostream & stream,
                                  RandomDistributionType type) {
   switch (type) {
     BOOST_PP_SEQ_FOR_EACH(AKANTU_RANDOM_DISTRIBUTION_TYPE_PRINT_CASE, _,
                           AKANTU_RANDOM_DISTRIBUTION_TYPES)
   default:
     stream << UInt(type) << " not a RandomDistributionType";
     break;
   }
   return stream;
 }
 #undef AKANTU_RANDOM_DISTRIBUTION_TYPE_PRINT_CASE
 
 /* -------------------------------------------------------------------------- */
 /* Some Helper                                                                */
 /* -------------------------------------------------------------------------- */
 template <typename T, class Distribution> class RandomDistributionTypeHelper {
   enum { value = _rdt_not_defined };
 };
 
 /* -------------------------------------------------------------------------- */
 #define AKANTU_RANDOM_DISTRIBUTION_TYPE_GET_TYPE(r, data, elem)                \
   template <typename T>                                                        \
       struct RandomDistributionTypeHelper<T, BOOST_PP_TUPLE_ELEM(2, 1, elem) < \
                                                  T>> {                         \
     enum {                                                                     \
       value = AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(                         \
           BOOST_PP_TUPLE_ELEM(2, 0, elem))                                     \
     };                                                                         \
+                                                                               \
+    static void printself(std::ostream & stream) {                             \
+      stream << BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 0, elem));           \
+    }                                                                          \
   };
 
 BOOST_PP_SEQ_FOR_EACH(AKANTU_RANDOM_DISTRIBUTION_TYPE_GET_TYPE, _,
                       AKANTU_RANDOM_DISTRIBUTION_TYPES)
 
 #undef AKANTU_RANDOM_DISTRIBUTION_TYPE_GET_TYPE
 
 /* -------------------------------------------------------------------------- */
 template <class T> class RandomDistribution {
 public:
   virtual T operator()(RandomGenerator<UInt> & gen) = 0;
   virtual std::unique_ptr<RandomDistribution<T>> make_unique() const = 0;
   virtual void printself(std::ostream & stream, int = 0) const = 0;
 };
 
 template <class T, class Distribution>
 class RandomDistributionProxy : public RandomDistribution<T> {
 public:
   explicit RandomDistributionProxy(Distribution dist)
       : distribution(std::move(dist)) {}
 
   T operator()(RandomGenerator<UInt> & gen) override {
     return distribution(gen);
   }
 
   std::unique_ptr<RandomDistribution<T>> make_unique() const override {
     return std::make_unique<RandomDistributionProxy<T, Distribution>>(
         distribution);
   }
 
   void printself(std::ostream & stream, int = 0) const override {
-    stream << RandomDistributionTypeHelper<T, Distribution>::value << " [ "
-           << distribution << " ]";
+    RandomDistributionTypeHelper<T, Distribution>::printself(stream);
+    stream << " [ " << distribution << " ]";
   }
 
 private:
   Distribution distribution;
 };
 
 /* -------------------------------------------------------------------------- */
 /* RandomParameter                                                            */
 /* -------------------------------------------------------------------------- */
 template <typename T> class RandomParameter {
 public:
   template <class Distribution>
   explicit RandomParameter(T base_value, Distribution dist)
       : base_value(base_value),
         type(RandomDistributionType(
             RandomDistributionTypeHelper<T, Distribution>::value)),
         distribution_proxy(
             std::make_unique<RandomDistributionProxy<T, Distribution>>(
                 std::move(dist))) {}
 
   explicit RandomParameter(T base_value)
       : base_value(base_value),
         type(RandomDistributionType(
             RandomDistributionTypeHelper<
                 T, std::uniform_real_distribution<T>>::value)),
         distribution_proxy(
             std::make_unique<
                 RandomDistributionProxy<T, std::uniform_real_distribution<T>>>(
                 std::uniform_real_distribution<T>(0., 0.))) {}
 
   RandomParameter(const RandomParameter & other)
       : base_value(other.base_value), type(other.type),
         distribution_proxy(other.distribution_proxy->make_unique()) {}
 
   RandomParameter & operator=(const RandomParameter & other) {
     distribution_proxy = other.distribution_proxy->make_unique();
     base_value = other.base_value;
     type = other.type;
     return *this;
   }
 
   inline void setBaseValue(const T & value) { this->base_value = value; }
   inline T getBaseValue() const { return this->base_value; }
 
   template <template <typename> class Generator, class iterator>
   void setValues(iterator it, iterator end) {
     RandomGenerator<UInt> gen;
     for (; it != end; ++it)
       *it = this->base_value + (*distribution_proxy)(gen);
   }
 
   virtual void printself(std::ostream & stream,
                          __attribute__((unused)) int indent = 0) const {
     stream << base_value;
-    stream << " " << *distribution_proxy;
+    stream << " + " << *distribution_proxy;
   }
 
 private:
   /// Value with no random variations
   T base_value;
 
   /// Random distribution type
   RandomDistributionType type;
 
   /// Proxy to store a std random distribution
   std::unique_ptr<RandomDistribution<T>> distribution_proxy;
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline std::ostream & operator<<(std::ostream & stream,
                                  RandomDistribution<T> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline std::ostream & operator<<(std::ostream & stream,
                                  RandomParameter<T> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // akantu
 
 #endif /* __AKANTU_AKA_RANDOM_GENERATOR_HH__ */
diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh
index 34fdf8041..dbe1e1846 100644
--- a/src/common/aka_types.hh
+++ b/src/common/aka_types.hh
@@ -1,1236 +1,1250 @@
 /**
  * @file   aka_types.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Feb 17 2011
  * @date last modification: Fri Jan 22 2016
  *
  * @brief  description of the "simple" types
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  2015 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu 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  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_error.hh"
 #include "aka_fwd.hh"
 #include "aka_math.hh"
 /* -------------------------------------------------------------------------- */
 #include <initializer_list>
 #include <iomanip>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_AKA_TYPES_HH__
 #define __AKANTU_AKA_TYPES_HH__
 
 namespace akantu {
 
 enum NormType { L_1 = 1, L_2 = 2, L_inf = UInt(-1) };
 
 /**
  * DimHelper is a class to generalize the setup of a dim array from 3
  * values. This gives a common interface in the TensorStorage class
  * independently of its derived inheritance (Vector, Matrix, Tensor3)
  * @tparam dim
  */
 template <UInt dim> struct DimHelper {
   static inline void setDims(UInt m, UInt n, UInt p, UInt dims[dim]);
 };
 
 /* -------------------------------------------------------------------------- */
 template <> struct DimHelper<1> {
   static inline void setDims(UInt m, __attribute__((unused)) UInt n,
                              __attribute__((unused)) UInt p, UInt dims[1]) {
     dims[0] = m;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <> struct DimHelper<2> {
   static inline void setDims(UInt m, UInt n, __attribute__((unused)) UInt p,
                              UInt dims[2]) {
     dims[0] = m;
     dims[1] = n;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <> struct DimHelper<3> {
   static inline void setDims(UInt m, UInt n, UInt p, UInt dims[3]) {
     dims[0] = m;
     dims[1] = n;
     dims[2] = p;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T, UInt ndim, class RetType> class TensorStorage;
 
 /* -------------------------------------------------------------------------- */
 /* Proxy classes                                                              */
 /* -------------------------------------------------------------------------- */
 namespace tensors {
 template <class A, class B> struct is_copyable {
   enum : bool { value = false };
 };
 
 template <class A> struct is_copyable<A, A> {
   enum : bool { value = true };
 };
 
 template <class A> struct is_copyable<A, typename A::RetType> {
   enum : bool { value = true };
 };
 
 template <class A> struct is_copyable<A, typename A::RetType::proxy> {
   enum : bool { value = true };
 };
 
 } // namespace tensors
 /**
  * @class TensorProxy aka_types.hh
  * @desc The TensorProxy class is a proxy class to the TensorStorage it handles
  * the
  * wrapped case. That is to say if an accessor should give access to a Tensor
  * wrapped on some data, like the Array<T>::iterator they can return a
  * TensorProxy that will be automatically transformed as a TensorStorage wrapped
  * on the same data
  * @tparam T stored type
  * @tparam ndim order of the tensor
  * @tparam RetType real derived type
  */
 template <typename T, UInt ndim, class _RetType> class TensorProxy {
 protected:
   using RetTypeProxy = typename _RetType::proxy;
 
   TensorProxy(T * data, UInt m, UInt n, UInt p) {
     DimHelper<ndim>::setDims(m, n, p, this->n);
     this->values = data;
   }
 
   template <class Other, typename = std::enable_if_t<
                              tensors::is_copyable<TensorProxy, Other>::value>>
   explicit TensorProxy(const Other & other) {
     this->values = other.storage();
     for (UInt i = 0; i < ndim; ++i)
       this->n[i] = other.size(i);
   }
 
 public:
   using RetType = _RetType;
 
   operator RetType() { return RetType(static_cast<RetTypeProxy &>(*this)); }
 
   UInt size(UInt i) const {
     AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim
                                                           << " dimensions, not "
                                                           << (i + 1));
     return n[i];
   }
 
   inline UInt size() const {
     UInt _size = 1;
     for (UInt d = 0; d < ndim; ++d)
       _size *= this->n[d];
     return _size;
   }
 
   T * storage() const { return values; }
 
   template <class Other, typename = std::enable_if_t<
                              tensors::is_copyable<TensorProxy, Other>::value>>
   inline TensorProxy & operator=(const Other & other) {
     AKANTU_DEBUG_ASSERT(
         other.size() == this->size(),
         "You are trying to copy two tensors with different sizes");
     memcpy(this->values, other.storage(), this->size() * sizeof(T));
     return *this;
   }
 
   // template <class Other, typename = std::enable_if_t<
   //                          tensors::is_copyable<TensorProxy, Other>::value>>
   // inline TensorProxy & operator=(const Other && other) {
   //   AKANTU_DEBUG_ASSERT(
   //       other.size() == this->size(),
   //       "You are trying to copy two tensors with different sizes");
   //   memcpy(this->values, other.storage(), this->size() * sizeof(T));
   //   return *this;
   // }
 
   template <typename O> inline RetTypeProxy & operator*=(const O & o) {
     RetType(*this) *= o;
     return static_cast<RetTypeProxy &>(*this);
   }
 
   template <typename O> inline RetTypeProxy & operator/=(const O & o) {
     RetType(*this) /= o;
     return static_cast<RetTypeProxy &>(*this);
   }
 
 protected:
   T * values;
   UInt n[ndim];
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T> class VectorProxy : public TensorProxy<T, 1, Vector<T>> {
   using parent = TensorProxy<T, 1, Vector<T>>;
   using type = Vector<T>;
 
 public:
   VectorProxy(T * data, UInt n) : parent(data, n, 0, 0) {}
   template <class Other> explicit VectorProxy(Other & src) : parent(src) {}
 
   /* ---------------------------------------------------------------------- */
   template <class Other>
   inline VectorProxy<T> & operator=(const Other & other) {
     parent::operator=(other);
     return *this;
   }
 
   // inline VectorProxy<T> & operator=(const VectorProxy && other) {
   //   parent::operator=(other);
   //   return *this;
   // }
 
   /* ------------------------------------------------------------------------ */
   T & operator()(UInt index) { return this->values[index]; };
   const T & operator()(UInt index) const { return this->values[index]; };
 };
 
 template <typename T> class MatrixProxy : public TensorProxy<T, 2, Matrix<T>> {
   using parent = TensorProxy<T, 2, Matrix<T>>;
   using type = Matrix<T>;
 
 public:
   MatrixProxy(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
   template <class Other> explicit MatrixProxy(Other & src) : parent(src) {}
 
   /* ---------------------------------------------------------------------- */
   template <class Other>
   inline MatrixProxy<T> & operator=(const Other & other) {
     parent::operator=(other);
     return *this;
   }
 };
 
 template <typename T>
 class Tensor3Proxy : public TensorProxy<T, 3, Tensor3<T>> {
   typedef TensorProxy<T, 3, Tensor3<T>> parent;
   using type = Tensor3<T>;
 
 public:
   Tensor3Proxy(T * data, UInt m, UInt n, UInt k) : parent(data, m, n, k) {}
   Tensor3Proxy(const Tensor3Proxy & src) : parent(src) {}
   Tensor3Proxy(const Tensor3<T> & src) : parent(src) {}
 
   /* ---------------------------------------------------------------------- */
   template <class Other>
   inline Tensor3Proxy<T> & operator=(const Other & other) {
     parent::operator=(other);
     return *this;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Tensor base class                                                          */
 /* -------------------------------------------------------------------------- */
 template <typename T, UInt ndim, class RetType> class TensorStorage {
 public:
   using value_type = T;
 
 protected:
   template <class TensorType> void copySize(const TensorType & src) {
     for (UInt d = 0; d < ndim; ++d)
       this->n[d] = src.size(d);
     this->_size = src.size();
   }
 
   TensorStorage() : values(nullptr) {
     for (UInt d = 0; d < ndim; ++d)
       this->n[d] = 0;
     _size = 0;
   }
 
   TensorStorage(const TensorProxy<T, ndim, RetType> & proxy) {
     this->copySize(proxy);
     this->values = proxy.storage();
     this->wrapped = true;
   }
 
 protected:
   TensorStorage(const TensorStorage & src) = delete;
 
 public:
   TensorStorage(const TensorStorage & src, bool deep_copy)
       : values(nullptr), wrapped(false) {
     if (deep_copy)
       this->deepCopy(src);
     else
       this->shallowCopy(src);
   }
 
 protected:
   TensorStorage(UInt m, UInt n, UInt p, const T & def) {
     DimHelper<ndim>::setDims(m, n, p, this->n);
 
     this->computeSize();
     this->values = new T[this->_size];
     this->set(def);
     this->wrapped = false;
   }
 
   TensorStorage(T * data, UInt m, UInt n, UInt p) {
     DimHelper<ndim>::setDims(m, n, p, this->n);
 
     this->computeSize();
     this->values = data;
     this->wrapped = true;
   }
 
 public:
   /* ------------------------------------------------------------------------ */
   template <class TensorType> inline void shallowCopy(const TensorType & src) {
     this->copySize(src);
     if (!this->wrapped)
       delete[] this->values;
     this->values = src.storage();
     this->wrapped = true;
   }
 
   /* ------------------------------------------------------------------------ */
   template <class TensorType> inline void deepCopy(const TensorType & src) {
     this->copySize(src);
     if (!this->wrapped)
       delete[] this->values;
     this->values = new T[this->_size];
     memcpy((void *)this->values, (void *)src.storage(),
            this->_size * sizeof(T));
     this->wrapped = false;
   }
 
   virtual ~TensorStorage() {
     if (!this->wrapped)
       delete[] this->values;
   }
 
   /* ------------------------------------------------------------------------ */
   inline TensorStorage & operator=(const RetType & src) {
     if (this != &src) {
       if (this->wrapped) {
         // this test is not sufficient for Tensor of order higher than 1
         AKANTU_DEBUG_ASSERT(this->_size == src.size(),
                             "Tensors of different size");
         memcpy((void *)this->values, (void *)src.storage(),
                this->_size * sizeof(T));
       } else {
         deepCopy(src);
       }
     }
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   template <class R>
   inline RetType & operator+=(const TensorStorage<T, ndim, R> & other) {
     T * a = this->storage();
     T * b = other.storage();
     AKANTU_DEBUG_ASSERT(
         _size == other.size(),
         "The two tensors do not have the same size, they cannot be subtracted");
     for (UInt i = 0; i < _size; ++i)
       *(a++) += *(b++);
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   template <class R>
   inline RetType & operator-=(const TensorStorage<T, ndim, R> & other) {
     T * a = this->storage();
     T * b = other.storage();
     AKANTU_DEBUG_ASSERT(
         _size == other.size(),
         "The two tensors do not have the same size, they cannot be subtracted");
     for (UInt i = 0; i < _size; ++i)
       *(a++) -= *(b++);
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   inline RetType & operator+=(const T & x) {
     T * a = this->values;
     for (UInt i = 0; i < _size; ++i)
       *(a++) += x;
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   inline RetType & operator-=(const T & x) {
     T * a = this->values;
     for (UInt i = 0; i < _size; ++i)
       *(a++) -= x;
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   inline RetType & operator*=(const T & x) {
     T * a = this->storage();
     for (UInt i = 0; i < _size; ++i)
       *(a++) *= x;
     return *(static_cast<RetType *>(this));
   }
 
   /* ---------------------------------------------------------------------- */
   inline RetType & operator/=(const T & x) {
     T * a = this->values;
     for (UInt i = 0; i < _size; ++i)
       *(a++) /= x;
     return *(static_cast<RetType *>(this));
   }
 
   /// Y = \alpha X + Y
   inline RetType & aXplusY(const TensorStorage & other, const T & alpha = 1.) {
     AKANTU_DEBUG_ASSERT(
         _size == other.size(),
         "The two tensors do not have the same size, they cannot be subtracted");
 
     Math::aXplusY(this->_size, alpha, other.storage(), this->storage());
     return *(static_cast<RetType *>(this));
   }
 
   /* ------------------------------------------------------------------------ */
   T * storage() const { return values; }
   UInt size() const { return _size; }
   UInt size(UInt i) const {
     AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim
                                                           << " dimensions, not "
                                                           << (i + 1));
     return n[i];
   };
   /* ------------------------------------------------------------------------ */
   inline void clear() { memset(values, 0, _size * sizeof(T)); };
   inline void set(const T & t) { std::fill_n(values, _size, t); };
 
   template <class TensorType> inline void copy(const TensorType & other) {
     AKANTU_DEBUG_ASSERT(
         _size == other.size(),
         "The two tensors do not have the same size, they cannot be copied");
     memcpy(values, other.storage(), _size * sizeof(T));
   }
 
   bool isWrapped() const { return this->wrapped; }
 
 protected:
   friend class Array<T>;
 
   inline void computeSize() {
     _size = 1;
     for (UInt d = 0; d < ndim; ++d)
       _size *= this->n[d];
   }
 
 protected:
   template <typename R, NormType norm_type> struct NormHelper {
     template <class Ten> static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it)
         _norm += std::pow(std::abs(*it), norm_type);
       return std::pow(_norm, 1. / norm_type);
     }
   };
 
   template <typename R> struct NormHelper<R, L_1> {
     template <class Ten> static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it)
         _norm += std::abs(*it);
       return _norm;
     }
   };
 
   template <typename R> struct NormHelper<R, L_2> {
     template <class Ten> static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it)
         _norm += *it * *it;
       return sqrt(_norm);
     }
   };
 
   template <typename R> struct NormHelper<R, L_inf> {
     template <class Ten> static R norm(const Ten & ten) {
       R _norm = 0.;
       R * it = ten.storage();
       R * end = ten.storage() + ten.size();
       for (; it < end; ++it)
         _norm = std::max(std::abs(*it), _norm);
       return _norm;
     }
   };
 
 public:
   /*----------------------------------------------------------------------- */
   /// "Entrywise" norm norm<L_p> @f[ \|\boldsymbol{T}\|_p = \left(
   /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}}
   /// @f]
   template <NormType norm_type> inline T norm() const {
     return NormHelper<T, norm_type>::norm(*this);
   }
 
 protected:
   UInt n[ndim];
   UInt _size;
   T * values;
   bool wrapped{false};
 };
 
 // template <typename T, UInt ndim, class RetType>
 // inline TensorProxy<T, ndim, RetType>::TensorProxy(
 //     const TensorStorage<T, ndim, RetType> & other) {
 //   this->values = other.storage();
 //   for (UInt i = 0; i < ndim; ++i)
 //     this->n[i] = other.size(i);
 // }
 
 /* -------------------------------------------------------------------------- */
 /* Vector                                                                     */
 /* -------------------------------------------------------------------------- */
 template <typename T> class Vector : public TensorStorage<T, 1, Vector<T>> {
   typedef TensorStorage<T, 1, Vector<T>> parent;
 
 public:
   using value_type = typename parent::value_type;
   using proxy = VectorProxy<T>;
 
 public:
   Vector() : parent() {}
   explicit Vector(UInt n, const T & def = T()) : parent(n, 0, 0, def) {}
   Vector(T * data, UInt n) : parent(data, n, 0, 0) {}
   Vector(const Vector & src, bool deep_copy = true) : parent(src, deep_copy) {}
   Vector(const VectorProxy<T> & src) : parent(src) {}
 
   Vector(std::initializer_list<T> list) : parent(list.size(), 0, 0, T()) {
     UInt i = 0;
     for (auto val : list) {
       operator()(i++) = val;
     }
   }
 
 public:
   ~Vector() override = default;
 
   /* ------------------------------------------------------------------------ */
   inline Vector & operator=(const Vector & src) {
     parent::operator=(src);
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   inline T & operator()(UInt i) {
     AKANTU_DEBUG_ASSERT((i < this->n[0]),
                         "Access out of the vector! "
                             << "Index (" << i
                             << ") is out of the vector of size (" << this->n[0]
                             << ")");
     return *(this->values + i);
   }
 
   inline const T & operator()(UInt i) const {
     AKANTU_DEBUG_ASSERT((i < this->n[0]),
                         "Access out of the vector! "
                             << "Index (" << i
                             << ") is out of the vector of size (" << this->n[0]
                             << ")");
     return *(this->values + i);
   }
 
   inline T & operator[](UInt i) { return this->operator()(i); }
   inline const T & operator[](UInt i) const { return this->operator()(i); }
 
   /* ------------------------------------------------------------------------ */
   inline Vector<T> & operator*=(Real x) { return parent::operator*=(x); }
   inline Vector<T> & operator/=(Real x) { return parent::operator/=(x); }
   /* ------------------------------------------------------------------------ */
   inline Vector<T> & operator*=(const Vector<T> & vect) {
     T * a = this->storage();
     T * b = vect.storage();
     for (UInt i = 0; i < this->_size; ++i)
       *(a++) *= *(b++);
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   inline Real dot(const Vector<T> & vect) const {
     return Math::vectorDot(this->values, vect.storage(), this->_size);
   }
 
   /* ------------------------------------------------------------------------ */
   inline Real mean() const {
     Real mean = 0;
     T * a = this->storage();
     for (UInt i = 0; i < this->_size; ++i)
       mean += *(a++);
     return mean / this->_size;
   }
 
   /* ------------------------------------------------------------------------ */
   inline Vector & crossProduct(const Vector<T> & v1, const Vector<T> & v2) {
     AKANTU_DEBUG_ASSERT(this->size() == 3,
                         "crossProduct is only defined in 3D (n=" << this->size()
                                                                  << ")");
     AKANTU_DEBUG_ASSERT(
         this->size() == v1.size() && this->size() == v2.size(),
         "crossProduct is not a valid operation non matching size vectors");
     Math::vectorProduct3(v1.storage(), v2.storage(), this->values);
     return *this;
   }
 
   /* ------------------------------------------------------------------------ */
   inline void solve(const Matrix<T> & A, const Vector<T> & b) {
     AKANTU_DEBUG_ASSERT(
         this->size() == A.rows() && this->_size == A.cols(),
         "The size of the solution vector mismatches the size of the matrix");
     AKANTU_DEBUG_ASSERT(
         this->_size == b._size,
         "The rhs vector has a mismatch in size with the matrix");
     Math::solve(this->_size, A.storage(), this->values, b.storage());
   }
 
   /* ------------------------------------------------------------------------ */
   template <bool tr_A>
   inline void mul(const Matrix<T> & A, const Vector<T> & x, Real alpha = 1.0);
   /* ------------------------------------------------------------------------ */
   inline Real norm() const { return parent::template norm<L_2>(); }
 
   template <NormType nt> inline Real norm() const {
     return parent::template norm<nt>();
   }
 
   /* ------------------------------------------------------------------------ */
   inline void normalize() {
     Real n = norm();
     operator/=(n);
   }
 
   /* ------------------------------------------------------------------------ */
   /// norm of (*this - x)
   inline Real distance(const Vector<T> & y) const {
     Real * vx = this->values;
     Real * vy = y.storage();
     Real sum_2 = 0;
     for (UInt i = 0; i < this->_size; ++i, ++vx, ++vy)
       sum_2 += (*vx - *vy) * (*vx - *vy);
     return sqrt(sum_2);
   }
 
   /* ------------------------------------------------------------------------ */
   inline bool equal(const Vector<T> & v,
                     Real tolerance = Math::getTolerance()) const {
     T * a = this->storage();
     T * b = v.storage();
     UInt i = 0;
     while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance))
       ++i;
     return i == this->_size;
   }
 
   /* ------------------------------------------------------------------------ */
   inline short compare(const Vector<T> & v,
                        Real tolerance = Math::getTolerance()) const {
     T * a = this->storage();
     T * b = v.storage();
     for (UInt i(0); i < this->_size; ++i, ++a, ++b) {
       if (std::abs(*a - *b) > tolerance)
         return (((*a - *b) > tolerance) ? 1 : -1);
     }
     return 0;
   }
 
   /* ------------------------------------------------------------------------ */
   inline bool operator==(const Vector<T> & v) const { return equal(v); }
   inline bool operator!=(const Vector<T> & v) const { return !operator==(v); }
   inline bool operator<(const Vector<T> & v) const { return compare(v) == -1; }
   inline bool operator>(const Vector<T> & v) const { return compare(v) == 1; }
 
   /* ------------------------------------------------------------------------ */
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const {
     std::string space;
     for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
       ;
 
     stream << "[";
     for (UInt i = 0; i < this->_size; ++i) {
       if (i != 0)
         stream << ", ";
       stream << this->values[i];
     }
     stream << "]";
   }
 
   //  friend class ::akantu::Array<T>;
 };
 
 using RVector = Vector<Real>;
 
 /* ------------------------------------------------------------------------ */
 template <>
 inline bool Vector<UInt>::equal(const Vector<UInt> & v,
                                 __attribute__((unused)) Real tolerance) const {
   UInt * a = this->storage();
   UInt * b = v.storage();
   UInt i = 0;
   while (i < this->_size && (*(a++) == *(b++)))
     ++i;
   return i == this->_size;
 }
 
 /* ------------------------------------------------------------------------ */
 /* Matrix                                                                   */
 /* ------------------------------------------------------------------------ */
 template <typename T> class Matrix : public TensorStorage<T, 2, Matrix<T>> {
   typedef TensorStorage<T, 2, Matrix<T>> parent;
 
 public:
   using value_type = typename parent::value_type;
   using proxy = MatrixProxy<T>;
 
 public:
   Matrix() : parent() {}
   Matrix(UInt m, UInt n, const T & def = T()) : parent(m, n, 0, def) {}
   Matrix(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
   Matrix(const Matrix & src, bool deep_copy = true) : parent(src, deep_copy) {}
   Matrix(const MatrixProxy<T> & src) : parent(src) {}
   Matrix(std::initializer_list<std::initializer_list<T>> list) {
     std::size_t n = 0;
     std::size_t m = list.size();
     for (auto row : list) {
       n = std::max(n, row.size());
     }
 
     DimHelper<2>::setDims(m, n, 0, this->n);
     this->computeSize();
     this->values = new T[this->_size];
     this->set(0);
 
     UInt i = 0, j = 0;
     for (auto & row : list) {
       for (auto & val : row) {
         at(i, j++) = val;
       }
       ++i;
       j = 0;
     }
   }
 
   virtual ~Matrix() = default;
   /* ------------------------------------------------------------------------ */
   inline Matrix & operator=(const Matrix & src) {
     parent::operator=(src);
     return *this;
   }
 
 public:
   /* ---------------------------------------------------------------------- */
   UInt rows() const { return this->n[0]; }
   UInt cols() const { return this->n[1]; }
 
   /* ---------------------------------------------------------------------- */
   inline T & at(UInt i, UInt j) {
     AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])),
                         "Access out of the matrix! "
                             << "Index (" << i << ", " << j
                             << ") is out of the matrix of size (" << this->n[0]
                             << ", " << this->n[1] << ")");
     return *(this->values + i + j * this->n[0]);
   }
 
   inline const T & at(UInt i, UInt j) const {
     AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])),
                         "Access out of the matrix! "
                             << "Index (" << i << ", " << j
                             << ") is out of the matrix of size (" << this->n[0]
                             << ", " << this->n[1] << ")");
     return *(this->values + i + j * this->n[0]);
   }
 
   /* ------------------------------------------------------------------------ */
   inline T & operator()(UInt i, UInt j) { return this->at(i, j); }
   inline const T & operator()(UInt i, UInt j) const { return this->at(i, j); }
 
   /// give a line vector wrapped on the column i
   inline VectorProxy<T> operator()(UInt j) {
     AKANTU_DEBUG_ASSERT(j < this->n[1],
                         "Access out of the matrix! "
                             << "You are trying to access the column vector "
                             << j << " in a matrix of size (" << this->n[0]
                             << ", " << this->n[1] << ")");
     return VectorProxy<T>(this->values + j * this->n[0], this->n[0]);
   }
 
   inline const VectorProxy<T> operator()(UInt j) const {
     AKANTU_DEBUG_ASSERT(j < this->n[1],
                         "Access out of the matrix! "
                             << "You are trying to access the column vector "
                             << j << " in a matrix of size (" << this->n[0]
                             << ", " << this->n[1] << ")");
     return VectorProxy<T>(this->values + j * this->n[0], this->n[0]);
   }
 
   inline T & operator[](UInt idx) { return *(this->values + idx); };
   inline const T & operator[](UInt idx) const { return *(this->values + idx); };
 
   /* ---------------------------------------------------------------------- */
   inline Matrix operator*(const Matrix & B) {
     Matrix C(this->rows(), B.cols());
     C.mul<false, false>(*this, B);
     return C;
   }
 
   /* ----------------------------------------------------------------------- */
   inline Matrix & operator*=(const T & x) { return parent::operator*=(x); }
 
   inline Matrix & operator*=(const Matrix & B) {
     Matrix C(*this);
     this->mul<false, false>(C, B);
     return *this;
   }
 
   /* ---------------------------------------------------------------------- */
   template <bool tr_A, bool tr_B>
   inline void mul(const Matrix & A, const Matrix & B, T alpha = 1.0) {
     UInt k = A.cols();
     if (tr_A)
       k = A.rows();
 
 #ifndef AKANTU_NDEBUG
     if (tr_B) {
       AKANTU_DEBUG_ASSERT(k == B.cols(),
                           "matrices to multiply have no fit dimensions");
       AKANTU_DEBUG_ASSERT(this->cols() == B.rows(),
                           "matrices to multiply have no fit dimensions");
     } else {
       AKANTU_DEBUG_ASSERT(k == B.rows(),
                           "matrices to multiply have no fit dimensions");
       AKANTU_DEBUG_ASSERT(this->cols() == B.cols(),
                           "matrices to multiply have no fit dimensions");
     }
     if (tr_A) {
       AKANTU_DEBUG_ASSERT(this->rows() == A.cols(),
                           "matrices to multiply have no fit dimensions");
     } else {
       AKANTU_DEBUG_ASSERT(this->rows() == A.rows(),
                           "matrices to multiply have no fit dimensions");
     }
 #endif // AKANTU_NDEBUG
 
     Math::matMul<tr_A, tr_B>(this->rows(), this->cols(), k, alpha, A.storage(),
                              B.storage(), 0., this->storage());
   }
 
   /* ---------------------------------------------------------------------- */
   inline void outerProduct(const Vector<T> & A, const Vector<T> & B) {
     AKANTU_DEBUG_ASSERT(
         A.size() == this->rows() && B.size() == this->cols(),
         "A and B are not compatible with the size of the matrix");
     for (UInt i = 0; i < this->rows(); ++i) {
       for (UInt j = 0; j < this->cols(); ++j) {
         this->values[i + j * this->rows()] += A[i] * B[j];
       }
     }
   }
 
 private:
   class EigenSorter {
   public:
     EigenSorter(const Vector<T> & eigs) : eigs(eigs) {}
 
     bool operator()(const UInt & a, const UInt & b) const {
       return (eigs(a) > eigs(b));
     }
 
   private:
     const Vector<T> & eigs;
   };
 
 public:
   /* ---------------------------------------------------------------------- */
   inline void eig(Vector<T> & eigenvalues, Matrix<T> & eigenvectors) const {
     AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
                         "eig is not a valid operation on a rectangular matrix");
     AKANTU_DEBUG_ASSERT(eigenvalues.size() == this->cols(),
                         "eigenvalues should be of size " << this->cols()
                                                          << ".");
 #ifndef AKANTU_NDEBUG
     if (eigenvectors.storage() != NULL)
       AKANTU_DEBUG_ASSERT((eigenvectors.rows() == eigenvectors.cols()) &&
                               (eigenvectors.rows() == this->cols()),
                           "Eigenvectors needs to be a square matrix of size "
                               << this->cols() << " x " << this->cols() << ".");
 #endif
 
     Matrix<T> tmp = *this;
     Vector<T> tmp_eigs(eigenvalues.size());
     Matrix<T> tmp_eig_vects(eigenvectors.rows(), eigenvectors.cols());
 
     if (tmp_eig_vects.rows() == 0 || tmp_eig_vects.cols() == 0)
       Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage());
     else
       Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage(),
                       tmp_eig_vects.storage());
 
     Vector<UInt> perm(eigenvalues.size());
     for (UInt i = 0; i < perm.size(); ++i)
       perm(i) = i;
 
     std::sort(perm.storage(), perm.storage() + perm.size(),
               EigenSorter(tmp_eigs));
 
     for (UInt i = 0; i < perm.size(); ++i)
       eigenvalues(i) = tmp_eigs(perm(i));
 
     if (tmp_eig_vects.rows() != 0 && tmp_eig_vects.cols() != 0)
       for (UInt i = 0; i < perm.size(); ++i) {
         for (UInt j = 0; j < eigenvectors.rows(); ++j) {
           eigenvectors(j, i) = tmp_eig_vects(j, perm(i));
         }
       }
   }
 
   /* ---------------------------------------------------------------------- */
   inline void eig(Vector<T> & eigenvalues) const {
     Matrix<T> empty;
     eig(eigenvalues, empty);
   }
 
   /* ---------------------------------------------------------------------- */
   inline void eye(T alpha = 1.) {
     AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
                         "eye is not a valid operation on a rectangular matrix");
     this->clear();
     for (UInt i = 0; i < this->cols(); ++i) {
       this->values[i + i * this->rows()] = alpha;
     }
   }
 
   /* ---------------------------------------------------------------------- */
   static inline Matrix<T> eye(UInt m, T alpha = 1.) {
     Matrix<T> tmp(m, m);
     tmp.eye(alpha);
     return tmp;
   }
 
   /* ---------------------------------------------------------------------- */
   inline T trace() const {
     AKANTU_DEBUG_ASSERT(
         this->cols() == this->rows(),
         "trace is not a valid operation on a rectangular matrix");
     T trace = 0.;
     for (UInt i = 0; i < this->rows(); ++i) {
       trace += this->values[i + i * this->rows()];
     }
     return trace;
   }
 
   /* ---------------------------------------------------------------------- */
   inline Matrix transpose() const {
     Matrix tmp(this->cols(), this->rows());
     for (UInt i = 0; i < this->rows(); ++i) {
       for (UInt j = 0; j < this->cols(); ++j) {
         tmp(j, i) = operator()(i, j);
       }
     }
     return tmp;
   }
 
   /* ---------------------------------------------------------------------- */
   inline void inverse(const Matrix & A) {
     AKANTU_DEBUG_ASSERT(A.cols() == A.rows(),
                         "inv is not a valid operation on a rectangular matrix");
     AKANTU_DEBUG_ASSERT(this->cols() == A.cols(),
                         "the matrix should have the same size as its inverse");
 
     if (this->cols() == 1)
       *this->values = 1. / *A.storage();
     else if (this->cols() == 2)
       Math::inv2(A.storage(), this->values);
     else if (this->cols() == 3)
       Math::inv3(A.storage(), this->values);
     else
       Math::inv(this->cols(), A.storage(), this->values);
   }
 
   /* --------------------------------------------------------------------- */
   inline T det() const {
     AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
                         "inv is not a valid operation on a rectangular matrix");
     if (this->cols() == 1)
       return *(this->values);
     else if (this->cols() == 2)
       return Math::det2(this->values);
     else if (this->cols() == 3)
       return Math::det3(this->values);
     else
       return Math::det(this->cols(), this->values);
   }
 
   /* --------------------------------------------------------------------- */
   inline T doubleDot(const Matrix<T> & other) const {
     AKANTU_DEBUG_ASSERT(
         this->cols() == this->rows(),
         "doubleDot is not a valid operation on a rectangular matrix");
     if (this->cols() == 1)
       return *(this->values) * *(other.storage());
     else if (this->cols() == 2)
       return Math::matrixDoubleDot22(this->values, other.storage());
     else if (this->cols() == 3)
       return Math::matrixDoubleDot33(this->values, other.storage());
     else
       AKANTU_DEBUG_ERROR("doubleDot is not defined for other spatial dimensions"
                          << " than 1, 2 or 3.");
     return T();
   }
 
   /* ---------------------------------------------------------------------- */
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const {
     std::string space;
     for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
       ;
 
     stream << "[";
     for (UInt i = 0; i < this->n[0]; ++i) {
       if (i != 0)
         stream << ", ";
       stream << "[";
       for (UInt j = 0; j < this->n[1]; ++j) {
         if (j != 0)
           stream << ", ";
         stream << operator()(i, j);
       }
       stream << "]";
     }
     stream << "]";
   };
 };
 
 /* ------------------------------------------------------------------------ */
 template <typename T>
 template <bool tr_A>
 inline void Vector<T>::mul(const Matrix<T> & A, const Vector<T> & x,
                            Real alpha) {
 #ifndef AKANTU_NDEBUG
   UInt n = x.size();
   if (tr_A) {
     AKANTU_DEBUG_ASSERT(n == A.rows(),
                         "matrix and vector to multiply have no fit dimensions");
     AKANTU_DEBUG_ASSERT(this->size() == A.cols(),
                         "matrix and vector to multiply have no fit dimensions");
   } else {
     AKANTU_DEBUG_ASSERT(n == A.cols(),
                         "matrix and vector to multiply have no fit dimensions");
     AKANTU_DEBUG_ASSERT(this->size() == A.rows(),
                         "matrix and vector to multiply have no fit dimensions");
   }
 #endif
   Math::matVectMul<tr_A>(A.rows(), A.cols(), alpha, A.storage(), x.storage(),
                          0., this->storage());
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline std::ostream & operator<<(std::ostream & stream,
                                  const Matrix<T> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline std::ostream & operator<<(std::ostream & stream,
                                  const Vector<T> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* ------------------------------------------------------------------------ */
 /* Tensor3                                                                  */
 /* ------------------------------------------------------------------------ */
 template <typename T> class Tensor3 : public TensorStorage<T, 3, Tensor3<T>> {
   typedef TensorStorage<T, 3, Tensor3<T>> parent;
 
 public:
   using value_type = typename parent::value_type;
   using proxy = Tensor3Proxy<T>;
 
 public:
   Tensor3() : parent(){};
 
   Tensor3(UInt m, UInt n, UInt p, const T & def = T()) : parent(m, n, p, def) {}
 
   Tensor3(T * data, UInt m, UInt n, UInt p) : parent(data, m, n, p) {}
 
   Tensor3(const Tensor3 & src, bool deep_copy = true)
       : parent(src, deep_copy) {}
 
 public:
   /* ------------------------------------------------------------------------ */
   inline Tensor3 & operator=(const Tensor3 & src) {
     parent::operator=(src);
     return *this;
   }
 
   /* ---------------------------------------------------------------------- */
   inline T & operator()(UInt i, UInt j, UInt k) {
     AKANTU_DEBUG_ASSERT(
         (i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]),
         "Access out of the tensor3! "
             << "You are trying to access the element "
             << "(" << i << ", " << j << ", " << k << ") in a tensor of size ("
             << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")");
     return *(this->values + (k * this->n[0] + i) * this->n[1] + j);
   }
 
   inline const T & operator()(UInt i, UInt j, UInt k) const {
     AKANTU_DEBUG_ASSERT(
         (i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]),
         "Access out of the tensor3! "
             << "You are trying to access the element "
             << "(" << i << ", " << j << ", " << k << ") in a tensor of size ("
             << this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")");
     return *(this->values + (k * this->n[0] + i) * this->n[1] + j);
   }
 
   inline MatrixProxy<T> operator()(UInt k) {
     AKANTU_DEBUG_ASSERT((k < this->n[2]),
                         "Access out of the tensor3! "
                             << "You are trying to access the slice " << k
                             << " in a tensor3 of size (" << this->n[0] << ", "
                             << this->n[1] << ", " << this->n[2] << ")");
     return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
                           this->n[0], this->n[1]);
   }
 
   inline const MatrixProxy<T> operator()(UInt k) const {
     AKANTU_DEBUG_ASSERT((k < this->n[2]),
                         "Access out of the tensor3! "
                             << "You are trying to access the slice " << k
                             << " in a tensor3 of size (" << this->n[0] << ", "
                             << this->n[1] << ", " << this->n[2] << ")");
     return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
                           this->n[0], this->n[1]);
   }
 
   inline MatrixProxy<T> operator[](UInt k) {
     return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
                           this->n[0], this->n[1]);
   }
 
   inline const MatrixProxy<T> operator[](UInt k) const {
     return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
                           this->n[0], this->n[1]);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 // support operations for the creation of other vectors
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Vector<T> operator*(const T & scalar, const Vector<T> & a) {
   Vector<T> r(a);
   r *= scalar;
   return r;
 }
 
 template <typename T>
 Vector<T> operator*(const Vector<T> & a, const T & scalar) {
   Vector<T> r(a);
   r *= scalar;
   return r;
 }
 
 template <typename T>
 Vector<T> operator/(const Vector<T> & a, const T & scalar) {
   Vector<T> r(a);
   r /= scalar;
   return r;
 }
 
+template <typename T>
+Vector<T> operator*(const Vector<T> & a, const Vector<T> & b) {
+  Vector<T> r(a);
+  r *= b;
+  return r;
+}
+
 template <typename T>
 Vector<T> operator+(const Vector<T> & a, const Vector<T> & b) {
   Vector<T> r(a);
   r += b;
   return r;
 }
 
 template <typename T>
 Vector<T> operator-(const Vector<T> & a, const Vector<T> & b) {
   Vector<T> r(a);
   r -= b;
   return r;
 }
 
+template <typename T>
+Vector<T> operator*(const Matrix<T> & A, const Vector<T> & b) {
+  Vector<T> r(b.size());
+  r.template mul<false>(A, b);
+  return r;
+}
+
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Matrix<T> operator*(const T & scalar, const Matrix<T> & a) {
   Matrix<T> r(a);
   r *= scalar;
   return r;
 }
 
 template <typename T>
 Matrix<T> operator*(const Matrix<T> & a, const T & scalar) {
   Matrix<T> r(a);
   r *= scalar;
   return r;
 }
 
 template <typename T>
 Matrix<T> operator/(const Matrix<T> & a, const T & scalar) {
   Matrix<T> r(a);
   r /= scalar;
   return r;
 }
 
 template <typename T>
 Matrix<T> operator+(const Matrix<T> & a, const Matrix<T> & b) {
   Matrix<T> r(a);
   r += b;
   return r;
 }
 
 template <typename T>
 Matrix<T> operator-(const Matrix<T> & a, const Matrix<T> & b) {
   Matrix<T> r(a);
   r -= b;
   return r;
 }
 
 } // akantu
 
 #endif /* __AKANTU_AKA_TYPES_HH__ */