From: Maximilian Bergbauer Date: Mon, 14 Aug 2023 10:49:40 +0000 (+0200) Subject: Implement ARM Neon intrinsics X-Git-Tag: relicensing~552^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a7e953982c1b3e4840521beba5b4700a501bcd01;p=dealii.git Implement ARM Neon intrinsics --- diff --git a/cmake/checks/check_01_cpu_features.cmake b/cmake/checks/check_01_cpu_features.cmake index dbd7a395af..eff830dcd3 100644 --- a/cmake/checks/check_01_cpu_features.cmake +++ b/cmake/checks/check_01_cpu_features.cmake @@ -28,6 +28,7 @@ # DEAL_II_HAVE_AVX (*) # DEAL_II_HAVE_AVX512 (*) # DEAL_II_HAVE_ALTIVEC (*) +# DEAL_II_HAVE_ARM_NEON (*) # DEAL_II_HAVE_OPENMP_SIMD (*) # DEAL_II_VECTORIZATION_WIDTH_IN_BITS # DEAL_II_OPENMP_SIMD_PRAGMA @@ -72,7 +73,7 @@ if(DEAL_II_ALLOW_PLATFORM_INTROSPECTION) # CMAKE_REQUIRED_FLAGS changes.. # unset_if_changed(CHECK_CPU_FEATURES_FLAGS_SAVED "${CMAKE_REQUIRED_FLAGS}" - DEAL_II_HAVE_SSE2 DEAL_II_HAVE_AVX DEAL_II_HAVE_AVX512 DEAL_II_HAVE_ALTIVEC + DEAL_II_HAVE_SSE2 DEAL_II_HAVE_AVX DEAL_II_HAVE_AVX512 DEAL_II_HAVE_ALTIVEC DEAL_II_HAVE_ARM_NEON ) CHECK_CXX_SOURCE_RUNS( @@ -245,6 +246,40 @@ if(DEAL_II_ALLOW_PLATFORM_INTROSPECTION) " DEAL_II_HAVE_ALTIVEC) + CHECK_CXX_SOURCE_RUNS( + " + #ifndef __ARM_NEON + #error Preprocessor flag not found + #endif + #include + #include + int main() + { + float64x2_t a, b; + const unsigned int vector_bytes = sizeof(float64x2_t); + const int n_vectors = vector_bytes/sizeof(double); + float64x2_t * data = + reinterpret_cast(aligned_alloc (vector_bytes, 2*vector_bytes)); + double * ptr = reinterpret_cast(&a); + ptr[0] = static_cast(1.0); + for (int i=1; i(2.25)); + data[0] = vaddq_f64 (a, b); + data[1] = vmulq_f64 (b, data[0]); + ptr = reinterpret_cast(&data[1]); + int return_value = 0; + if (ptr[0] != 7.3125) + return_value = 1; + for (int i=1; i +(Maximilian Bergbauer, Peter Munch, 2023/08/24) diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index 9203616264..b320580d2d 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -138,6 +138,8 @@ #define DEAL_II_COMPILER_VECTORIZATION_LEVEL 0 #endif +#define DEAL_II_HAVE_ARM_NEON @DEAL_II_HAVE_ARM_NEON@ + #define DEAL_II_OPENMP_SIMD_PRAGMA @DEAL_II_OPENMP_SIMD_PRAGMA@ diff --git a/include/deal.II/base/numbers.h b/include/deal.II/base/numbers.h index 360b7cda76..896fb4c375 100644 --- a/include/deal.II/base/numbers.h +++ b/include/deal.II/base/numbers.h @@ -130,6 +130,8 @@ namespace internal 8; #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__) 4; +#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON) + 4; #else 1; #endif diff --git a/include/deal.II/base/vectorization.h b/include/deal.II/base/vectorization.h index 20059d4f88..aaeead517f 100644 --- a/include/deal.II/base/vectorization.h +++ b/include/deal.II/base/vectorization.h @@ -69,7 +69,9 @@ # undef vector # undef pixel # undef bool -# else +# elif defined(__ARM_NEON) +# include +# elif defined(__x86_64__) # include # endif @@ -966,6 +968,565 @@ vectorized_transpose_and_store(const bool add_into, #ifndef DOXYGEN +# if defined(DEAL_II_HAVE_ARM_NEON) && defined(__ARM_NEON) + +/** + * Specialization for double and ARM Neon. + */ +template <> +class VectorizedArray + : public VectorizedArrayBase, 2> +{ +public: + /** + * This gives the type of the array elements. + */ + using value_type = double; + + /** + * Default empty constructor, leaving the data in an uninitialized state + * similar to float/double. + */ + VectorizedArray() = default; + + /** + * Construct an array with the given scalar broadcast to all lanes. + */ + VectorizedArray(const double scalar) + { + this->operator=(scalar); + } + + /** + * Construct an array with the given initializer list. + */ + template + VectorizedArray(const std::initializer_list &list) + : VectorizedArrayBase, 2>(list) + {} + + /** + * This function can be used to set all data fields to a given scalar. + */ + VectorizedArray & + operator=(const double x) & + { + data = vdupq_n_f64(x); + return *this; + } + + /** + * Assign a scalar to the current object. This overload is used for + * rvalue references; because it does not make sense to assign + * something to a temporary, the function is deleted. + */ + VectorizedArray & + operator=(const double scalar) && = delete; + + /** + * Access operator. + */ + double & + operator[](const unsigned int comp) + { + return *(reinterpret_cast(&data) + comp); + } + + /** + * Constant access operator. + */ + const double & + operator[](const unsigned int comp) const + { + return *(reinterpret_cast(&data) + comp); + } + + /** + * Addition. + */ + VectorizedArray & + operator+=(const VectorizedArray &vec) + { + data = vaddq_f64(data, vec.data); + return *this; + } + + /** + * Subtraction. + */ + VectorizedArray & + operator-=(const VectorizedArray &vec) + { + data = vsubq_f64(data, vec.data); + return *this; + } + + /** + * Multiplication. + */ + VectorizedArray & + operator*=(const VectorizedArray &vec) + { + data = vmulq_f64(data, vec.data); + return *this; + } + + /** + * Division. + */ + VectorizedArray & + operator/=(const VectorizedArray &vec) + { + data = vdivq_f64(data, vec.data); + return *this; + } + + /** + * Load @p size() from memory into the calling class, starting at + * the given address. The memory need not be aligned by 16 bytes, as opposed + * to casting a double address to VectorizedArray*. + */ + void + load(const double *ptr) + { + data = vld1q_f64(ptr); + } + + DEAL_II_ALWAYS_INLINE + void + load(const float *ptr) + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < 2; ++i) + data[i] = ptr[i]; + } + + /** + * Write the content of the calling class into memory in form of @p + * size() to the given address. The memory need not be aligned by + * 16 bytes, as opposed to casting a double address to + * VectorizedArray*. + */ + void + store(double *ptr) const + { + vst1q_f64(ptr, data); + } + + DEAL_II_ALWAYS_INLINE + void + store(float *ptr) const + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < 2; ++i) + ptr[i] = data[i]; + } + + /** + * @copydoc VectorizedArray::streaming_store() + * @note Memory must be aligned by 16 bytes. + */ + DEAL_II_ALWAYS_INLINE + void + streaming_store(double *ptr) const + { + Assert(reinterpret_cast(ptr) % 16 == 0, + ExcMessage("Memory not aligned")); + vst1q_f64(ptr, data); + } + + /** + * Load @p size() from memory into the calling class, starting at + * the given address and with given offsets, each entry from the offset + * providing one element of the vectorized array. + * + * This operation corresponds to the following code (but uses a more + * efficient implementation in case the hardware allows for that): + * @code + * for (unsigned int v=0; v::size(); ++v) + * this->operator[](v) = base_ptr[offsets[v]]; + * @endcode + */ + void + gather(const double *base_ptr, const unsigned int *offsets) + { + for (unsigned int i = 0; i < 2; ++i) + *(reinterpret_cast(&data) + i) = base_ptr[offsets[i]]; + } + + /** + * Write the content of the calling class into memory in form of @p + * size() to the given address and the given offsets, filling the + * elements of the vectorized array into each offset. + * + * This operation corresponds to the following code (but uses a more + * efficient implementation in case the hardware allows for that): + * @code + * for (unsigned int v=0; v::size(); ++v) + * base_ptr[offsets[v]] = this->operator[](v); + * @endcode + */ + void + scatter(const unsigned int *offsets, double *base_ptr) const + { + for (unsigned int i = 0; i < 2; ++i) + base_ptr[offsets[i]] = *(reinterpret_cast(&data) + i); + } + + /** + * Returns sum over entries of the data field, $\sum_{i=1}^{\text{size}()} + * this->data[i]$. + */ + double + sum() const + { + return vaddvq_f64(data); + } + + /** + * Actual data field. To be consistent with the standard layout type and to + * enable interaction with external SIMD functionality, this member is + * declared public. + */ + mutable float64x2_t data; + +private: + /** + * Return the square root of this field. Not for use in user code. Use + * sqrt(x) instead. + */ + VectorizedArray + get_sqrt() const + { + VectorizedArray res; + res.data = vsqrtq_f64(data); + return res; + } + + /** + * Return the absolute value of this field. Not for use in user code. Use + * abs(x) instead. + */ + VectorizedArray + get_abs() const + { + VectorizedArray res; + res.data = vabsq_f64(data); + return res; + } + + /** + * Return the component-wise maximum of this field and another one. Not for + * use in user code. Use max(x,y) instead. + */ + VectorizedArray + get_max(const VectorizedArray &other) const + { + VectorizedArray res; + res.data = vmaxq_f64(data, other.data); + return res; + } + + /** + * Return the component-wise minimum of this field and another one. Not for + * use in user code. Use min(x,y) instead. + */ + VectorizedArray + get_min(const VectorizedArray &other) const + { + VectorizedArray res; + res.data = vminq_f64(data, other.data); + return res; + } + + // Make a few functions friends. + template + friend VectorizedArray + std::sqrt(const VectorizedArray &); + template + friend VectorizedArray + std::abs(const VectorizedArray &); + template + friend VectorizedArray + std::max(const VectorizedArray &, + const VectorizedArray &); + template + friend VectorizedArray + std::min(const VectorizedArray &, + const VectorizedArray &); +}; + +/** + * Specialization for float and ARM Neon. + */ +template <> +class VectorizedArray + : public VectorizedArrayBase, 4> +{ +public: + /** + * This gives the type of the array elements. + */ + using value_type = float; + + /** + * Default empty constructor, leaving the data in an uninitialized state + * similar to float/double. + */ + VectorizedArray() = default; + + /** + * Construct an array with the given scalar broadcast to all lanes. + */ + VectorizedArray(const float scalar) + { + this->operator=(scalar); + } + + /** + * Construct an array with the given initializer list. + */ + template + VectorizedArray(const std::initializer_list &list) + : VectorizedArrayBase, 4>(list) + {} + + /** + * This function can be used to set all data fields to a given scalar. + */ + VectorizedArray & + operator=(const float x) & + { + data = vdupq_n_f32(x); + return *this; + } + + /** + * Assign a scalar to the current object. This overload is used for + * rvalue references; because it does not make sense to assign + * something to a temporary, the function is deleted. + */ + VectorizedArray & + operator=(const float scalar) && = delete; + + /** + * Access operator. + */ + value_type & + operator[](const unsigned int comp) + { + return *(reinterpret_cast(&data) + comp); + } + + /** + * Constant access operator. + */ + const value_type & + operator[](const unsigned int comp) const + { + return *(reinterpret_cast(&data) + comp); + } + + /** + * Addition. + */ + VectorizedArray & + operator+=(const VectorizedArray &vec) + { + data = vaddq_f32(data, vec.data); + return *this; + } + + /** + * Subtraction. + */ + VectorizedArray & + operator-=(const VectorizedArray &vec) + { + data = vsubq_f32(data, vec.data); + return *this; + } + + /** + * Multiplication. + */ + VectorizedArray & + operator*=(const VectorizedArray &vec) + { + data = vmulq_f32(data, vec.data); + return *this; + } + + /** + * Division. + */ + VectorizedArray & + operator/=(const VectorizedArray &vec) + { + data = vdivq_f32(data, vec.data); + return *this; + } + + /** + * Load @p size() from memory into the calling class, starting at + * the given address. The memory need not be aligned by 16 bytes, as opposed + * to casting a float address to VectorizedArray*. + */ + void + load(const float *ptr) + { + data = vld1q_f32(ptr); + } + + /** + * Write the content of the calling class into memory in form of @p + * size() to the given address. The memory need not be aligned by + * 16 bytes, as opposed to casting a float address to + * VectorizedArray*. + */ + void + store(float *ptr) const + { + vst1q_f32(ptr, data); + } + + /** + * @copydoc VectorizedArray::streaming_store() + * @note Memory must be aligned by 16 bytes. + */ + DEAL_II_ALWAYS_INLINE + void + streaming_store(float *ptr) const + { + Assert(reinterpret_cast(ptr) % 16 == 0, + ExcMessage("Memory not aligned")); + vst1q_f32(ptr, data); + } + + /** + * Load @p size() from memory into the calling class, starting at + * the given address and with given offsets, each entry from the offset + * providing one element of the vectorized array. + * + * This operation corresponds to the following code (but uses a more + * efficient implementation in case the hardware allows for that): + * @code + * for (unsigned int v=0; v::size(); ++v) + * this->operator[](v) = base_ptr[offsets[v]]; + * @endcode + */ + void + gather(const float *base_ptr, const unsigned int *offsets) + { + for (unsigned int i = 0; i < 4; ++i) + *(reinterpret_cast(&data) + i) = base_ptr[offsets[i]]; + } + + /** + * Write the content of the calling class into memory in form of @p + * size() to the given address and the given offsets, filling the + * elements of the vectorized array into each offset. + * + * This operation corresponds to the following code (but uses a more + * efficient implementation in case the hardware allows for that): + * @code + * for (unsigned int v=0; v::size(); ++v) + * base_ptr[offsets[v]] = this->operator[](v); + * @endcode + */ + void + scatter(const unsigned int *offsets, float *base_ptr) const + { + for (unsigned int i = 0; i < 4; ++i) + base_ptr[offsets[i]] = *(reinterpret_cast(&data) + i); + } + + /** + * Returns sum over entries of the data field, $\sum_{i=1}^{\text{size}()} + * this->data[i]$. + */ + float + sum() const + { + return vaddvq_f32(data); + } + + /** + * Actual data field. To be consistent with the standard layout type and to + * enable interaction with external SIMD functionality, this member is + * declared public. + */ + mutable float32x4_t data; + +private: + /** + * Return the square root of this field. Not for use in user code. Use + * sqrt(x) instead. + */ + VectorizedArray + get_sqrt() const + { + VectorizedArray res; + res.data = vsqrtq_f32(data); + return res; + } + + /** + * Return the absolute value of this field. Not for use in user code. Use + * abs(x) instead. + */ + VectorizedArray + get_abs() const + { + VectorizedArray res; + res.data = vabsq_f32(data); + return res; + } + + /** + * Return the component-wise maximum of this field and another one. Not for + * use in user code. Use max(x,y) instead. + */ + VectorizedArray + get_max(const VectorizedArray &other) const + { + VectorizedArray res; + res.data = vmaxq_f32(data, other.data); + return res; + } + + /** + * Return the component-wise minimum of this field and another one. Not for + * use in user code. Use min(x,y) instead. + */ + VectorizedArray + get_min(const VectorizedArray &other) const + { + VectorizedArray res; + res.data = vminq_f32(data, other.data); + return res; + } + + // Make a few functions friends. + template + friend VectorizedArray + std::sqrt(const VectorizedArray &); + template + friend VectorizedArray + std::abs(const VectorizedArray &); + template + friend VectorizedArray + std::max(const VectorizedArray &, + const VectorizedArray &); + template + friend VectorizedArray + std::min(const VectorizedArray &, + const VectorizedArray &); +}; + + +# endif + # if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__) /** @@ -1466,10 +2027,6 @@ public: */ using value_type = float; - /** - * This function can be used to set all data fields to a given scalar. - */ - /** * Default empty constructor, leaving the data in an uninitialized state * similar to float/double. @@ -1492,6 +2049,9 @@ public: : VectorizedArrayBase, 4>(list) {} + /** + * This function can be used to set all data fields to a given scalar. + */ DEAL_II_ALWAYS_INLINE VectorizedArray & operator=(const float x) & @@ -5636,6 +6196,89 @@ compare_and_apply_mask(const VectorizedArray &left, return result; } +# endif + +# if defined(DEAL_II_HAVE_ARM_NEON) && defined(__ARM_NEON) + +template +DEAL_II_ALWAYS_INLINE inline VectorizedArray +compare_and_apply_mask(const VectorizedArray &left, + const VectorizedArray &right, + const VectorizedArray &true_values, + const VectorizedArray &false_values) +{ + uint32x4_t mask; + switch (predicate) + { + case SIMDComparison::equal: + mask = vceqq_f32(left.data, right.data); + break; + case SIMDComparison::not_equal: + mask = vmvnq_u32(vceqq_f32(left.data, right.data)); + break; + case SIMDComparison::less_than: + mask = vcltq_f32(left.data, right.data); + break; + case SIMDComparison::less_than_or_equal: + mask = vcleq_f32(left.data, right.data); + break; + case SIMDComparison::greater_than: + mask = vcgtq_f32(left.data, right.data); + break; + case SIMDComparison::greater_than_or_equal: + mask = vcgeq_f32(left.data, right.data); + break; + } + + VectorizedArray result; + result.data = vreinterpretq_f32_u32(vorrq_u32( + vandq_u32(mask, vreinterpretq_u32_f32(true_values.data)), + vandq_u32(vmvnq_u32(mask), vreinterpretq_u32_f32(false_values.data)))); + + return result; +} + + +template +DEAL_II_ALWAYS_INLINE inline VectorizedArray +compare_and_apply_mask(const VectorizedArray &left, + const VectorizedArray &right, + const VectorizedArray &true_values, + const VectorizedArray &false_values) +{ + uint64x2_t mask; + switch (predicate) + { + case SIMDComparison::equal: + mask = vceqq_f64(left.data, right.data); + break; + case SIMDComparison::not_equal: + mask = vreinterpretq_u64_u32( + vmvnq_u32(vreinterpretq_u32_u64(vceqq_f64(left.data, right.data)))); + break; + case SIMDComparison::less_than: + mask = vcltq_f64(left.data, right.data); + break; + case SIMDComparison::less_than_or_equal: + mask = vcleq_f64(left.data, right.data); + break; + case SIMDComparison::greater_than: + mask = vcgtq_f64(left.data, right.data); + break; + case SIMDComparison::greater_than_or_equal: + mask = vcgeq_f64(left.data, right.data); + break; + } + + VectorizedArray result; + result.data = vreinterpretq_f64_u64(vorrq_u64( + vandq_u64(mask, vreinterpretq_u64_f64(true_values.data)), + vandq_u64(vreinterpretq_u64_u32(vmvnq_u32(vreinterpretq_u32_u64(mask))), + vreinterpretq_u64_f64(false_values.data)))); + + return result; +} + # endif #endif // DOXYGEN