From ba1f915e3a70288cc05dd86d72b38ae77cdcca04 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 25 Oct 2020 21:21:54 +0100 Subject: [PATCH] Different fixes --- .../deal.II/matrix_free/evaluation_kernels.h | 43 +++++- include/deal.II/matrix_free/fe_evaluation.h | 136 +++++++++++++----- 2 files changed, 136 insertions(+), 43 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index a45c1b2d22..67e437b511 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -2496,6 +2496,7 @@ namespace internal auto n_components = proc.n_components; auto integrate = proc.integrate; auto global_vector_ptr = proc.global_vector_ptr; + auto &sm_ptr = proc.sm_ptr; auto &data = proc.data; auto &dof_info = proc.dof_info; auto values_quad = proc.values_quad; @@ -2516,6 +2517,7 @@ namespace internal static const int fe_degree = Processor::fe_degree_; using VectorizedArrayType = typename Processor::VectorizedArrayType_; + using Number = typename Processor::Number_; using Number2_ = typename Processor::Number2_; const unsigned int cell = cells[0]; @@ -2999,7 +3001,8 @@ namespace internal const bool vectorization_possible = (n_face_orientations == 1) && - (n_filled_lanes == VectorizedArrayType::size()); + (n_filled_lanes == VectorizedArrayType::size()) && + (sm_ptr != nullptr); std::array vector_ptrs = {}; @@ -3009,7 +3012,20 @@ namespace internal if (n_face_orientations == 1) { for (unsigned int v = 0; v < n_filled_lanes; ++v) - vector_ptrs[v] = vector_ptr + indices[v]; + if (sm_ptr == nullptr) + { + vector_ptrs[v] = vector_ptr + indices[v]; + } + else + { + const auto &temp = + dof_info.dof_indices_contiguous_sm + [dof_access_index] + [cell * VectorizedArrayType::size() + v]; + vector_ptrs[v] = const_cast( + sm_ptr->operator[](temp.first).data() + + temp.second); + } } else if (n_face_orientations == VectorizedArrayType::size()) { @@ -3017,10 +3033,25 @@ namespace internal v < VectorizedArrayType::size(); ++v) if (cells[v] != numbers::invalid_unsigned_int) - vector_ptrs[v] = - vector_ptr + - dof_info.dof_indices_contiguous[dof_access_index] + { + if (sm_ptr == nullptr) + { + vector_ptrs[v] = + vector_ptr + + dof_info + .dof_indices_contiguous[dof_access_index] [cells[v]]; + } + else + { + const auto &temp = + dof_info.dof_indices_contiguous_sm + [dof_access_index][cells[v]]; + vector_ptrs[v] = const_cast( + sm_ptr->operator[](temp.first).data() + + temp.second); + } + } } else { @@ -3255,6 +3286,7 @@ namespace internal static const int fe_degree_ = fe_degree; static const int n_q_points_1d_ = n_q_points_1d; using VectorizedArrayType_ = VectorizedArrayType; + using Number_ = Number; using Number2_ = const Number2; Processor( @@ -3530,6 +3562,7 @@ namespace internal static const int fe_degree_ = fe_degree; static const int n_q_points_1d_ = n_q_points_1d; using VectorizedArrayType_ = VectorizedArrayType; + using Number_ = Number; using Number2_ = Number2; diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 759c64f7fe..b4790b7a05 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -4670,25 +4670,33 @@ FEEvaluationBase:: VectorizedArrayType::size()> vector_ptrs = {}; - (void)comp; - - /* for (unsigned int v = 0; v < n_filled_lanes; ++v) { + Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(), + ExcIndexRange( + ind, 0, this->dof_info->dof_indices_contiguous_sm.size())); + Assert(this->cell * VectorizedArrayType::size() + v < + this->dof_info->dof_indices_contiguous_sm[ind].size(), + ExcIndexRange( + this->cell * VectorizedArrayType::size() + v, + 0, + this->dof_info->dof_indices_contiguous_sm[ind].size())); + const auto &temp = this->dof_info->dof_indices_contiguous_sm [ind][this->cell * VectorizedArrayType::size() + v]; if (temp.first != numbers::invalid_unsigned_int) vector_ptrs[v] = const_cast( - vectors_sm[comp]->operator[](temp.first).data() + temp.second); + vectors_sm[comp]->operator[](temp.first).data() + temp.second + + this->dof_info->component_dof_indices_offset + [this->active_fe_index][this->first_selected_component]); else vector_ptrs[v] = nullptr; } for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v) vector_ptrs[v] = nullptr; - */ return vector_ptrs; }; @@ -4697,16 +4705,17 @@ FEEvaluationBase:: n_lanes == VectorizedArrayType::size()) { if (n_components == 1 || n_fe_components == 1) - for (unsigned int comp = 0; comp < n_components; ++comp) - { - auto vector_ptrs = compute_vector_ptrs(comp); - - operation.process_dofs_vectorized_transpose( - this->data->dofs_per_component_on_cell, - vector_ptrs, - values_dofs[comp], - vector_selector); - } + { + for (unsigned int comp = 0; comp < n_components; ++comp) + { + auto vector_ptrs = compute_vector_ptrs(comp); + operation.process_dofs_vectorized_transpose( + this->data->dofs_per_component_on_cell, + vector_ptrs, + values_dofs[comp], + vector_selector); + } + } else { auto vector_ptrs = compute_vector_ptrs(0); @@ -4720,7 +4729,8 @@ FEEvaluationBase:: else for (unsigned int comp = 0; comp < n_components; ++comp) { - auto vector_ptrs = compute_vector_ptrs(comp); + auto vector_ptrs = compute_vector_ptrs( + (n_components == 1 || n_fe_components == 1) ? comp : 0); for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell; ++i) @@ -4945,18 +4955,29 @@ namespace internal typename std::enable_if::value, VectorType>::type * = nullptr> const std::vector> * - get_shared_vector_data(VectorType &vec /*other parameters?*/) + get_shared_vector_data(VectorType & vec, + const bool is_valid_mode_for_sm, + const unsigned int active_fe_index, + const internal::MatrixFreeFunctions::DoFInfo *dof_info) { - // return &vec.shared_vector_data(); - (void)vec; - return nullptr; + // note: no hp is supported + if (is_valid_mode_for_sm && + dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/] + .size() > 0 && + active_fe_index == 0) + return &vec.shared_vector_data(); + else + return nullptr; } template ::value, VectorType>::type * = nullptr> const std::vector> * - get_shared_vector_data(VectorType &) + get_shared_vector_data(VectorType &, + const bool, + const unsigned int, + const internal::MatrixFreeFunctions::DoFInfo *) { return nullptr; } @@ -4974,7 +4995,11 @@ namespace internal IsBlockVector::type>::value>:: BaseVectorType::value_type>> *, n_components>> - get_vector_data(VectorType &src, const unsigned int first_index) + get_vector_data(VectorType & src, + const unsigned int first_index, + const bool is_valid_mode_for_sm, + const unsigned int active_fe_index, + const internal::MatrixFreeFunctions::DoFInfo *dof_info) { // select between block vectors and non-block vectors. Note that the number // of components is checked in the internal data @@ -5002,7 +5027,10 @@ namespace internal d + first_index); for (unsigned int d = 0; d < n_components; ++d) - src_data.second[d] = get_shared_vector_data(*src_data.first[d]); + src_data.second[d] = get_shared_vector_data(*src_data.first[d], + is_valid_mode_for_sm, + active_fe_index, + dof_info); return src_data; } @@ -5020,8 +5048,13 @@ inline void FEEvaluationBase:: read_dof_values(const VectorType &src, const unsigned int first_index) { - const auto src_data = - internal::get_vector_data(src, first_index); + const auto src_data = internal::get_vector_data( + src, + first_index, + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); internal::VectorReader reader; read_write_operation(reader, @@ -5047,8 +5080,13 @@ inline void FEEvaluationBase:: read_dof_values_plain(const VectorType &src, const unsigned int first_index) { - const auto src_data = - internal::get_vector_data(src, first_index); + const auto src_data = internal::get_vector_data( + src, + first_index, + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); internal::VectorReader reader; read_write_operation(reader, @@ -5082,8 +5120,13 @@ FEEvaluationBase:: internal::ExcAccessToUninitializedField()); # endif - const auto dst_data = - internal::get_vector_data(dst, first_index); + const auto dst_data = internal::get_vector_data( + dst, + first_index, + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); internal::VectorDistributorLocalToGlobal distributor; @@ -5109,8 +5152,13 @@ FEEvaluationBase:: internal::ExcAccessToUninitializedField()); # endif - const auto dst_data = - internal::get_vector_data(dst, first_index); + const auto dst_data = internal::get_vector_data( + dst, + first_index, + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); internal::VectorSetter setter; read_write_operation(setter, dst_data.first, dst_data.second, mask); @@ -5136,8 +5184,13 @@ FEEvaluationBase:: internal::ExcAccessToUninitializedField()); # endif - const auto dst_data = - internal::get_vector_data(dst, first_index); + const auto dst_data = internal::get_vector_data( + dst, + first_index, + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); internal::VectorSetter setter; read_write_operation(setter, dst_data.first, dst_data.second, mask, false); @@ -8617,9 +8670,12 @@ FEFaceEvaluationdof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); if (this->dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && @@ -8789,8 +8845,12 @@ FEFaceEvaluationis_interior_face == false) == false, ExcNotImplemented()); - // TODO - const auto shared_vector_data = internal::get_shared_vector_data(destination); + const auto shared_vector_data = internal::get_shared_vector_data( + destination, + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell, + this->active_fe_index, + this->dof_info); // TODO: this copying should not be necessary once we have introduced // an internal-data structure -- 2.39.5