diff --git a/src/samplefunctional.hh b/src/samplefunctional.hh
index 6f2d2de496efa95f5670cb0b865739e6807061e6..7b06cbabdd8543e3a3f4731ec7cf7082761049b4 100644
--- a/src/samplefunctional.hh
+++ b/src/samplefunctional.hh
@@ -29,7 +29,7 @@ class SampleFunctional {
     return y * v + func_(v.two_norm()); // <1/2 Av - b,v> + H(|v|)
   }
 
-  SmallVector descentDirection(const SmallVector x) const {
+  void descentDirection(const SmallVector x, SmallVector &ret) const {
     if (x == SmallVector(0.0)) {
       SmallVector d = smoothGradient(x);
       // Decline of the smooth part in the negative gradient direction
@@ -38,16 +38,17 @@ class SampleFunctional {
           func_.rightDifferential(0.0) * d.two_norm(); // TODO: is this correct?
       double combinedDecline = smoothDecline + nonlinearDecline;
 
-      return (combinedDecline < 0) ? d : SmallVector(0.0);
+      ret = (combinedDecline < 0) ? d : SmallVector(0.0);
+      return;
     }
 
     SmallVector const pg = upperGradient(x);
     SmallVector const mg = lowerGradient(x);
-    SmallVector ret;
     // TODO: collinearity checks suck
     if (pg * x == pg.two_norm() * x.two_norm() &&
         -(mg * x) == mg.two_norm() * x.two_norm()) {
-      return SmallVector(0);
+      ret = SmallVector(0);
+      return;
     } else if (pg * x >= 0 && mg * x >= 0) {
       ret = pg;
     } else if (pg * x <= 0 && mg * x <= 0) {
@@ -56,7 +57,6 @@ class SampleFunctional {
       ret = project(smoothGradient(x), x);
     }
     ret *= -1;
-    return ret;
   }
 
   SmallMatrix A;
@@ -97,7 +97,8 @@ template <class Functional>
 void minimise(const Functional J, const typename Functional::SmallVector x,
               typename Functional::SmallVector &corr) {
   typedef typename Functional::SmallVector SmallVector;
-  SmallVector descDir = J.descentDirection(x);
+  SmallVector descDir;
+  J.descentDirection(x, descDir);
 
   if (descDir == SmallVector(0.0)) {
     corr = SmallVector(0.0);
diff --git a/src/test-gradient-method.cc b/src/test-gradient-method.cc
index 5ff673375849a0d5b1ed4ae97796da1c5f3602b3..3e6c2d2c3a707a601cba3918f15e19648dcaa227 100644
--- a/src/test-gradient-method.cc
+++ b/src/test-gradient-method.cc
@@ -55,7 +55,8 @@ void testSampleFunction() {
   SampleFunctional::SmallVector error;
   error[0] = -(102 - 1 + 2 / sqrt(5));
   error[1] = -(161.5 - 2 + 4 / sqrt(5));
-  SampleFunctional::SmallVector returned = J.descentDirection(start);
+  SampleFunctional::SmallVector returned;
+  J.descentDirection(start, returned);
   error -= returned;
   assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
 
@@ -97,7 +98,8 @@ void testSampleFunctionNonsmooth() {
     start /= (b.two_norm() + 1e-12);
     assert(start.two_norm() < 1);
 
-    SampleFunctional::SmallVector returned = J.descentDirection(start);
+    SampleFunctional::SmallVector returned;
+    J.descentDirection(start, returned);
     error[0] = -(7 / sqrt(5) - 1);
     error[1] = -(11.5 / sqrt(5) - 2);
     error -= returned;
@@ -111,7 +113,8 @@ void testSampleFunctionNonsmooth() {
     start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
     assert(start.two_norm() > 1);
 
-    SampleFunctional::SmallVector returned = J.descentDirection(start);
+    SampleFunctional::SmallVector returned;
+    J.descentDirection(start, returned);
     error[0] = -(8 / sqrt(5) - 1);
     error[1] = -(13.5 / sqrt(5) - 2);
     error -= returned;
@@ -148,7 +151,8 @@ void testTrivialFunction() {
   SampleFunctional::SmallVector error;
   error[0] = -101;
   error[1] = -159.5;
-  SampleFunctional::SmallVector returned = J.descentDirection(start);
+  SampleFunctional::SmallVector returned;
+  J.descentDirection(start, returned);
   error -= returned;
   assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?