From 84b69dca389e97c61801f5952f68a194ef6088a5 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 26 Apr 2018 13:35:01 +0200 Subject: [PATCH] Use an enum rather than a plain number to index into dof index arrays. --- include/deal.II/matrix_free/dof_info.h | 28 ++++++- .../deal.II/matrix_free/dof_info.templates.h | 75 ++++++++++--------- include/deal.II/matrix_free/fe_evaluation.h | 58 ++++++++------ .../matrix_free/matrix_free.templates.h | 14 ++-- 4 files changed, 105 insertions(+), 70 deletions(-) diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index 914b045e0a..55a36b3e39 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -246,6 +246,26 @@ namespace internal contiguous }; + /** + * Enum used to distinguish the data arrays for the vectorization type + * in cells and faces. + */ + enum DoFAccessIndex : unsigned char + { + /** + * The data index for the faces designated as interior + */ + dof_access_face_interior = 0, + /** + * The data index for the faces designated as exterior + */ + dof_access_face_exterior = 1, + /** + * The data index for the cells + */ + dof_access_cell = 2 + }; + /** * Stores the dimension of the underlying DoFHandler. Since the indices * are not templated, this is the variable that makes the dimension @@ -264,11 +284,11 @@ namespace internal unsigned int vectorization_length; /** - * Stores the index storage variant of all cell batches. + * Stores the index storage variant of all cell and face batches. * * The three arrays given here address the types for the faces decorated * as interior (0), the faces decorated with as exterior (1), and the - * cells (2). + * cells (2) according to CellOrFaceAccess. */ std::vector index_storage_variants[3]; @@ -319,7 +339,7 @@ namespace internal * * The three arrays given here address the types for the faces decorated * as interior (0), the faces decorated with as exterior (1), and the - * cells (2). + * cells (2) according to CellOrFaceAccess. */ std::vector dof_indices_contiguous[3]; @@ -330,7 +350,7 @@ namespace internal * * The three arrays given here address the types for the faces decorated * as interior (0), the faces decorated with as exterior (1), and the - * cells (2). + * cells (2) according to CellOrFaceAccess. */ std::vector n_vectorization_lanes_filled[3]; diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 518b93eb7d..997b0a5c19 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -614,17 +614,17 @@ no_constraint: if (vectorization_length > 1) AssertDimension(row_starts.size()/vectorization_length/n_components, irregular_cells.size()); - index_storage_variants[2].resize(irregular_cells.size(), - IndexStorageVariants::full); - n_vectorization_lanes_filled[2].resize(irregular_cells.size()); + index_storage_variants[dof_access_cell].resize(irregular_cells.size(), + IndexStorageVariants::full); + n_vectorization_lanes_filled[dof_access_cell].resize(irregular_cells.size()); for (unsigned int i=0; i 0) - n_vectorization_lanes_filled[2][i] = irregular_cells[i]; + n_vectorization_lanes_filled[dof_access_cell][i] = irregular_cells[i]; else - n_vectorization_lanes_filled[2][i] = vectorization_length; + n_vectorization_lanes_filled[dof_access_cell][i] = vectorization_length; - dof_indices_contiguous[2].resize(irregular_cells.size()*vectorization_length, - numbers::invalid_unsigned_int); + dof_indices_contiguous[dof_access_cell].resize(irregular_cells.size()*vectorization_length, + numbers::invalid_unsigned_int); dof_indices_interleaved.resize(dof_indices.size(), numbers::invalid_unsigned_int); std::vector index_kinds(static_cast(IndexStorageVariants::contiguous)+1); @@ -632,7 +632,7 @@ no_constraint: for (unsigned int i=0; i no compression possible bool has_constraints = false; @@ -646,7 +646,7 @@ no_constraint: } } if (has_constraints) - index_storage_variants[2][i] = IndexStorageVariants::full; + index_storage_variants[dof_access_cell][i] = IndexStorageVariants::full; else { bool indices_are_contiguous = true; @@ -680,30 +680,33 @@ no_constraint: indices_are_interleaved_and_contiguous) { for (unsigned int j=0; jdof_indices[row_starts[(i*vectorization_length+j)*n_components].first]; } if (indices_are_contiguous) { - index_storage_variants[2][i] = IndexStorageVariants::contiguous; + index_storage_variants[dof_access_cell][i] = IndexStorageVariants::contiguous; } else { const unsigned int *dof_indices = &this->dof_indices[row_starts[i*vectorization_length*n_components].first]; if (n_comp == vectorization_length) - index_storage_variants[2][i] = IndexStorageVariants::interleaved; + index_storage_variants[dof_access_cell][i] = IndexStorageVariants::interleaved; else - index_storage_variants[2][i] = IndexStorageVariants::full; + index_storage_variants[dof_access_cell][i] = IndexStorageVariants::full; + + // do not use interleaved storage if two vectorized + // components point to the same field (scatter not possible) for (unsigned int k=0; kdof_indices_interleaved[row_starts[i*vectorization_length*n_components].first]; @@ -713,7 +716,7 @@ no_constraint: } } } - index_kinds[static_cast(index_storage_variants[2][i])]++; + index_kinds[static_cast(index_storage_variants[dof_access_cell][i])]++; } } @@ -726,22 +729,24 @@ no_constraint: { AssertDimension(length, vectorization_length); - index_storage_variants[0].resize(faces.size(), IndexStorageVariants::full); - dof_indices_contiguous[0].resize(faces.size()*length, numbers::invalid_unsigned_int); - n_vectorization_lanes_filled[0].resize(faces.size()); + index_storage_variants[dof_access_face_interior] + .resize(faces.size(), IndexStorageVariants::full); + dof_indices_contiguous[dof_access_face_interior] + .resize(faces.size()*length, numbers::invalid_unsigned_int); + n_vectorization_lanes_filled[dof_access_face_interior].resize(faces.size()); // all interior faces come before the boundary faces - unsigned int n_plus_faces = 0; - for (; n_plus_faces < faces.size(); ++n_plus_faces) - if (faces[n_plus_faces].cells_exterior[0] == numbers::invalid_unsigned_int) + unsigned int n_exterior_faces = 0; + for (; n_exterior_faces < faces.size(); ++n_exterior_faces) + if (faces[n_exterior_faces].cells_exterior[0] == numbers::invalid_unsigned_int) break; - index_storage_variants[1].resize(n_plus_faces, IndexStorageVariants::full); - dof_indices_contiguous[1].resize(n_plus_faces*length, numbers::invalid_unsigned_int); - n_vectorization_lanes_filled[1].resize(n_plus_faces); + index_storage_variants[dof_access_face_exterior].resize(n_exterior_faces, IndexStorageVariants::full); + dof_indices_contiguous[dof_access_face_exterior].resize(n_exterior_faces*length, numbers::invalid_unsigned_int); + n_vectorization_lanes_filled[dof_access_face_exterior].resize(n_exterior_faces); for (unsigned int face=0; face quadrature_weights (mapping_data->descriptor[active_quad_index].quadrature_weights.begin()), cell (numbers::invalid_unsigned_int), is_interior_face (is_interior_face), - face_vector_access_index (is_face ? (is_interior_face ? 0 : 1) : 2), + dof_access_index (is_face ? + (is_interior_face ? + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior + : + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) + : + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell), cell_type (internal::MatrixFreeFunctions::general), dof_values_initialized (false), values_quad_initialized (false), @@ -2829,7 +2835,7 @@ FEEvaluationBase cell (0), cell_type (internal::MatrixFreeFunctions::general), is_interior_face (true), - face_vector_access_index (numbers::invalid_unsigned_int), + dof_access_index (internal::MatrixFreeFunctions::DoFInfo::dof_access_cell), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2895,7 +2901,7 @@ FEEvaluationBase cell (numbers::invalid_unsigned_int), cell_type (internal::MatrixFreeFunctions::general), is_interior_face (other.is_interior_face), - face_vector_access_index (other.face_vector_access_index), + dof_access_index (other.dof_access_index), dof_values_initialized (false), values_quad_initialized (false), gradients_quad_initialized(false), @@ -2969,7 +2975,7 @@ FEEvaluationBase cell = numbers::invalid_unsigned_int; cell_type = internal::MatrixFreeFunctions::general; is_interior_face = other.is_interior_face; - face_vector_access_index = other.face_vector_access_index; + dof_access_index = other.dof_access_index; // Create deep copy of mapped geometry for use in parallel... if (other.mapped_geometry.get() != nullptr) @@ -3649,8 +3655,9 @@ FEEvaluationBase // Case 2: contiguous indices which use reduced storage of indices and can // use vectorized load/store operations -> go to separate function AssertIndexRange(cell, - dof_info->index_storage_variants[face_vector_access_index].size()); - if (dof_info->index_storage_variants[is_face ? face_vector_access_index : 2][cell] >= + dof_info->index_storage_variants[dof_access_index].size()); + if (dof_info->index_storage_variants[is_face ? dof_access_index : + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] >= internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) { read_write_operation_contiguous(operation, src); @@ -3662,11 +3669,13 @@ FEEvaluationBase constexpr unsigned int n_vectorization = VectorizedArray::n_array_elements; const unsigned int dofs_per_component = this->data->dofs_per_component_on_cell; - if (dof_info->index_storage_variants[is_face ? face_vector_access_index : 2][cell] == + if (dof_info->index_storage_variants[is_face ? dof_access_index : + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved) { const unsigned int *dof_indices = - &dof_info->dof_indices_interleaved[dof_info->row_starts[cell*n_vectorization*n_fe_components+first_selected_component].first]; + &dof_info->dof_indices_interleaved[dof_info->row_starts[cell*n_fe_components*n_vectorization].first]+ + dof_info->component_dof_indices_offset[active_fe_index][first_selected_component]*n_vectorization; if (n_components == 1 || n_fe_components == 1) for (unsigned int i=0; i unsigned int cells_copied[VectorizedArray::n_array_elements]; const unsigned int *cells; unsigned int n_vectorization_actual = - dof_info->n_vectorization_lanes_filled[face_vector_access_index][cell]; + dof_info->n_vectorization_lanes_filled[dof_access_index][cell]; bool has_constraints = false; if (is_face) { - if (face_vector_access_index == 2) + if (dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) for (unsigned int v=0; v::n_array_elements+v; cells = - face_vector_access_index == 2 ? + dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ? &cells_copied[0] : (is_interior_face ? @@ -3953,7 +3962,8 @@ FEEvaluationBase std::integral_constant::value> vector_selector; - const unsigned int ind = is_face ? face_vector_access_index : 2; + const internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex ind = + is_face ? dof_access_index : internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; const std::vector &dof_indices_cont = dof_info->dof_indices_contiguous[ind]; @@ -6362,9 +6372,9 @@ FEFaceEvaluation return; this->cell = face_index; - this->face_vector_access_index = this->is_interior_face ? 0 : 1; - if (face_index >= this->matrix_info->get_task_info().refinement_edge_face_partition_data[0]) - this->face_vector_access_index = 0; + this->dof_access_index = this->is_interior_face ? + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior : + internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior; Assert (this->mapping_data != nullptr, ExcNotInitialized()); const unsigned int n_vectors = VectorizedArray::n_array_elements; const internal::MatrixFreeFunctions::FaceToCellTopology &faces = @@ -6436,7 +6446,7 @@ FEFaceEvaluation this->face_orientation = 0; this->subface_index = GeometryInfo::max_children_per_cell; this->face_no = face_number; - this->face_vector_access_index = 2; + this->dof_access_index = internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; const unsigned int offsets = this->matrix_info->get_mapping_info().face_data_by_cells[this->quad_no]. @@ -6639,10 +6649,10 @@ FEFaceEvaluation internal::VectorReader reader; - if (this->dof_info->index_storage_variants[this->face_vector_access_index][this->cell] == + if (this->dof_info->index_storage_variants[this->dof_access_index][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous && - this->dof_info->n_vectorization_lanes_filled[this->face_vector_access_index][this->cell] == + this->dof_info->n_vectorization_lanes_filled[this->dof_access_index][this->cell] == VectorizedArray::n_array_elements && ((evaluate_gradients == false && this->data->nodal_at_cell_boundaries == true) || @@ -6650,7 +6660,7 @@ FEFaceEvaluation fe_degree > 1))) { const unsigned int *indices = &this->dof_info->dof_indices_contiguous - [this->face_vector_access_index][this->cell*VectorizedArray::n_array_elements]; + [this->dof_access_index][this->cell*VectorizedArray::n_array_elements]; if (evaluate_gradients == true && this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite) { @@ -6789,10 +6799,10 @@ FEFaceEvaluation internal::VectorDistributorLocalToGlobal writer; - if (this->dof_info->index_storage_variants[this->face_vector_access_index][this->cell] == + if (this->dof_info->index_storage_variants[this->dof_access_index][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous && - this->dof_info->n_vectorization_lanes_filled[this->face_vector_access_index][this->cell] == + this->dof_info->n_vectorization_lanes_filled[this->dof_access_index][this->cell] == VectorizedArray::n_array_elements && ((integrate_gradients == false && this->data->nodal_at_cell_boundaries == true) || @@ -6800,7 +6810,7 @@ FEFaceEvaluation fe_degree > 1))) { const unsigned int *indices = &this->dof_info->dof_indices_contiguous - [this->face_vector_access_index][this->cell*VectorizedArray::n_array_elements]; + [this->dof_access_index][this->cell*VectorizedArray::n_array_elements]; if (integrate_gradients == true && this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite) @@ -6914,7 +6924,7 @@ FEFaceEvaluation ::quadrature_point (const unsigned int q) const { AssertIndexRange (q, n_q_points); - if (this->face_vector_access_index < 2) + if (this->dof_access_index < 2) { Assert(this->mapping_data->quadrature_point_offsets.empty() == false, ExcNotImplemented()); diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 8fb6a6244f..2518d3cd29 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1333,10 +1333,10 @@ void MatrixFree::initialize_indices AssertIndexRange(face_info.faces[f].cells_interior[v], n_macro_cells_before* VectorizedArray::n_array_elements); - if (dof_info[no].index_storage_variants[1][f] >= + if (dof_info[no].index_storage_variants[internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior][f] >= internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous && - dof_info[no].dof_indices_contiguous[2][face_info.faces[f].cells_exterior[v]] >= + dof_info[no].dof_indices_contiguous[internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][face_info.faces[f].cells_exterior[v]] >= part.local_size()) { const unsigned int p = face_info.faces[f].cells_exterior[v]; @@ -1349,7 +1349,7 @@ void MatrixFree::initialize_indices = shape_info(dof_info[no].global_base_element_offset+e,0,0,0); for (unsigned int j=0; j::initialize_indices } AssertDimension(i, dof_info[no].dofs_per_cell[0]*stride); } - else if (dof_info[no].index_storage_variants[1][f] < + else if (dof_info[no].index_storage_variants[internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior][f] < internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) has_noncontiguous_cell = true; } @@ -1411,10 +1411,10 @@ void MatrixFree::initialize_indices { AssertIndexRange(face_info.faces[f].cells_interior[v], n_macro_cells_before*VectorizedArray::n_array_elements); - if (dof_info[no].index_storage_variants[1][f] >= + if (dof_info[no].index_storage_variants[internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior][f] >= internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous && - dof_info[no].dof_indices_contiguous[2][face_info.faces[f].cells_exterior[v]] >= + dof_info[no].dof_indices_contiguous[internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][face_info.faces[f].cells_exterior[v]] >= part.local_size()) { const unsigned int p = face_info.faces[f].cells_exterior[v]; @@ -1427,7 +1427,7 @@ void MatrixFree::initialize_indices = shape_info(dof_info[no].global_base_element_offset+e,0,0,0); for (unsigned int j=0; j<2*shape.dofs_per_component_on_face; ++j) ghost_indices. - push_back(part.local_to_global(dof_info[no].dof_indices_contiguous[2][p] + push_back(part.local_to_global(dof_info[no].dof_indices_contiguous[internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][p] +i+ shape.face_to_cell_index_hermite (face_info.faces[f].exterior_face_no,j)*stride)); -- 2.39.5