diff --git a/src/samplefunctional.hh b/src/samplefunctional.hh
index 1f5c409ed5af33a35791c428bd404b2351362164..a32e96645407c98d8ac43d9387cccef888b568c8 100644
--- a/src/samplefunctional.hh
+++ b/src/samplefunctional.hh
@@ -26,7 +26,7 @@ template <int dim> class SampleFunctional {
                    MyNonlinearity<dim> const &phi)
       : A(A), b(b), phi(phi) {}
 
-  double operator()(const SmallVector v) const {
+  double operator()(SmallVector const v) const {
     SmallVector y;
     A.mv(v, y);            // y = Av
     y /= 2;                // y = 1/2 Av
@@ -34,7 +34,7 @@ template <int dim> class SampleFunctional {
     return y * v + phi(v); // <1/2 Av - b,v> + H(|v|)
   }
 
-  void descentDirection(const SmallVector x, SmallVector &ret) const {
+  void descentDirection(SmallVector const x, SmallVector &ret) const {
     // Check the squared norm rather than each component because
     // negative_projection() divides by it
     if (x.two_norm2() == 0.0) {
@@ -87,21 +87,21 @@ template <int dim> class SampleFunctional {
 
 private:
   // Gradient of the smooth part
-  SmallVector smoothGradient(const SmallVector x) const {
+  SmallVector smoothGradient(SmallVector const x) const {
     SmallVector y;
     A.mv(x, y); // y = Av
     y -= b;     // y = Av - b
     return y;
   }
 
-  SmallVector upperGradient(const SmallVector x) const {
+  SmallVector upperGradient(SmallVector const x) const {
     SmallVector y;
     phi.upperGradient(x, y);
     y += smoothGradient(x);
     return y;
   }
 
-  SmallVector lowerGradient(const SmallVector x) const {
+  SmallVector lowerGradient(SmallVector const x) const {
     SmallVector y;
     phi.lowerGradient(x, y);
     y += smoothGradient(x);
@@ -109,8 +109,8 @@ template <int dim> class SampleFunctional {
   }
 
   // No normalising is done!
-  SmallVector negative_projection(const SmallVector z,
-                                  const SmallVector x) const {
+  SmallVector negative_projection(SmallVector const z,
+                                  SmallVector const x) const {
     SmallVector y = z;
     y.axpy(-(z * x) / x.two_norm2(), x);
     return y;
@@ -118,7 +118,7 @@ template <int dim> class SampleFunctional {
 };
 
 template <class Functional>
-void minimise(const Functional J, typename Functional::SmallVector &x,
+void minimise(Functional const J, typename Functional::SmallVector &x,
               size_t steps = 1,
               Bisection const &bisection =
                   Bisection(0.0, // acceptError: Stop if the search interval has