From d72348060bbaf8a7309dd47de4c4a8a4a9728ff5 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 9 Sep 2011 17:13:35 +0200
Subject: [PATCH] v need not be normalised

---
 src/mynonlinearity.cc | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/src/mynonlinearity.cc b/src/mynonlinearity.cc
index 4efdd603..3e9f9028 100644
--- a/src/mynonlinearity.cc
+++ b/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);
     }
   }
 
-- 
GitLab