SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/RootFinding - QuadraticEquation.cpp 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             : #include "NumericalAlgorithms/RootFinding/QuadraticEquation.hpp"
       5             : 
       6             : #include <gsl/gsl_poly.h>
       7             : #include <limits>
       8             : 
       9             : #include "DataStructures/DataVector.hpp"
      10             : #include "Utilities/EqualWithinRoundoff.hpp"
      11             : #include "Utilities/ErrorHandling/Assert.hpp"
      12             : #include "Utilities/ErrorHandling/Error.hpp"
      13             : #include "Utilities/GenerateInstantiations.hpp"
      14             : 
      15           1 : double positive_root(const double a, const double b, const double c) noexcept {
      16             :   const auto roots = real_roots(a, b, c);
      17             :   ASSERT(roots[0] <= 0.0 and roots[1] >= 0.0,
      18             :          "There are two positive roots, " << roots[0] << " and " << roots[1]
      19             :          << ", with a=" << a << " b=" << b << " c=" << c);
      20             :   return roots[1];
      21             : }
      22             : 
      23           1 : std::array<double, 2> real_roots(const double a, const double b,
      24             :                                  const double c) noexcept {
      25             :   double x0 = std::numeric_limits<double>::signaling_NaN();
      26             :   double x1 = std::numeric_limits<double>::signaling_NaN();
      27             :   // clang-tidy: value stored ... never read (true if in Release Build)
      28             :   // NOLINTNEXTLINE
      29             :   const int num_real_roots = gsl_poly_solve_quadratic(a, b, c, &x0, &x1);
      30             :   ASSERT(num_real_roots == 2,
      31             :          "There are only " << num_real_roots << " real roots with a=" << a
      32             :          << " b=" << b << " c=" << c);
      33             :   return {{x0, x1}};
      34             : }
      35             : 
      36             : namespace detail {
      37             : template <typename T>
      38             : struct root_between_values_impl;
      39             : enum class RootToChoose { min, max };
      40             : }  // namespace detail
      41             : 
      42             : template <typename T>
      43           1 : T smallest_root_greater_than_value_within_roundoff(
      44             :     const T& a, const T& b, const T& c, const double value) noexcept {
      45             :   return detail::root_between_values_impl<T>::function(
      46             :       a, b, c, value, std::numeric_limits<double>::max(),
      47             :       detail::RootToChoose::min);
      48             : }
      49             : 
      50             : template <typename T>
      51           1 : T largest_root_between_values_within_roundoff(const T& a, const T& b,
      52             :                                               const T& c,
      53             :                                               const double min_value,
      54             :                                               const double max_value) noexcept {
      55             :   return detail::root_between_values_impl<T>::function(
      56             :       a, b, c, min_value, max_value, detail::RootToChoose::max);
      57             : }
      58             : 
      59             : namespace detail {
      60             : template <>
      61             : struct root_between_values_impl<double> {
      62             :   static double function(const double a, const double b, const double c,
      63             :                          const double min_value, const double max_value,
      64             :                          const RootToChoose min_or_max) noexcept {
      65             :     // Roots are returned in increasing order.
      66             :     const auto roots = real_roots(a, b, c);
      67             : 
      68             :     const std::array<bool, 2> roots_out_of_bounds_low{
      69             :         {roots[0] < min_value and
      70             :              not equal_within_roundoff(roots[0], min_value),
      71             :          roots[1] < min_value and
      72             :              not equal_within_roundoff(roots[1], min_value)}};
      73             :     const std::array<bool, 2> roots_out_of_bounds_high{
      74             :         {roots[0] > max_value and
      75             :              not equal_within_roundoff(roots[0], max_value),
      76             :          roots[1] > max_value and
      77             :              not equal_within_roundoff(roots[1], max_value)}};
      78             : 
      79             :     double return_value = std::numeric_limits<double>::signaling_NaN();
      80             :     const auto error_message = [&a, &b, &c,
      81             :                                 &roots](const std::string& message) noexcept {
      82             :       ERROR(message << " Roots are " << roots[0] << " and " << roots[1]
      83             :                     << ", with a=" << a << " b=" << b << " c=" << c);
      84             :     };
      85             : 
      86             :     if (min_or_max == RootToChoose::min) {
      87             :       // Check roots[0] first because it is the smallest
      88             :       if (roots_out_of_bounds_low[0]) {
      89             :         if (roots_out_of_bounds_low[1]) {
      90             :           error_message("No root >= (within roundoff) min_value.");
      91             :         }
      92             :         if (roots_out_of_bounds_high[1]) {
      93             :           error_message(
      94             :               "No root between min_value and max_value (within roundoff).");
      95             :         }
      96             :         return_value = roots[1];
      97             :       } else {
      98             :         if (roots_out_of_bounds_high[0]) {
      99             :           error_message("No root <= (within roundoff) max_value.");
     100             :         }
     101             :         return_value = roots[0];
     102             :       }
     103             :     } else {
     104             :       // Check roots[1] first because it is the largest
     105             :       if (roots_out_of_bounds_high[1]) {
     106             :         if (roots_out_of_bounds_high[0]) {
     107             :           error_message("No root <= (within roundoff) max_value.");
     108             :         }
     109             :         if (roots_out_of_bounds_low[0]) {
     110             :           error_message(
     111             :               "No root between min_value and max_value (within "
     112             :               "roundoff).");
     113             :         }
     114             :         return_value = roots[0];
     115             :       } else {
     116             :         if (roots_out_of_bounds_low[1]) {
     117             :           error_message("No root >= (within roundoff) min_value.");
     118             :         }
     119             :         return_value = roots[1];
     120             :       }
     121             :     }
     122             :     return return_value;
     123             :   }
     124             : };
     125             : 
     126             : template <>
     127             : struct root_between_values_impl<DataVector> {
     128             :   static DataVector function(const DataVector& a, const DataVector& b,
     129             :                              const DataVector& c, const double min_value,
     130             :                              const double max_value,
     131             :                              const RootToChoose min_or_max) noexcept {
     132             :     ASSERT(a.size() == b.size(),
     133             :            "Size mismatch a vs b: " << a.size() << " " << b.size());
     134             :     ASSERT(a.size() == c.size(),
     135             :            "Size mismatch a vs c: " << a.size() << " " << c.size());
     136             :     DataVector result(a.size());
     137             :     for (size_t i = 0; i < a.size(); ++i) {
     138             :       result[i] = root_between_values_impl<double>::function(
     139             :           a[i], b[i], c[i], min_value, max_value, min_or_max);
     140             :     }
     141             :     return result;
     142             :   }
     143             : };
     144             : }  // namespace detail
     145             : 
     146             : // Explicit instantiations
     147             : /// \cond
     148             : #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
     149             : 
     150             : #define INSTANTIATE(_, data)                                               \
     151             :   template DTYPE(data) smallest_root_greater_than_value_within_roundoff(   \
     152             :       const DTYPE(data) & a, const DTYPE(data) & b, const DTYPE(data) & c, \
     153             :       double value) noexcept;                                              \
     154             :   template DTYPE(data) largest_root_between_values_within_roundoff(        \
     155             :       const DTYPE(data) & a, const DTYPE(data) & b, const DTYPE(data) & c, \
     156             :       double min_value, double max_value) noexcept;
     157             : 
     158             : GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector))
     159             : 
     160             : #undef DTYPE
     161             : #undef INSTANTIATE
     162             : /// \endcond

Generated by: LCOV version 1.14