Set minimise() free

src/bisection-example-new.cc

@@ -18,6 +18,7 @@ class SampleFunctional {
 public:
   typedef Dune::FieldVector<double, dimension> SmallVector;
   typedef Dune::FieldMatrix<double, dimension, dimension> SmallMatrix;
+  typedef Function FunctionType;
 
   SampleFunctional(SmallMatrix A, SmallVector b) : A_(A), b_(b), func_() {}
 
@@ -59,50 +60,10 @@ class SampleFunctional {
     return ret;
   }
 
-  SmallVector minimise(const SmallVector x) const {
-    SmallVector descDir = descentDirection(x);
-
-    if (descDir == SmallVector(0.0))
-      return SmallVector(0.0);
-
-    /* We have
-
-       1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>
-
-       since A is symmetric.
-     */
-    SmallVector tmp;
-
-    A_.mv(descDir, tmp);                 // Av
-    double const rest_A = tmp * descDir; // <Av,v>
-
-    A_.mv(x, tmp);                              // Au
-    double const rest_b = (b_ - tmp) * descDir; // <b-Au,v>
-
-    typedef MyNonlinearity<dimension, Function> MyNonlinearityType;
-    MyNonlinearityType phi;
-    typedef DirectionalConvexFunction<MyNonlinearityType>
-    MyDirectionalConvexFunctionType;
-    MyDirectionalConvexFunctionType rest(rest_A, rest_b, phi, x, descDir);
-
-    // FIXME: default values are used
-    Bisection bisection;
-    int count;
-    // FIXME: does x_old = 1 make any sense?!
-    double const m = bisection.minimize(rest, 0.0, 1.0, count);
-    Dune::dverb << "Number of iterations in the bisection method: " << count
-                << std::endl;
-    ;
-
-    SmallVector middle = descDir;
-    middle *= m;
-    return middle;
-  }
-
-private:
   SmallMatrix A_;
   SmallVector b_;
 
+private:
   Function func_;
 
   // Gradient of the smooth part
@@ -133,6 +94,51 @@ class SampleFunctional {
   }
 };
 
+// TODO: take a reference to a correction
+
+template <class Functional>
+typename Functional::SmallVector minimise(
+    const Functional J, const typename Functional::SmallVector x) {
+  typename Functional::SmallVector descDir = J.descentDirection(x);
+
+  if (descDir == typename Functional::SmallVector(0.0))
+    return typename Functional::SmallVector(0.0);
+
+  /* We have
+
+     1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>
+
+     since A is symmetric.
+  */
+  typename Functional::SmallVector tmp;
+
+  J.A_.mv(descDir, tmp);               // Av
+  double const rest_A = tmp * descDir; // <Av,v>
+
+  J.A_.mv(x, tmp);                              // Au
+  double const rest_b = (J.b_ - tmp) * descDir; // <b-Au,v>
+
+  typedef MyNonlinearity<Functional::SmallVector::dimension,
+                         typename Functional::FunctionType> MyNonlinearityType;
+  MyNonlinearityType phi;
+  typedef DirectionalConvexFunction<MyNonlinearityType>
+  MyDirectionalConvexFunctionType;
+  MyDirectionalConvexFunctionType rest(rest_A, rest_b, phi, x, descDir);
+
+  // FIXME: default values are used
+  Bisection bisection;
+  int count;
+  // FIXME: does x_old = 1 make any sense?!
+  double const m = bisection.minimize(rest, 0.0, 1.0, count);
+  Dune::dverb << "Number of iterations in the bisection method: " << count
+              << std::endl;
+  ;
+
+  typename Functional::SmallVector middle = descDir;
+  middle *= m;
+  return middle;
+}
+
 template <int dim, class Function>
 void functionTester(
     SampleFunctional<dim, Function> J,
@@ -140,7 +146,7 @@ void functionTester(
   typename SampleFunctional<dim, Function>::SmallVector correction;
   std::cout << "Old value: J(...) = " << J(start) << std::endl;
   for (size_t i = 1; i <= runs; ++i) {
-    correction = J.minimise(start);
+    correction = minimise(J, start);
     start += correction;
     if (i != runs)
       std::cout << "New value: J(...) = " << J(start) << std::endl;