From 5a07d338fea3d08cbecd4d24c29f0120d4900796 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 18 Jan 2022 08:26:14 +0100 Subject: [PATCH] FEEval: replace usages of n_vectorization_actual --- include/deal.II/matrix_free/fe_evaluation.h | 129 ++++++++++++-------- 1 file changed, 79 insertions(+), 50 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 8f6703f041..f8964cbbb3 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3441,17 +3441,20 @@ FEEvaluationBase:: // Case 2: contiguous indices which use reduced storage of indices and can // use vectorized load/store operations -> go to separate function - AssertIndexRange( - this->cell, - this->dof_info->index_storage_variants[this->dof_access_index].size()); - if (this->dof_info->index_storage_variants - [is_face ? this->dof_access_index : - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell] >= - internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) + if (this->cell != numbers::invalid_unsigned_int) { - read_write_operation_contiguous(operation, src, src_sm, mask); - return; + AssertIndexRange( + this->cell, + this->dof_info->index_storage_variants[this->dof_access_index].size()); + if (this->dof_info->index_storage_variants + [is_face ? this->dof_access_index : + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell] >= internal::MatrixFreeFunctions::DoFInfo:: + IndexStorageVariants::contiguous) + { + read_write_operation_contiguous(operation, src, src_sm, mask); + return; + } } // Case 3: standard operation with one index per degree of freedom -> go on @@ -3461,10 +3464,6 @@ FEEvaluationBase:: ExcNotImplemented("Masking currently not implemented for " "non-contiguous DoF storage")); - unsigned int n_vectorization_actual = - this->dof_info - ->n_vectorization_lanes_filled[this->dof_access_index][this->cell]; - const std::array &cells = this->get_cell_ids(); @@ -3472,8 +3471,9 @@ FEEvaluationBase:: if (is_face == false) { - for (unsigned int v = 0; v < n_vectorization_actual; ++v) - if (this->dof_info->hanging_node_constraint_masks.size() > 0 && + for (unsigned int v = 0; v < n_lanes; ++v) + if (cells[v] != numbers::invalid_unsigned_int && + this->dof_info->hanging_node_constraint_masks.size() > 0 && this->dof_info->hanging_node_constraint_masks [cells[v] * this->n_fe_components + this->first_selected_component] != @@ -3538,8 +3538,11 @@ FEEvaluationBase:: if (is_face) { - for (unsigned int v = 0; v < n_vectorization_actual; ++v) + for (unsigned int v = 0; v < n_lanes; ++v) { + if (cells[v] == numbers::invalid_unsigned_int) + continue; + Assert(cells[v] < this->dof_info->row_starts.size() - 1, ExcInternalError()); const std::pair *my_index_start = @@ -3559,8 +3562,11 @@ FEEvaluationBase:: } else { - for (unsigned int v = 0; v < n_vectorization_actual; ++v) + for (unsigned int v = 0; v < n_lanes; ++v) { + if (cells[v] == numbers::invalid_unsigned_int) + continue; + const std::pair *my_index_start = &this->dof_info->row_starts[cells[v] * this->n_fe_components + this->first_selected_component]; @@ -3586,32 +3592,45 @@ FEEvaluationBase:: } } + if (std::count_if(cells.begin(), cells.end(), [](const auto i) { + return i != numbers::invalid_unsigned_int; + }) < n_lanes) + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_empty(values_dofs[comp][i]); + // Case where we have no constraints throughout the whole cell: Can go // through the list of DoFs directly if (!has_constraints && apply_constraints) { - if (n_vectorization_actual < n_lanes) - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_component; ++i) - operation.process_empty(values_dofs[comp][i]); if (n_components == 1 || this->n_fe_components == 1) { - for (unsigned int v = 0; v < n_vectorization_actual; ++v) - for (unsigned int i = 0; i < dofs_per_component; ++i) - for (unsigned int comp = 0; comp < n_components; ++comp) - operation.process_dof(dof_indices[v][i], - *src[comp], - values_dofs[comp][i][v]); + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (cells[v] == numbers::invalid_unsigned_int) + continue; + + for (unsigned int i = 0; i < dofs_per_component; ++i) + for (unsigned int comp = 0; comp < n_components; ++comp) + operation.process_dof(dof_indices[v][i], + *src[comp], + values_dofs[comp][i][v]); + } } else { for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int v = 0; v < n_vectorization_actual; ++v) - for (unsigned int i = 0; i < dofs_per_component; ++i) - operation.process_dof( - dof_indices[v][comp * dofs_per_component + i], - *src[0], - values_dofs[comp][i][v]); + for (unsigned int v = 0; v < n_lanes; ++v) + { + if (cells[v] == numbers::invalid_unsigned_int) + continue; + + for (unsigned int i = 0; i < dofs_per_component; ++i) + operation.process_dof( + dof_indices[v][comp * dofs_per_component + i], + *src[0], + values_dofs[comp][i][v]); + } } return; } @@ -3621,12 +3640,12 @@ FEEvaluationBase:: // holds local number on cell, index iterates over the elements of // index_local_to_global and dof_indices points to the global indices stored // in index_local_to_global - if (n_vectorization_actual < n_lanes) - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < dofs_per_component; ++i) - operation.process_empty(values_dofs[comp][i]); - for (unsigned int v = 0; v < n_vectorization_actual; ++v) + + for (unsigned int v = 0; v < n_lanes; ++v) { + if (cells[v] == numbers::invalid_unsigned_int) + continue; + const unsigned int cell_index = cells[v]; const unsigned int cell_dof_index = cell_index * this->n_fe_components + this->first_selected_component; @@ -3850,6 +3869,8 @@ FEEvaluationBase:: values_dofs[c] = const_cast(this->values_dofs) + c * dofs_per_component; + Assert(this->cell != numbers::invalid_unsigned_int, ExcNotImplemented()); + // Simple case: We have contiguous storage, so we can simply copy out the // data if ((this->dof_info->index_storage_variants[ind][this->cell] == @@ -4282,10 +4303,6 @@ FEEvaluationBase:: this->dof_info->hanging_node_constraint_masks.size() == 0) return; // nothing to do with faces - const unsigned int n_vectorization_actual = - this->dof_info - ->n_vectorization_lanes_filled[this->dof_access_index][this->cell]; - constexpr unsigned int n_lanes = VectorizedArrayType::size(); std::array constraint_mask; @@ -4295,8 +4312,15 @@ FEEvaluationBase:: const std::array &cells = this->get_cell_ids(); - for (unsigned int v = 0; v < n_vectorization_actual; ++v) + for (unsigned int v = 0; v < n_lanes; ++v) { + if (cells[v] == numbers::invalid_unsigned_int) + { + constraint_mask[v] = + internal::MatrixFreeFunctions::ConstraintKinds::unconstrained; + continue; + } + const unsigned int cell_index = cells[v]; const unsigned int cell_dof_index = cell_index * this->n_fe_components + this->first_selected_component; @@ -4312,10 +4336,6 @@ FEEvaluationBase:: if (hn_available == false) return; // no hanging node on cell batch -> nothing to do - for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v) - constraint_mask[v] = - internal::MatrixFreeFunctions::ConstraintKinds::unconstrained; - internal::FEEvaluationHangingNodesFactory:: apply(n_components, this->data->data.front().fe_degree, @@ -6874,8 +6894,12 @@ FEEvaluationjacobian = &this->mapping_data->jacobians[0][offsets]; this->J_value = &this->mapping_data->JxW_values[offsets]; - for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i) + unsigned int i = 0; + for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell); + ++i) this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i; + for (; i < VectorizedArrayType::size(); ++i) + this->cell_ids[i] = numbers::invalid_unsigned_int; # ifdef DEBUG this->dof_values_initialized = false; @@ -7596,8 +7620,13 @@ FEFaceEvaluationreinit_face(this->matrix_free->get_face_info(face_index)); - for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i) + + unsigned int i = 0; + for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell); + ++i) this->cell_or_face_ids[i] = face_index * VectorizedArrayType::size() + i; + for (; i < VectorizedArrayType::size(); ++i) + this->cell_or_face_ids[i] = numbers::invalid_unsigned_int; this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index]; const unsigned int offsets = -- 2.39.5