Test directions of descent

src/test-gradient-method.cc

@@ -60,6 +60,65 @@ void testSampleFunction() {
   functionTester(J, start, 6);
 }
 
+void testSampleFunctionNonsmooth() {
+  int const dim = 2;
+  typedef Dune::SampleFunctional<dim, Dune::SampleFunction> SampleFunctional;
+
+  SampleFunctional::SmallMatrix A;
+  A[0][0] = 3;
+  A[0][1] = 1.5;
+  A[1][0] = 1.5;
+  A[1][1] = 4;
+  SampleFunctional::SmallVector b;
+  b[0] = 1;
+  b[1] = 2;
+
+  SampleFunctional J(A, b);
+
+  SampleFunctional::SmallVector start;
+  SampleFunctional::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);
+
+    SampleFunctional::SmallVector returned = J.descentDirection(start);
+    error[0] = -(7 / sqrt(5) - 1);
+    error[1] = -(11.5 / sqrt(5) - 2);
+    error -= returned;
+    assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
+
+    functionTester(J, start, 6);
+  }
+  std::cout << std::endl;
+  {
+    start = b;
+    start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;
+    assert(start.two_norm() > 1);
+
+    SampleFunctional::SmallVector returned = J.descentDirection(start);
+    error[0] = -(8 / sqrt(5) - 1);
+    error[1] = -(13.5 / sqrt(5) - 2);
+    error -= returned;
+    assert(error.two_norm() < 1e-10); // FIXME: 1e-10 sounds reasonable. Is it?
+
+    functionTester(J, start, 6);
+  }
+}
+
 void testTrivialFunction() {
   int const dim = 2;
   typedef Dune::SampleFunctional<dim> SampleFunctional;
@@ -98,6 +157,8 @@ int main() {
   try {
     testSampleFunction();
     std::cout << std::endl;
+    testSampleFunctionNonsmooth();
+    std::cout << std::endl;
     testTrivialFunction();
     return 0;
   }