From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 16 Jan 2024 18:49:59 +0000 (-0700)
Subject: Use a better initialization scheme.
X-Git-Tag: relicensing~145^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16480%2Fhead;p=dealii.git

Use a better initialization scheme.
---

diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index b951af8709..a452ce5411 100644
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -1320,18 +1320,119 @@ Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
 }
 
 
+#  ifdef DEAL_II_HAVE_CXX20
 
 template <int rank_, int dim, typename Number>
 constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE
 Tensor<rank_, dim, Number>::Tensor()
+  : values(
+      // In order to initialize the std::array<Number,dim>, we would need a
+      // brace-enclosed list of length 'dim'. There is no way in C++ to create
+      // such a list in-place, but we can come up with a lambda function that
+      // expands such a list via template-pack expansion, and then uses this
+      // list to initialize a std::array which it then returns.
+      //
+      // The trick to come up with such a lambda function is to have a function
+      // that takes an argument that depends on a template-pack of integers.
+      // We will call the function with an integer list of length 'dim', and
+      // in the function itself expand that pack in a way that it serves as
+      // a brace-enclosed list of initializers for a std::array.
+      //
+      // Of course, we do not want to initialize the array with the integers,
+      // but with  zeros. (Or, more correctly, a zero of the element type.)
+      // The canonical way to do this would be using the comma operator:
+      //     (sequence_element, 0.0)
+      // returns zero, and
+      //     (sequence, 0.0)...
+      // returns a list of zeros of the right length. Unfortunately, some
+      // compilers then warn that the left side of the comma expression has
+      // no effect -- well, bummer, that was of course exactly the idea.
+      // We could work around this by using
+      //     (sequence_element * 0.0)
+      // instead, assuming that the compiler will optimize (known) integer
+      // times zero to zero, and similarly for (known) integer times times
+      // default-initialized tensor.
+      //
+      // But, instead of relying on compiler optimizations, a better way is
+      // to simply have another (nested) lambda function that takes the
+      // integer sequence element as an argument and ignores it, just
+      // returning a zero instead.
+      []<std::size_t... I>(
+        const std::index_sequence<I...> &) constexpr -> decltype(values) {
+        if constexpr (rank_ == 1)
+          {
+            auto get_zero_and_ignore_argument = [](int) {
+              return internal::NumberType<Number>::value(0.0);
+            };
+            return {{(get_zero_and_ignore_argument(I))...}};
+          }
+        else
+          {
+            auto get_zero_and_ignore_argument = [](int) {
+              return Tensor<rank_ - 1, dim, Number>();
+            };
+            return {{(get_zero_and_ignore_argument(I))...}};
+          }
+      }(std::make_index_sequence<dim>()))
+{}
+
+#  else
+
+// The C++17 case works in essence the same, except that we can't use a
+// lambda function with explicit template parameters, i.e., we can't do
+//   []<std::size_t... I>(const std::index_sequence<I...> &)
+// as above because that's a C++20 feature. Lambda functions in C++17 can
+// have template packs as arguments, but we need the ability to *name*
+// that template pack (the 'I' above) and that's not possible in C++17.
+//
+// We work around this by moving the lambda function to a global function
+// and using the traditional template syntax on it.
+namespace internal
 {
-  // The default constructor of std::array does not initialize its elements.
-  // So we have to do it by hand.
-  std::fill(values.begin(),
-            values.end(),
-            internal::NumberType<Number>::value(0.0));
-}
+  namespace TensorInitialization
+  {
+    template <int rank, int dim, typename Number, std::size_t... I>
+    constexpr std::array<typename Tensor<rank, dim, Number>::value_type, dim>
+    make_zero_array(const std::index_sequence<I...> &)
+    {
+      static_assert(sizeof...(I) == dim, "This is bad.");
+
+      // First peel off the case dim==0. If we don't, some compilers
+      // will warn below that we define these lambda functions but
+      // never use them (because the expanded list has zero elements,
+      // and the get_zero_and_ignore_argument() function is not used...)
+      if constexpr (dim == 0)
+        {
+          return {};
+        }
+      else if constexpr (rank == 1)
+        {
+          auto get_zero_and_ignore_argument = [](int) {
+            return internal::NumberType<Number>::value(0.0);
+          };
+          return {{(get_zero_and_ignore_argument(I))...}};
+        }
+      else
+        {
+          auto get_zero_and_ignore_argument = [](int) {
+            return Tensor<rank - 1, dim, Number>();
+          };
+          return {{(get_zero_and_ignore_argument(I))...}};
+        }
+    }
+  } // namespace TensorInitialization
+} // namespace internal
+
 
+template <int rank_, int dim, typename Number>
+constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE
+Tensor<rank_, dim, Number>::Tensor()
+  : values(internal::TensorInitialization::make_zero_array<rank_, dim, Number>(
+      std::make_index_sequence<dim>()))
+{}
+
+
+#  endif
 
 
 template <int rank_, int dim, typename Number>
@@ -1557,8 +1658,16 @@ constexpr inline bool
 Tensor<rank_, dim, Number>::operator==(
   const Tensor<rank_, dim, OtherNumber> &p) const
 {
+#  ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
+  Assert(!(std::is_same_v<Number, adouble> ||
+           std::is_same_v<OtherNumber, adouble>),
+         ExcMessage(
+           "The Tensor equality operator for ADOL-C taped numbers has not yet "
+           "been extended to support advanced branching."));
+#  endif
+
   for (unsigned int i = 0; i < dim; ++i)
-    if (values[i] != p.values[i])
+    if (numbers::values_are_not_equal(values[i], p.values[i]))
       return false;
   return true;
 }