Rescale for numerical stability

dune/tectonic/myblockproblem.hh

@@ -74,7 +74,11 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
 
   double computeDampingParameter(VectorType const &u,
                                  VectorType const &projected_v) const {
-    VectorType const v = projected_v;
+    VectorType v = projected_v;
+
+    double const vnorm = v.two_norm();
+
+    v /= vnorm; // Rescale for numerical stability
 
     VectorType tmp = problem.f;
     problem.A.mmv(u, tmp);
@@ -94,8 +98,8 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem {
         localA, localb, problem.phi, u, v);
 
     int bisectionsteps = 0;
-    Bisection bisection(0.0, 1.0, 1e-12, true, 0);            // TODO
-    return bisection.minimize(psi, 1.0, 0.0, bisectionsteps); // TODO
+    Bisection bisection(0.0, 1.0, 1e-12, true, 0);                      // TODO
+    return bisection.minimize(psi, vnorm, 0.0, bisectionsteps) / vnorm; // TODO
   }
 
   void assembleTruncate(VectorType const &u, Linearization &linearization,

src/TODO.org

-* check if v should be normalised by default
 * use nested iteration to obtain better iteration to start with
 * fix up octave bindings (low-pri)