Modularise minimise()

dune/tectonic/samplefunctional.hh

@@ -133,6 +133,87 @@ template <int dim> class SampleFunctional {
   }
 };
 
+template <class Functional>
+void descentMinimisation(Functional const &J,
+                         typename Functional::SmallVector &x,
+                         typename Functional::SmallVector const &descDir,
+                         Bisection const &bisection) {
+  typedef typename Functional::SmallVector SmallVector;
+  typedef typename Functional::NonlinearityType LocalNonlinearityType;
+
+  // {{{ Construct a restriction of J to the line x + t * descDir
+
+  /* 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 = J.b;               //  b
+  J.A.mmv(x, tmp);                     //  b-Au
+  double const JRestb = tmp * descDir; // <b-Au,v>
+
+  J.A.mv(descDir, tmp);                //  Av
+  double const JRestA = tmp * descDir; // <Av,v>
+
+  MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
+      JRestA, JRestb, *J.phi, x, descDir);
+  // }}}
+
+  { // Debug
+    Interval<double> D;
+    JRest.subDiff(0, D);
+
+    dverb
+        << "## Directional derivative (as per subdifferential of restriction): "
+        << D[1] << " (coordinates of the restriction)" << std::endl;
+    /*
+      It is possible that this differs quite a lot from the
+      directional derivative computed in the descentDirection()
+      method:
+
+      If phi is nonsmooth at x, so that the directional
+      derivatives jump, and |x| is computed to be too small or too
+      large globally or locally, the locally computed
+      subdifferential and the globally computed subdifferential
+      will no longer coincide!
+
+      The assertion D[1] <= 0 may thus fail.
+    */
+  }
+
+  int count;
+  double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
+  dverb << "Number of iterations in the bisection method: " << count
+        << std::endl;
+  ;
+
+  x.axpy(stepsize, descDir);
+}
+
+template <class Functional>
+void tangentialMinimisation(Functional const &J,
+                            typename Functional::SmallVector &x,
+                            typename Functional::SmallVector const &descDir,
+                            Bisection const &bisection) {
+  typedef typename Functional::NonlinearityType LocalNonlinearityType;
+  typedef typename Functional::SmallVector SmallVector;
+
+  CircularConvexFunction<LocalNonlinearityType> const JRest(J.A, J.b, *J.phi, x,
+                                                            descDir);
+
+  int count;
+  double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
+  dverb << "Number of iterations in the bisection method: " << count
+        << std::endl;
+  ;
+
+  // Since x is used in the computation of the rhs, do not write to it directly
+  SmallVector tmp;
+  JRest.cartesian(stepsize, tmp);
+  x = tmp;
+}
+
 template <class Functional>
 void minimise(Functional const &J, typename Functional::SmallVector &x,
               size_t steps, Bisection const &bisection) {
@@ -145,75 +226,13 @@ void minimise(Functional const &J, typename Functional::SmallVector &x,
     if (descDir == SmallVector(0.0))
       return;
 
-    typedef typename Functional::NonlinearityType LocalNonlinearityType;
     if (linesearchp) {
-      // {{{ Construct a restriction of J to the line x + t * descDir
-
-      /* 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 = J.b;               //  b
-      J.A.mmv(x, tmp);                     //  b-Au
-      double const JRestb = tmp * descDir; // <b-Au,v>
-
-      J.A.mv(descDir, tmp);                //  Av
-      double const JRestA = tmp * descDir; // <Av,v>
-
-      MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
-          JRestA, JRestb, *J.phi, x, descDir);
-      // }}}
-
-      { // Debug
-        Interval<double> D;
-        JRest.subDiff(0, D);
-
-        dverb << "## Directional derivative (as per subdifferential of "
-                 "restriction): " << D[1] << " (coordinates of the restriction)"
-              << std::endl;
-        /*
-          It is possible that this differs quite a lot from the
-          directional derivative computed in the descentDirection()
-          method:
-
-          If phi is nonsmooth at x, so that the directional
-          derivatives jump, and |x| is computed to be too small or too
-          large globally or locally, the locally computed
-          subdifferential and the globally computed subdifferential
-          will no longer coincide!
-
-          The assertion D[1] <= 0 may thus fail.
-        */
-      }
-
-      int count;
-      double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
-      dverb << "Number of iterations in the bisection method: " << count
-            << std::endl;
-      ;
-
-      x.axpy(stepsize, descDir);
+      descentMinimisation(J, x, descDir, bisection);
     } else {
       Bisection slowBisection(bisection);
       slowBisection.setFastQuadratic(false);
 
-      CircularConvexFunction<LocalNonlinearityType> const JRest(
-          J.A, J.b, *J.phi, x, descDir);
-
-      int count;
-      double const stepsize = slowBisection.minimize(JRest, 0.0, 1.0, count);
-      dverb << "Number of iterations in the bisection method: " << count
-            << std::endl;
-      ;
-
-      // Since x is used in the computation of the rhs, do not write to it
-      // directly
-      SmallVector tmp;
-      JRest.cartesian(stepsize, tmp);
-      x = tmp;
+      tangentialMinimisation(J, x, descDir, slowBisection);
     }
   }
 }