From 2d0b9deec7649bebc6fadb05bf9778599624e7c0 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 14 Dec 2006 10:14:10 +0000
Subject: [PATCH] bugfix

[[Imported from SVN: r1090]]
---
 src/rodsolver.cc | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/rodsolver.cc b/src/rodsolver.cc
index 58382c5e1..76377e17a 100644
--- a/src/rodsolver.cc
+++ b/src/rodsolver.cc
@@ -217,7 +217,7 @@ void RodSolver<GridType>::solve()
         std::cout << "Correction: " << std::endl << corr << std::endl;
         
         printf("infinity norm of the correction: %g\n", corr.infinity_norm());
-        if (corr.infinity_norm() < 1e-7) {
+        if (corr.infinity_norm() < 1e-5) {
             std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
             break;
         }
@@ -256,7 +256,7 @@ void RodSolver<GridType>::solve()
         // Note that rhs = -g
         CorrectionType tmp(corr.size());
         tmp = 0;
-        hessianMatrix_->mmv(corr, tmp);
+        hessianMatrix_->umv(corr, tmp);
         double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
         
         std::cout << "Model decrease: " << modelDecrease 
@@ -269,7 +269,10 @@ void RodSolver<GridType>::solve()
 //             std::cout << "old energy: " << oldEnergy << "   new energy: " << energy << std::endl;
 //             exit(0);
         }
-        
+
+        if (std::abs(oldEnergy-energy) < 1e-12)
+            std::cout << "Suspecting rounding problems" << std::endl;
+
         // //////////////////////////////////////////////
         //   Check for acceptance of the step
         // //////////////////////////////////////////////
@@ -279,7 +282,8 @@ void RodSolver<GridType>::solve()
             x_ = newIterate;
             trustRegionRadius *= 2;
             
-        } else if ( (oldEnergy-energy) / modelDecrease > 0.01) {
+        } else if ( (oldEnergy-energy) / modelDecrease > 0.01
+                    || std::abs(oldEnergy-energy) < 1e-12) {
             // successful iteration
             x_ = newIterate;
             
-- 
GitLab