diff --git a/dune/tectonic/spatial-solving/tnnmg/linesearchsolver.hh b/dune/tectonic/spatial-solving/tnnmg/linesearchsolver.hh index 587c81567d57ba7f49e90243018ba692b18c02e2..2b9e86dc06ea38c304d38d67c58679c87a485b37 100644 --- a/dune/tectonic/spatial-solving/tnnmg/linesearchsolver.hh +++ b/dune/tectonic/spatial-solving/tnnmg/linesearchsolver.hh @@ -16,25 +16,37 @@ class LineSearchSolver public: template<class Vector, class Functional, class BitVector> void operator()(Vector& x, const Functional& f, const BitVector& ignore) const { - x = 0; + x = 1.0; if (ignore) return; - /*Dune::Solvers::Interval<double> D; + Dune::Solvers::Interval<double> D; D = f.subDifferential(0); + + std::cout << "f.A " << f.quadraticPart() << " f.b " << f.linearPart() << std::endl; + + std::cout << D[0] << " " << D[1] << std::endl; + std::cout << "domain: " << f.domain()[0] << " " << f.domain()[1] << std::endl; + + if (D[1] > 0) // NOTE: Numerical instability can actually get us here return; -*/ + + if (almost_equal(f.domain()[0], f.domain()[1], 2)) { x = f.domain()[0]; + std::cout << "no interval: " << x << std::endl; return; } int bisectionsteps = 0; - const Bisection globalBisection(0.0, 1.0, 0.0, 0.0);; + const Bisection globalBisection; //(0.0, 1.0, 0.0, 0.0); - x = globalBisection.minimize(f, 0.0, 0.0, bisectionsteps); + x = globalBisection.minimize(f, f.scaling(), 0.0, bisectionsteps); + std::cout << "x: " << x << "scaling: " << f.scaling(); + x /= f.scaling(); + std::cout << "final x: " << x << std::endl; //x = f.domain().projectIn(x); } };