yat/utility/RootFinder.cc

Code
Comments
Other
Rev Date Author Line
4237 17 Sep 22 peter 1 // $Id$
4237 17 Sep 22 peter 2
4237 17 Sep 22 peter 3 /*
4237 17 Sep 22 peter 4   Copyright (C) 2022 Peter Johansson
4237 17 Sep 22 peter 5
4237 17 Sep 22 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4237 17 Sep 22 peter 7
4237 17 Sep 22 peter 8   The yat library is free software; you can redistribute it and/or
4237 17 Sep 22 peter 9   modify it under the terms of the GNU General Public License as
4237 17 Sep 22 peter 10   published by the Free Software Foundation; either version 3 of the
4237 17 Sep 22 peter 11   License, or (at your option) any later version.
4237 17 Sep 22 peter 12
4237 17 Sep 22 peter 13   The yat library is distributed in the hope that it will be useful,
4237 17 Sep 22 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4237 17 Sep 22 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4237 17 Sep 22 peter 16   General Public License for more details.
4237 17 Sep 22 peter 17
4237 17 Sep 22 peter 18   You should have received a copy of the GNU General Public License
4237 17 Sep 22 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4237 17 Sep 22 peter 20 */
4237 17 Sep 22 peter 21
4237 17 Sep 22 peter 22 #include <config.h>
4237 17 Sep 22 peter 23
4237 17 Sep 22 peter 24 #include "RootFinder.h"
4237 17 Sep 22 peter 25
4237 17 Sep 22 peter 26 #include <gsl/gsl_errno.h>
4237 17 Sep 22 peter 27
4237 17 Sep 22 peter 28 #include <cassert>
4237 17 Sep 22 peter 29 #include <memory>
4237 17 Sep 22 peter 30
4237 17 Sep 22 peter 31 namespace theplu {
4237 17 Sep 22 peter 32 namespace yat {
4237 17 Sep 22 peter 33 namespace utility {
4237 17 Sep 22 peter 34
4237 17 Sep 22 peter 35   RootFinder::RootFinder(const gsl_root_fsolver_type* t)
4237 17 Sep 22 peter 36     : type_(t), solver_(gsl_root_fsolver_alloc(t))
4237 17 Sep 22 peter 37   {}
4237 17 Sep 22 peter 38
4237 17 Sep 22 peter 39
4237 17 Sep 22 peter 40   const char* RootFinder::name(void) const
4237 17 Sep 22 peter 41   {
4237 17 Sep 22 peter 42     assert(solver_.get());
4237 17 Sep 22 peter 43     return gsl_root_fsolver_name(solver_.get());
4237 17 Sep 22 peter 44   }
4237 17 Sep 22 peter 45
4237 17 Sep 22 peter 46
4237 17 Sep 22 peter 47   RootFinder::Interval::Interval(double epsabs, double epsrel)
4237 17 Sep 22 peter 48     : epsabs_(epsabs), epsrel_(epsrel)
4237 17 Sep 22 peter 49   {}
4237 17 Sep 22 peter 50
4237 17 Sep 22 peter 51
4237 17 Sep 22 peter 52   bool RootFinder::Interval::operator()(const gsl_root_fsolver* s)
4237 17 Sep 22 peter 53   {
4237 17 Sep 22 peter 54     double lower = gsl_root_fsolver_x_lower(s);
4237 17 Sep 22 peter 55     double upper = gsl_root_fsolver_x_upper(s);
4237 17 Sep 22 peter 56     return gsl_root_test_interval(lower, upper, epsabs_, epsrel_) ==
4237 17 Sep 22 peter 57       GSL_SUCCESS;
4237 17 Sep 22 peter 58   }
4237 17 Sep 22 peter 59
4237 17 Sep 22 peter 60
4237 17 Sep 22 peter 61   RootFinder::Delta::Delta(double epsabs, double epsrel)
4237 17 Sep 22 peter 62     : epsabs_(epsabs), epsrel_(epsrel)
4237 17 Sep 22 peter 63   {}
4237 17 Sep 22 peter 64
4237 17 Sep 22 peter 65
4237 17 Sep 22 peter 66   bool RootFinder::Delta::operator()(const gsl_root_fsolver* s)
4237 17 Sep 22 peter 67   {
4237 17 Sep 22 peter 68     double x = gsl_root_fsolver_root(s);
4237 17 Sep 22 peter 69     if (!prev_) {
4237 17 Sep 22 peter 70       prev_.reset(new double(x));
4237 17 Sep 22 peter 71       return false;
4237 17 Sep 22 peter 72     }
4237 17 Sep 22 peter 73
4237 17 Sep 22 peter 74     int status = gsl_root_test_delta(*prev_, x, epsabs_, epsrel_);
4237 17 Sep 22 peter 75     *prev_ = x;
4237 17 Sep 22 peter 76     return status == GSL_SUCCESS;
4237 17 Sep 22 peter 77   }
4237 17 Sep 22 peter 78
4237 17 Sep 22 peter 79
4237 17 Sep 22 peter 80   void RootFinder::GslFree::operator()(gsl_root_fsolver* solver) const
4237 17 Sep 22 peter 81   {
4237 17 Sep 22 peter 82     gsl_root_fsolver_free(solver);
4237 17 Sep 22 peter 83   }
4237 17 Sep 22 peter 84 }}}

Generated on Sat Oct 14 02:25:13 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)