From 2e1d937b6f556ae6d339238ece9eb3a4db92bfdf Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 6 Sep 2011 18:41:36 +0200
Subject: [PATCH] Handle case x = 0 in the modified gradient method

---
 src/bisection-example-flexible.cc | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/src/bisection-example-flexible.cc b/src/bisection-example-flexible.cc
index f4d40890..e8add65b 100644
--- a/src/bisection-example-flexible.cc
+++ b/src/bisection-example-flexible.cc
@@ -49,6 +49,8 @@ class SampleFunctional {
 
   SmallVector minimise(const SmallVector x, unsigned int iterations) const {
     SmallVector descDir = ModifiedGradient(x);
+    if (descDir == SmallVector(0.0))
+      return SmallVector(0.0);
 
     Dune::dverb << "Starting at x with J(x) = " << operator()(x) << std::endl;
     Dune::dverb << "Minimizing in direction w with dJ(x,w) = "
@@ -129,9 +131,16 @@ class SampleFunctional {
   }
 
   SmallVector ModifiedGradient(const SmallVector x) const {
-    if (x == SmallVector(0.0))
-      // TODO
-      DUNE_THROW(Dune::Exception, "The case x = 0 is not yet handled.");
+    if (x == SmallVector(0.0)) {
+      SmallVector d = SmoothGrad(x);
+      // Decline of the smooth part in the negative gradient direction
+      double smoothDecline = -(d * d);
+      double nonlinearDecline =
+          func_.rightDifferential(0.0) * d.two_norm(); // TODO: is this correct?
+      double combinedDecline = smoothDecline + nonlinearDecline;
+
+      return (combinedDecline < 0) ? d : SmallVector(0.0);
+    }
 
     SmallVector const pg = PlusGrad(x);
     SmallVector const mg = MinusGrad(x);
-- 
GitLab