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

Generated by: LCOV version 1.14