diff --git a/src/surface/isopowerlaw.hh b/src/surface/isopowerlaw.hh
index b0dfab1..42ff39d 100644
--- a/src/surface/isopowerlaw.hh
+++ b/src/surface/isopowerlaw.hh
@@ -1,104 +1,105 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 ISOPOWERLAW_H
 #define ISOPOWERLAW_H
 /* -------------------------------------------------------------------------- */
 #include "filter.hh"
 #include "grid_hermitian.hh"
 #include "tamaas.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /// Class representing an isotropic power law spectrum
 template <UInt dim>
 class Isopowerlaw : public Filter<dim> {
 public:
   /// Default constructor
   Isopowerlaw() = default;
   /// Destructor
   virtual ~Isopowerlaw() = default;
 
 public:
   /// Compute filter coefficients
   virtual void
   computeFilter(GridHermitian<Real, dim>& filter_coefficients) const;
 
   /// Compute a point of the PSD
-  inline Real operator()(const VectorProxy<Real, dim>& q_vec) const;
+  __device__ __host__ inline Real
+  operator()(const VectorProxy<Real, dim>& q_vec) const;
 
   /// Analytical rms of heights
   Real rmsHeights() const;
 
   /// Analytical moments
   std::vector<Real> moments() const;
 
   /// Analytical Nayak's parameter
   Real alpha() const;
 
   /// Analytical RMS slopes
   Real rmsSlopes() const;
 
 public:
   TAMAAS_ACCESSOR(q0, UInt, Q0);
   TAMAAS_ACCESSOR(q1, UInt, Q1);
   TAMAAS_ACCESSOR(q2, UInt, Q2);
   TAMAAS_ACCESSOR(hurst, Real, Hurst);
 
 protected:
   UInt q0, q1, q2;
   Real hurst;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline implementations                                                     */
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 inline Real Isopowerlaw<dim>::operator()(const VectorProxy<Real, dim>& q_vec) const {
   const Real C = 1.;
   const Real q = q_vec.l2norm();
 
   Real val;
 
   if (q < q0)
     val = 0;
   else if (q > q2)
     val = 0;
   else if (q < q1)
     val = C;
   else
     val = C * std::pow((q / q1), -2 * (hurst + 1));
 
   return std::sqrt(val);
 }
 
 __END_TAMAAS__
 
 #endif  // ISOPOWERLAW_H
diff --git a/tests/test_grid.cpp b/tests/test_grid.cpp
index eb6b202..5dc9261 100644
--- a/tests/test_grid.cpp
+++ b/tests/test_grid.cpp
@@ -1,199 +1,201 @@
 /**
  *
  * @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 "grid_view.hh"
 
 using namespace tamaas;
 
 /* -------------------------------------------------------------------------- */
 
 // Testing 1D creation
 TEST(TestGridCreation, Creation1d) {
   Grid<Real, 1> grid{{5}, 2};
   std::array<UInt, 1> correct_size{{5}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 2> correct_strides{{2, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 // Testing 2D creation
 TEST(TestGridCreation, Creation2d) {
   Grid<Real, 2> grid({5, 2}, 3);
   std::array<UInt, 2> correct_size{{5, 2}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 3> correct_strides{{6, 3, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 // Testing 3D creation
 TEST(TestGridCreation, Creation3d) {
   Grid<Real, 3> grid({8, 5, 2}, 3);
   std::array<UInt, 3> correct_size{{8, 5, 2}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 4> correct_strides{{30, 15, 3, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 /* -------------------------------------------------------------------------- */
 
 // Testing iterators in STL function iota
 TEST(TestGridIterators, Iota) {
   constexpr UInt N = 20;
   Grid<UInt, 1> grid({N}, 1);
   std::iota(grid.begin(), grid.end(), 0);
 
   const UInt * p = grid.getInternalData();
   for (UInt i = 0 ; i < N ; ++i)
     ASSERT_TRUE(p[i] == i) << "Iota fill failed";
 }
 
 // Testing filling grid with OpenMP loop
 TEST(TestGridIterators, OpenMPLoops) {
   Grid<UInt, 1> grid({20}, 1);
 
 #ifndef USE_CUDA
 #pragma omp parallel for
 #endif
   for (auto it = grid.begin() ; it < grid.end() ; ++it) {
     UInt i = it - grid.begin();
     *it = i;
   }
 
   std::vector<UInt> correct(20);
   std::iota(correct.begin(), correct.end(), 0);
 
   ASSERT_TRUE(compare(grid, correct)) << "Fill using OpenMP loop failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 // Testing scalar simple loop-based operators
 TEST(TestGridOperators, LoopOperators) {
   Grid<UInt, 1> grid({20}, 1);
   grid = 1;
   EXPECT_TRUE(std::all_of(grid.begin(), grid.end(), [](UInt x) {
     return x == 1;
   })) << "Assignement operator failed";
   grid += 1;
   EXPECT_TRUE(std::all_of(grid.begin(), grid.end(), [](UInt x) {
     return x == 2;
   })) << "Plus-equal operator failed";
 }
 
 // Testing loop-based functions with reductions
 TEST(TestGridOperators, Stats) {
   constexpr UInt N = 20;
   Grid<Real, 1> grid({N}, 1);
   std::iota(grid.begin(), grid.end(), 0);
 
   auto b = grid.begin(), e = grid.end();
   Real min = *std::min_element(b, e);
   Real max = *std::max_element(b, e);
   Real acc = std::accumulate(b, e, 0);
   Real mea = acc / N;
   Real var = 35;
 
   EXPECT_DOUBLE_EQ(grid.min(), min) << "Minimum function failed";
   EXPECT_DOUBLE_EQ(grid.max(), max) << "Maximum function failed";
   EXPECT_DOUBLE_EQ(grid.sum(), acc) <<"Sum function failed";
   EXPECT_DOUBLE_EQ(grid.mean(), mea) << "Mean function failed";
   EXPECT_DOUBLE_EQ(grid.var(), var) << "Var function failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 TEST(TestGridView, View2D) {
   Grid<UInt, 3> grid3d({10, 10, 10}, 3);
   std::iota(grid3d.begin(), grid3d.end(), 0);
 
   for (UInt i : Loop::range(10)) {
     auto grid_view = make_view(grid3d, i);
     std::vector<UInt> reference(300);
     std::iota(reference.begin(), reference.end(), 300 * i);
     EXPECT_TRUE(
         std::equal(grid_view.begin(), grid_view.end(), reference.begin()));
   }
 }
 
 /* -------------------------------------------------------------------------- */
+
+struct BroadcastSet123 {
+  CUDA_LAMBDA inline void operator()(VectorProxy<UInt, 3>&& v) {
+    v(0) = 1;
+    v(1) = 2;
+    v(2) = 3;
+  }
+};
+
 TEST(TestGridOperators, BroadcastOperators) {
   Grid<UInt, 2> grid({10, 10}, 3), solution({10, 10}, 3);
+  auto set = BroadcastSet123();
 
   // Filling up solution
-  Loop::stridedLoop([](VectorProxy<UInt, 3>&& v) {
-      v(0) = 1;
-      v(1) = 2;
-      v(2) = 3;
-    }, solution);
+  Loop::stridedLoop(set, solution);
 
   Vector<UInt, 3> v;
   v(0) = 1;
   v(1) = 2;
   v(2) = 3;
 
   grid = v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast assignement failure";
 
   v(1) = 1;
   v(2) = 1;
 
   solution += 1;
   grid += v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast += failure";
 
   solution -= 1;
   grid -= v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast -= failure";
 
   v(1) = 2;
   v(2) = 3;
 
   solution *= solution;
   grid *= v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast *= failure";
 
-  Loop::stridedLoop([](VectorProxy<UInt, 3>&& v) {
-      v(0) = 1;
-      v(1) = 2;
-      v(2) = 3;
-    }, solution);
+  Loop::stridedLoop(set, solution);
   grid /= v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast /= failure";
 }
diff --git a/tests/test_loop.cpp b/tests/test_loop.cpp
index 4124a30..66cb4eb 100644
--- a/tests/test_loop.cpp
+++ b/tests/test_loop.cpp
@@ -1,254 +1,257 @@
 /**
  *
  * @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";
 }
 
+struct VectorReduction {
+  CUDA_LAMBDA Vector<UInt, 3> operator()(VectorProxy<UInt, 3>&& v) { return v; }
+};
+
 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>(VectorReduction(), grid);
 
   ASSERT_EQ(res(0), 100);
   ASSERT_EQ(res(1), 200);
   ASSERT_EQ(res(2), 300);
 }