From 9b43f48f85ae84d819c71cd235c3692f7ee0d69d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 16 Jan 2024 11:49:59 -0700 Subject: [PATCH] Use a better initialization scheme. --- include/deal.II/base/tensor.h | 123 ++++++++++++++++++++++++++++++++-- 1 file changed, 116 insertions(+), 7 deletions(-) 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::Tensor(const ArrayLike &initializer, } +# ifdef DEAL_II_HAVE_CXX20 template constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE Tensor::Tensor() + : values( + // In order to initialize the std::array, 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. + []( + const std::index_sequence &) constexpr -> decltype(values) { + if constexpr (rank_ == 1) + { + auto get_zero_and_ignore_argument = [](int) { + return internal::NumberType::value(0.0); + }; + return {{(get_zero_and_ignore_argument(I))...}}; + } + else + { + auto get_zero_and_ignore_argument = [](int) { + return Tensor(); + }; + return {{(get_zero_and_ignore_argument(I))...}}; + } + }(std::make_index_sequence())) +{} + +# 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 +// [](const std::index_sequence &) +// 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::value(0.0)); -} + namespace TensorInitialization + { + template + constexpr std::array::value_type, dim> + make_zero_array(const std::index_sequence &) + { + 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::value(0.0); + }; + return {{(get_zero_and_ignore_argument(I))...}}; + } + else + { + auto get_zero_and_ignore_argument = [](int) { + return Tensor(); + }; + return {{(get_zero_and_ignore_argument(I))...}}; + } + } + } // namespace TensorInitialization +} // namespace internal + +template +constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE +Tensor::Tensor() + : values(internal::TensorInitialization::make_zero_array( + std::make_index_sequence())) +{} + + +# endif template @@ -1557,8 +1658,16 @@ constexpr inline bool Tensor::operator==( const Tensor &p) const { +# ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING + Assert(!(std::is_same_v || + std::is_same_v), + 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; } -- 2.39.5