From: Peter Munch Date: Tue, 14 Dec 2021 22:30:04 +0000 (+0100) Subject: Global coarsening: compress weights X-Git-Tag: v9.4.0-rc1~720^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F13099%2Fhead;p=dealii.git Global coarsening: compress weights --- diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index ae7c1c81e3..5948e7c355 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -2586,6 +2586,40 @@ namespace internal } + template + inline void + weight_fe_q_dofs_by_entity(const VectorizedArray *weights, + const unsigned int n_components, + const int loop_length_non_template, + VectorizedArray *data) + { + const int loop_length = loop_length_template != -1 ? + loop_length_template : + loop_length_non_template; + + Assert(loop_length > 0, ExcNotImplemented()); + Assert(loop_length < 100, ExcNotImplemented()); + unsigned int degree_to_3[100]; + degree_to_3[0] = 0; + for (int i = 1; i < loop_length - 1; ++i) + degree_to_3[i] = 1; + degree_to_3[loop_length - 1] = 2; + for (unsigned int c = 0; c < n_components; ++c) + for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k) + for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j) + { + const unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j]; + data[0] *= weights[shift]; + // loop bound as int avoids compiler warnings in case loop_length + // == 1 (polynomial degree 0) + for (int i = 1; i < loop_length - 1; ++i) + data[i] *= weights[shift + 1]; + data[loop_length - 1] *= weights[shift + 2]; + data += loop_length; + } + } + + } // end of namespace internal diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h index e1a42451cb..35767cde1a 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h @@ -419,6 +419,13 @@ private: */ std::vector weights; + /** + * Weights for continuous elements, compressed into 3^dim doubles per + * cell if possible. + */ + std::vector, Utilities::pow(3, dim)>> + weights_compressed; + /** * DoF indices of the coarse cells, expressed in indices local to the MPI * rank. 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 ed233633cd..b3dfc48ecc 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -40,6 +40,7 @@ #include #include +#include #include #include @@ -1079,6 +1080,121 @@ namespace internal + /** + * Try to compress weights to Utilities::pow(3, dim) doubles. + * Weights of a cell can be compressed if all components have the + * same weights. + */ + template + static void + compress_weights( + MGTwoLevelTransfer> + &transfer) + { + const Number *weights = transfer.weights.data(); + + // Helper function that extracts the weight. Before doing so, + // it checks the values already set. + const auto set = [](auto &mask, const auto &weight) { + if (mask == -1.0 || mask == weight) + { + mask = weight; + return true; + } + + return false; + }; + + std::vector> masks; + + // loop over all schemes + for (const auto &scheme : transfer.schemes) + { + const int loop_length = scheme.degree_fine + 1; + unsigned int degree_to_3[100]; + degree_to_3[0] = 0; + for (int i = 1; i < loop_length - 1; ++i) + degree_to_3[i] = 1; + degree_to_3[loop_length - 1] = 2; + + // loop over all cells + for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell) + { + std::vector> mask( + transfer.n_components); + + // loop over all components + for (unsigned int c = 0; c < transfer.n_components; ++c) + mask[c].fill(-1.0); + + for (unsigned int c = 0; c < transfer.n_components; ++c) + for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k) + for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j) + { + const unsigned int shift = + 9 * degree_to_3[k] + 3 * degree_to_3[j]; + + if (!set(mask[c][shift], weights[0])) + return; // early return, since the entry already + // had a different weight + + for (int i = 1; i < loop_length - 1; ++i) + if (!set(mask[c][shift + 1], weights[i])) + return; + + if (!set(mask[c][shift + 2], weights[loop_length - 1])) + return; + + weights += loop_length; + } + + if (std::find_if(mask.begin(), mask.end(), [&](const auto &a) { + return a != mask[0]; + }) != mask.end()) + return; + + masks.push_back(mask[0]); + } + } + + // vectorize weights + std::vector, Utilities::pow(3, dim)>> + masks_vectorized; + + const unsigned int n_lanes = VectorizedArray::size(); + + unsigned int c = 0; + + for (const auto &scheme : transfer.schemes) + { + for (unsigned int cell = 0; cell < scheme.n_coarse_cells; + cell += n_lanes) + { + const unsigned int n_lanes_filled = + (cell + n_lanes > scheme.n_coarse_cells) ? + (scheme.n_coarse_cells - cell) : + n_lanes; + + std::array, Utilities::pow(3, dim)> mask; + + for (unsigned int v = 0; v < n_lanes_filled; ++v, ++c) + { + for (unsigned int i = 0; + i < Utilities::pow(3, dim); + ++i) + mask[i][v] = masks[c][i]; + } + + masks_vectorized.push_back(mask); + } + } + + // copy result + transfer.weights_compressed = masks_vectorized; + } + + + public: template static void @@ -1592,6 +1708,9 @@ namespace internal weights_1 += transfer.schemes[1].n_dofs_per_cell_fine; } }); + + if (is_feq) + compress_weights(transfer); } } @@ -1669,6 +1788,16 @@ namespace internal ExcNotImplemented()); #endif + const bool is_feq = + std::all_of(dof_handler_fine.get_fe_collection().begin(), + dof_handler_fine.get_fe_collection().end(), + [](const auto &fe) { + return fe.n_base_elements() == 1 && + (dynamic_cast *>( + &fe.base_element(0)) != nullptr); + ; + }); + const auto process_cells = [&](const auto &fu) { loop_over_active_or_level_cells( dof_handler_coarse, mg_level_coarse, [&](const auto &cell_coarse) { @@ -2096,6 +2225,9 @@ namespace internal weights_[fe_pair_no] += transfer.schemes[fe_pair_no].n_dofs_per_cell_fine; }); + + if (is_feq) + compress_weights(transfer); } } }; @@ -2360,9 +2492,14 @@ MGTwoLevelTransfer>:: const unsigned int *indices_coarse = level_dof_indices_coarse.data(); const unsigned int *indices_fine = level_dof_indices_fine.data(); const Number * weights = nullptr; + const std::array, Utilities::pow(3, dim)> + *weights_compressed = nullptr; if (fine_element_is_continuous) - weights = this->weights.data(); + { + weights = this->weights.data(); + weights_compressed = this->weights_compressed.data(); + } for (const auto &scheme : schemes) { @@ -2370,29 +2507,44 @@ MGTwoLevelTransfer>:: if (scheme.prolongation_matrix.size() == 0 && scheme.prolongation_matrix_1d.size() == 0) { - for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell) + for (unsigned int cell = 0; cell < scheme.n_coarse_cells; + cell += n_lanes) { - if ((scheme.n_dofs_per_cell_fine != 0) && - (scheme.n_dofs_per_cell_coarse != 0)) + const unsigned int n_lanes_filled = + (cell + n_lanes > scheme.n_coarse_cells) ? + (scheme.n_coarse_cells - cell) : + n_lanes; + + // read from source vector + for (unsigned int v = 0; v < n_lanes_filled; ++v) { + if ((scheme.n_dofs_per_cell_fine != 0) && + (scheme.n_dofs_per_cell_coarse != 0)) + { + if (fine_element_is_continuous) + for (unsigned int i = 0; + i < scheme.n_dofs_per_cell_fine; + ++i) + vec_fine_ptr->local_element(indices_fine[i]) += + read_dof_values(indices_coarse[i], vec_coarse) * + weights[i]; + else + for (unsigned int i = 0; + i < scheme.n_dofs_per_cell_fine; + ++i) + vec_fine_ptr->local_element(indices_fine[i]) += + read_dof_values(indices_coarse[i], vec_coarse); + } + + indices_fine += scheme.n_dofs_per_cell_fine; + indices_coarse += scheme.n_dofs_per_cell_coarse; + if (fine_element_is_continuous) - for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; - ++i) - vec_fine_ptr->local_element(indices_fine[i]) += - read_dof_values(indices_coarse[i], vec_coarse) * - weights[i]; - else - for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; - ++i) - vec_fine_ptr->local_element(indices_fine[i]) += - read_dof_values(indices_coarse[i], vec_coarse); + weights += scheme.n_dofs_per_cell_fine; } - indices_fine += scheme.n_dofs_per_cell_fine; - indices_coarse += scheme.n_dofs_per_cell_coarse; - if (fine_element_is_continuous) - weights += scheme.n_dofs_per_cell_fine; + weights_compressed += 1; } continue; @@ -2444,9 +2596,20 @@ MGTwoLevelTransfer>:: // ------------------------------ fine ------------------------------- // weight and write into dst vector + if (fine_element_is_continuous && this->weights_compressed.size() > 0) + { + internal::weight_fe_q_dofs_by_entity( + weights_compressed->data(), + n_components, + scheme.degree_fine + 1, + evaluation_data_fine.begin()); + weights_compressed += 1; + } + for (unsigned int v = 0; v < n_lanes_filled; ++v) { - if (fine_element_is_continuous) + if (fine_element_is_continuous && + this->weights_compressed.size() == 0) for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i) vec_fine_ptr->local_element(indices_fine[i]) += evaluation_data_fine[i][v] * weights[i]; @@ -2519,9 +2682,14 @@ MGTwoLevelTransfer>:: const unsigned int *indices_coarse = level_dof_indices_coarse.data(); const unsigned int *indices_fine = level_dof_indices_fine.data(); const Number * weights = nullptr; + const std::array, Utilities::pow(3, dim)> + *weights_compressed = nullptr; if (fine_element_is_continuous) - weights = this->weights.data(); + { + weights = this->weights.data(); + weights_compressed = this->weights_compressed.data(); + } for (const auto &scheme : schemes) { @@ -2529,33 +2697,47 @@ MGTwoLevelTransfer>:: if (scheme.prolongation_matrix.size() == 0 && scheme.prolongation_matrix_1d.size() == 0) { - for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell) + for (unsigned int cell = 0; cell < scheme.n_coarse_cells; + cell += n_lanes) { - if ((scheme.n_dofs_per_cell_fine != 0) && - (scheme.n_dofs_per_cell_coarse != 0)) + const unsigned int n_lanes_filled = + (cell + n_lanes > scheme.n_coarse_cells) ? + (scheme.n_coarse_cells - cell) : + n_lanes; + + for (unsigned int v = 0; v < n_lanes_filled; ++v) { + if ((scheme.n_dofs_per_cell_fine != 0) && + (scheme.n_dofs_per_cell_coarse != 0)) + { + if (fine_element_is_continuous) + for (unsigned int i = 0; + i < scheme.n_dofs_per_cell_fine; + ++i) + distribute_local_to_global( + indices_coarse[i], + vec_fine_ptr->local_element(indices_fine[i]) * + weights[i], + this->vec_coarse); + else + for (unsigned int i = 0; + i < scheme.n_dofs_per_cell_fine; + ++i) + distribute_local_to_global( + indices_coarse[i], + vec_fine_ptr->local_element(indices_fine[i]), + this->vec_coarse); + } + + indices_fine += scheme.n_dofs_per_cell_fine; + indices_coarse += scheme.n_dofs_per_cell_coarse; + if (fine_element_is_continuous) - for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; - ++i) - distribute_local_to_global(indices_coarse[i], - vec_fine_ptr->local_element( - indices_fine[i]) * - weights[i], - this->vec_coarse); - else - for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; - ++i) - distribute_local_to_global(indices_coarse[i], - vec_fine_ptr->local_element( - indices_fine[i]), - this->vec_coarse); + weights += scheme.n_dofs_per_cell_fine; } - indices_fine += scheme.n_dofs_per_cell_fine; - indices_coarse += scheme.n_dofs_per_cell_coarse; - if (fine_element_is_continuous) - weights += scheme.n_dofs_per_cell_fine; + weights_compressed += 1; } continue; @@ -2583,7 +2765,8 @@ MGTwoLevelTransfer>:: // read from source vector and weight for (unsigned int v = 0; v < n_lanes_filled; ++v) { - if (fine_element_is_continuous) + if (fine_element_is_continuous && + this->weights_compressed.size() == 0) for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i) evaluation_data_fine[i][v] = vec_fine_ptr->local_element(indices_fine[i]) * weights[i]; @@ -2598,6 +2781,16 @@ MGTwoLevelTransfer>:: weights += scheme.n_dofs_per_cell_fine; } + if (fine_element_is_continuous && this->weights_compressed.size() > 0) + { + internal::weight_fe_q_dofs_by_entity( + weights_compressed->data(), + n_components, + scheme.degree_fine + 1, + evaluation_data_fine.begin()); + weights_compressed += 1; + } + // ------------------------------ fine ------------------------------- for (int c = n_components - 1; c >= 0; --c) { diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc index 9762e2eea7..8010b0fa6f 100644 --- a/source/multigrid/mg_transfer_matrix_free.cc +++ b/source/multigrid/mg_transfer_matrix_free.cc @@ -29,6 +29,7 @@ #include #include +#include #include #include @@ -355,41 +356,6 @@ MGTransferMatrixFree::restrict_and_add( -namespace -{ - template - void - weight_dofs_on_child(const VectorizedArray *weights, - const unsigned int n_components, - const unsigned int fe_degree, - VectorizedArray * data) - { - Assert(fe_degree > 0, ExcNotImplemented()); - Assert(fe_degree < 100, ExcNotImplemented()); - const int loop_length = degree != -1 ? 2 * degree + 1 : 2 * fe_degree + 1; - unsigned int degree_to_3[100]; - degree_to_3[0] = 0; - for (int i = 1; i < loop_length - 1; ++i) - degree_to_3[i] = 1; - degree_to_3[loop_length - 1] = 2; - for (unsigned int c = 0; c < n_components; ++c) - for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k) - for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j) - { - const unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j]; - data[0] *= weights[shift]; - // loop bound as int avoids compiler warnings in case loop_length - // == 1 (polynomial degree 0) - for (int i = 1; i < loop_length - 1; ++i) - data[i] *= weights[shift + 1]; - data[loop_length - 1] *= weights[shift + 2]; - data += loop_length; - } - } -} // namespace - - - template template void @@ -494,10 +460,13 @@ MGTransferMatrixFree::do_prolongate_add( c * n_scalar_cell_dofs, fe_degree + 1, 2 * fe_degree + 1); - weight_dofs_on_child( + internal::weight_fe_q_dofs_by_entity( &weights_on_refined[to_level - 1][(cell / vec_size) * three_to_dim], n_components, - fe_degree, + 2 * fe_degree + 1, evaluation_data.begin()); } else @@ -575,12 +544,14 @@ MGTransferMatrixFree::do_restrict_add( // perform tensorized operation if (element_is_continuous) { - weight_dofs_on_child( - &weights_on_refined[from_level - 1] - [(cell / vec_size) * three_to_dim], - n_components, - fe_degree, - evaluation_data.data()); + internal::weight_fe_q_dofs_by_entity< + dim, + degree != -1 ? 2 * degree + 1 : -1, + Number>(&weights_on_refined[from_level - 1] + [(cell / vec_size) * three_to_dim], + n_components, + 2 * fe_degree + 1, + evaluation_data.data()); for (unsigned int c = 0; c < n_components; ++c) internal::FEEvaluationImplBasisChange< internal::evaluate_general,