From ba75576e4a89ddfbfe238783fe54432298a9f041 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 23 Oct 2020 08:50:52 +0200 Subject: [PATCH] MatrixFree: pass SM vectors to functions --- .../deal.II/matrix_free/evaluation_kernels.h | 52 +++++--- .../matrix_free/evaluation_template_factory.h | 14 +- .../evaluation_template_factory.templates.h | 16 ++- include/deal.II/matrix_free/fe_evaluation.h | 125 ++++++++++++++---- include/deal.II/matrix_free/type_traits.h | 23 ++++ 5 files changed, 170 insertions(+), 60 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index a27a0e3a9d..ce85c837ba 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3186,9 +3186,10 @@ namespace internal { template static bool - run(const unsigned int n_components, - const unsigned int n_face_orientations, - const Number2 * src_ptr, + run(const unsigned int n_components, + const unsigned int n_face_orientations, + const Number2 * src_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -3211,9 +3212,12 @@ namespace internal return false; } + (void)sm_ptr; + Processor p(n_components, false, src_ptr, + sm_ptr, data, dof_info, values_quad, @@ -3248,9 +3252,10 @@ namespace internal using Number2_ = const Number2; Processor( - const unsigned int n_components, - const bool integrate, - const Number2 * global_vector_ptr, + const unsigned int n_components, + const bool integrate, + const Number2 * global_vector_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -3270,6 +3275,7 @@ namespace internal : n_components(n_components) , integrate(integrate) , global_vector_ptr(global_vector_ptr) + , sm_ptr(sm_ptr) , data(data) , dof_info(dof_info) , values_quad(values_quad) @@ -3402,9 +3408,10 @@ namespace internal subface_index); } - const unsigned int n_components; - const bool integrate; - const Number2 * global_vector_ptr; + const unsigned int n_components; + const bool integrate; + const Number2 * global_vector_ptr; + const std::vector> *sm_ptr; const MatrixFreeFunctions::ShapeInfo &data; const MatrixFreeFunctions::DoFInfo & dof_info; VectorizedArrayType * values_quad; @@ -3432,9 +3439,10 @@ namespace internal { template static bool - run(const unsigned int n_components, - const unsigned int n_face_orientations, - Number2 * dst_ptr, + run(const unsigned int n_components, + const unsigned int n_face_orientations, + Number2 * dst_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_array, @@ -3453,6 +3461,8 @@ namespace internal face_orientations, const Table<2, unsigned int> &orientation_map) { + (void)sm_ptr; + if (dst_ptr == nullptr) { AssertDimension(n_face_orientations, 1); @@ -3482,6 +3492,7 @@ namespace internal n_components, true, dst_ptr, + sm_ptr, data, dof_info, values_quad, @@ -3517,10 +3528,11 @@ namespace internal Processor( - VectorizedArrayType *values_array, - const unsigned int n_components, - const bool integrate, - Number2 * global_vector_ptr, + VectorizedArrayType * values_array, + const unsigned int n_components, + const bool integrate, + Number2 * global_vector_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -3541,6 +3553,7 @@ namespace internal , n_components(n_components) , integrate(integrate) , global_vector_ptr(global_vector_ptr) + , sm_ptr(sm_ptr) , data(data) , dof_info(dof_info) , values_quad(values_quad) @@ -3695,9 +3708,10 @@ namespace internal VectorizedArrayType *values_array; - const unsigned int n_components; - const bool integrate; - Number2 * global_vector_ptr; + const unsigned int n_components; + const bool integrate; + Number2 * global_vector_ptr; + const std::vector> *sm_ptr; const MatrixFreeFunctions::ShapeInfo &data; const MatrixFreeFunctions::DoFInfo & dof_info; VectorizedArrayType * values_quad; diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index dbb4e57a42..3052af2dff 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -99,9 +99,10 @@ namespace internal static bool gather_evaluate( - const unsigned int n_components, - const std::size_t n_face_orientations, - const Number * src_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + const Number * src_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -121,9 +122,10 @@ namespace internal static bool integrate_scatter( - const unsigned int n_components, - const std::size_t n_face_orientations, - Number * dst_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + Number * dst_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_array, diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index 168cf1c82c..1a5c254899 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -200,9 +200,10 @@ namespace internal template bool FEFaceEvaluationFactory::gather_evaluate( - const unsigned int n_components, - const std::size_t n_face_orientations, - const Number * src_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + const Number * src_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -230,6 +231,7 @@ namespace internal n_components, n_face_orientations, src_ptr, + sm_ptr, data, dof_info, values_quad, @@ -252,9 +254,10 @@ namespace internal template bool FEFaceEvaluationFactory::integrate_scatter( - const unsigned int n_components, - const std::size_t n_face_orientations, - Number * dst_ptr, + const unsigned int n_components, + const std::size_t n_face_orientations, + Number * dst_ptr, + const std::vector> *sm_ptr, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_array, @@ -283,6 +286,7 @@ namespace internal n_components, n_face_orientations, dst_ptr, + sm_ptr, data, dof_info, values_array, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index a1be617d7b..9d65dc211f 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -1170,10 +1170,14 @@ protected: */ template void - read_write_operation(const VectorOperation & operation, - const std::array &vectors, - const std::bitset &mask, - const bool apply_constraints = true) const; + read_write_operation( + const VectorOperation & operation, + const std::array &vectors, + const std::array< + const std::vector> *, + n_components_> & vectors_sm, + const std::bitset &mask, + const bool apply_constraints = true) const; /** * A unified function to read from and write into vectors based on the given @@ -1185,8 +1189,11 @@ protected: template void read_write_operation_contiguous( - const VectorOperation & operation, - const std::array & vectors, + const VectorOperation & operation, + const std::array &vectors, + const std::array< + const std::vector> *, + n_components_> & vectors_sm, const std::bitset &mask) const; /** @@ -4181,10 +4188,14 @@ template inline void FEEvaluationBase:: - read_write_operation(const VectorOperation & operation, - const std::array &src, - const std::bitset &mask, - const bool apply_constraints) const + read_write_operation( + const VectorOperation & operation, + const std::array &src, + const std::array< + const std::vector> *, + n_components_> & src_sm, + const std::bitset &mask, + const bool apply_constraints) const { // Case 1: No MatrixFree object given, simple case because we do not need to // process constraints and need not care about vectorization -> go to @@ -4226,7 +4237,7 @@ FEEvaluationBase:: [this->cell] >= internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) { - read_write_operation_contiguous(operation, src, mask); + read_write_operation_contiguous(operation, src, src_sm, mask); return; } @@ -4579,10 +4590,15 @@ template inline void FEEvaluationBase:: read_write_operation_contiguous( - const VectorOperation & operation, - const std::array & src, + const VectorOperation & operation, + const std::array &src, + const std::array< + const std::vector> *, + n_components_> & vectors_sm, const std::bitset &mask) const { + (void)vectors_sm; // TODO: use it + // This functions processes the functions read_dof_values, // distribute_local_to_global, and set_dof_values with the same code for // contiguous cell indices (DG case). The distinction between these three @@ -4837,30 +4853,67 @@ namespace internal return nullptr; } - template - std::array::type, - IsBlockVector::type>:: - value>::BaseVectorType *, - n_components> - get_vector_data(VectorType &src, const unsigned int first_index) + template ::value, + VectorType>::type * = nullptr> + const std::vector> * + get_shared_vector_data(VectorType &vec /*other parameters?*/) { - // select between block vectors and non-block vectors. Note that the number - // of components is checked in the internal data + return vec.shared_vector_data(); + } + + template ::value, + VectorType>::type * = nullptr> + const std::vector> * + get_shared_vector_data(VectorType &) + { + return nullptr; + } + + template + std::pair< std::array::type, IsBlockVector::type>:: value>::BaseVectorType *, - n_components> + n_components>, + std::array< + const std::vector::type, + IsBlockVector::type>::value>:: + BaseVectorType::value_type>> *, + n_components>> + get_vector_data(VectorType &src, const unsigned int first_index) + { + // select between block vectors and non-block vectors. Note that the number + // of components is checked in the internal data + std::pair< + std::array::type, + IsBlockVector::type>:: + value>::BaseVectorType *, + n_components>, + std::array< + const std::vector< + ArrayView::type, + IsBlockVector::type>:: + value>::BaseVectorType::value_type>> *, + n_components>> src_data; + for (unsigned int d = 0; d < n_components; ++d) - src_data[d] = internal::BlockVectorSelector< + src_data.first[d] = internal::BlockVectorSelector< typename std::remove_const::type, IsBlockVector::type>::value>:: get_vector_component( const_cast::type &>(src), d + first_index); + for (unsigned int d = 0; d < n_components; ++d) + src_data.second[d] = get_shared_vector_data(*src_data.first[d]); + return src_data; } } // namespace internal @@ -4882,7 +4935,8 @@ FEEvaluationBase:: internal::VectorReader reader; read_write_operation(reader, - src_data, + src_data.first, + src_data.second, std::bitset().flip(), true); @@ -4908,7 +4962,8 @@ FEEvaluationBase:: internal::VectorReader reader; read_write_operation(reader, - src_data, + src_data.first, + src_data.second, std::bitset().flip(), false); @@ -4942,7 +4997,7 @@ FEEvaluationBase:: internal::VectorDistributorLocalToGlobal distributor; - read_write_operation(distributor, dst_data, mask); + read_write_operation(distributor, dst_data.first, dst_data.second, mask); } @@ -4968,7 +5023,7 @@ FEEvaluationBase:: internal::get_vector_data(dst, first_index); internal::VectorSetter setter; - read_write_operation(setter, dst_data, mask); + read_write_operation(setter, dst_data.first, dst_data.second, mask); } @@ -4995,7 +5050,7 @@ FEEvaluationBase:: internal::get_vector_data(dst, first_index); internal::VectorSetter setter; - read_write_operation(setter, dst_data, mask, false); + read_write_operation(setter, dst_data.first, dst_data.second, mask, false); } @@ -8472,6 +8527,10 @@ FEFaceEvaluationdof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && this->is_interior_face == false) @@ -8488,6 +8547,7 @@ FEFaceEvaluation(input_vector), + shared_vector_data, *this->data, *this->dof_info, this->begin_values(), @@ -8528,6 +8588,7 @@ FEFaceEvaluation(input_vector), + shared_vector_data, *this->data, *this->dof_info, this->begin_values(), @@ -8551,6 +8612,7 @@ FEFaceEvaluation(input_vector), + shared_vector_data, *this->data, *this->dof_info, this->begin_values(), @@ -8637,6 +8699,9 @@ FEFaceEvaluationis_interior_face == false) == false, ExcNotImplemented()); + // TODO + const auto shared_vector_data = internal::get_shared_vector_data(destination); + // TODO: this copying should not be necessary once we have introduced // an internal-data structure std::array cells_ = {}; @@ -8657,6 +8722,7 @@ FEFaceEvaluation(destination), + shared_vector_data, *this->data, *this->dof_info, this->begin_dof_values(), @@ -8688,6 +8754,7 @@ FEFaceEvaluation(destination), + shared_vector_data, *this->data, *this->dof_info, this->begin_dof_values(), diff --git a/include/deal.II/matrix_free/type_traits.h b/include/deal.II/matrix_free/type_traits.h index cbf49bd2bf..b57c05049a 100644 --- a/include/deal.II/matrix_free/type_traits.h +++ b/include/deal.II/matrix_free/type_traits.h @@ -166,6 +166,29 @@ namespace internal const bool has_begin::value; + // same as above to check + // ... T::shared_vector_data() const + template + struct has_shared_vector_data + { + private: + static void + detect(...); + + template + static decltype(std::declval().shared_vector_data()) + detect(const U &); + + public: + static const bool value = + !std::is_same()))>::value; + }; + + // We need to have a separate declaration for static const members + template + const bool has_shared_vector_data::value; + + // type trait for vector T and Number to see if // we can do vectorized load/save. // for VectorReader and VectorDistributorLocalToGlobal we assume that -- 2.39.5