linesearchsolver.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_TECTONIC_LINESEARCHSOLVER_HH
#define DUNE_TECTONIC_LINESEARCHSOLVER_HH

#include <dune/tnnmg/problem-classes/bisection.hh>

#include "../../utils/almostequal.hh"

/**
 * \brief A local solver for scalar quadratic obstacle problems with nonlinearity
 * using bisection
 */
class LineSearchSolver
{
public:
  template<class Vector, class Functional, class BitVector>
  void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
    x = 1.0;

    if (ignore)
        return;

    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);

    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);
  }
};

#endif