diff --git a/configure.ac b/configure.ac
index d04cabd11e5e24d4f5ab718c3b3c6f7ab3b3dc39..7416cf732ce292b78e9f2986b74a5f86605891e2 100644
--- a/configure.ac
+++ b/configure.ac
@@ -4,7 +4,7 @@ AC_PREREQ(2.50)
 DUNE_AC_INIT # gets module version from dune.module file
 AM_INIT_AUTOMAKE
 AM_SILENT_RULES
-AC_CONFIG_SRCDIR([src/test-gradient-method.cc])
+AC_CONFIG_SRCDIR([src/one-body-sample.cc])
 AM_CONFIG_HEADER([config.h])
 
 
diff --git a/src/Makefile.am b/src/Makefile.am
index 2e267a9767ec7a62e98b3f8cf49de2d54edf6e31..e0ade588c821f1aa1618dceb955baf179df0cd11 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,7 +1,17 @@
 SUBDIRS = octave
 
 check_PROGRAMS = \
-	test-gradient-method
+	test-gradient-horrible \
+	test-gradient-horrible-logarithmic \
+	test-gradient-identity \
+	test-gradient-sample \
+	test-gradient-sample-3d \
+	test-gradient-sample-nonsmooth \
+	test-gradient-sample-steep \
+	test-gradient-sample-steep2 \
+	test-gradient-sample-verysteep \
+	test-gradient-sample2 \
+	test-gradient-trivial
 
 bin_PROGRAMS = \
 	one-body-sample-2D \
@@ -63,10 +73,7 @@ one_body_sample_3D_SOURCES = \
 one_body_sample_3D_CPPFLAGS = \
 	$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\" -DDIM=3
 
-test_gradient_method_SOURCES = \
-	test-gradient-method.cc
-
-TESTS= test-gradient-method
+TESTS = $(check_PROGRAMS)
 
 # Some are for clang, others are for gcc
 AM_CXXFLAGS = \
diff --git a/src/test-gradient-horrible-logarithmic.cc b/src/test-gradient-horrible-logarithmic.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c6aaa9ff4917ed574678f91e2a96eddb0d56eef3
--- /dev/null
+++ b/src/test-gradient-horrible-logarithmic.cc
@@ -0,0 +1,50 @@
+/* Checks if the algorithm converges regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 3;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 4;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto f = Dune::make_shared<Dune::HorribleFunctionLogarithmic const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start = b;
+  start *= 17;
+
+  double const ret1 = functionTester(J, start, 6);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 12);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 4);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-horrible.cc b/src/test-gradient-horrible.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f864e891ca0c7a79ee78445936860771035d3fd7
--- /dev/null
+++ b/src/test-gradient-horrible.cc
@@ -0,0 +1,50 @@
+/* Checks if the algorithm converges regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 3;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 4;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto f = Dune::make_shared<Dune::HorribleFunction const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start = b;
+  start *= 17;
+
+  double const ret1 = functionTester(J, start, 7);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 8);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 4);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-identity.cc b/src/test-gradient-identity.cc
new file mode 100644
index 0000000000000000000000000000000000000000..596b30086d7b5ad97a088ce099c031f23896ed10
--- /dev/null
+++ b/src/test-gradient-identity.cc
@@ -0,0 +1,50 @@
+/* Checks if the algorithm converges regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 3;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 4;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto f = Dune::make_shared<Dune::LinearFunction const>(1);
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start = b;
+  start *= 17;
+
+  double const ret1 = functionTester(J, start, 6);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 10);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 3);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-method.cc b/src/test-gradient-method.cc
deleted file mode 100644
index 3348aea7a0acab5afdc3e4cad834e717fad9de08..0000000000000000000000000000000000000000
--- a/src/test-gradient-method.cc
+++ /dev/null
@@ -1,497 +0,0 @@
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <cassert>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/shared_ptr.hh>
-
-#include <dune/tectonic/samplefunctional.hh>
-
-#include "test-gradient-method-nicefunction.hh"
-#include "test-gradient-method-helper.hh"
-
-void testIdentity() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 3;
-  A[0][1] = A[1][0] = 1.5;
-  A[1][1] = 4;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-
-  auto f = Dune::make_shared<Dune::LinearFunction const>(1);
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start = b;
-  start *= 17;
-
-  double const ret1 = functionTester(J, start, 6);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 10);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 3);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-void testSampleFunction() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 3;
-  A[0][1] = A[1][0] = 1.5;
-  A[1][1] = 4;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-
-  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start = b;
-  start *= 17;
-
-  /*
-    j(x)
-    = Ax - b + 2/|x| x
-    = 17*(6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
-    = (102 - 1 + 2/sqrt(5), 161.5 - 2 + 4/sqrt(5))
-  */
-  Functional::SmallVector error;
-  error[0] = -(102 - 1 + 2 / sqrt(5));
-  error[1] = -(161.5 - 2 + 4 / sqrt(5));
-  Functional::SmallVector returned;
-  J.descentDirection(start, returned);
-  error -= returned;
-  assert(error.two_norm() < 1e-10);
-
-  double const ret1 = functionTester(J, start, 6);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 10);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 3);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-void testSampleFunctionNonsmooth() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 3;
-  A[0][1] = A[1][0] = 1.5;
-  A[1][1] = 4;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-
-  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start;
-  Functional::SmallVector error;
-  /*
-    for x = b/|b|:
-
-    j_(x)
-    = Ax - b + 1/|x| x
-    = 1/sqrt(5) (6, 9.5) - (1, 2) + 1/sqrt(5) (1, 2)
-    = (7/sqrt(5) - 1, 11.5/sqrt(5) - 2)
-
-    j+(x)
-    = Ax - b + 2/|x| x
-    = 1/sqrt(5) (6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
-    = (8/sqrt(5) - 1, 13.5/sqrt(5) - 2)
-  */
-  {
-    start = b;
-    start /= (b.two_norm() + 1e-12);
-    assert(start.two_norm() < 1);
-
-    Functional::SmallVector returned;
-    J.descentDirection(start, returned);
-    error[0] = -(7 / sqrt(5) - 1);
-    error[1] = -(11.5 / sqrt(5) - 2);
-    error -= returned;
-    assert(error.two_norm() < 1e-10);
-
-    functionTester(J, start, 6);
-  }
-  {
-    start = b;
-    start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
-    assert(start.two_norm() > 1);
-
-    Functional::SmallVector returned;
-    J.descentDirection(start, returned);
-    error[0] = -(8 / sqrt(5) - 1);
-    error[1] = -(13.5 / sqrt(5) - 2);
-    error -= returned;
-    assert(error.two_norm() < 1e-10);
-
-    functionTester(J, start, 6);
-  }
-}
-
-void testTrivialFunction() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 3;
-  A[0][1] = A[1][0] = 1.5;
-  A[1][1] = 4;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-
-  auto f = Dune::make_shared<Dune::TrivialFunction const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start = b;
-  start *= 17;
-
-  /*
-    j(x)
-    = Ax - b
-    = 17*(6, 9.5) - (1, 2)
-    = (102 - 1, 161.5 - 2)
-  */
-  Functional::SmallVector error;
-  error[0] = -101;
-  error[1] = -159.5;
-  Functional::SmallVector returned;
-  J.descentDirection(start, returned);
-  error -= returned;
-  assert(error.two_norm() < 1e-10);
-
-  double const ret1 = functionTester(J, start, 6);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 16);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-}
-
-void testHorribleFunction() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 3;
-  A[0][1] = A[1][0] = 1.5;
-  A[1][1] = 4;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-
-  auto f = Dune::make_shared<Dune::HorribleFunction const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start = b;
-  start *= 17;
-
-  double const ret1 = functionTester(J, start, 7);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 8);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 4);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-void testHorribleFunctionLogarithmic() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 3;
-  A[0][1] = A[1][0] = 1.5;
-  A[1][1] = 4;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-
-  auto f = Dune::make_shared<Dune::HorribleFunctionLogarithmic const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start = b;
-  start *= 17;
-
-  double const ret1 = functionTester(J, start, 6);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 12);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 4);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-void testSampleFunction3D() {
-  int const dim = 3;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A(0);
-  A[0][0] = 3;
-  A[0][1] = A[1][0] = 1.5;
-  A[1][1] = 4;
-  A[1][2] = A[2][1] = 1.5;
-  A[2][2] = 5;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-  b[2] = 3;
-
-  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start = b;
-  start *= 17;
-
-  double const ret1 = functionTester(J, start, 9);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-  start[2] = -15;
-
-  double const ret2 = functionTester(J, start, 15);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-  start[2] = 0;
-
-  double const ret3 = functionTester(J, start, 5);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-// Checks if reaching the minimum in one step leads to problems
-void testSampleFunction2() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 1;
-  A[0][1] = A[1][0] = 0;
-  A[1][1] = 1;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 1;
-
-  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start = b;
-
-  double const ret1 = functionTester(J, start, 2);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 7);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 2);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-void testSampleFunctionSteep1() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 1;
-  A[0][1] = A[1][0] = 0;
-  A[1][1] = 1;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2;
-
-  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start;
-
-  start[0] = 0;
-  start[1] = 1;
-
-  double const ret1 = functionTester(J, start, 3);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 3);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 1);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-void testSampleFunctionSteep2() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 1;
-  A[0][1] = A[1][0] = 0;
-  A[1][1] = 1;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2.5;
-
-  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start;
-
-  start[0] = 0;
-  start[1] = 1;
-
-  double const ret1 = functionTester(J, start, 1);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 3);
-  assert(std::abs(ret1 - ret2) < 1e-5);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 1);
-  assert(std::abs(ret1 - ret3) < 1e-5);
-}
-
-void testSteepFunction() {
-  int const dim = 2;
-  typedef Dune::SampleFunctional<dim> Functional;
-
-  Functional::SmallMatrix A;
-  A[0][0] = 1;
-  A[0][1] = A[1][0] = 0;
-  A[1][1] = 1;
-  Functional::SmallVector b;
-  b[0] = 1;
-  b[1] = 2.5;
-
-  auto f = Dune::make_shared<Dune::SampleFunction<100> const>();
-  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
-  Functional J(A, b, phi);
-
-  Functional::SmallVector start;
-
-  start[0] = 0;
-  start[1] = 1;
-
-  double const ret1 = functionTester(J, start, 1);
-
-  // Something random
-  start[0] = 279;
-  start[1] = -96;
-
-  double const ret2 = functionTester(J, start, 4);
-  assert(std::abs(ret1 - ret2) < 1e-8);
-
-  start[0] = 0;
-  start[1] = 0;
-
-  double const ret3 = functionTester(J, start, 1);
-  assert(std::abs(ret2 - ret3) < 1e-5);
-}
-
-int main() {
-  try {
-    std::cout << "testIdentity:" << std::endl;
-    testIdentity();
-
-    std::cout << std::endl << "testSampleFunction:" << std::endl;
-    testSampleFunction();
-
-    std::cout << std::endl << "testSampleFunctionNonsmooth:" << std::endl;
-    testSampleFunctionNonsmooth();
-
-    std::cout << std::endl << "testTrivialFunction:" << std::endl;
-    testTrivialFunction();
-
-    std::cout << std::endl << "testHorribleFunction:" << std::endl;
-    testHorribleFunction();
-
-    std::cout << std::endl << "testHorribleFunctionLogarithmic:" << std::endl;
-    testHorribleFunctionLogarithmic();
-
-    std::cout << std::endl << "testSampleFunction2:" << std::endl;
-    testSampleFunction2();
-
-    std::cout << std::endl << "testSampleFunction3D:" << std::endl;
-    testSampleFunction3D();
-
-    std::cout << std::endl << "testSampleFunctionSteep1:" << std::endl;
-    testSampleFunctionSteep1();
-
-    std::cout << std::endl << "testSampleFunctionSteep2:" << std::endl;
-    testSampleFunctionSteep2();
-
-    std::cout << std::endl << "testSteepFunction:" << std::endl;
-    testSteepFunction();
-  }
-  catch (Dune::Exception &e) {
-    Dune::derr << "Dune reported error: " << e << std::endl;
-  }
-}
diff --git a/src/test-gradient-sample-3d.cc b/src/test-gradient-sample-3d.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1ed6a05d07896341ceec7e718b20cf7b8f1345dd
--- /dev/null
+++ b/src/test-gradient-sample-3d.cc
@@ -0,0 +1,55 @@
+/* Checks if the algorithm converges regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 3;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A(0);
+  A[0][0] = 3;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 4;
+  A[1][2] = A[2][1] = 1.5;
+  A[2][2] = 5;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+  b[2] = 3;
+
+  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start = b;
+  start *= 17;
+
+  double const ret1 = functionTester(J, start, 9);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+  start[2] = -15;
+
+  double const ret2 = functionTester(J, start, 15);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+  start[2] = 0;
+
+  double const ret3 = functionTester(J, start, 5);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-sample-nonsmooth.cc b/src/test-gradient-sample-nonsmooth.cc
new file mode 100644
index 0000000000000000000000000000000000000000..6fb3861567cd382f76ce1a2c30487886112fdd41
--- /dev/null
+++ b/src/test-gradient-sample-nonsmooth.cc
@@ -0,0 +1,77 @@
+/* Checks if the descent direction is computed correctly
+   using the analytic solution */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 3;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 4;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start;
+  Functional::SmallVector error;
+  /*
+    for x = b/|b|:
+
+    j_(x)
+    = Ax - b + 1/|x| x
+    = 1/sqrt(5) (6, 9.5) - (1, 2) + 1/sqrt(5) (1, 2)
+    = (7/sqrt(5) - 1, 11.5/sqrt(5) - 2)
+
+    j+(x)
+    = Ax - b + 2/|x| x
+    = 1/sqrt(5) (6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
+    = (8/sqrt(5) - 1, 13.5/sqrt(5) - 2)
+  */
+  {
+    start = b;
+    start /= (b.two_norm() + 1e-12);
+    assert(start.two_norm() < 1);
+
+    Functional::SmallVector returned;
+    J.descentDirection(start, returned);
+    error[0] = -(7 / sqrt(5) - 1);
+    error[1] = -(11.5 / sqrt(5) - 2);
+    error -= returned;
+    assert(error.two_norm() < 1e-10);
+
+    functionTester(J, start, 6);
+  }
+  {
+    start = b;
+    start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
+    assert(start.two_norm() > 1);
+
+    Functional::SmallVector returned;
+    J.descentDirection(start, returned);
+    error[0] = -(8 / sqrt(5) - 1);
+    error[1] = -(13.5 / sqrt(5) - 2);
+    error -= returned;
+    assert(error.two_norm() < 1e-10);
+
+    functionTester(J, start, 6);
+  }
+}
diff --git a/src/test-gradient-sample-steep.cc b/src/test-gradient-sample-steep.cc
new file mode 100644
index 0000000000000000000000000000000000000000..dfe498f2c018d50237f24afaaa9c5b02dc09e79b
--- /dev/null
+++ b/src/test-gradient-sample-steep.cc
@@ -0,0 +1,52 @@
+/* Checks if the algorithm converges regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 1;
+  A[0][1] = A[1][0] = 0;
+  A[1][1] = 1;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start;
+
+  start[0] = 0;
+  start[1] = 1;
+
+  double const ret1 = functionTester(J, start, 3);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 3);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 1);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-sample-steep2.cc b/src/test-gradient-sample-steep2.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f5892cd5f848438a27a04f9449a396ce20137736
--- /dev/null
+++ b/src/test-gradient-sample-steep2.cc
@@ -0,0 +1,52 @@
+/* Checks if the algorithm converges regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 1;
+  A[0][1] = A[1][0] = 0;
+  A[1][1] = 1;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2.5;
+
+  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start;
+
+  start[0] = 0;
+  start[1] = 1;
+
+  double const ret1 = functionTester(J, start, 1);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 3);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 1);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-sample-verysteep.cc b/src/test-gradient-sample-verysteep.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0cc230e6ab02431e069d6487e1b13b6c0431e467
--- /dev/null
+++ b/src/test-gradient-sample-verysteep.cc
@@ -0,0 +1,52 @@
+/* Checks if the algorithm converges regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 1;
+  A[0][1] = A[1][0] = 0;
+  A[1][1] = 1;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2.5;
+
+  auto f = Dune::make_shared<Dune::SampleFunction<100> const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start;
+
+  start[0] = 0;
+  start[1] = 1;
+
+  double const ret1 = functionTester(J, start, 1);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 4);
+  assert(std::abs(ret1 - ret2) < 1e-8);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 1);
+  assert(std::abs(ret2 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-sample.cc b/src/test-gradient-sample.cc
new file mode 100644
index 0000000000000000000000000000000000000000..80f0916281e9c50e7f569242e7cce4cfdb3c694c
--- /dev/null
+++ b/src/test-gradient-sample.cc
@@ -0,0 +1,66 @@
+/* Checks if the descent direction is computed correctly using the
+   analytic solution; also checks if the algorithm converges
+   regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 3;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 4;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start = b;
+  start *= 17;
+
+  /*
+    j(x)
+    = Ax - b + 2/|x| x
+    = 17*(6, 9.5) - (1, 2) + 2/sqrt(5) (1, 2)
+    = (102 - 1 + 2/sqrt(5), 161.5 - 2 + 4/sqrt(5))
+  */
+  Functional::SmallVector error;
+  error[0] = -(102 - 1 + 2 / sqrt(5));
+  error[1] = -(161.5 - 2 + 4 / sqrt(5));
+  Functional::SmallVector returned;
+  J.descentDirection(start, returned);
+  error -= returned;
+  assert(error.two_norm() < 1e-10);
+
+  double const ret1 = functionTester(J, start, 6);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 10);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 3);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-sample2.cc b/src/test-gradient-sample2.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ce6382ff9f94fd4cb9e043c9ddc0c5de4a640991
--- /dev/null
+++ b/src/test-gradient-sample2.cc
@@ -0,0 +1,51 @@
+/* Checks if the algorithm converges regardless of where it starts;
+   in particular, that converging in one step does not cause us to
+   diverge again */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 1;
+  A[0][1] = A[1][0] = 0;
+  A[1][1] = 1;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 1;
+
+  auto f = Dune::make_shared<Dune::SampleFunction<2> const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start = b;
+
+  double const ret1 = functionTester(J, start, 2);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 7);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+
+  start[0] = 0;
+  start[1] = 0;
+
+  double const ret3 = functionTester(J, start, 2);
+  assert(std::abs(ret1 - ret3) < 1e-5);
+}
diff --git a/src/test-gradient-trivial.cc b/src/test-gradient-trivial.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ed16c4f7725009b8d5c5f74bed16e8acda791178
--- /dev/null
+++ b/src/test-gradient-trivial.cc
@@ -0,0 +1,60 @@
+/* Checks if the descent direction is computed correctly using the
+   analytic solution; also checks if the algorithm converges
+   regardless of where it starts */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/samplefunctional.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim> Functional;
+
+  Functional::SmallMatrix A;
+  A[0][0] = 3;
+  A[0][1] = A[1][0] = 1.5;
+  A[1][1] = 4;
+  Functional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  auto f = Dune::make_shared<Dune::TrivialFunction const>();
+  auto phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+  Functional J(A, b, phi);
+
+  Functional::SmallVector start = b;
+  start *= 17;
+
+  /*
+    j(x)
+    = Ax - b
+    = 17*(6, 9.5) - (1, 2)
+    = (102 - 1, 161.5 - 2)
+  */
+  Functional::SmallVector error;
+  error[0] = -101;
+  error[1] = -159.5;
+  Functional::SmallVector returned;
+  J.descentDirection(start, returned);
+  error -= returned;
+  assert(error.two_norm() < 1e-10);
+
+  double const ret1 = functionTester(J, start, 6);
+
+  // Something random
+  start[0] = 279;
+  start[1] = -96;
+
+  double const ret2 = functionTester(J, start, 16);
+  assert(std::abs(ret1 - ret2) < 1e-5);
+}