From 220513e8cf84acb1e682bdacbf02d2f6a8f74bf3 Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Tue, 7 May 2019 19:05:33 -0700 Subject: [PATCH] Tensor is now constexpr CONSEXPR macro added to Tensor functions. New constructors added to VectorizedArray Tensor free functions all made constexpr constexpr added to tensor_accessors.h default initialization for VectorizedArray constexpr feature added to Tensor class adding DEAL_II_ALWAYS_INLINE back to functions in tensor_accessors.h Add assignement operator ans default constructor back Apply suggestions from code review Typos fixed Co-Authored-By: Daniel Arndt Indentation Improve check_01_cxx_features.cmake Fix changelog entry static constexpr converted to back static const for VectorizedArray one more constexpr tensor test test failure fixed added DEAL_II_CONSTEXPR to three more functions removed member initialization in VectorizedArray Avoid warning constructor added to VectorizedArray reverted the changes in vectorization.h asserting VectorizedArray is POD --- cmake/checks/check_01_cxx_features.cmake | 13 + doc/news/changes/minor/20190526RezaRastak | 4 + include/deal.II/base/config.h.in | 1 + include/deal.II/base/tensor.h | 398 ++++++++++++---------- include/deal.II/base/tensor_accessors.h | 61 ++-- source/base/vectorization.cc | 13 + tests/base/vectorization_07.cc | 36 -- tests/base/vectorization_07.output | 2 - tests/tensors/constexpr_tensor.cc | 95 ++++++ tests/tensors/constexpr_tensor.output | 78 +++++ 10 files changed, 451 insertions(+), 250 deletions(-) create mode 100644 doc/news/changes/minor/20190526RezaRastak delete mode 100644 tests/base/vectorization_07.cc delete mode 100644 tests/base/vectorization_07.output create mode 100644 tests/tensors/constexpr_tensor.cc create mode 100644 tests/tensors/constexpr_tensor.output diff --git a/cmake/checks/check_01_cxx_features.cmake b/cmake/checks/check_01_cxx_features.cmake index edd3f7996c..2285e206fa 100644 --- a/cmake/checks/check_01_cxx_features.cmake +++ b/cmake/checks/check_01_cxx_features.cmake @@ -28,6 +28,7 @@ # DEAL_II_HAVE_FP_EXCEPTIONS # DEAL_II_HAVE_COMPLEX_OPERATOR_OVERLOADS # +# DEAL_II_CONSTEXPR # DEAL_II_FALLTHROUGH # @@ -753,6 +754,18 @@ CHECK_CXX_SOURCE_COMPILES( " DEAL_II_HAVE_CXX14_CONSTEXPR_CAN_CALL_NONCONSTEXPR) +# +# The macro DEAL_II_CONSTEXPR allows using c++ constexpr features in a portable way. +# Here we enable it only when a constexpr function can call simple non-constexpr +# functions. This requirement is probabely very conservative in most cases, but +# it will prevent breaking builds with certain compilers. +# +IF (DEAL_II_HAVE_CXX14_CONSTEXPR_CAN_CALL_NONCONSTEXPR) + SET(DEAL_II_CONSTEXPR "constexpr") +ELSE() + SET(DEAL_II_CONSTEXPR " ") +ENDIF() + # # Not all compilers with C++17 support include the new special math # functions. Check this separately so that we can use C++17 compilers that don't diff --git a/doc/news/changes/minor/20190526RezaRastak b/doc/news/changes/minor/20190526RezaRastak new file mode 100644 index 0000000000..475e038763 --- /dev/null +++ b/doc/news/changes/minor/20190526RezaRastak @@ -0,0 +1,4 @@ +New: Tensor can be constructed and manipulated in constexpr +setting (if supported by the compiler). +
+(Reza Rastak, 2019/05/26) diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index d4e79ef91c..a86966b272 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -99,6 +99,7 @@ #cmakedefine DEAL_II_ALWAYS_INLINE @DEAL_II_ALWAYS_INLINE@ #cmakedefine DEAL_II_RESTRICT @DEAL_II_RESTRICT@ #cmakedefine DEAL_II_COMPILER_HAS_DIAGNOSTIC_PRAGMA +#cmakedefine DEAL_II_CONSTEXPR @DEAL_II_CONSTEXPR@ /* * A variable to tell if the compiler used in the current compilation process diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index d0a350fc8d..d9f8f4c676 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -142,7 +143,7 @@ public: * * @note This function can also be used in CUDA device code. */ - DEAL_II_CUDA_HOST_DEV + constexpr DEAL_II_CUDA_HOST_DEV Tensor(); /** @@ -151,19 +152,19 @@ public: * Number. */ template - Tensor(const Tensor<0, dim, OtherNumber> &initializer); + constexpr Tensor(const Tensor<0, dim, OtherNumber> &initializer); /** * Constructor, where the data is copied from a C-style array. */ template - Tensor(const OtherNumber &initializer); + constexpr Tensor(const OtherNumber &initializer); /** * Return a pointer to the first element of the underlying storage. */ - Number * - begin_raw(); + DEAL_II_CONSTEXPR Number * + begin_raw(); /** * Return a const pointer to the first element of the underlying storage. @@ -174,8 +175,8 @@ public: /** * Return a pointer to the element past the end of the underlying storage. */ - Number * - end_raw(); + DEAL_II_CONSTEXPR Number * + end_raw(); /** * Return a const pointer to the element past the end of the underlying @@ -193,7 +194,7 @@ public: * * @note This function can also be used in CUDA device code. */ - DEAL_II_CUDA_HOST_DEV operator Number &(); + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator Number &(); /** * Return a reference to the encapsulated Number object. Since rank-0 @@ -203,7 +204,7 @@ public: * * @note This function can also be used in CUDA device code. */ - DEAL_II_CUDA_HOST_DEV operator const Number &() const; + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator const Number &() const; /** * Assignment from tensors with different underlying scalar type. This @@ -211,8 +212,8 @@ public: * Number. */ template - Tensor & - operator=(const Tensor<0, dim, OtherNumber> &rhs); + DEAL_II_CONSTEXPR Tensor & + operator=(const Tensor<0, dim, OtherNumber> &rhs); #ifdef __INTEL_COMPILER /** @@ -221,8 +222,8 @@ public: * copy constructor for Sacado::Rad::ADvar types automatically. * See https://github.com/dealii/dealii/pull/5865. */ - Tensor & - operator=(const Tensor<0, dim, Number> &rhs); + DEAL_II_CONSTEXPR Tensor & + operator=(const Tensor<0, dim, Number> &rhs); #endif /** @@ -230,14 +231,14 @@ public: * that the @p OtherNumber type is convertible to @p Number. */ template - Tensor & - operator=(const OtherNumber &d); + DEAL_II_CONSTEXPR Tensor & + operator=(const OtherNumber &d); /** * Test for equality of two tensors. */ template - bool + DEAL_II_CONSTEXPR bool operator==(const Tensor<0, dim, OtherNumber> &rhs) const; /** @@ -251,15 +252,15 @@ public: * Add another scalar */ template - Tensor & - operator+=(const Tensor<0, dim, OtherNumber> &rhs); + DEAL_II_CONSTEXPR Tensor & + operator+=(const Tensor<0, dim, OtherNumber> &rhs); /** * Subtract another scalar. */ template - Tensor & - operator-=(const Tensor<0, dim, OtherNumber> &rhs); + DEAL_II_CONSTEXPR Tensor & + operator-=(const Tensor<0, dim, OtherNumber> &rhs); /** * Multiply the scalar with a factor. @@ -267,15 +268,15 @@ public: * @note This function can also be used in CUDA device code. */ template - DEAL_II_CUDA_HOST_DEV Tensor & - operator*=(const OtherNumber &factor); + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor & + operator*=(const OtherNumber &factor); /** * Divide the scalar by factor. */ template - Tensor & - operator/=(const OtherNumber &factor); + DEAL_II_CONSTEXPR Tensor & + operator/=(const OtherNumber &factor); /** * Tensor with inverted entries. @@ -295,7 +296,7 @@ public: * and indeed the state where all elements have a zero value is the state * right after construction of such an object. */ - void + DEAL_II_CONSTEXPR void clear(); /** @@ -303,8 +304,8 @@ public: * the absolute squares of all entries. For the present case of rank-1 * tensors, this equals the usual l2 norm of the vector. */ - real_type - norm() const; + DEAL_II_CONSTEXPR real_type + norm() const; /** * Return the square of the Frobenius-norm of a tensor, i.e. the sum of the @@ -339,7 +340,7 @@ private: * Internal helper function for unroll. */ template - void + DEAL_II_CONSTEXPR void unroll_recursion(Vector &result, unsigned int & start_index) const; @@ -444,7 +445,7 @@ public: /** * Constructor, where the data is copied from a C-style array. */ - explicit Tensor(const array_type &initializer); + constexpr explicit Tensor(const array_type &initializer); /** * Constructor from tensors with different underlying scalar type. This @@ -452,13 +453,13 @@ public: * Number. */ template - Tensor(const Tensor &initializer); + constexpr Tensor(const Tensor &initializer); /** * Constructor that converts from a "tensor of tensors". */ template - Tensor( + constexpr Tensor( const Tensor<1, dim, Tensor> &initializer); /** @@ -473,7 +474,8 @@ public: * * @note This function can also be used in CUDA device code. */ - DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i); + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV value_type & + operator[](const unsigned int i); /** * Read-only access operator. @@ -486,18 +488,19 @@ public: /** * Read access using TableIndices indices */ - const Number &operator[](const TableIndices &indices) const; + DEAL_II_CONSTEXPR const Number & + operator[](const TableIndices &indices) const; /** * Read and write access using TableIndices indices */ - Number &operator[](const TableIndices &indices); + DEAL_II_CONSTEXPR Number &operator[](const TableIndices &indices); /** * Return a pointer to the first element of the underlying storage. */ - Number * - begin_raw(); + DEAL_II_CONSTEXPR Number * + begin_raw(); /** * Return a const pointer to the first element of the underlying storage. @@ -508,8 +511,8 @@ public: /** * Return a pointer to the element past the end of the underlying storage. */ - Number * - end_raw(); + DEAL_II_CONSTEXPR Number * + end_raw(); /** * Return a pointer to the element past the end of the underlying storage. @@ -523,8 +526,8 @@ public: * Number. */ template - Tensor & - operator=(const Tensor &rhs); + DEAL_II_CONSTEXPR Tensor & + operator=(const Tensor &rhs); /** * This operator assigns a scalar to a tensor. To avoid confusion with what @@ -532,14 +535,14 @@ public: * value allowed for d, allowing the intuitive notation * t=0 to reset all elements of the tensor to zero. */ - Tensor & - operator=(const Number &d); + DEAL_II_CONSTEXPR Tensor & + operator=(const Number &d); /** * Test for equality of two tensors. */ template - bool + DEAL_II_CONSTEXPR bool operator==(const Tensor &) const; /** @@ -553,15 +556,15 @@ public: * Add another tensor. */ template - Tensor & - operator+=(const Tensor &); + DEAL_II_CONSTEXPR Tensor & + operator+=(const Tensor &); /** * Subtract another tensor. */ template - Tensor & - operator-=(const Tensor &); + DEAL_II_CONSTEXPR Tensor & + operator-=(const Tensor &); /** * Scale the tensor by factor, i.e. multiply all components by @@ -570,21 +573,21 @@ public: * @note This function can also be used in CUDA device code. */ template - DEAL_II_CUDA_HOST_DEV Tensor & - operator*=(const OtherNumber &factor); + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor & + operator*=(const OtherNumber &factor); /** * Scale the vector by 1/factor. */ template - Tensor & - operator/=(const OtherNumber &factor); + DEAL_II_CONSTEXPR Tensor & + operator/=(const OtherNumber &factor); /** * Unary minus operator. Negate all entries of a tensor. */ - Tensor - operator-() const; + DEAL_II_CONSTEXPR Tensor + operator-() const; /** * Reset all values to zero. @@ -598,7 +601,7 @@ public: * and indeed the state where all elements have a zero value is the state * right after construction of such an object. */ - void + DEAL_II_CONSTEXPR void clear(); /** @@ -607,7 +610,7 @@ public: * tensors, this equals the usual l2 norm of the vector. */ - typename numbers::NumberTraits::real_type + DEAL_II_CONSTEXPR typename numbers::NumberTraits::real_type norm() const; /** @@ -616,8 +619,9 @@ public: * * @note This function can also be used in CUDA device code. */ - DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type - norm_square() const; + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV + typename numbers::NumberTraits::real_type + norm_square() const; /** * Fill a vector with all tensor elements. @@ -627,22 +631,22 @@ public: * fastest. */ template - void + DEAL_II_CONSTEXPR void unroll(Vector &result) const; /** * Return an unrolled index in the range [0,dim^rank-1] for the element of * the tensor indexed by the argument to the function. */ - static unsigned int + static DEAL_II_CONSTEXPR unsigned int component_to_unrolled_index(const TableIndices &indices); /** * Opposite of component_to_unrolled_index: For an index in the range * [0,dim^rank-1], return which set of indices it would correspond to. */ - static TableIndices - unrolled_to_component_indices(const unsigned int i); + static DEAL_II_CONSTEXPR TableIndices + unrolled_to_component_indices(const unsigned int i); /** * Determine an estimate for the memory consumption (in bytes) of this @@ -677,10 +681,18 @@ private: * Internal helper function for unroll. */ template - void + DEAL_II_CONSTEXPR void unroll_recursion(Vector &result, unsigned int & start_index) const; + /** + * This constructor is for internal use. It provides a way + * to create constexpr constructors for Tensor + */ + template + constexpr Tensor(const ArrayLike &initializer, + std_cxx14::index_sequence); + /** * Allow an arbitrary Tensor to access the underlying values. */ @@ -719,8 +731,8 @@ namespace internal return t; } - static Tensor - value(const T &t) + static DEAL_II_CONSTEXPR Tensor + value(const T &t) { Tensor tmp; tmp = t; @@ -737,16 +749,16 @@ namespace internal return t; } - static Tensor> - value(const T &t) + static DEAL_II_CONSTEXPR Tensor> + value(const T &t) { Tensor> tmp; tmp = internal::NumberType>::value(t); return tmp; } - static Tensor> - value(const VectorizedArray &t) + static DEAL_II_CONSTEXPR Tensor> + value(const VectorizedArray &t) { Tensor> tmp; tmp = t; @@ -760,34 +772,33 @@ namespace internal template -DEAL_II_CUDA_HOST_DEV inline Tensor<0, dim, Number>::Tensor() +constexpr DEAL_II_CUDA_HOST_DEV +Tensor<0, dim, Number>::Tensor() // Some auto-differentiable numbers need explicit - // zero initialization. - : value(internal::NumberType::value(0.0)) + // zero initialization such as adtl::adouble. + : Tensor{0.0} {} template template -inline Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer) -{ - value = internal::NumberType::value(initializer); -} +constexpr Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer) + : value(internal::NumberType::value(initializer)) +{} template template -inline Tensor<0, dim, Number>::Tensor(const Tensor<0, dim, OtherNumber> &p) -{ - value = p.value; -} +constexpr Tensor<0, dim, Number>::Tensor(const Tensor<0, dim, OtherNumber> &p) + : Tensor{p.value} +{} template -inline Number * +DEAL_II_CONSTEXPR inline Number * Tensor<0, dim, Number>::begin_raw() { return std::addressof(value); @@ -805,7 +816,7 @@ Tensor<0, dim, Number>::begin_raw() const template -inline Number * +DEAL_II_CONSTEXPR inline Number * Tensor<0, dim, Number>::end_raw() { return begin_raw() + n_independent_components; @@ -823,7 +834,8 @@ Tensor<0, dim, Number>::end_raw() const template -inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &() +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>:: + operator Number &() { // We cannot use Assert inside a CUDA kernel #ifndef __CUDA_ARCH__ @@ -835,8 +847,8 @@ inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &() template -inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>:: - operator const Number &() const +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>:: + operator const Number &() const { // We cannot use Assert inside a CUDA kernel #ifndef __CUDA_ARCH__ @@ -849,7 +861,7 @@ inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>:: template template -inline Tensor<0, dim, Number> & +DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p) { value = internal::NumberType::value(p); @@ -859,7 +871,7 @@ Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p) #ifdef __INTEL_COMPILER template -inline Tensor<0, dim, Number> & +DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator=(const Tensor<0, dim, Number> &p) { value = p.value; @@ -870,7 +882,7 @@ Tensor<0, dim, Number>::operator=(const Tensor<0, dim, Number> &p) template template -inline Tensor<0, dim, Number> & +DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator=(const OtherNumber &d) { value = internal::NumberType::value(d); @@ -880,7 +892,7 @@ Tensor<0, dim, Number>::operator=(const OtherNumber &d) template template -inline bool +DEAL_II_CONSTEXPR inline bool Tensor<0, dim, Number>::operator==(const Tensor<0, dim, OtherNumber> &p) const { #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING @@ -906,7 +918,7 @@ Tensor<0, dim, Number>::operator!=(const Tensor<0, dim, OtherNumber> &p) const template template -inline Tensor<0, dim, Number> & +DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator+=(const Tensor<0, dim, OtherNumber> &p) { value += p.value; @@ -916,7 +928,7 @@ Tensor<0, dim, Number>::operator+=(const Tensor<0, dim, OtherNumber> &p) template template -inline Tensor<0, dim, Number> & +DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator-=(const Tensor<0, dim, OtherNumber> &p) { value -= p.value; @@ -930,7 +942,7 @@ namespace internal namespace ComplexWorkaround { template - inline DEAL_II_CUDA_HOST_DEV void + DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV void multiply_assign_scalar(Number &val, const OtherNumber &s) { val *= s; @@ -938,7 +950,7 @@ namespace internal #ifdef __CUDA_ARCH__ template - inline DEAL_II_CUDA_HOST_DEV void + DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV void multiply_assign_scalar(std::complex &, const OtherNumber &) { printf("This function is not implemented for std::complex!\n"); @@ -951,7 +963,7 @@ namespace internal template template -inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> & +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator*=(const OtherNumber &s) { internal::ComplexWorkaround::multiply_assign_scalar(value, s); @@ -962,7 +974,7 @@ Tensor<0, dim, Number>::operator*=(const OtherNumber &s) template template -inline Tensor<0, dim, Number> & +DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> & Tensor<0, dim, Number>::operator/=(const OtherNumber &s) { value /= s; @@ -979,7 +991,7 @@ Tensor<0, dim, Number>::operator-() const template -inline typename Tensor<0, dim, Number>::real_type +DEAL_II_CONSTEXPR inline typename Tensor<0, dim, Number>::real_type Tensor<0, dim, Number>::norm() const { Assert(dim != 0, @@ -1003,7 +1015,7 @@ inline typename Tensor<0, dim, Number>::real_type DEAL_II_CUDA_HOST_DEV template template -inline void +DEAL_II_CONSTEXPR inline void Tensor<0, dim, Number>::unroll_recursion(Vector &result, unsigned int & index) const { @@ -1015,7 +1027,7 @@ Tensor<0, dim, Number>::unroll_recursion(Vector &result, template -inline void +DEAL_II_CONSTEXPR inline void Tensor<0, dim, Number>::clear() { // Some auto-differentiable numbers need explicit @@ -1035,42 +1047,45 @@ Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int) /*-------------------- Inline functions: Tensor --------------------*/ - template -inline DEAL_II_ALWAYS_INLINE -Tensor::Tensor(const array_type &initializer) +template +DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( + const ArrayLike &initializer, + std_cxx14::index_sequence) + : values{Tensor(initializer[indices])...} { - for (unsigned int i = 0; i < dim; ++i) - values[i] = Tensor(initializer[i]); + static_assert(sizeof...(indices) == dim, + "dim should match the number of indices"); } +template +DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( + const array_type &initializer) + : Tensor(initializer, std_cxx14::make_index_sequence{}) +{} + + template template -inline DEAL_II_ALWAYS_INLINE -Tensor::Tensor( +DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( const Tensor &initializer) -{ - for (unsigned int i = 0; i != dim; ++i) - values[i] = Tensor(initializer[i]); -} + : Tensor(initializer, std_cxx14::make_index_sequence{}) +{} template template -inline DEAL_II_ALWAYS_INLINE -Tensor::Tensor( +DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( const Tensor<1, dim, Tensor> &initializer) -{ - for (unsigned int i = 0; i < dim; ++i) - values[i] = initializer[i]; -} + : Tensor(initializer, std_cxx14::make_index_sequence{}) +{} template template -constexpr DEAL_II_ALWAYS_INLINE Tensor:: - operator Tensor<1, dim, Tensor>() const +DEAL_II_ALWAYS_INLINE constexpr Tensor:: +operator Tensor<1, dim, Tensor>() const { return Tensor<1, dim, Tensor>(values); } @@ -1082,10 +1097,11 @@ namespace internal namespace TensorSubscriptor { template - inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV ArrayElementType & - subscript(ArrayElementType * values, - const unsigned int i, - std::integral_constant) + DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE + DEAL_II_CUDA_HOST_DEV ArrayElementType & + subscript(ArrayElementType * values, + const unsigned int i, + std::integral_constant) { // We cannot use Assert in a CUDA kernel #ifndef __CUDA_ARCH__ @@ -1094,9 +1110,19 @@ namespace internal return values[i]; } + // The variables within this struct will be referenced in the next function. + // It is a workaround that allows returning a reference to a static variable + // while allowing constexpr evaluation of the function. + // It has to be defined outside the function because constexpr functions + // cannot define static variables + template + struct Uninitialized + { + static ArrayElementType value; + }; template - ArrayElementType & + DEAL_II_CONSTEXPR inline ArrayElementType & subscript(ArrayElementType *, const unsigned int, std::integral_constant) @@ -1105,15 +1131,14 @@ namespace internal false, ExcMessage( "Cannot access elements of an object of type Tensor.")); - static ArrayElementType t; - return t; + return Uninitialized::value; } } // namespace TensorSubscriptor } // namespace internal template -inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV typename Tensor::value_type &Tensor:: operator[](const unsigned int i) { @@ -1123,9 +1148,11 @@ inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV template -constexpr DEAL_II_ALWAYS_INLINE - DEAL_II_CUDA_HOST_DEV const typename Tensor::value_type & - Tensor::operator[](const unsigned int i) const +DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV const typename Tensor< + rank_, + dim, + Number>::value_type &Tensor:: + operator[](const unsigned int i) const { return dealii::internal::TensorSubscriptor::subscript( values, i, std::integral_constant()); @@ -1133,8 +1160,8 @@ constexpr DEAL_II_ALWAYS_INLINE template -inline const Number &Tensor:: - operator[](const TableIndices &indices) const +DEAL_II_CONSTEXPR inline const Number &Tensor:: + operator[](const TableIndices &indices) const { Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor")); @@ -1145,8 +1172,8 @@ inline const Number &Tensor:: template -inline Number &Tensor:: - operator[](const TableIndices &indices) +DEAL_II_CONSTEXPR inline Number &Tensor:: + operator[](const TableIndices &indices) { Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor")); @@ -1157,7 +1184,7 @@ inline Number &Tensor:: template -inline Number * +DEAL_II_CONSTEXPR inline Number * Tensor::begin_raw() { return std::addressof( @@ -1177,7 +1204,7 @@ Tensor::begin_raw() const template -inline Number * +DEAL_II_CONSTEXPR inline Number * Tensor::end_raw() { return begin_raw() + n_independent_components; @@ -1196,7 +1223,7 @@ Tensor::end_raw() const template template -inline DEAL_II_ALWAYS_INLINE Tensor & +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor & Tensor::operator=(const Tensor &t) { if (dim > 0) @@ -1206,7 +1233,7 @@ Tensor::operator=(const Tensor &t) template -inline DEAL_II_ALWAYS_INLINE Tensor & +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor & Tensor::operator=(const Number &d) { Assert(numbers::value_is_zero(d), @@ -1221,7 +1248,7 @@ Tensor::operator=(const Number &d) template template -inline bool +DEAL_II_CONSTEXPR inline bool Tensor:: operator==(const Tensor &p) const { @@ -1239,7 +1266,7 @@ operator==(const Tensor &p) const // implement this function here template <> template <> -inline bool +DEAL_II_CONSTEXPR inline bool Tensor<1, 0, double>::operator==(const Tensor<1, 0, double> &) const { return true; @@ -1258,7 +1285,7 @@ operator!=(const Tensor &p) const template template -inline Tensor & +DEAL_II_CONSTEXPR inline Tensor & Tensor::operator+=(const Tensor &p) { for (unsigned int i = 0; i < dim; ++i) @@ -1269,7 +1296,7 @@ Tensor::operator+=(const Tensor &p) template template -inline Tensor & +DEAL_II_CONSTEXPR inline Tensor & Tensor::operator-=(const Tensor &p) { for (unsigned int i = 0; i < dim; ++i) @@ -1280,7 +1307,7 @@ Tensor::operator-=(const Tensor &p) template template -inline DEAL_II_CUDA_HOST_DEV Tensor & +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor & Tensor::operator*=(const OtherNumber &s) { for (unsigned int i = 0; i < dim; ++i) @@ -1291,7 +1318,7 @@ Tensor::operator*=(const OtherNumber &s) template template -inline Tensor & +DEAL_II_CONSTEXPR inline Tensor & Tensor::operator/=(const OtherNumber &s) { for (unsigned int i = 0; i < dim; ++i) @@ -1301,7 +1328,7 @@ Tensor::operator/=(const OtherNumber &s) template -inline Tensor +DEAL_II_CONSTEXPR inline Tensor Tensor::operator-() const { Tensor tmp; @@ -1314,7 +1341,7 @@ Tensor::operator-() const template -inline typename numbers::NumberTraits::real_type +DEAL_II_CONSTEXPR inline typename numbers::NumberTraits::real_type Tensor::norm() const { return std::sqrt(norm_square()); @@ -1322,8 +1349,9 @@ Tensor::norm() const template -inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type -Tensor::norm_square() const +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV + typename numbers::NumberTraits::real_type + Tensor::norm_square() const { typename numbers::NumberTraits::real_type s = internal::NumberType< typename numbers::NumberTraits::real_type>::value(0.0); @@ -1336,7 +1364,7 @@ Tensor::norm_square() const template template -inline void +DEAL_II_CONSTEXPR inline void Tensor::unroll(Vector &result) const { AssertDimension(result.size(), @@ -1349,7 +1377,7 @@ Tensor::unroll(Vector &result) const template template -inline void +DEAL_II_CONSTEXPR inline void Tensor::unroll_recursion(Vector &result, unsigned int & index) const { @@ -1359,7 +1387,7 @@ Tensor::unroll_recursion(Vector &result, template -inline unsigned int +DEAL_II_CONSTEXPR inline unsigned int Tensor::component_to_unrolled_index( const TableIndices &indices) { @@ -1412,7 +1440,7 @@ namespace internal template -inline TableIndices +DEAL_II_CONSTEXPR inline TableIndices Tensor::unrolled_to_component_indices(const unsigned int i) { Assert(i < n_independent_components, @@ -1433,7 +1461,7 @@ Tensor::unrolled_to_component_indices(const unsigned int i) template -inline void +DEAL_II_CONSTEXPR inline void Tensor::clear() { for (unsigned int i = 0; i < dim; ++i) @@ -1619,7 +1647,7 @@ constexpr DEAL_II_ALWAYS_INLINE * @relatesalso Tensor */ template -inline DEAL_II_ALWAYS_INLINE +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor -inline Tensor< - rank, - dim, - typename ProductType::type>::type> -operator/(const Tensor &t, const OtherNumber &factor) +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE + Tensor::type>::type> + operator/(const Tensor &t, const OtherNumber &factor) { // recurse over the base objects Tensor::type> tt; @@ -1688,7 +1716,7 @@ operator/(const Tensor &t, const OtherNumber &factor) * @relatesalso Tensor */ template -inline DEAL_II_ALWAYS_INLINE +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor::type> operator+(const Tensor & p, const Tensor &q) @@ -1710,7 +1738,7 @@ inline DEAL_II_ALWAYS_INLINE * @relatesalso Tensor */ template -inline DEAL_II_ALWAYS_INLINE +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor::type> operator-(const Tensor & p, const Tensor &q) @@ -1759,7 +1787,7 @@ template -inline DEAL_II_ALWAYS_INLINE +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE typename Tensor::type>::tensor_type @@ -1816,7 +1844,7 @@ template -inline DEAL_II_ALWAYS_INLINE +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE typename Tensor::type>::tensor_type @@ -1890,7 +1918,7 @@ template -inline +DEAL_II_CONSTEXPR inline typename Tensor::type>::tensor_type @@ -1972,9 +2000,10 @@ inline * @author Matthias Maier, 2015 */ template -inline DEAL_II_ALWAYS_INLINE typename ProductType::type -scalar_product(const Tensor & left, - const Tensor &right) +DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE + typename ProductType::type + scalar_product(const Tensor & left, + const Tensor &right) { typename ProductType::type result; TensorAccessors::contract(result, left, right); @@ -2009,10 +2038,11 @@ template