diff --git a/src/core/static_types.hh b/src/core/static_types.hh
index 303093f..2aad759 100644
--- a/src/core/static_types.hh
+++ b/src/core/static_types.hh
@@ -1,328 +1,328 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __STATIC_TYPES_HH__
 #define __STATIC_TYPES_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 
 __BEGIN_TAMAAS__
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 template <UInt acc, UInt n, UInt... ns>
 struct product_tail_rec : product_tail_rec<acc * n, ns...> {};
 
 template <UInt acc, UInt n>
 struct product_tail_rec<acc, n> : std::integral_constant<UInt, acc * n> {};
 
 template <UInt N, UInt n, UInt... ns>
 struct get_rec : get_rec<N - 1, ns...> {};
 
 template <UInt n, UInt... ns>
 struct get_rec<0, n, ns...> : std::integral_constant<UInt, n> {};
 }
 
 template <UInt... ns>
 struct product : detail::product_tail_rec<1, ns...> {};
 
 template <UInt N, UInt... ns>
 struct get : detail::get_rec<N, ns...> {};
 
 template <typename T>
 struct is_arithmetic : std::is_arithmetic<T> {};
 
 template <typename T>
 struct is_arithmetic<thrust::complex<T>> : std::true_type {};
 
 /**
  * @brief Static Array
  *
  * This class is meant to be a small and fast object for intermediate
  * calculations, possibly on wrapped memory belonging to a grid. Support type
  * show be either a pointer or a C array. It should not contain any virtual
  * method.
  */
 template <typename DataType, typename SupportType, UInt size>
 class StaticArray {
   static_assert(std::is_array<SupportType>::value ||
                     std::is_pointer<SupportType>::value,
                 "the support type of StaticArray should be either a pointer or "
                 "a C-array");
 
   using T = DataType;
 
 public:
   /// Access operator
   __device__ __host__ T& operator()(UInt i) {
     // TAMAAS_ASSERT(i < n, "Access out of bounds");
     return _mem[i];
   }
 
   /// Access operator
   __device__ __host__ const T& operator()(UInt i) const {
     // TAMAAS_ASSERT(i < n, "Access out of bounds");
     return _mem[i];
   }
 
   /// L2 norm squared
   __device__ __host__ T l2squared() const {
     T res = 0;
     for (UInt i = 0; i < size; ++i)
-      res += _mem[i] * _mem[i];
+      res += (*this)(i) * (*this)(i);
     return res;
   }
 
   /// L2 norm
   __device__ __host__ T l2norm() const { return std::sqrt(l2squared()); }
 
 #define VECTOR_OP(op)                                                          \
   template <typename DT, typename ST>                                          \
   __device__ __host__ void operator op(const StaticArray<DT, ST, size>& o) {   \
     for (UInt i = 0; i < size; ++i)                                            \
       (*this)(i) op o(i);                                                      \
   }
 
   VECTOR_OP(+=)
   VECTOR_OP(-=)
   VECTOR_OP(*=)
   VECTOR_OP(/=)
 #undef VECTOR_OP
 
 #define SCALAR_OP(op)                                                          \
   template <typename T1>                                                       \
   __device__ __host__                                                          \
       typename std::enable_if<is_arithmetic<T1>::value, StaticArray&>::type    \
       operator op(const T1& x) {                                               \
     for (UInt i = 0; i < size; ++i)                                            \
       (*this)(i) op x;                                                         \
     return *this;                                                              \
   }
 
   SCALAR_OP(+=)
   SCALAR_OP(-=)
   SCALAR_OP(*=)
   SCALAR_OP(/=)
   SCALAR_OP(=)
 #undef SCALAR_OP
 
   /// Overriding the implit copy operator
   StaticArray& operator=(const StaticArray& o) { return this->copy(o); }
 
   template <typename DT, typename ST>
   StaticArray& copy(const StaticArray<DT, ST, size>& o) {
     for (UInt i = 0; i < size; ++i)
       (*this)(i) = o(i);
     return *this;
   }
 
 protected:
   SupportType _mem;
 };
 
 /**
  * @brief Static Tensor
  *
  * This class implements a multi-dimensional tensor behavior.
  */
 template <typename DataType, typename SupportType = DataType*, UInt... dims>
 class StaticTensor
     : public StaticArray<DataType, SupportType, product<dims...>::value> {
   using parent = StaticArray<DataType, SupportType, product<dims...>::value>;
   using T = DataType;
 
 public:
   static constexpr UInt dim = sizeof...(dims);
   using parent::operator=;
 
 private:
   template <typename... Idx>
   static UInt unpackOffset(UInt offset, UInt index, Idx... rest) {
     constexpr UInt size = sizeof...(rest);
     offset += index;
     offset *= get<dim - size, dims...>::value;
     return unpackOffset(offset, rest...);
   }
 
   template <typename... Idx>
   static UInt unpackOffset(UInt offset, UInt index) {
     return offset + index;
   }
 
 public:
   template <typename... Idx>
   const T& operator()(Idx... idx) const {
-    return this->_mem[unpackOffset(0, idx...)];
+    return parent::operator()(unpackOffset(0, idx...));
   }
 
   template <typename... Idx>
   T& operator()(Idx... idx) {
-    return this->_mem[unpackOffset(0, idx...)];
+    return parent::operator()(unpackOffset(0, idx...));
   }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Common Static Types */
 /* -------------------------------------------------------------------------- */
 template <typename DataType, typename SupportType, UInt n, UInt m>
 using StaticMatrix = StaticTensor<DataType, SupportType, n, m>;
 
 /// Vector class with size determined at compile-time
 template <typename DataType, typename SupportType, UInt n>
 class StaticVector : public StaticTensor<DataType, SupportType, n> {
   using T = DataType;
 
 public:
   using StaticTensor<DataType, SupportType, n>::operator=;
   /// Matrix-vector product
   template <typename DT1, typename ST1, typename DT2, typename ST2, UInt m>
   __device__ __host__ void mul(const StaticMatrix<DT1, ST1, n, m>& mat,
                                const StaticVector<DT2, ST2, m>& vec) {
     *this = T(0);
     for (UInt i = 0; i < n; ++i)
       for (UInt j = 0; j < m; ++j)
         (*this)(i) += mat(i, j) * vec(j);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 /* Proxy Static Types */
 /* -------------------------------------------------------------------------- */
 
 /// Proxy type for tensor
 template <template <typename, typename, UInt...> class StaticParent, typename T,
           UInt... dims>
 class TensorProxy : public StaticParent<T, T*, dims...> {
 public:
   /// Implicit wrap constructor
   __device__ __host__ TensorProxy(T& spot) { this->_mem = &spot; }
 
   using StaticTensor<T, T*, dims...>::operator=;
 };
 
 template <typename T, UInt n, UInt m>
 using MatrixProxy = TensorProxy<StaticMatrix, T, n, m>;
 
 template <typename T, UInt n>
 using VectorProxy = TensorProxy<StaticVector, T, n>;
 
 /* -------------------------------------------------------------------------- */
 /* On the stack static types */
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, typename, UInt...> class StaticParent, typename T,
           UInt... dims>
 class Tensor : public StaticParent<T, T[product<dims...>::value], dims...> {
   using parent = StaticParent<T, T[product<dims...>::value], dims...>;
 
 public:
   using parent::operator=;
   using parent::copy;
 
   /// Default constructor
   Tensor() = default;
   /// Construct with default value
   Tensor(T val) { *this = val; }
 
   /// Construct by copy from static tensor
   template <typename DT, typename ST>
   Tensor(const StaticParent<DT, ST, dims...>& o) {
     this->copy(o);
   }
 };
 
 template <typename T, UInt n, UInt m>
 using Matrix = Tensor<StaticMatrix, T, n, m>;
 
 template <typename T, UInt n>
 using Vector = Tensor<StaticVector, T, n>;
 
 /* -------------------------------------------------------------------------- */
 /* Arithmetic operators creating temporaries */
 /* -------------------------------------------------------------------------- */
 
 // template <template <typename, typename, UInt...> class StaticParent,
 //           typename DT1, typename ST1, typename DT2, typename ST2, UInt...
 //           dims>
 // Tensor<StaticParent, decltype(DT1(0) + DT2(0)), dims...>
 // operator+(const StaticParent<DT1, ST1, dims...>& t1,
 //           const StaticParent<DT2, ST2, dims...>& t2) {
 //   Tensor<StaticParent, decltype(DT1(0) + DT2(0)), dims...> a(t1);
 //   return a += t2;
 // }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
 Vector<decltype(DT1(0) + DT2(0)), dim>
 operator+(const StaticVector<DT1, ST1, dim>& a,
           const StaticVector<DT2, ST2, dim>& b) {
   Vector<decltype(DT1(0) + DT2(0)), dim> res(a);
   res += b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt dim>
 Vector<decltype(DT1(0) - DT2(0)), dim>
 operator-(const StaticVector<DT1, ST1, dim>& a,
           const StaticVector<DT2, ST2, dim>& b) {
   Vector<decltype(DT1(0) - DT2(0)), dim> res(a);
   res -= b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 Matrix<decltype(DT1(0) + DT2(0)), n, m>
 operator+(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticMatrix<DT2, ST2, n, m>& b) {
   Matrix<decltype(DT1(0) + DT2(0)), n, m> res(a);
   res += b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 Matrix<decltype(DT1(0) - DT2(0)), n, m>
 operator-(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticMatrix<DT2, ST2, n, m>& b) {
   Matrix<decltype(DT1(0) - DT2(0)), n, m> res(a);
   res -= b;
   return res;
 }
 
 template <typename DT1, typename ST1, typename DT2, typename ST2, UInt n,
           UInt m>
 Vector<decltype(DT1(0) * DT2(0)), n>
 operator*(const StaticMatrix<DT1, ST1, n, m>& a,
           const StaticVector<DT2, ST2, m>& b) {
   Vector<decltype(DT1(0) * DT2(0)), n> res;
   res.mul(a, b);
   return res;
 }
 
 __END_TAMAAS__
 
 #endif  // __STATIC_TYPES_HH__
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index 71fff04..f29d8a5 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,221 +1,219 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "westergaard.hh"
 #include "fft_plan_manager.hh"
 #include "influence.hh"
 #include "loop.hh"
 #include "model.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 Westergaard<type>::Westergaard(Model* model) : BEEngine(model) {
   // Copy sizes
   std::array<UInt, trait::dimension> sizes;
   std::copy(model->getDiscretization().begin(),
             model->getDiscretization().end(), sizes.begin());
   auto hermitian_sizes =
       GridHermitian<Real, trait::dimension>::hermitianDimensions(sizes);
 
   // Copy boundary sizes
   auto sizes_it = sizes.begin();
   std::array<UInt, trait::boundary_dimension> boundary_hermitian_sizes;
   // Ignoring first dimension if model is volumetric
   if (trait::dimension > trait::boundary_dimension)
     ++sizes_it;
   std::copy(sizes_it, sizes.end(), boundary_hermitian_sizes.begin());
   boundary_hermitian_sizes =
       GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
           boundary_hermitian_sizes);
 
   constexpr UInt nb_components = trait::components;
 
   buffer.setNbComponents(nb_components);
   buffer.resize(boundary_hermitian_sizes);
   influence.setNbComponents(nb_components * nb_components);
   influence.resize(hermitian_sizes);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::initInfluence() {
   TAMAAS_EXCEPTION("Not implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define INFLUENCE_BASIC_MACRO                                                  \
   do {                                                                         \
     constexpr UInt bdim = trait::boundary_dimension;                           \
     auto wavevectors =                                                         \
         FFTransform<Real, bdim>::computeFrequencies<true>(influence.sizes());  \
     Real E_star = model->getHertzModulus();                                    \
-    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);          \
+    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);           \
     Loop::stridedLoop(influence::Basic<bdim>(E_star, domain), wavevectors,     \
                       influence);                                              \
     influence(0) = 0;                                                          \
   } while (0)
 
 template <>
 void Westergaard<model_type::basic_2d>::initInfluence() {
   INFLUENCE_BASIC_MACRO;
 }
 
 template <>
 void Westergaard<model_type::basic_1d>::initInfluence() {
   INFLUENCE_BASIC_MACRO;
 }
 
 #undef INFLUENCE_BASIC_MACRO
 
 #define INFLUENCE_SURFACE_MACRO                                                \
   do {                                                                         \
     constexpr UInt bdim = trait::boundary_dimension;                           \
     auto wavevectors =                                                         \
         FFTransform<Real, bdim>::computeFrequencies<true>(influence.sizes());  \
     Real E = model->getYoungModulus();                                         \
     Real nu = model->getPoissonRatio();                                        \
-    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);          \
+    VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);           \
     Loop::stridedLoop(influence::Surface<bdim>(E, nu, domain), wavevectors,    \
                       influence);                                              \
+    MatrixProxy<Complex, bdim + 1, bdim + 1> inf_0(influence(0));              \
+    inf_0 = 0;                                                                 \
   } while (0)
 
 template <>
 void Westergaard<model_type::surface_1d>::initInfluence() {
   INFLUENCE_SURFACE_MACRO;
 }
 
 template <>
 void Westergaard<model_type::surface_2d>::initInfluence() {
   INFLUENCE_SURFACE_MACRO;
 }
 
 #undef INFLUENCE_SURFACE_MACRO
 
 /* ------------------------------------------------------------------------ */
 
 #define NEUMANN_BASIC(type)                                                    \
   template <>                                                                  \
   void Westergaard<type>::solveNeumann(GridBase<Real>& in,                     \
                                        GridBase<Real>& out) const {            \
     auto apply = [](decltype(buffer)& buffer,                                  \
                     const decltype(influence)& influence) {                    \
       buffer *= influence;                                                     \
       buffer(0) = 0;                                                           \
     };                                                                         \
                                                                                \
     fourierApply(apply, in, out);                                              \
   }
 
 #define DIRICHLET_BASIC(type)                                                  \
   template <>                                                                  \
   void Westergaard<type>::solveDirichlet(GridBase<Real>& in,                   \
                                          GridBase<Real>& out) const {          \
     auto apply = [](decltype(buffer)& buffer,                                  \
                     const decltype(influence)& influence) {                    \
       buffer /= influence;                                                     \
       buffer(0) = 0;                                                           \
     };                                                                         \
                                                                                \
     fourierApply(apply, in, out);                                              \
   }
 
 NEUMANN_BASIC(model_type::basic_1d);
 NEUMANN_BASIC(model_type::basic_2d);
 DIRICHLET_BASIC(model_type::basic_1d);
 DIRICHLET_BASIC(model_type::basic_2d);
 
 #undef NEUMANN_BASIC
 #undef DIRICHLET_BASIC
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 /// Class needed for cuda
 template <UInt d>
 struct SurfaceApply {
   void operator()(GridHermitian<Real, d>& buffer_,
                   const GridHermitian<Real, d>& influence) {
     Loop::stridedLoop(
         [] CUDA_LAMBDA(VectorProxy<Complex, d + 1> && buf,
                        MatrixProxy<const Complex, d + 1, d + 1> && inf) {
-          Vector<Complex, d + 1> copy(buf);
-          buf.mul(inf, copy);
+	  buf.copy(inf * buf);
         },
         buffer_, influence);
-
-    VectorProxy<Complex, d + 1> buff_0(buffer_(0));
-    buff_0 = 0;
   }
 };
 }  // namespace detail
 
 #define NEUMANN_SURFACE(type)                                                  \
   template <>                                                                  \
   void Westergaard<type>::solveNeumann(GridBase<Real>& in,                     \
                                        GridBase<Real>& out) const {            \
     fourierApply(detail::SurfaceApply<trait::boundary_dimension>(), in, out);  \
   }
 
 /// \TODO Find a way to make it work with cuda
 NEUMANN_SURFACE(model_type::surface_1d);
 NEUMANN_SURFACE(model_type::surface_2d);
 
 #undef NEUMANN_SURFACE
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::solveNeumann(GridBase<Real>& /*neumann*/,
                                      GridBase<Real>& /*dirichlet*/) const {
   TAMAAS_EXCEPTION("Not yet implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::solveDirichlet(GridBase<Real>& /*neumann*/,
                                        GridBase<Real>& /*dirichlet*/) const {
   TAMAAS_EXCEPTION("Not yet implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template class Westergaard<model_type::basic_1d>;
 template class Westergaard<model_type::basic_2d>;
 template class Westergaard<model_type::surface_1d>;
 template class Westergaard<model_type::surface_2d>;
 template class Westergaard<model_type::volume_1d>;
 template class Westergaard<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
diff --git a/tests/test_loop.cpp b/tests/test_loop.cpp
index df1a9d4..4124a30 100644
--- a/tests/test_loop.cpp
+++ b/tests/test_loop.cpp
@@ -1,255 +1,254 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "test.hh"
 #include "grid.hh"
 #include "static_types.hh"
 #include "loop.hh"
 /* -------------------------------------------------------------------------- */
 /* WARNING: here we cannot use lambdas for tests because GoogleTest declares  */
 /* test functions as private members of classes which is incompatible with    */
 /* cuda's extended lambdas. I know... it's f*cking stupid                     */
 /* -------------------------------------------------------------------------- */
 
 using namespace tamaas;
 
 template <typename T>
 struct AddOneInplace {
   CUDA_LAMBDA void operator()(T& x) { x += 1; }
 };
 
 // Testing loops on one grid
 TEST(TestLoops, OneArgument) {
   Grid<Real, 1> grid({20}, 1);
   Grid<Real, 1> solution({20}, 1);
 
   auto add_one = [] (Real & x) { return x + 1; };
 
   std::iota(grid.begin(), grid.end(), 1);
 
   // Makeing solution
   std::transform(grid.begin(), grid.end(), solution.begin(), add_one);
 
   auto add_one_inplace = AddOneInplace<Real>();
 
   Loop::loop(add_one_inplace, grid);
 
   ASSERT_TRUE(compare(grid, solution, AreFloatEqual())) << "One argument loop failed";
 }
 
 struct PrimalTest {
   CUDA_LAMBDA void operator()(Int& primal, Int& val) {
     val = (primal > 0)? -1:1;
   }
 };
 
 // Testing loops on two grids
 TEST(TestLoops, TwoArguments) {
   // Why no ints?
   Grid<Int, 2> grid({20, 20}, 1);
   Grid<Int, 2> primal({20, 20}, 1);
   Grid<Int, 2> solution({20, 20}, 1);
 
   primal(0, 0) = 1; primal(0, 1) = 1;
   primal(1, 0) = 1; primal(1, 1) = 1;
 
   std::transform(primal.begin(), primal.end(), solution.begin(),
 		 [] (Int & primal) {
 		   return (primal > 0)? -1:1;
 		 });
 
   auto primal_test = PrimalTest();
 
   Loop::loop(primal_test, primal, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "Two argument loop failed";
 }
 
 struct AssignUInt {
   CUDA_LAMBDA void operator()(UInt& x, UInt i) {
     x = i;
   }
 };
 
 // Testing an enumeration
 TEST(TestLoops, Enumeration) {
   Grid<UInt, 1> grid({100}, 1);
   Grid<UInt, 1> solution({100}, 1);
   std::iota(solution.begin(), solution.end(), 0);
 
   auto assign_uint = AssignUInt();
 
   Loop::loop(assign_uint, grid, Loop::arange<UInt>(100));
 
   ASSERT_TRUE(compare(solution, grid)) << "Enumeration loop failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 struct Identity {
   CUDA_LAMBDA UInt operator()(UInt& x) const { return x; }
 };
 
 // Testing one grid reductions
 TEST(TestReductions, OneArgument) {
   Grid<UInt, 1> grid({6}, 1);
   std::iota(grid.begin(), grid.end(), 1);
 
   const auto id = Identity();
 
   // Sum reduction
   UInt sol = std::accumulate(grid.begin(), grid.end(), 0, std::plus<UInt>());
   UInt red = Loop::reduce<operation::plus>(id, grid);
   ASSERT_TRUE(sol == red) << "Addition reduction failed on one argument";
 
   // Product reduction
   sol = std::accumulate(grid.begin(), grid.end(), 1, std::multiplies<UInt>());
   red = Loop::reduce<operation::times>(id, grid);
   ASSERT_TRUE(sol == red) << "Multiplication reduction failed on one argument";
 
   // Min reduction
   sol = *std::min_element(grid.begin(), grid.end());
   red = Loop::reduce<operation::min>(id, grid);
   ASSERT_TRUE(sol == red) << "Min reduction failed on one argument";
 
   // Max reduction
   sol = *std::max_element(grid.begin(), grid.end());
   red = Loop::reduce<operation::max>(id, grid);
   ASSERT_TRUE(sol == red) << "Max reduction failed on one argument";
 }
 
 struct PrimalReduce {
   CUDA_LAMBDA UInt operator()(UInt& p, UInt& val) {
     return (p > 0)? val:0;
   }
 };
 
 TEST(TestReductions, TwoArguments) {
   Grid<UInt, 1> grid({20}, 1);
   Grid<UInt, 1> primal({20}, 1);
   grid = 1;
   primal(0) = 1;
   primal(1) = 1;
 
   auto primal_reduce = PrimalReduce();
 
   // Reduce on values where primal > 0
   UInt red = Loop::reduce<operation::plus>(primal_reduce, primal, grid);
   ASSERT_TRUE(red == 2) << "Two args reduction failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 TEST(TestStridedLoops, OneArgument) {
   Grid<UInt, 2> grid({10, 10}, 2);
   std::iota(grid.begin(), grid.end(), 1);
 
   Grid<UInt, 2> solution({10, 10}, 2);
   solution = grid;
 
   std::for_each(solution.begin(), solution.end(), [] (UInt & x) {
       if (x % 2 == 1)
 	x += 1;
     });
 
   auto add_one_inplace = AddOneInplace<UInt>();
 
   Loop::stridedLoop(add_one_inplace, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "One argument strided loop failed";
 }
 
 template <typename T> using WrapVector = VectorProxy<T, 2>;
 
 struct AddOneVector {
   CUDA_LAMBDA void operator()(WrapVector<UInt>&& x) {
     x(0) += 1;
   }
 };
 
 TEST(TestStridedLoops, VectorStride) {
   Grid<UInt, 2> grid({10, 10}, 2);
   std::iota(grid.begin(), grid.end(), 1);
 
   Grid<UInt, 2> solution({10, 10}, 2);
   solution = grid;
 
   std::for_each(solution.begin(), solution.end(), [] (UInt & x) {
       if (x % 2 == 1)
 	x += 1;
     });
 
   auto add_one_inplace = AddOneVector();
 
   Loop::stridedLoop(add_one_inplace, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "Static vector strided loop failed";
 }
 
 template <typename T> using WrapMatrix = MatrixProxy<T, 2, 2>;
 
 struct SetOneMatrix {
   CUDA_LAMBDA void operator()(WrapMatrix<UInt>&& x) {
     x(0, 0) = 1;
     x(1, 1) = 1;
   }
 };
 
 TEST(TestStridedLoops, MatrixStride) {
   Grid<UInt, 2> grid({10, 10}, 4);
   Grid<UInt, 2> solution({10, 10}, 4);
   std::iota(solution.begin(), solution.end(), 0);
 
   std::for_each(solution.begin(), solution.end(), [] (UInt & x) {
       if (x % 4 == 0 || x % 4 == 3)
 	x = 1;
       else
 	x = 0;
     });
 
   auto set_one = SetOneMatrix();
 
   Loop::stridedLoop(set_one, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "Static matrix strided loop failed";
 }
 
 TEST(TestStridedReduction, VectorReduce) {
   Grid<UInt, 2> grid({10, 10}, 3);
 
   Loop::stridedLoop([](VectorProxy<UInt, 3>&& v) {
       v(0) = 1;
       v(1) = 2;
       v(2) = 3;
     }, grid);
 
-  auto res = Loop::stridedReduce<operation::plus>([](VectorProxy<UInt, 3>&& v) -> Vector<UInt, 3> {
-      return v;
-    }, grid);
+  auto res = Loop::stridedReduce<operation::plus>(
+      [](VectorProxy<UInt, 3>&& v) -> Vector<UInt, 3> { return v; }, grid);
 
   ASSERT_EQ(res(0), 100);
   ASSERT_EQ(res(1), 200);
   ASSERT_EQ(res(2), 300);
 }