SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/RootFinding - GslMultiRoot.hpp Hit Total Coverage
Commit: 2db722c93a8e9b106e406b439b79c8e05c2057fb Lines: 4 5 80.0 %
Date: 2021-03-03 22:01:00
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cmath>
       8             : #include <cstddef>
       9             : #include <gsl/gsl_errno.h>
      10             : #include <gsl/gsl_matrix_double.h>
      11             : #include <gsl/gsl_multiroots.h>
      12             : #include <gsl/gsl_vector_double.h>
      13             : #include <ostream>
      14             : #include <string>
      15             : 
      16             : #include "Parallel/Printf.hpp"
      17             : #include "Utilities/ErrorHandling/Error.hpp"
      18             : #include "Utilities/ErrorHandling/Exceptions.hpp"
      19             : #include "Utilities/Gsl.hpp"
      20             : #include "Utilities/Requires.hpp"
      21             : #include "Utilities/TypeTraits/CreateIsCallable.hpp"
      22             : 
      23             : namespace RootFinder {
      24             : namespace gsl_multiroot_detail {
      25             : template <size_t Dim, typename Solver>
      26             : void print_state(const size_t iteration_number, const Solver* const solver,
      27             :                  const bool print_header = false) noexcept {
      28             :   if (print_header) {
      29             :     Parallel::printf("Iter\t");
      30             :     for (size_t i = 0; i < Dim; ++i) {
      31             :       Parallel::printf(" x[%u]\t", i);
      32             :     }
      33             :     for (size_t i = 0; i < Dim; ++i) {
      34             :       Parallel::printf(" f[%u]\t", i);
      35             :     }
      36             :     Parallel::printf("\n");
      37             :   }
      38             : 
      39             :   Parallel::printf("%u\t", iteration_number);
      40             :   for (size_t i = 0; i < Dim; ++i) {
      41             :     Parallel::printf("%3.4f  ", gsl_vector_get(solver->x, i));
      42             :   }
      43             :   for (size_t i = 0; i < Dim; ++i) {
      44             :     Parallel::printf("%1.3e  ", gsl_vector_get(solver->f, i));
      45             :   }
      46             :   Parallel::printf("\n");
      47             : }
      48             : 
      49             : template <size_t Dim>
      50             : std::array<double, Dim> gsl_to_std_array(const gsl_vector* const x) noexcept {
      51             :   std::array<double, Dim> input_as_std_array{};
      52             :   for (size_t i = 0; i < Dim; i++) {
      53             :     gsl::at(input_as_std_array, i) = gsl_vector_get(x, i);
      54             :   }
      55             :   return input_as_std_array;
      56             : }
      57             : 
      58             : template <size_t Dim>
      59             : void gsl_vector_set_with_std_array(
      60             :     gsl_vector* const func,
      61             :     const std::array<double, Dim>& result_as_std_array) noexcept {
      62             :   for (size_t i = 0; i < Dim; i++) {
      63             :     gsl_vector_set(func, i, gsl::at(result_as_std_array, i));
      64             :   }
      65             : }
      66             : 
      67             : template <size_t Dim>
      68             : void gsl_matrix_set_with_std_array(
      69             :     gsl_matrix* const matrix,
      70             :     const std::array<std::array<double, Dim>, Dim>& matrix_array) noexcept {
      71             :   for (size_t i = 0; i < Dim; i++) {
      72             :     for (size_t j = 0; j < Dim; j++) {
      73             :       gsl_matrix_set(matrix, i, j, gsl::at(gsl::at(matrix_array, i), j));
      74             :     }
      75             :   }
      76             : }
      77             : 
      78             : // The gsl_multiroot_function_fdf expects its functions to be of the form
      79             : // int (* f) (const gsl_vector * x, void * params, gsl_vector * f).
      80             : // However, we would like to be able to perform rootfinding on functions
      81             : // of the form std::array<double, Dim> f(const std::array<double, Dim>& x).
      82             : // So, we pass the function wrapper below to gsl_multiroot_function_fdf.
      83             : // In the gsl documentation the third parameter is refered to as "void* params",
      84             : // referring to the parameters that select out a particular function out of a
      85             : // family of possible functions described in a class, but we instead pass the
      86             : // pointer to the entire function object itself here. The type of the function
      87             : // object is passed through with the Function template parameter.
      88             : template <size_t Dim, typename Function>
      89             : int gsl_multirootfunctionfdf_wrapper_f(const gsl_vector* const x,
      90             :                                        void* const untyped_function_object,
      91             :                                        gsl_vector* const f_of_x) noexcept {
      92             :   const auto function_object =
      93             :       static_cast<const Function*>(untyped_function_object);
      94             :   gsl_vector_set_with_std_array(
      95             :       f_of_x, function_object->operator()(gsl_to_std_array<Dim>(x)));
      96             :   return GSL_SUCCESS;
      97             : }
      98             : 
      99             : template <size_t Dim, typename Function>
     100             : int gsl_multirootfunctionfdf_wrapper_df(const gsl_vector* const x,
     101             :                                         void* const untyped_function_object,
     102             :                                         gsl_matrix* const jacobian) noexcept {
     103             :   const auto function_object =
     104             :       static_cast<const Function*>(untyped_function_object);
     105             :   gsl_matrix_set_with_std_array(
     106             :       jacobian, function_object->jacobian(gsl_to_std_array<Dim>(x)));
     107             : 
     108             :   return GSL_SUCCESS;
     109             : }
     110             : 
     111             : template <size_t Dim, typename Function>
     112             : int gsl_multirootfunctionfdf_wrapper_fdf(const gsl_vector* const x,
     113             :                                          void* const untyped_function_object,
     114             :                                          gsl_vector* const f_of_x,
     115             :                                          gsl_matrix* const jacobian) noexcept {
     116             :   const auto function_object =
     117             :       static_cast<const Function*>(untyped_function_object);
     118             :   const std::array<double, Dim> x_as_std_array = gsl_to_std_array<Dim>(x);
     119             :   gsl_vector_set_with_std_array(f_of_x,
     120             :                                 function_object->operator()(x_as_std_array));
     121             :   gsl_matrix_set_with_std_array(jacobian,
     122             :                                 function_object->jacobian(x_as_std_array));
     123             :   return GSL_SUCCESS;
     124             : }
     125             : 
     126             : CREATE_IS_CALLABLE(jacobian)
     127             : CREATE_IS_CALLABLE_V(jacobian)
     128             : }  // namespace gsl_multiroot_detail
     129             : 
     130             : /*!
     131             :  *  \ingroup NumericalAlgorithmsGroup
     132             :  *  \brief The different options for the verbosity of gsl_multiroot.
     133             :  */
     134           1 : enum class Verbosity {
     135             :   /// Do not print anything
     136             :   Silent,
     137             :   /// Print only "success" or "failed" on termination
     138             :   Quiet,
     139             :   /// Print final functions values on termination
     140             :   Verbose,
     141             :   /// Print function values on every iteration
     142             :   Debug
     143             : };
     144             : 
     145             : std::ostream& operator<<(std::ostream& /*os*/,
     146             :                          const Verbosity& /*verbosity*/) noexcept;
     147             : 
     148             : /*!
     149             :  *  \ingroup NumericalAlgorithmsGroup
     150             :  *  \brief The different options for the rootfinding method of gsl_multiroot.
     151             :  *
     152             :  *  This enum is for setting the method used the rootfinder.
     153             :  *  The precise method used by the gsl rootfinder depends on whether or not the
     154             :  *  function passed to it has a callable `jacobian` member function. In the
     155             :  *  case where it doesn't, the jacobian is approximated with a finite
     156             :  *  difference. For example, if the Method specified is Hybrid, gsl will use
     157             :  *  the gsl_multiroot_fdfsolver_hybridj method in the case where a `jacobian`
     158             :  *  is provided, and gsl_multiroot_fsolver_hybrid in the case where one isn't.
     159             :  *  See
     160             :  *  [GSL's documentation for multidimensional
     161             :  *  rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
     162             :  *  for information on the different methods.
     163             :  *  \note gsl does not provide a finite difference version for the modified
     164             :  *  Newton method (gsl_multiroot_fdfsolver_gnewton). In the case where a
     165             :  *  jacobian is not provided the method used will be a non-modified Newton
     166             :  *  method.
     167             :  *
     168             :  */
     169           1 : enum class Method {
     170             :   /// Hybrid of Newton's method along with following the gradient direction.
     171             :   /// \note Sometimes Hybrids works only with the Absolute stopping condition.
     172             :   Hybrids,
     173             :   /// "Unscaled version of Hybrids that uses a spherical trust region," see
     174             :   /// GSL documentation for more details.
     175             :   Hybrid,
     176             :   /// If an analytic jacobian is provided, gsl uses a modification of Newton's
     177             :   /// method to improve global convergence. Uses vanilla Newton's method if no
     178             :   /// jacobian is provided.
     179             :   Newton
     180             : };
     181             : 
     182             : std::ostream& operator<<(std::ostream& /*os*/,
     183             :                          const Method& /*method*/) noexcept;
     184             : 
     185             : /*!
     186             :  *  \ingroup NumericalAlgorithmsGroup
     187             :  *  \brief The different options for the convergence criterion of gsl_multiroot.
     188             :  *
     189             :  *  See
     190             :  *  [GSL's documentation for multidimensional
     191             :  *  rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
     192             :  *  for information on the different stopping conditions.
     193             :  */
     194           1 : enum class StoppingCondition {
     195             :   /// See GSL documentation for gsl_multiroot_test_delta.
     196             :   AbsoluteAndRelative,
     197             :   /// See GSL documentation for gsl_multiroot_test_residual.
     198             :   Absolute
     199             : };
     200             : 
     201             : std::ostream& operator<<(std::ostream& /*os*/,
     202             :                          const StoppingCondition& /*condition*/) noexcept;
     203             : 
     204             : namespace gsl_multiroot_detail {
     205             : template <typename SolverType, typename SolverAlloc, typename SolverSet,
     206             :           typename SolverIterate, typename SolverFree, size_t Dim,
     207             :           typename Function>
     208             : std::array<double, Dim> gsl_multiroot_impl(
     209             :     Function& f, const std::array<double, Dim>& initial_guess,
     210             :     const double absolute_tolerance, const size_t maximum_iterations,
     211             :     const double relative_tolerance, const Verbosity verbosity,
     212             :     const double maximum_absolute_tolerance, const Method method,
     213             :     const SolverType solver_type, const StoppingCondition condition,
     214             :     const SolverAlloc solver_alloc, const SolverSet solver_set,
     215             :     const SolverIterate solver_iterate, const SolverFree solver_free) {
     216             :   // Check for valid stopping condition:
     217             :   if (UNLIKELY(condition != StoppingCondition::AbsoluteAndRelative and
     218             :                condition != StoppingCondition::Absolute)) {
     219             :     ERROR(
     220             :         "Invalid stopping condition. Has to be either AbsoluteAndRelative"
     221             :         "or Absolute.");
     222             :   }
     223             : 
     224             :   // Supply gsl_root with the initial guess:
     225             :   gsl_vector* const gsl_root = gsl_vector_alloc(Dim);
     226             :   gsl_vector_set_with_std_array(gsl_root, initial_guess);
     227             :   auto* const solver = solver_alloc(solver_type, Dim);
     228             :   solver_set(solver, &f, gsl_root);
     229             : 
     230             :   // Take iterations:
     231             :   int status;
     232             :   size_t iteration_number = 0;
     233             :   do {
     234             :     if (UNLIKELY(verbosity == Verbosity::Debug)) {
     235             :       print_state<Dim>(iteration_number, solver, iteration_number == 0);
     236             :     }
     237             :     iteration_number++;
     238             :     status = solver_iterate(solver);
     239             :     // Check if solver is stuck
     240             :     if (UNLIKELY(status == GSL_ENOPROG)) {
     241             :       if (UNLIKELY(verbosity == Verbosity::Debug)) {
     242             :         Parallel::printf(
     243             :             "The iteration is not making any progress, preventing the "
     244             :             "algorithm from continuing.");
     245             :       }
     246             :       break;
     247             :     }
     248             :     if (condition == StoppingCondition::AbsoluteAndRelative) {
     249             :       status = gsl_multiroot_test_delta(solver->dx, solver->x,
     250             :                                         absolute_tolerance, relative_tolerance);
     251             :     } else {  // condition is StoppingCondition::Absolute
     252             :       // NOTE: Sometimes hybridsj works only with the test_residual condition
     253             :       status = gsl_multiroot_test_residual(solver->f, absolute_tolerance);
     254             :     }
     255             :   } while (status == GSL_CONTINUE and iteration_number < maximum_iterations);
     256             :   if (UNLIKELY(verbosity == Verbosity::Verbose or
     257             :                verbosity == Verbosity::Debug)) {
     258             :     Parallel::printf("Finished iterating:\n");
     259             :     print_state<Dim>(iteration_number, solver, verbosity == Verbosity::Verbose);
     260             :   }
     261             :   bool success = (status == GSL_SUCCESS);
     262             :   if (UNLIKELY(verbosity != Verbosity::Silent)) {
     263             :     Parallel::printf("\n");
     264             :     if (not success) {
     265             :       const std::string ascii_divider = std::string(70, '#');
     266             :       const std::string failure_message =
     267             :           ascii_divider + "\n\t\tWARNING: Root Finding FAILED\n" +
     268             :           ascii_divider;
     269             :       Parallel::printf("%s\n", failure_message);
     270             :     } else {
     271             :       Parallel::printf("Root finder converged.\n");
     272             :     }
     273             :   }
     274             :   // If maximum_absolute_tolerance is given, return success = true
     275             :   // as long as maximum_absolute_tolerance is achieved even if the
     276             :   // root finder doesn't converge.
     277             :   bool success_with_tolerance = true;
     278             :   bool failed_root_is_forgiven = false;
     279             :   if (not success and maximum_absolute_tolerance > 0.0) {
     280             :     for (size_t i = 0; i < Dim; ++i) {
     281             :       if (fabs(gsl_vector_get(solver->f, i)) > maximum_absolute_tolerance) {
     282             :         success_with_tolerance = false;
     283             :       }
     284             :     }
     285             :     failed_root_is_forgiven = success_with_tolerance;
     286             :   }
     287             :   if (success_with_tolerance and maximum_absolute_tolerance > 0.0) {
     288             :     success = true;
     289             :   }
     290             :   if (UNLIKELY(failed_root_is_forgiven and (verbosity == Verbosity::Verbose or
     291             :                                             verbosity == Verbosity::Debug))) {
     292             :     Parallel::printf(
     293             :         "The failed root was forgiven as each component was found to be under "
     294             :         "maximum_absolute_tolerance %f",
     295             :         maximum_absolute_tolerance);
     296             :   }
     297             : 
     298             :   if (UNLIKELY(not success)) {
     299             :     std::stringstream error_message;
     300             :     error_message << "The root find failed and was not forgiven. An exception "
     301             :                      "has been thrown.\n"
     302             :                   << "The gsl error returned is: " << gsl_strerror(status)
     303             :                   << "\n"
     304             :                   << "Verbosity: " << verbosity << "\n"
     305             :                   << "Method: " << method << "\n"
     306             :                   << "StoppingCondition: " << condition << "\n"
     307             :                   << "Maximum absolute tolerance: "
     308             :                   << maximum_absolute_tolerance << "\n"
     309             :                   << "Absolute tolerance: " << absolute_tolerance << "\n"
     310             :                   << "Relative tolerance: " << relative_tolerance << "\n"
     311             :                   << "Maximum number of iterations: " << maximum_iterations
     312             :                   << "\n"
     313             :                   << "Number of iterations reached: " << iteration_number
     314             :                   << "\n"
     315             :                   << "The last value of f in the root solver is:\n";
     316             :     for (size_t i = 0; i < Dim; i++) {
     317             :       error_message << gsl_vector_get(solver->f, i) << "\n";
     318             :     }
     319             :     error_message << "The last value of x in the root solver is:\n";
     320             :     for (size_t i = 0; i < Dim; i++) {
     321             :       error_message << gsl_vector_get(solver->x, i) << "\n";
     322             :     }
     323             :     error_message << "The last value of dx in the root solver is:\n";
     324             :     for (size_t i = 0; i < Dim; i++) {
     325             :       error_message << gsl_vector_get(solver->dx, i) << "\n";
     326             :     }
     327             : 
     328             :     if (UNLIKELY(verbosity == Verbosity::Debug)) {
     329             :       Parallel::printf("Error: %s\n", gsl_strerror(status));
     330             :     if (iteration_number >= maximum_iterations) {
     331             :       Parallel::printf(
     332             :           "The number of iterations (%d) has reached the maximum number of "
     333             :           "iterations (%d)\n",
     334             :           iteration_number, maximum_iterations);
     335             :     } else {
     336             :       Parallel::printf(
     337             :           "The number of iterations (%d) failed to reach the maximum number of "
     338             :           "iterations (%d)\n",
     339             :           iteration_number, maximum_iterations);
     340             :     }
     341             :     }
     342             :     throw convergence_error(error_message.str());
     343             :   }
     344             : 
     345             :   // Store the converged root in result
     346             :   std::array<double, Dim> result = gsl_to_std_array<Dim>(solver->x);
     347             :   solver_free(solver);
     348             :   gsl_vector_free(gsl_root);
     349             : 
     350             :   return result;
     351             : }
     352             : 
     353             : void print_rootfinding_parameters(Method method, double absolute_tolerance,
     354             :                                   double relative_tolerance,
     355             :                                   double maximum_absolute_tolerance,
     356             :                                   StoppingCondition condition) noexcept;
     357             : }  // namespace gsl_multiroot_detail
     358             : 
     359             : // @{
     360             : /*!
     361             :  * \ingroup NumericalAlgorithmsGroup
     362             :  * \brief A multidimensional root finder supporting Newton and Hybrid
     363             :  * methods, as well as modified methods based on these.
     364             :  *
     365             :  * This root finder accepts function objects with and without a callable
     366             :  * `jacobian` member function. This member function both accepts and
     367             :  * returns a `std::array<double, Dim>`, the dimension of the domain and range
     368             :  * of the function the root find is being performed on. Whether the jacobian
     369             :  * is provided determines the details of the implementation of the
     370             :  * root-finding method that is selected by the user using the Method enum.
     371             :  * That is, whether the jacobian is computed analytically via the `jacobian`
     372             :  * member function, or whether the jacobian is computed numerically via a
     373             :  * finite difference approximation.
     374             :  * \note GSL does not provide a finite difference version of its modified
     375             :  * Newton method, so the unmodified one is used instead when the user
     376             :  * uses the Method::Newton method.
     377             :  *
     378             :  * The user can select one of two possible criteria for convergence,
     379             :  * StoppingCondition::Absolute, where the sum of the absolute values of the
     380             :  * components of the residual vector f are compared against the value
     381             :  * provided to `absolute_tolerance`, and
     382             :  * StoppingCondition::AbsoluteAndRelative, where the size of the most recent
     383             :  * step taken in the root-finding iteration is compared against
     384             :  * `absolute_tolerance` + `relative_tolerance` * |x_i|, for each component.
     385             :  * In either case, a `maximum_absolute_tolerance` may be specified if the user
     386             :  * anticipates that the convergence criterion specified with StoppingCondition
     387             :  * will be too strict for a few points out of a population of points found with
     388             :  * a sequence of root finds.
     389             :  *
     390             :  * See
     391             :  * [GSL's documentation for multidimensional
     392             :  * rootfinding](https://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html)
     393             :  * for reference.
     394             :  *
     395             :  * \param func Function whose root is to be found.
     396             :  * \param initial_guess Contains initial guess.
     397             :  * \param absolute_tolerance The absolute tolerance.
     398             :  * \param maximum_iterations The maximum number of iterations.
     399             :  * \param relative_tolerance The relative tolerance.
     400             :  * \param verbosity Whether to print diagnostic messages.
     401             :  * \param maximum_absolute_tolerance Acceptable absolute tolerance when
     402             :  *                                   root finder doesn't converge.
     403             :  *                                   You may wish to use this if there
     404             :  *                                   are only a few "problematic" points where
     405             :  *                                   it is difficult to do a precise root find.
     406             :  * \param method The method to use. See the documentation for the Method enum.
     407             :  * \param condition The convergence condition to use. See the documentation
     408             :  *                                   for the StoppingCondition enum.
     409             :  */
     410             : template <size_t Dim, typename Function,
     411             :           Requires<gsl_multiroot_detail::is_jacobian_callable_v<
     412             :               Function, std::array<double, Dim>>> = nullptr>
     413           1 : std::array<double, Dim> gsl_multiroot(
     414             :     const Function& func, const std::array<double, Dim>& initial_guess,
     415             :     const double absolute_tolerance, const size_t maximum_iterations,
     416             :     const double relative_tolerance = 0.0,
     417             :     const Verbosity verbosity = Verbosity::Silent,
     418             :     const double maximum_absolute_tolerance = 0.0,
     419             :     const Method method = Method::Newton,
     420             :     const StoppingCondition condition = StoppingCondition::Absolute) {
     421             :   gsl_multiroot_function_fdf gsl_func = {
     422             :       &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
     423             :       &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_df<Dim, Function>,
     424             :       &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_fdf<Dim,
     425             :                                                                   Function>,
     426             :       Dim, const_cast<Function*>(&func)};  //NOLINT
     427             : 
     428             :   // Set up method for solver:
     429             :   const gsl_multiroot_fdfsolver_type* solver_type;
     430             :   if (method == Method::Newton) {
     431             :     solver_type = gsl_multiroot_fdfsolver_gnewton;
     432             :   } else if (method == Method::Hybrids) {
     433             :     solver_type = gsl_multiroot_fdfsolver_hybridsj;
     434             :   } else if (method == Method::Hybrid) {
     435             :     solver_type = gsl_multiroot_fdfsolver_hybridj;
     436             :   } else {
     437             :     ERROR(
     438             :         "Invalid method. Has to be one of Newton, Hybrids or "
     439             :         "Hybrid.");
     440             :   }
     441             :   // Print initial parameters
     442             :   if (UNLIKELY(verbosity == Verbosity::Verbose or
     443             :                verbosity == Verbosity::Debug)) {
     444             :     gsl_multiroot_detail::print_rootfinding_parameters(
     445             :         method, absolute_tolerance, relative_tolerance,
     446             :         maximum_absolute_tolerance, condition);
     447             :   }
     448             :   return gsl_multiroot_detail::gsl_multiroot_impl(
     449             :       gsl_func, initial_guess, absolute_tolerance, maximum_iterations,
     450             :       relative_tolerance, verbosity, maximum_absolute_tolerance, method,
     451             :       solver_type, condition, &gsl_multiroot_fdfsolver_alloc,
     452             :       &gsl_multiroot_fdfsolver_set, &gsl_multiroot_fdfsolver_iterate,
     453             :       &gsl_multiroot_fdfsolver_free);
     454             : }
     455             : 
     456             : template <size_t Dim, typename Function,
     457             :           Requires<not gsl_multiroot_detail::is_jacobian_callable_v<
     458             :               Function, std::array<double, Dim>>> = nullptr>
     459             : std::array<double, Dim> gsl_multiroot(
     460             :     const Function& func, const std::array<double, Dim>& initial_guess,
     461             :     const double absolute_tolerance, const size_t maximum_iterations,
     462             :     const double relative_tolerance = 0.0,
     463             :     const Verbosity verbosity = Verbosity::Silent,
     464             :     const double maximum_absolute_tolerance = 0.0,
     465             :     const Method method = Method::Newton,
     466             :     const StoppingCondition condition = StoppingCondition::Absolute) {
     467             :   gsl_multiroot_function gsl_func = {
     468             :       &gsl_multiroot_detail::gsl_multirootfunctionfdf_wrapper_f<Dim, Function>,
     469             :       Dim, const_cast<Function*>(&func)};  // NOLINT
     470             : 
     471             :   // Set up method for solver:
     472             :   const gsl_multiroot_fsolver_type* solver_type;
     473             :   if (method == Method::Newton) {
     474             :     solver_type = gsl_multiroot_fsolver_dnewton;
     475             :   } else if (method == Method::Hybrids) {
     476             :     solver_type = gsl_multiroot_fsolver_hybrids;
     477             :   } else if (method == Method::Hybrid) {
     478             :     solver_type = gsl_multiroot_fsolver_hybrid;
     479             :   } else {
     480             :     ERROR(
     481             :         "Invalid method. Has to be one of Newton, Hybrids or "
     482             :         "Hybrid.");
     483             :   }
     484             :   // Print initial parameters
     485             :   if (UNLIKELY(verbosity == Verbosity::Verbose or
     486             :                verbosity == Verbosity::Debug)) {
     487             :     gsl_multiroot_detail::print_rootfinding_parameters(
     488             :         method, absolute_tolerance, relative_tolerance,
     489             :         maximum_absolute_tolerance, condition);
     490             :   }
     491             :   return gsl_multiroot_detail::gsl_multiroot_impl(
     492             :       gsl_func, initial_guess, absolute_tolerance, maximum_iterations,
     493             :       relative_tolerance, verbosity, maximum_absolute_tolerance, method,
     494             :       solver_type, condition, &gsl_multiroot_fsolver_alloc,
     495             :       &gsl_multiroot_fsolver_set, &gsl_multiroot_fsolver_iterate,
     496             :       &gsl_multiroot_fsolver_free);
     497             : }
     498             : // @}
     499             : }  // namespace RootFinder

Generated by: LCOV version 1.14