Have minimise() make multiple steps

src/myblockproblem.hh

@@ -32,6 +32,7 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
   class IterateObject;
 
   MyBlockProblem(MyConvexProblemType& problem) : problem(problem) {
+    // TODO: Is it clever to create a bisection here?
     bisection = Bisection(0.0, 1.0, 1e-15, true, 1e-14);
   };
 
@@ -118,10 +119,7 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject {
       Dune::SampleFunctional<block_size> localJ(*localA, localb, phi);
 
       LocalVectorType correction;
-      for (size_t i = 1; i <= 10; ++i) { // FIXME: hardcoded value
-        Dune::minimise(localJ, ui, correction);
-        ui += correction;
-      }
+      Dune::minimise(localJ, ui, 10); // FIXME: hardcoded value
     }
   }
 

src/samplefunctional.hh

@@ -118,8 +118,8 @@ template <int dimension> class SampleFunctional {
 };
 
 template <class Functional>
-void minimise(const Functional J, const typename Functional::SmallVector x,
-              typename Functional::SmallVector &corr,
+void minimise(const Functional J, typename Functional::SmallVector &x,
+              size_t steps = 1,
               Bisection const &bisection =
                   Bisection(0.0, // acceptError: Stop if the search interval has
                                  // become smaller than this number
@@ -130,59 +130,60 @@ void minimise(const Functional J, const typename Functional::SmallVector x,
                                    // minimization
 {
   typedef typename Functional::SmallVector SmallVector;
-  SmallVector descDir;
-  J.descentDirection(x, descDir);
 
-  if (descDir == SmallVector(0.0)) {
-    corr = SmallVector(0.0);
-    return;
-  }
+  for (size_t step = 0; step < steps; ++step) {
+    SmallVector descDir;
+    J.descentDirection(x, descDir);
 
-  // {{{ Construct a restriction of J to the line x + t * descDir
+    if (descDir == SmallVector(0.0))
+      return;
 
-  /* We have
+    // {{{ Construct a restriction of J to the line x + t * descDir
 
-     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>
+    /* We have
 
-     since A is symmetric.
-  */
-  SmallVector tmp;
+       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>
 
-  J.A.mv(descDir, tmp);                // Av
-  double const JRestA = tmp * descDir; // <Av,v>
+       since A is symmetric.
+    */
+    SmallVector tmp;
 
-  J.A.mv(x, tmp);                              // Au
-  double const JRestb = (J.b - tmp) * descDir; // <b-Au,v>
+    J.A.mv(descDir, tmp);                // Av
+    double const JRestA = tmp * descDir; // <Av,v>
 
-  typedef typename Functional::NonlinearityType MyNonlinearityType;
-  MyNonlinearityType phi = J.phi;
-  typedef DirectionalConvexFunction<MyNonlinearityType>
-  MyDirectionalConvexFunctionType;
-  // FIXME: We cannot pass J.phi directly because the constructor
-  // does not allow for constant arguments
-  MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir);
-  // }}}
+    J.A.mv(x, tmp);                              // Au
+    double const JRestb = (J.b - tmp) * descDir; // <b-Au,v>
 
-  { // Debug
-    Interval<double> D;
-    JRest.subDiff(0, D);
+    typedef typename Functional::NonlinearityType MyNonlinearityType;
+    MyNonlinearityType phi = J.phi;
+    typedef DirectionalConvexFunction<MyNonlinearityType>
+    MyDirectionalConvexFunctionType;
+    // FIXME: We cannot pass J.phi directly because the constructor
+    // does not allow for constant arguments
+    MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir);
+    // }}}
 
-    dverb
-        << "## Directional derivative (as per subdifferential of restriction): "
-        << D[1] << " (coordinates of the restriction)" << std::endl;
-    assert(D[1] <=
-           0); // We should not be minimising in this direction otherwise
-  }
+    { // Debug
+      Interval<double> D;
+      JRest.subDiff(0, D);
 
-  int count;
-  // FIXME: The value of x_old should not matter if the factor is 1.0, correct?
-  double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
-  dverb << "Number of iterations in the bisection method: " << count
-        << std::endl;
-  ;
+      dverb << "## Directional derivative (as per subdifferential of "
+               "restriction): " << D[1] << " (coordinates of the restriction)"
+            << std::endl;
+      assert(D[1] <=
+             0); // We should not be minimising in this direction otherwise
+    }
+
+    int count;
+    // FIXME: The value of x_old should not matter if the factor is 1.0,
+    // correct?
+    double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
+    dverb << "Number of iterations in the bisection method: " << count
+          << std::endl;
+    ;
 
-  corr = descDir;
-  corr *= stepsize;
+    x.axpy(stepsize, descDir);
+  }
 }
 }
 #endif

src/test-gradient-method.cc

@@ -15,14 +15,8 @@ template <int dim>
 double functionTester(Dune::SampleFunctional<dim> J,
                       typename Dune::SampleFunctional<dim>::SmallVector &start,
                       size_t runs) {
-  typename Dune::SampleFunctional<dim>::SmallVector correction;
   std::cout << "Old value: J(...) = " << J(start) << std::endl;
-  for (size_t i = 1; i <= runs; ++i) {
-    Dune::minimise(J, start, correction);
-    start += correction;
-    if (i != runs)
-      std::cout << "New value: J(...) = " << J(start) << std::endl;
-  }
+  Dune::minimise(J, start, runs);
   double const final = J(start);
   std::cout << "Final value J(...) = " << final << std::endl;
   return final;