v need not be normalised

src/mynonlinearity.cc

@@ -24,15 +24,22 @@ class MyNonlinearity {
   typedef SmallVector VectorType;
   typedef SmallMatrix MatrixType;
 
+  // directional subdifferential: at u on the line u + t*v
   void directionalSubDiff(VectorType u, VectorType v, Interval<double>& D) {
     if (u == SmallVector(0.0)) {
       D[0] = D[1] = func_.rightDifferential(0);
-    } else if (u * v > 0) {
-      D[1] = (v * u) * func_.rightDifferential(u.two_norm()) / u.two_norm();
-      D[0] = (v * u) * func_.leftDifferential(u.two_norm()) / u.two_norm();
+      return;
+    }
+    double const un = u.two_norm();
+    double const ndotp = (u * v) / (un * v.two_norm());
+    if (ndotp > 0) {
+      // If we had |v| = 1, this would be
+      // <f'_pm(|u|)*u/|u|,v>
+      D[1] = ndotp * func_.rightDifferential(un);
+      D[0] = ndotp * func_.leftDifferential(un);
     } else {
-      D[1] = (v * u) * func_.leftDifferential(u.two_norm()) / u.two_norm();
-      D[0] = (v * u) * func_.rightDifferential(u.two_norm()) / u.two_norm();
+      D[1] = ndotp * func_.leftDifferential(un);
+      D[0] = ndotp * func_.rightDifferential(un);
     }
   }