]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Don't include Sacado.hpp in base/numbers.h.
authorDavid Wells <drwells@email.unc.edu>
Wed, 3 Aug 2022 15:39:02 +0000 (11:39 -0400)
committerDavid Wells <drwells@email.unc.edu>
Thu, 4 Aug 2022 22:03:45 +0000 (18:03 -0400)
Sacado.hpp, when preprocessed, is over 200,000 lines. We can get things working
in an equivalent way by using argument dependent lookup instead.

include/deal.II/base/numbers.h
include/deal.II/base/symmetric_tensor.h
include/deal.II/base/tensor.h
include/deal.II/differentiation/ad/sacado_math.h

index 74d8976ee2692b3a83017ac507029f606fc86eac..292816d2dd042065554fb146e2c53d592f6a8fac 100644 (file)
@@ -137,15 +137,6 @@ DEAL_II_NAMESPACE_CLOSE
 
 #  include <adolc/adouble.h> // Taped double
 #endif
-// Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
-// but header indirectly references numbers.h. We therefore simply
-// import the whole Sacado header at this point to get the math
-// functions imported into the standard namespace.
-#ifdef DEAL_II_TRILINOS_WITH_SACADO
-DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
-#  include <Sacado.hpp>
-DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
-#endif
 
 namespace std
 {
@@ -608,7 +599,9 @@ namespace numbers
   typename NumberTraits<number>::real_type
   NumberTraits<number>::abs(const number &x)
   {
-    return std::abs(x);
+    // Make things work with AD types
+    using std::abs;
+    return abs(x);
   }
 
 
@@ -626,7 +619,9 @@ namespace numbers
   typename NumberTraits<std::complex<number>>::real_type
   NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
   {
-    return std::abs(x);
+    // Make things work with AD types
+    using std::abs;
+    return abs(x);
   }
 
 
index d69798c9a6c12a5411ed096647012f4f6cb7c338..6396f05d58fc95c132e3c6cef8634202f564d79b 100644 (file)
@@ -1563,7 +1563,7 @@ namespace internal
         // whether we have to fear that the matrix is not regular.
         Number diagonal_sum = internal::NumberType<Number>::value(0.0);
         for (unsigned int i = 0; i < N; ++i)
-          diagonal_sum += std::fabs(tmp.data[i][i]);
+          diagonal_sum += numbers::NumberTraits<Number>::abs(tmp.data[i][i]);
         const Number typical_diagonal_element =
           diagonal_sum / static_cast<double>(N);
         (void)typical_diagonal_element;
@@ -1576,12 +1576,12 @@ namespace internal
           {
             // Pivot search: search that part of the line on and right of the
             // diagonal for the largest element.
-            Number       max = std::fabs(tmp.data[j][j]);
-            unsigned int r   = j;
+            Number max     = numbers::NumberTraits<Number>::abs(tmp.data[j][j]);
+            unsigned int r = j;
             for (unsigned int i = j + 1; i < N; ++i)
-              if (std::fabs(tmp.data[i][j]) > max)
+              if (numbers::NumberTraits<Number>::abs(tmp.data[i][j]) > max)
                 {
-                  max = std::fabs(tmp.data[i][j]);
+                  max = numbers::NumberTraits<Number>::abs(tmp.data[i][j]);
                   r   = i;
                 }
 
@@ -2343,25 +2343,25 @@ namespace internal
   compute_norm(const typename SymmetricTensorAccessors::
                  StorageType<2, dim, Number>::base_tensor_type &data)
   {
+    // Make things work with AD types
+    using std::sqrt;
     switch (dim)
       {
         case 1:
           return numbers::NumberTraits<Number>::abs(data[0]);
 
         case 2:
-          return std::sqrt(
-            numbers::NumberTraits<Number>::abs_square(data[0]) +
-            numbers::NumberTraits<Number>::abs_square(data[1]) +
-            2. * numbers::NumberTraits<Number>::abs_square(data[2]));
+          return sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
+                      numbers::NumberTraits<Number>::abs_square(data[1]) +
+                      2. * numbers::NumberTraits<Number>::abs_square(data[2]));
 
         case 3:
-          return std::sqrt(
-            numbers::NumberTraits<Number>::abs_square(data[0]) +
-            numbers::NumberTraits<Number>::abs_square(data[1]) +
-            numbers::NumberTraits<Number>::abs_square(data[2]) +
-            2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
-            2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
-            2. * numbers::NumberTraits<Number>::abs_square(data[5]));
+          return sqrt(numbers::NumberTraits<Number>::abs_square(data[0]) +
+                      numbers::NumberTraits<Number>::abs_square(data[1]) +
+                      numbers::NumberTraits<Number>::abs_square(data[2]) +
+                      2. * numbers::NumberTraits<Number>::abs_square(data[3]) +
+                      2. * numbers::NumberTraits<Number>::abs_square(data[4]) +
+                      2. * numbers::NumberTraits<Number>::abs_square(data[5]));
 
         default:
           {
@@ -2375,7 +2375,7 @@ namespace internal
               return_value +=
                 2. * numbers::NumberTraits<Number>::abs_square(data[d]);
 
-            return std::sqrt(return_value);
+            return sqrt(return_value);
           }
       }
   }
@@ -2387,6 +2387,8 @@ namespace internal
   compute_norm(const typename SymmetricTensorAccessors::
                  StorageType<4, dim, Number>::base_tensor_type &data)
   {
+    // Make things work with AD types
+    using std::sqrt;
     switch (dim)
       {
         case 1:
@@ -2416,7 +2418,7 @@ namespace internal
                 return_value +=
                   4. * numbers::NumberTraits<Number>::abs_square(data[i][j]);
 
-            return std::sqrt(return_value);
+            return sqrt(return_value);
           }
       }
   }
index 4cde9bde1834bae1e47fefba75ccfe0af5193dd7..b382b2b401b3acbc69c4c388cfbb36f83edd5f19 100644 (file)
@@ -1727,7 +1727,9 @@ template <int rank_, int dim, typename Number>
 inline typename numbers::NumberTraits<Number>::real_type
 Tensor<rank_, dim, Number>::norm() const
 {
-  return std::sqrt(norm_square());
+  // Make things work with AD types
+  using std::sqrt;
+  return sqrt(norm_square());
 }
 
 
@@ -3009,7 +3011,7 @@ l1_norm(const Tensor<2, dim, Number> &t)
     {
       Number sum = internal::NumberType<Number>::value(0.0);
       for (unsigned int i = 0; i < dim; ++i)
-        sum += std::fabs(t[i][j]);
+        sum += numbers::NumberTraits<Number>::abs(t[i][j]);
 
       if (sum > max)
         max = sum;
@@ -3035,7 +3037,7 @@ linfty_norm(const Tensor<2, dim, Number> &t)
     {
       Number sum = internal::NumberType<Number>::value(0.0);
       for (unsigned int j = 0; j < dim; ++j)
-        sum += std::fabs(t[i][j]);
+        sum += numbers::NumberTraits<Number>::abs(t[i][j]);
 
       if (sum > max)
         max = sum;
index 655d281ce73921ad329061275a841b0c4f185860..bd9d8cde9378561a1da7aa81e8d06fbdd5f9c330 100644 (file)
 
 #  ifndef DOXYGEN
 
+DEAL_II_NAMESPACE_OPEN
 /**
- * Define some missing fundamental math functions
+ * Implementation of the error function for real-valued Sacado numbers.
  */
-namespace std
+template <
+  typename ADNumberType,
+  typename = std::enable_if_t<
+    dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value &&
+    dealii::Differentiation::AD::is_real_valued_ad_number<ADNumberType>::value>>
+inline ADNumberType
+erf(ADNumberType x)
 {
-  /**
-   * Implementation of the error function for real-valued Sacado numbers.
-   */
-  template <
-    typename ADNumberType,
-    typename = std::enable_if_t<
-      dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value &&
-      dealii::Differentiation::AD::is_real_valued_ad_number<
-        ADNumberType>::value>>
-  inline ADNumberType
-  erf(ADNumberType x)
-  {
-    // Reference:
-    // Handbook of Mathematical Functions: with Formulas, Graphs, and
-    // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on
-    // Mathematics. 1972.
-    //
-    // Current implementation source:
-    // https://www.johndcook.com/blog/cpp_erf/
-    // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
-    //
-    // Note: This implementation has a reported maximum round-off error
-    // of 1.5e-7.
-
-    // Constants
-    const double a1 = 0.254829592;
-    const double a2 = -0.284496736;
-    const double a3 = 1.421413741;
-    const double a4 = -1.453152027;
-    const double a5 = 1.061405429;
-    const double p  = 0.3275911;
-
-    // Save the sign of x
-    const bool neg_val =
-      (x < dealii::internal::NumberType<ADNumberType>::value(0.0) ? true :
-                                                                    false);
-    x = std::fabs(x);
-
-    // Abramowitz1972a equation 7.1.26
-    const ADNumberType t = 1.0 / (1.0 + p * x);
-    const ADNumberType y =
-      1.0 -
-      (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
-
-    if (!neg_val)
-      return y;
-    else
-      return -y;
-  }
-
-
-  /**
-   * Implementation of the complementary error function for Sacado numbers.
-   */
-  template <
-    typename ADNumberType,
-    typename = std::enable_if_t<
-      dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value>>
-  inline ADNumberType
-  erfc(const ADNumberType &x)
-  {
-    return 1.0 - std::erf(x);
-  }
-
-} // namespace std
+  // Reference:
+  // Handbook of Mathematical Functions: with Formulas, Graphs, and
+  // Mathematical Tables Abramowitz, M. and Stegun, I. Dover Books on
+  // Mathematics. 1972.
+  //
+  // Current implementation source:
+  // https://www.johndcook.com/blog/cpp_erf/
+  // https://www.johndcook.com/blog/2009/01/19/stand-alone-error-function-erf/
+  //
+  // Note: This implementation has a reported maximum round-off error
+  // of 1.5e-7.
+
+  // Constants
+  const double a1 = 0.254829592;
+  const double a2 = -0.284496736;
+  const double a3 = 1.421413741;
+  const double a4 = -1.453152027;
+  const double a5 = 1.061405429;
+  const double p  = 0.3275911;
+
+  // Save the sign of x
+  const bool neg_val =
+    (x < dealii::internal::NumberType<ADNumberType>::value(0.0) ? true : false);
+  x = std::fabs(x);
+
+  // Abramowitz1972a equation 7.1.26
+  const ADNumberType t = 1.0 / (1.0 + p * x);
+  const ADNumberType y =
+    1.0 -
+    (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * std::exp(-x * x);
+
+  if (!neg_val)
+    return y;
+  else
+    return -y;
+}
+
+
+/**
+ * Implementation of the complementary error function for Sacado numbers.
+ */
+template <typename ADNumberType,
+          typename = std::enable_if_t<
+            dealii::Differentiation::AD::is_sacado_number<ADNumberType>::value>>
+inline ADNumberType
+erfc(const ADNumberType &x)
+{
+  return 1.0 - std::erf(x);
+}
+
+DEAL_II_NAMESPACE_CLOSE
 
 #  endif // DOXYGEN
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.