make solver type configurable

src/root_solver.cpp

@@ -40,6 +40,7 @@ RGFlow<Root>::RGFlow()
    , convergence_tester(NULL)
    , initial_guesser(NULL)
    , precision_goal(1e-4)
+   , solver_type(Hybrid)
 {
 }
 
@@ -110,7 +111,7 @@ int RGFlow<Root>::find_root()
    gsl_multiroot_function f = {func, dimension, this};
 
    const gsl_multiroot_fsolver_type* solver_type
-      = gsl_multiroot_fsolver_hybrid; // @todo make configurable
+      = get_gsl_solver_type();
 
    gsl_multiroot_fsolver* solver
       = gsl_multiroot_fsolver_alloc(solver_type, dimension);
@@ -310,6 +311,21 @@ void RGFlow<Root>::set_model_parameters(const Eigen::ArrayXd& p)
    assert(offset == get_dimension());
 }
 
+const gsl_multiroot_fsolver_type* RGFlow<Root>::get_gsl_solver_type() const
+{
+   switch (solver_type) {
+   case Hybrid:  return gsl_multiroot_fsolver_hybrid;
+   case HybridS: return gsl_multiroot_fsolver_hybrids;
+   case Broyden: return gsl_multiroot_fsolver_broyden;
+   case Newton:  return gsl_multiroot_fsolver_dnewton;
+   default: {
+      std::stringstream sstr;
+      sstr << "RGFlow<Root>: Unknow solver type: " << solver_type;
+      throw SetupError(sstr.str());
+   }
+   }
+}
+
 /**
  * Print state of the root solver
  *

src/root_solver.hpp

@@ -51,6 +51,8 @@ class Root_model;
 template<>
 class RGFlow<Root> {
 public:
+   enum ERootFinder_t { Hybrid, HybridS, Broyden, Newton };
+
    RGFlow();
    ~RGFlow();
 
@@ -67,6 +69,8 @@ class RGFlow<Root> {
    void set_precision_goal(double);
    /// set initial guesser
    void set_initial_guesser(Initial_guesser<Root>*);
+   /// set solver type
+   void set_solver_type(ERootFinder_t);
    /// solves the boundary value problem
    void solve();
 
@@ -95,6 +99,7 @@ class RGFlow<Root> {
    Convergence_tester<Root>* convergence_tester; ///< the convergence tester
    Initial_guesser<Root>* initial_guesser;       ///< does initial guess
    double precision_goal;                ///< precision goal
+   ERootFinder_t solver_type;          ///< solver type
 
    static int func(const gsl_vector*, void*, gsl_vector*);
    static Eigen::ArrayXd to_eigen_array(const gsl_vector*);
@@ -109,6 +114,7 @@ class RGFlow<Root> {
    int find_root();                    ///< find root
    std::size_t get_dimension() const;  ///< calculate dimension
    unsigned int get_max_iterations() const; ///< returns max. number of iterations
+   const gsl_multiroot_fsolver_type* get_gsl_solver_type() const;
    void initial_guess();               ///< initial guess
    void print_state(const gsl_multiroot_fsolver*, std::size_t) const;
    void set_model_parameters(const Eigen::ArrayXd&); ///< set model parameters