From a86c9ca873bef90edad6909186d8fd935414c5dd Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 17 Apr 2024 18:48:32 +0200 Subject: [PATCH] Non-nested MG transfer: Avoid repeated access of arrays --- .../mg_transfer_global_coarsening.templates.h | 101 ++++++++++-------- 1 file changed, 55 insertions(+), 46 deletions(-) diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index a0e12f38a1..f58c6d5c14 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -4900,9 +4900,6 @@ MGTwoLevelTransferNonNested::prolongate_and_add_internal_comp( EvaluatorTypeTraits; using value_type = typename Traits::value_type; - std::vector evaluation_point_results; - std::vector buffer; - const auto evaluation_function = [&](auto &values, const auto &cell_data) { this->signals_non_nested.prolongation_cell_loop(true); std::vector solution_values; @@ -4938,38 +4935,48 @@ MGTwoLevelTransferNonNested::prolongate_and_add_internal_comp( }; this->signals_non_nested.prolongation(true); + + std::vector evaluation_point_results; + std::vector buffer; rpe.template evaluate_and_process(evaluation_point_results, buffer, evaluation_function); this->signals_non_nested.prolongation(false); - // Weight operator in case some points are owned by multiple cells. - if (rpe.is_map_unique() == false) - { - const auto evaluation_point_results_temp = evaluation_point_results; - evaluation_point_results.clear(); - evaluation_point_results.reserve(rpe.get_point_ptrs().size() - 1); + // Keep a vector of typical inverse touch counts that avoid divisions, all + // other cases are handled by a division in the code below. + std::array typical_weights; + for (unsigned int i = 0; i < typical_weights.size(); ++i) + typical_weights[i] = Number(1) / Number(i + 1); - const auto &ptr = rpe.get_point_ptrs(); - - for (unsigned int i = 0; i < ptr.size() - 1; ++i) + const bool must_interpolate = (rpe.is_map_unique() == false); + const auto &ptr = rpe.get_point_ptrs(); + for (unsigned int j = 0; j < ptr.size() - 1; ++j) + { + value_type result; + // Weight operator in case some points are owned by multiple cells. + if (must_interpolate) { - const auto n_entries = ptr[i + 1] - ptr[i]; - - value_type result{}; + const unsigned int n_entries = ptr[j + 1] - ptr[j]; - if (n_entries > 0) + if (n_entries > 1) { - for (unsigned int j = 0; j < n_entries; ++j) - result += evaluation_point_results_temp[ptr[i] + j]; - result /= Number(n_entries); + result = {}; + for (unsigned int k = 0; k < n_entries; ++k) + result += evaluation_point_results[ptr[j] + k]; + if (n_entries <= typical_weights.size()) + result *= typical_weights[n_entries - 1]; + else + result /= Number(n_entries); } - evaluation_point_results.push_back(result); + else if (n_entries == 1) + result = evaluation_point_results[ptr[j]]; + else + result = {}; } - } + else + result = evaluation_point_results[j]; - for (unsigned int j = 0; j < evaluation_point_results.size(); ++j) - { if (level_dof_indices_fine_ptrs.empty()) { for (unsigned int c = 0; c < n_components; ++c) @@ -4978,7 +4985,7 @@ MGTwoLevelTransferNonNested::prolongate_and_add_internal_comp( this->level_dof_indices_fine.size()); dst.local_element( this->level_dof_indices_fine[n_components * j + c]) += - internal::access(evaluation_point_results[j], c); + internal::access(result, c); } } else @@ -4992,7 +4999,7 @@ MGTwoLevelTransferNonNested::prolongate_and_add_internal_comp( this->level_dof_indices_fine.size()); dst.local_element( this->level_dof_indices_fine[n_components * i + c]) += - internal::access(evaluation_point_results[j], c); + internal::access(result, c); } } } @@ -5030,10 +5037,17 @@ MGTwoLevelTransferNonNested::restrict_and_add_internal_comp( std::vector evaluation_point_results; std::vector buffer; - evaluation_point_results.resize(rpe.get_point_ptrs().size() - 1); + std::array typical_weights; + for (unsigned int i = 0; i < typical_weights.size(); ++i) + typical_weights[i] = Number(1) / Number(i + 1); + + const bool must_interpolate = (rpe.is_map_unique() == false); + const auto &ptr = rpe.get_point_ptrs(); + evaluation_point_results.resize(ptr.size() - 1); for (unsigned int j = 0; j < evaluation_point_results.size(); ++j) { + value_type result; if (level_dof_indices_fine_ptrs.empty()) { for (unsigned int c = 0; c < n_components; ++c) @@ -5041,14 +5055,13 @@ MGTwoLevelTransferNonNested::restrict_and_add_internal_comp( AssertIndexRange(n_components * j + c, this->level_dof_indices_fine.size()); - internal::access(evaluation_point_results[j], c) = - src.local_element( - this->level_dof_indices_fine[n_components * j + c]); + internal::access(result, c) = src.local_element( + this->level_dof_indices_fine[n_components * j + c]); } } else { - evaluation_point_results[j] = value_type(); + result = value_type(); for (unsigned int i = this->level_dof_indices_fine_ptrs[j]; i < this->level_dof_indices_fine_ptrs[j + 1]; @@ -5057,26 +5070,22 @@ MGTwoLevelTransferNonNested::restrict_and_add_internal_comp( { AssertIndexRange(n_components * i + c, this->level_dof_indices_fine.size()); - internal::access(evaluation_point_results[j], c) += - src.local_element( - this->level_dof_indices_fine[n_components * i + c]); + internal::access(result, c) += src.local_element( + this->level_dof_indices_fine[n_components * i + c]); } } - } - - // Weight operator in case some points are owned by multiple cells. - if (rpe.is_map_unique() == false) - { - const auto &ptr = rpe.get_point_ptrs(); - - for (unsigned int i = 0; i < ptr.size() - 1; ++i) + if (must_interpolate) { - const auto n_entries = ptr[i + 1] - ptr[i]; - if (n_entries == 0) - continue; - - evaluation_point_results[i] /= Number(n_entries); + const unsigned int n_entries = ptr[j + 1] - ptr[j]; + if (n_entries > 1) + { + if (n_entries <= typical_weights.size()) + result *= typical_weights[n_entries - 1]; + else + result /= Number(n_entries); + } } + evaluation_point_results[j] = result; } const auto evaluation_function = [&](const auto &values, -- 2.39.5