root_finder.hpp

Go to the documentation of this file.

// ====================================================================
// This file is part of FlexibleSUSY.
//
// FlexibleSUSY is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// FlexibleSUSY is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FlexibleSUSY.  If not, see
// <http://www.gnu.org/licenses/>.
// ====================================================================
 
#ifndef ROOT_FINDER_H
#define ROOT_FINDER_H
 
#include <algorithm>
#include <cmath>
#include <string>
#include <utility>
#include <Eigen/Core>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
 
#include "logger.hpp"
#include "error.hpp"
#include "ewsb_solver.hpp"
#include "gsl_multiroot_fsolver.hpp"
#include "gsl_utils.hpp"
#include "gsl_vector.hpp"
 
namespace flexiblesusy {
 
template <std::size_t dimension>
class Root_finder : public EWSB_solver {
public:
   using Vector_t = Eigen::Matrix<double,dimension,1>;
   using Function_t = std::function<Vector_t(const Vector_t&)>;
   enum Solver_type { GSLHybrid, GSLHybridS, GSLBroyden, GSLNewton };
 
   Root_finder() = default;
   Root_finder(const Function_t&, std::size_t, double, Solver_type solver_type_ = GSLHybrid);
   virtual ~Root_finder() = default;
 
   void set_function(const Function_t& f) { function = f; }
   void set_precision(double p) { precision = p; }
   void set_max_iterations(std::size_t n) { max_iterations = n; }
   void set_solver_type(Solver_type t) { solver_type = t; }
   int find_root(const Vector_t&);
 
   // EWSB_solver interface methods
   virtual std::string name() const override { return std::string("Root_finder<") + solver_type_name() + ">"; }
   virtual int solve(const Eigen::VectorXd&) override;
   virtual Eigen::VectorXd get_solution() const override { return root; }
 
private:
   std::size_t max_iterations{100};    
   double precision{1.e-2};            
   Vector_t root{Vector_t::Zero()};    
   Function_t function{nullptr};       
   Solver_type solver_type{GSLHybrid}; 
 
   const char* solver_type_name() const;
   const gsl_multiroot_fsolver_type* solver_type_to_gsl_pointer() const;
   static int gsl_function(const gsl_vector*, void*, gsl_vector*);
 
   static bool is_finite(const Vector_t& v) {
      return std::any_of(v.data(), v.data() + v.size(),
                         [](double x) { return std::isfinite(x); });
   }
};
 
template <std::size_t dimension>
Root_finder<dimension>::Root_finder(
   const Function_t& function_,
   std::size_t max_iterations_,
   double precision_,
   Solver_type solver_type_
)
   : max_iterations(max_iterations_)
   , precision(precision_)
   , function(function_)
   , solver_type(solver_type_)
{
}
 
template <std::size_t dimension>
int Root_finder<dimension>::find_root(const Vector_t& start)
{
   if (!function)
      throw SetupError("Root_finder: function not callable");
 
   int status;
   std::size_t iter = 0;
   void* parameters = &function;
   gsl_multiroot_function f = {gsl_function, dimension, parameters};
   GSL_multiroot_fsolver solver(solver_type_to_gsl_pointer(), dimension, &f, to_GSL_vector(start));
 
#ifndef ENABLE_DEBUG
   gsl_set_error_handler_off();
#endif
 
   do {
      iter++;
      status = solver.iterate();
      solver.print_state(iter);
 
      if (status)   // check if solver is stuck
         break;
 
      status = solver.test_residual(precision);
   } while (status == GSL_CONTINUE && iter < max_iterations);
 
   VERBOSE_MSG("\t\t\tRoot_finder status = " << gsl_strerror(status));
 
   root = to_eigen_vector<dimension>(solver.get_root().raw());
 
   return status;
}
 
template <std::size_t dimension>
int Root_finder<dimension>::gsl_function(const gsl_vector* x, void* params, gsl_vector* f)
{
   if (!flexiblesusy::is_finite(x)) {
      gsl_vector_set_all(f, std::numeric_limits<double>::max());
      return GSL_EDOM;
   }
 
   Function_t* fun = static_cast<Function_t*>(params);
   int status = GSL_SUCCESS;
   const Vector_t arg(to_eigen_vector<dimension>(x));
   Vector_t result;
   result.setConstant(std::numeric_limits<double>::max());
 
   try {
      result = (*fun)(arg);
      status = is_finite(result) ? GSL_SUCCESS : GSL_EDOM;
   } catch (const flexiblesusy::Error&) {
      status = GSL_EDOM;
   }
 
   // copy result -> f
   for (std::size_t i = 0; i < dimension; i++) {
      gsl_vector_set(f, i, result(i));
   }
 
   return status;
}
 
template <std::size_t dimension>
const char* Root_finder<dimension>::solver_type_name() const
{
   switch (solver_type) {
   case GSLHybrid : return "GSLHybrid";
   case GSLHybridS: return "GSLHybridS";
   case GSLBroyden: return "GSLBroyden";
   case GSLNewton : return "GSLNewton";
   default:
      throw SetupError("Unknown root solver type: "
                       + std::to_string(solver_type));
   }
 
   return "unknown";
}
 
template <std::size_t dimension>
const gsl_multiroot_fsolver_type* Root_finder<dimension>::solver_type_to_gsl_pointer() const
{
   switch (solver_type) {
   case GSLHybrid : return gsl_multiroot_fsolver_hybrid;
   case GSLHybridS: return gsl_multiroot_fsolver_hybrids;
   case GSLBroyden: return gsl_multiroot_fsolver_broyden;
   case GSLNewton : return gsl_multiroot_fsolver_dnewton;
   default:
      throw SetupError("Unknown root solver type: "
                       + std::to_string(solver_type));
   }
 
   return nullptr;
}
 
template <std::size_t dimension>
int Root_finder<dimension>::solve(const Eigen::VectorXd& start)
{
   return (find_root(start) == GSL_SUCCESS ?
           EWSB_solver::SUCCESS : EWSB_solver::FAIL);
}
 
} // namespace flexiblesusy
 
#endif