From a2646f58f20015952685b6d4c29ff399f1c7a428 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 24 Jun 2024 18:02:44 +0200 Subject: [PATCH] Enable matrix-free loops with ArrayView arguments --- include/deal.II/matrix_free/type_traits.h | 88 +++++++++++-------- .../matrix_free/vector_access_internal.h | 32 ++++--- .../matrix_free/fe_evaluation_type_traits.cc | 6 +- ...tion_type_traits.with_trilinos=true.output | 1 - 4 files changed, 70 insertions(+), 57 deletions(-) diff --git a/include/deal.II/matrix_free/type_traits.h b/include/deal.II/matrix_free/type_traits.h index 9ed8fadf55..7b7ebb8936 100644 --- a/include/deal.II/matrix_free/type_traits.h +++ b/include/deal.II/matrix_free/type_traits.h @@ -28,13 +28,14 @@ DEAL_II_NAMESPACE_OPEN #ifndef DOXYGEN + namespace internal { // // type traits for FEEvaluation // - // a helper type-trait that leverage SFINAE to figure out if type T has + // a helper type-trait leveraging SFINAE to figure out if type T has // ... T::local_element() const template using local_element_t = decltype(std::declval().local_element(0)); @@ -44,7 +45,7 @@ namespace internal - // a helper type-trait that leverage SFINAE to figure out if type T has + // a helper type-trait leveraging SFINAE to figure out if type T has // void T::add_local_element(const uint, const typename T::value_type) template using add_local_element_t = @@ -56,7 +57,7 @@ namespace internal - // a helper type-trait that leverage SFINAE to figure out if type T has + // a helper type-trait leveraging SFINAE to figure out if type T has // void T::set_local_element(const uint, const typename T::value_type) template using set_local_element_t = @@ -67,7 +68,6 @@ namespace internal is_supported_operation; - // same as above to check // bool T::partitioners_are_compatible(const Utilities::MPI::Partitioner &) // const @@ -103,27 +103,6 @@ namespace internal is_supported_operation; - - // type trait for vector T and Number to see if - // we can do vectorized load/save. - // for VectorReader and VectorDistributorLocalToGlobal we assume that - // if both begin() and local_element() - // exist, then begin() + offset == local_element(offset) - template - struct is_vectorizable - { - static const bool value = - has_begin && - (has_local_element || - is_serial_vector>::value) && - std::is_same_v; - }; - - // We need to have a separate declaration for static const members - template - const bool is_vectorizable::value; - - // // type-traits for Matrix-Free // @@ -176,17 +155,6 @@ namespace internal constexpr bool has_communication_block_size = is_supported_operation; - - - // type trait for vector T to see if - // we need to do any data exchange for this vector type at all. - // is_serial_vector<> would have been enough, but in some circumstances - // (like calculation of diagonals for matrix-free operators) - // a dummy InVector == unsigned int is provided. - // Thus we have to treat this case as well. - template - using not_parallel_vector_t = std::bool_constant::value>; - /** * A predicate stating whether something is a vector type. We test this * by seeing whether the `is_serial_vector` type is declared for the @@ -212,6 +180,54 @@ namespace internal constexpr bool is_not_parallel_vector = (is_supported_operation == false) || (is_supported_operation == true); + + + + /** + * A type trait for vector T indicating serial vectors with access through + * operator[], which are the deal.II serial vectors and ArrayView + */ + template + struct is_serial_vector_or_array + { + static const bool value = + is_supported_operation; + }; + + /** + * A type trait for vector T indicating serial vectors with access through + * operator[], which are the deal.II serial vectors and ArrayView + */ + template + struct is_serial_vector_or_array> + { + static const bool value = true; + }; + + // We need to have a separate declaration for static const members + template + const bool is_serial_vector_or_array::value; + + + + // type trait for vector T and Number to see if + // we can do vectorized load/save. + // for VectorReader and VectorDistributorLocalToGlobal we assume that + // if both begin() and local_element() + // exist, then begin() + offset == local_element(offset) + template + struct is_vectorizable + { + static const bool value = + has_begin && + (has_local_element || is_serial_vector_or_array::value) && + std::is_same_v; + }; + + // We need to have a separate declaration for static const members + template + const bool is_vectorizable::value; + } // namespace internal #endif diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h index b371c650e3..b156b25a2b 100644 --- a/include/deal.II/matrix_free/vector_access_internal.h +++ b/include/deal.II/matrix_free/vector_access_internal.h @@ -35,30 +35,28 @@ namespace internal { // below we use type-traits from matrix-free/type_traits.h - // access to generic const vectors that have operator (). - // FIXME: this is wrong for Trilinos/PETSc MPI vectors - // where we should first do Partitioner::local_to_global() - template < - typename VectorType, - std::enable_if_t, VectorType> * = nullptr> + + + // access to serial const vectors that have operator[]. + template ::value, + VectorType> * = nullptr> inline typename VectorType::value_type vector_access(const VectorType &vec, const unsigned int entry) { - return vec(entry); + return vec[entry]; } - // access to generic non-const vectors that have operator (). - // FIXME: this is wrong for Trilinos/PETSc MPI vectors - // where we should first do Partitioner::local_to_global() - template < - typename VectorType, - std::enable_if_t, VectorType> * = nullptr> + // access to serial non-const vectors that have operator[]. + template ::value, + VectorType> * = nullptr> inline typename VectorType::value_type & vector_access(VectorType &vec, const unsigned int entry) { - return vec(entry); + return vec[entry]; } @@ -136,7 +134,7 @@ namespace internal const types::global_dof_index entry, const typename VectorType::value_type &val) { - vec(entry) += val; + vec[entry] += val; } @@ -483,7 +481,7 @@ namespace internal const VectorType &vec, Number &res) const { - res = vec(index); + res = vec[index]; } @@ -953,7 +951,7 @@ namespace internal VectorType &vec, Number &res) const { - vec(index) = res; + vec[index] = res; } diff --git a/tests/matrix_free/fe_evaluation_type_traits.cc b/tests/matrix_free/fe_evaluation_type_traits.cc index 927804c12f..a217ab5a40 100644 --- a/tests/matrix_free/fe_evaluation_type_traits.cc +++ b/tests/matrix_free/fe_evaluation_type_traits.cc @@ -71,14 +71,14 @@ public: using value_type = Number; Number - operator()(const unsigned int) const + operator[](const unsigned int) const { deallog << "Dummy2::operator() const" << std::endl; return Number(); } Number & - operator()(const unsigned int) + operator[](const unsigned int) { deallog << "Dummy2::operator()" << std::endl; return dummy; @@ -112,7 +112,7 @@ main() // we expect local_element() to be called deallog << "internal::vector_access:" << std::endl; internal::vector_access(dummy, 0); - internal::vector_access(dummy2, 0); + // internal::vector_access(dummy2, 0); // now check has_partitioners_are_compatible: diff --git a/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output b/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output index b7e821fc99..ed4d98ded6 100644 --- a/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output +++ b/tests/matrix_free/fe_evaluation_type_traits.with_trilinos=true.output @@ -7,7 +7,6 @@ DEAL::Dummy2 = 0 DEAL::Vector = 0 DEAL::internal::vector_access: DEAL::Dummy::local_element() -DEAL::Dummy2::operator() DEAL::has_partitioners_are_compatible: DEAL::LinearAlgebra::distributed::Vector = 1 DEAL::TrilinosWrappers::MPI::Vector = 0 -- 2.39.5