Add more debugging output

src/samplefunctional.hh

@@ -51,8 +51,12 @@ class SampleFunctional {
       return;
     } else if (pg * x >= 0 && mg * x >= 0) {
       ret = pg;
+      Dune::dverb << "## Expected directional derivative: " << -(pg * mg)
+                  << " (coordinates of the restriction)" << std::endl;
     } else if (pg * x <= 0 && mg * x <= 0) {
       ret = mg;
+      Dune::dverb << "## Expected directional derivative: " << -(mg * pg)
+                  << " (coordinates of the restriction)" << std::endl;
     } else {
       ret = project(smoothGradient(x), x);
     }
@@ -129,6 +133,15 @@ void minimise(const Functional J, const typename Functional::SmallVector x,
   MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir);
   // }}}
 
+  { // Debug
+    Interval<double> D;
+    JRest.subDiff(0, D);
+    Dune::dverb << "## Descent according to JRest.subDiff: " << D[1]
+                << std::endl;
+    assert(D[1] <=
+           0); // We should not be minimising in this direction otherwise
+  }
+
   // FIXME: default values are used
   Bisection bisection;
   int count;

src/test-gradient-method.cc

@@ -108,6 +108,8 @@ void testSampleFunctionNonsmooth() {
     functionTester(J, start, 6);
   }
   std::cout << std::endl;
+  std::cout << std::endl;
+  std::cout << std::endl;
   {
     start = b;
     start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;