From: Martin Kronbichler Date: Mon, 6 Aug 2018 09:01:28 +0000 (+0200) Subject: Add contiguous-interleaved storage option in MatrixFree. X-Git-Tag: v9.1.0-rc1~806^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=823ee33c1dc9a8c8fdb96f6e0b6d3a35106d2855;p=dealii.git Add contiguous-interleaved storage option in MatrixFree. --- diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index 3c558f06d0..6839cd1d35 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -265,7 +265,56 @@ namespace internal * is directed to the array `dof_indices_contiguous` with the index * `cell_index*n_vectorization*n_components`. */ - contiguous + contiguous, + /** + * This value indicates that the indices with a cell are contiguous and + * interleaved for vectorization, i.e., the first DoF index on a cell + * to the four or eight cells in the vectorization batch come first, + * than the second DoF index, and so on. Furthermore, the interleaving + * between cells implies that only the batches for vectorization can be + * accessed efficiently, whereas there is a strided access for getting + * only some of the entries. + * + * The two additional categories `interleaved_contiguous_strided` and + * `interleaved_contiguous_mixed_strides` are a consequence of this + * storage type. The former is for faces where at least one of the two + * adjacent sides will break with the interleaved storage. We then have + * to make a strided access as described in the next category. The last + * category `interleaved_contiguous_mixed_strides` appears in the ghost + * layer, see the more detailed description of that category below. + * Again, this is something that cannot be avoided in general once we + * interleave the indices between cells. + * + * For a cell/face of this index type, the data access in + * FEEvaluationBase is directed to the array `dof_indices_contiguous` + * with the index `cell_index*n_vectorization*n_components`. + */ + interleaved_contiguous, + /** + * Similar to interleaved_contiguous storage, but for the case when the + * interleaved indices are only contiguous within the degrees of + * freedom, but not also over the components of a vectorized array. + * This happens typically on faces with DG where the cells have + * `interleaved_contiguous` storage but the faces' numbering is not the + * same as the cell's numbering. For a + * cell/face of this index type, the data access in FEEvaluationBase + * is directed to the array `dof_indices_contiguous` with the index + * `cell_index*n_vectorization*n_components`. + */ + interleaved_contiguous_strided, + /** + * Similar to interleaved_contiguous_separate storage, but for the case + * when the interleaved indices are not `n_vectorization apart`. This + * happens typically within the ghost layer of DG where the remote + * owner has applied an interleaved storage and the current processor + * only sees some of the cells. For a + * cell/face of this index type, the data access in FEEvaluationBase + * is directed to the array `dof_indices_contiguous` with the index + * `cell_index*n_vectorization*n_components`, including the array + * `dof_indices_interleave_strides` for the information about the + * actual stride. + */ + interleaved_contiguous_mixed_strides }; /** @@ -366,6 +415,16 @@ namespace internal */ std::vector dof_indices_contiguous[3]; + /** + * Compressed index storage for faster access than through @p + * dof_indices used according to the description in IndexStorageVariants. + * + * The three arrays given here address the types for the faces decorated + * as minus (0), the faces decorated with as plus (1), and the cells + * (2). + */ + std::vector dof_indices_interleave_strides[3]; + /** * Caches the number of indices filled when vectorizing. This * information can implicitly deduced from the row_starts data fields, diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 844994c95b..caaba870dd 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -167,6 +167,14 @@ namespace internal start_components.clear(); row_starts_plain_indices.clear(); plain_dof_indices.clear(); + dof_indices_interleaved.clear(); + for (unsigned int i = 0; i < 3; ++i) + { + index_storage_variants[i].clear(); + dof_indices_contiguous[i].clear(); + dof_indices_interleave_strides[i].clear(); + n_vectorization_lanes_filled[i].clear(); + } store_plain_indices = false; cell_active_fe_index.clear(); max_fe_index = 0; @@ -699,9 +707,14 @@ namespace internal numbers::invalid_unsigned_int); dof_indices_interleaved.resize(dof_indices.size(), numbers::invalid_unsigned_int); + dof_indices_interleave_strides[dof_access_cell].resize( + irregular_cells.size() * vectorization_length, + numbers::invalid_unsigned_int); std::vector index_kinds( - static_cast(IndexStorageVariants::contiguous) + 1); + static_cast( + IndexStorageVariants::interleaved_contiguous_mixed_strides) + + 1); std::vector offsets(vectorization_length); for (unsigned int i = 0; i < irregular_cells.size(); ++i) { @@ -745,8 +758,10 @@ namespace internal break; } } + bool indices_are_interleaved_and_contiguous = (ndofs > 1 && n_comp == vectorization_length); + { const unsigned int *dof_indices = this->dof_indices.data() + @@ -760,6 +775,7 @@ namespace internal break; } } + if (indices_are_contiguous || indices_are_interleaved_and_contiguous) { @@ -772,52 +788,209 @@ namespace internal .first]; } - if (indices_are_contiguous) + if (indices_are_interleaved_and_contiguous) + { + Assert(n_comp == vectorization_length, ExcInternalError()); + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants::interleaved_contiguous; + for (unsigned int j = 0; j < n_comp; ++j) + dof_indices_interleave_strides[2][i * vectorization_length + + j] = n_comp; + } + else if (indices_are_contiguous) { index_storage_variants[dof_access_cell][i] = IndexStorageVariants::contiguous; + for (unsigned int j = 0; j < n_comp; ++j) + dof_indices_interleave_strides[2][i * vectorization_length + + j] = 1; } else { + int indices_are_interleaved_and_mixed = 2; const unsigned int *dof_indices = - this->dof_indices.data() + - row_starts[i * vectorization_length * n_components].first; - if (n_comp == vectorization_length) - index_storage_variants[dof_access_cell][i] = - IndexStorageVariants::interleaved; - else - 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) + &this->dof_indices[row_starts[i * vectorization_length * + n_components] + .first]; + for (unsigned int j = 0; j < n_comp; ++j) + offsets[j] = + dof_indices[j * ndofs + 1] - dof_indices[j * ndofs]; for (unsigned int k = 0; k < ndofs; ++k) - for (unsigned int l = 0; l < n_comp; ++l) - for (unsigned int j = l + 1; j < n_comp; ++j) - if (dof_indices[j * ndofs + k] == - dof_indices[l * ndofs + k]) + for (unsigned int j = 0; j < n_comp; ++j) + if (dof_indices[j * ndofs + k] != + dof_indices[j * ndofs] + k * offsets[j]) + { + indices_are_interleaved_and_mixed = 0; + break; + } + if (indices_are_interleaved_and_mixed == 2) + { + for (unsigned int j = 0; j < n_comp; ++j) + dof_indices_interleave_strides + [dof_access_cell][i * vectorization_length + j] = + offsets[j]; + for (unsigned int j = 0; j < n_comp; ++j) + dof_indices_contiguous[dof_access_cell] + [i * vectorization_length + j] = + dof_indices[j * ndofs]; + for (unsigned int j = 0; j < n_comp; ++j) + if (offsets[j] != vectorization_length) { - index_storage_variants[dof_access_cell][i] = - IndexStorageVariants::full; + indices_are_interleaved_and_mixed = 1; break; } - if (index_storage_variants[dof_access_cell][i] != - IndexStorageVariants::full) + if (indices_are_interleaved_and_mixed == 1 || + n_comp != vectorization_length) + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants:: + interleaved_contiguous_mixed_strides; + else + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants::interleaved_contiguous_strided; + } + else { - unsigned int *interleaved_dof_indices = - this->dof_indices_interleaved.data() + + const unsigned int *dof_indices = + this->dof_indices.data() + row_starts[i * vectorization_length * n_components] .first; + if (n_comp == vectorization_length) + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants::interleaved; + else + 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; k < ndofs; ++k) - for (unsigned int j = 0; j < n_comp; ++j) - interleaved_dof_indices[k * n_comp + j] = - dof_indices[j * ndofs + k]; + for (unsigned int l = 0; l < n_comp; ++l) + for (unsigned int j = l + 1; j < n_comp; ++j) + if (dof_indices[j * ndofs + k] == + dof_indices[l * ndofs + k]) + { + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants::full; + break; + } } } } index_kinds[static_cast( index_storage_variants[dof_access_cell][i])]++; } + + // Cleanup phase: we want to avoid single cells with different properties + // than the bulk of the domain in order to avoid extra checks in the face + // identification. + + // Step 1: check whether the interleaved indices were only assigned to + // the single cell within a vectorized array. + auto fix_single_interleaved_indices = + [&](const IndexStorageVariants variant) { + if (index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous_mixed_strides)] > + 0 && + index_kinds[static_cast(variant)] > 0) + for (unsigned int i = 0; i < irregular_cells.size(); ++i) + { + if (index_storage_variants[dof_access_cell][i] == + IndexStorageVariants:: + interleaved_contiguous_mixed_strides && + n_vectorization_lanes_filled[dof_access_cell][i] == 1 && + (variant != IndexStorageVariants::contiguous || + dof_indices_interleave_strides[dof_access_cell] + [i * vectorization_length] == + 1)) + { + index_storage_variants[dof_access_cell][i] = variant; + index_kinds[static_cast( + IndexStorageVariants:: + interleaved_contiguous_mixed_strides)]--; + index_kinds[static_cast(variant)]++; + } + } + }; + + fix_single_interleaved_indices(IndexStorageVariants::full); + fix_single_interleaved_indices(IndexStorageVariants::contiguous); + fix_single_interleaved_indices(IndexStorageVariants::interleaved); + + unsigned int n_interleaved = + index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous)] + + index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous_strided)] + + index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous_mixed_strides)]; + + // Step 2: fix single contiguous cell among others with interleaved + // storage + if (n_interleaved > 0 && index_kinds[static_cast( + IndexStorageVariants::contiguous)]) + for (unsigned int i = 0; i < irregular_cells.size(); ++i) + if (index_storage_variants[dof_access_cell][i] == + IndexStorageVariants::contiguous) + { + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants::interleaved_contiguous_mixed_strides; + index_kinds[static_cast( + IndexStorageVariants::contiguous)]--; + index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous_mixed_strides)]++; + } + + n_interleaved = + index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous)] + + index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous_strided)] + + index_kinds[static_cast( + IndexStorageVariants::interleaved_contiguous_mixed_strides)]; + + // Step 3: Interleaved cells are left but also some non-contiguous ones + // -> revert all to full storage + if (n_interleaved > 0 && + index_kinds[static_cast(IndexStorageVariants::full)] + + index_kinds[static_cast( + IndexStorageVariants::interleaved)] > + 0) + for (unsigned int i = 0; i < irregular_cells.size(); ++i) + if (index_storage_variants[dof_access_cell][i] > + IndexStorageVariants::contiguous) + { + index_kinds[static_cast( + index_storage_variants[2][i])]--; + if (n_vectorization_lanes_filled[dof_access_cell][i] == + vectorization_length) + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants::interleaved; + else + index_storage_variants[dof_access_cell][i] = + IndexStorageVariants::full; + index_kinds[static_cast( + index_storage_variants[dof_access_cell][i])]++; + } + + // Step 4: Copy the interleaved indices into their own data structure + for (unsigned int i = 0; i < irregular_cells.size(); ++i) + if (index_storage_variants[dof_access_cell][i] == + IndexStorageVariants::interleaved) + { + const unsigned int ndofs = + dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0]; + const unsigned int *dof_indices = + &this->dof_indices + [row_starts[i * vectorization_length * n_components].first]; + unsigned int *interleaved_dof_indices = + &this->dof_indices_interleaved + [row_starts[i * vectorization_length * n_components].first]; + for (unsigned int k = 0; k < ndofs; ++k) + for (unsigned int j = 0; j < vectorization_length; ++j) + interleaved_dof_indices[k * vectorization_length + j] = + dof_indices[j * ndofs + k]; + } } @@ -833,6 +1006,8 @@ namespace internal faces.size(), IndexStorageVariants::full); dof_indices_contiguous[dof_access_face_interior].resize( faces.size() * length, numbers::invalid_unsigned_int); + dof_indices_interleave_strides[dof_access_face_interior].resize( + faces.size() * length, numbers::invalid_unsigned_int); n_vectorization_lanes_filled[dof_access_face_interior].resize( faces.size()); @@ -846,6 +1021,8 @@ namespace internal n_exterior_faces, IndexStorageVariants::full); dof_indices_contiguous[dof_access_face_exterior].resize( n_exterior_faces * length, numbers::invalid_unsigned_int); + dof_indices_interleave_strides[dof_access_face_exterior].resize( + faces.size() * length, numbers::invalid_unsigned_int); n_vectorization_lanes_filled[dof_access_face_exterior].resize( n_exterior_faces); @@ -854,6 +1031,7 @@ namespace internal auto face_computation = [&](const DoFAccessIndex face_index, const unsigned int * cell_indices_face) { bool is_contiguous = false; + bool is_interleaved = false; bool needs_full_storage = false; for (unsigned int v = 0; v < length && @@ -861,22 +1039,74 @@ namespace internal ++v) { n_vectorization_lanes_filled[face_index][face]++; + if (index_storage_variants[dof_access_cell] + [cell_indices_face[v] / length] >= + IndexStorageVariants::interleaved_contiguous) + is_interleaved = true; if (index_storage_variants[dof_access_cell] [cell_indices_face[v] / length] == IndexStorageVariants::contiguous) is_contiguous = true; + if (index_storage_variants[dof_access_cell] + [cell_indices_face[v] / length] >= + IndexStorageVariants::contiguous) + dof_indices_interleave_strides[face_index][face * length + + v] = + dof_indices_interleave_strides[dof_access_cell] + [cell_indices_face[v]]; if (index_storage_variants[dof_access_cell] [cell_indices_face[v] / length] < IndexStorageVariants::contiguous) needs_full_storage = true; } - if (is_contiguous) + Assert(!(is_interleaved && is_contiguous), + ExcMessage("Unsupported index compression found")); + + if (is_interleaved || is_contiguous) for (unsigned int v = 0; v < n_vectorization_lanes_filled[face_index][face]; ++v) dof_indices_contiguous[face_index][face * length + v] = dof_indices_contiguous[dof_access_cell][cell_indices_face[v]]; - if (is_contiguous && !needs_full_storage) + if (is_interleaved) + { + bool is_also_contiguous = + n_vectorization_lanes_filled[face_index][face] == length; + for (unsigned int v = 0; + v < n_vectorization_lanes_filled[face_index][face]; + ++v) + if (dof_indices_contiguous[face_index][face * length + v] != + dof_indices_contiguous[face_index][face * length] + v || + dof_indices_interleave_strides[dof_access_cell] + [cell_indices_face[v]] != + length) + is_also_contiguous = false; + + if (is_also_contiguous) + { + index_storage_variants[face_index][face] = + IndexStorageVariants::interleaved_contiguous; + } + else + { + bool all_indices_same_offset = + n_vectorization_lanes_filled[face_index][face] == length; + for (unsigned int v = 0; + v < n_vectorization_lanes_filled[face_index][face]; + ++v) + if (dof_indices_interleave_strides + [dof_access_cell][cell_indices_face[v]] != length) + all_indices_same_offset = false; + if (all_indices_same_offset) + index_storage_variants[face_index][face] = + IndexStorageVariants::interleaved_contiguous_strided; + else + index_storage_variants[face_index][face] = + IndexStorageVariants:: + interleaved_contiguous_mixed_strides; + } + } + else if (is_contiguous && !needs_full_storage) index_storage_variants[face_index][face] = IndexStorageVariants::contiguous; else diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h index b8cd0248cc..2fd382aa06 100644 --- a/include/deal.II/matrix_free/evaluation_selector.h +++ b/include/deal.II/matrix_free/evaluation_selector.h @@ -86,7 +86,8 @@ namespace internal Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients) + const bool integrate_gradients, + const bool sum_into_values_array = false) { internal::FEEvaluationImpl< internal::MatrixFreeFunctions::tensor_general, @@ -101,7 +102,7 @@ namespace internal scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } }; @@ -205,7 +206,8 @@ namespace internal Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients) + const bool integrate_gradients, + const bool sum_into_values_array = false) { const int runtime_degree = shape_info.fe_degree; constexpr unsigned int start_n_q_points = degree + 1; @@ -217,7 +219,8 @@ namespace internal gradients_quad, scratch_data, integrate_values, - integrate_gradients); + integrate_gradients, + sum_into_values_array); else Factory:: integrate(shape_info, @@ -226,7 +229,8 @@ namespace internal gradients_quad, scratch_data, integrate_values, - integrate_gradients); + integrate_gradients, + sum_into_values_array); } }; @@ -329,7 +333,8 @@ namespace internal Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients) + const bool integrate_gradients, + const bool sum_into_values_array) { const int runtime_n_q_points_1d = shape_info.n_q_points_1d; if (runtime_n_q_points_1d == n_q_points_1d) @@ -346,7 +351,7 @@ namespace internal scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); else if (degree < n_q_points_1d) internal::FEEvaluationImplTransformToCollocation< dim, @@ -360,7 +365,7 @@ namespace internal scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); else internal::FEEvaluationImpl< internal::MatrixFreeFunctions::tensor_symmetric, @@ -375,7 +380,7 @@ namespace internal scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else Factory:: @@ -385,7 +390,8 @@ namespace internal gradients_quad, scratch_data, integrate_values, - integrate_gradients); + integrate_gradients, + sum_into_values_array); } }; @@ -437,7 +443,8 @@ namespace internal Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients) + const bool integrate_gradients, + const bool sum_into_values_array = false) { Assert(shape_info.element_type <= internal::MatrixFreeFunctions::tensor_symmetric, @@ -448,7 +455,8 @@ namespace internal gradients_quad, scratch_data, integrate_values, - integrate_gradients); + integrate_gradients, + sum_into_values_array); } } // namespace EvaluationSelectorImplementation } // namespace internal @@ -505,7 +513,8 @@ struct SelectEvaluator Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients); + const bool integrate_gradients, + const bool sum_into_values_array = false); }; /** @@ -555,7 +564,8 @@ struct SelectEvaluator Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients); + const bool integrate_gradients, + const bool sum_into_values_array = false); }; //----------------------Implementation for SelectEvaluator--------------------- @@ -709,7 +719,8 @@ SelectEvaluator::integrate( Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients) + const bool integrate_gradients, + const bool sum_into_values_array) { Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError()); @@ -726,7 +737,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else if (fe_degree < n_q_points_1d && shape_info.element_type <= @@ -744,7 +755,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric) @@ -762,7 +773,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) @@ -780,7 +791,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor) @@ -798,7 +809,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general) @@ -815,7 +826,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else AssertThrow(false, ExcNotImplemented()); @@ -914,7 +925,8 @@ SelectEvaluator::integrate( Number * gradients_quad, Number * scratch_data, const bool integrate_values, - const bool integrate_gradients) + const bool integrate_gradients, + const bool sum_into_values_array) { if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0) @@ -932,7 +944,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor) @@ -950,7 +962,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); } else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general) @@ -966,7 +978,7 @@ SelectEvaluator::integrate( scratch_data, integrate_values, integrate_gradients, - false); + sum_into_values_array); else internal::EvaluationSelectorImplementation:: symmetric_selector_integrate( @@ -976,7 +988,8 @@ SelectEvaluator::integrate( gradients_quad, scratch_data, integrate_values, - integrate_gradients); + integrate_gradients, + sum_into_values_array); } #endif // DOXYGEN diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 1459d498d8..daafb5ad8f 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3454,6 +3454,36 @@ namespace internal res = vector_access(vec, index); } + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + const Number *vec_ptr = vec.begin() + dof_index; + for (unsigned int i = 0; i < dofs_per_cell; + ++i, vec_ptr += VectorizedArray::n_array_elements) + dof_values[i].load(vec_ptr); + } + + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; + ++v) + dof_values[i][v] = vector_access( + vec, dof_index + v + i * VectorizedArray::n_array_elements); + } + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -3560,6 +3590,42 @@ namespace internal vector_access(vec, index) += res; } + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + Number *vec_ptr = vec.begin() + dof_index; + for (unsigned int i = 0; i < dofs_per_cell; + ++i, vec_ptr += VectorizedArray::n_array_elements) + { + VectorizedArray tmp; + tmp.load(vec_ptr); + tmp += dof_values[i]; + tmp.store(vec_ptr); + } + } + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; + ++v) + vector_access(vec, + dof_index + v + + i * VectorizedArray::n_array_elements) += + dof_values[i][v]; + } + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -3670,6 +3736,37 @@ namespace internal vector_access(vec, index) = res; } + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + Number *vec_ptr = vec.begin() + dof_index; + for (unsigned int i = 0; i < dofs_per_cell; + ++i, vec_ptr += VectorizedArray::n_array_elements) + dof_values[i].store(vec_ptr); + } + + template + void + process_dofs_vectorized(const unsigned int dofs_per_cell, + const unsigned int dof_index, + VectorType & vec, + VectorizedArray *dof_values, + std::integral_constant) const + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; + ++v) + vector_access(vec, + dof_index + v + + i * VectorizedArray::n_array_elements) = + dof_values[i][v]; + } + template void process_dofs_vectorized_transpose(const unsigned int dofs_per_cell, @@ -4294,6 +4391,37 @@ FEEvaluationBase:: const std::vector &dof_indices_cont = dof_info->dof_indices_contiguous[ind]; + + // Simple case: We have contiguous storage, so we can simply copy out the + // data + if (dof_info->index_storage_variants[ind][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous) + { + const unsigned int dof_index = + dof_indices_cont[cell * VectorizedArray::n_array_elements] + + dof_info->component_dof_indices_offset[active_fe_index] + [first_selected_component] * + VectorizedArray::n_array_elements; + if (n_components == 1 || n_fe_components == 1) + for (unsigned int comp = 0; comp < n_components; ++comp) + operation.process_dofs_vectorized(data->dofs_per_component_on_cell, + dof_index, + *src[comp], + values_dofs[comp], + vector_selector); + else + operation.process_dofs_vectorized(data->dofs_per_component_on_cell * + n_components, + dof_index, + *src[0], + values_dofs[0], + vector_selector); + return; + } + + // More general case: Must go through the components one by one and apply + // some transformations const unsigned int vectorization_populated = dof_info->n_vectorization_lanes_filled[ind][this->cell]; unsigned int dof_indices[VectorizedArray::n_array_elements]; @@ -4301,7 +4429,9 @@ FEEvaluationBase:: dof_indices[v] = dof_indices_cont[cell * VectorizedArray::n_array_elements + v] + dof_info->component_dof_indices_offset[active_fe_index] - [first_selected_component]; + [first_selected_component] * + dof_info->dof_indices_interleave_strides + [ind][cell * VectorizedArray::n_array_elements + v]; for (unsigned int v = vectorization_populated; v < VectorizedArray::n_array_elements; ++v) @@ -4311,40 +4441,149 @@ FEEvaluationBase:: // constraints and that the indices within each element are contiguous if (vectorization_populated == VectorizedArray::n_array_elements) { - if (n_components == 1 || n_fe_components == 1) - for (unsigned int comp = 0; comp < n_components; ++comp) - operation.process_dofs_vectorized_transpose( - data->dofs_per_component_on_cell, - dof_indices, - *src[comp], - values_dofs[comp], - vector_selector); + if (dof_info->index_storage_variants[ind][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + contiguous) + { + if (n_components == 1 || n_fe_components == 1) + for (unsigned int comp = 0; comp < n_components; ++comp) + operation.process_dofs_vectorized_transpose( + data->dofs_per_component_on_cell, + dof_indices, + *src[comp], + values_dofs[comp], + vector_selector); + else + operation.process_dofs_vectorized_transpose( + data->dofs_per_component_on_cell * n_components, + dof_indices, + *src[0], + &values_dofs[0][0], + vector_selector); + } + else if (dof_info->index_storage_variants[ind][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_strided) + { + if (n_components == 1 || n_fe_components == 1) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) + { + for (unsigned int comp = 0; comp < n_components; ++comp) + operation.process_dof_gather( + dof_indices, + *src[comp], + i * VectorizedArray::n_array_elements, + values_dofs[comp][i], + vector_selector); + } + else + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + ++i) + { + operation.process_dof_gather( + dof_indices, + *src[0], + (comp * data->dofs_per_component_on_cell + i) * + VectorizedArray::n_array_elements, + values_dofs[comp][i], + vector_selector); + } + } else - operation.process_dofs_vectorized_transpose( - data->dofs_per_component_on_cell * n_components, - dof_indices, - *src[0], - &values_dofs[0][0], - vector_selector); + { + Assert(dof_info->index_storage_variants[ind][cell] == + internal::MatrixFreeFunctions::DoFInfo:: + IndexStorageVariants::interleaved_contiguous_mixed_strides, + ExcNotImplemented()); + const unsigned int *offsets = + &dof_info->dof_indices_interleave_strides + [ind][VectorizedArray::n_array_elements * cell]; + if (n_components == 1 || n_fe_components == 1) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) + { + for (unsigned int comp = 0; comp < n_components; ++comp) + operation.process_dof_gather(dof_indices, + *src[comp], + 0, + values_dofs[comp][i], + vector_selector); + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + dof_indices[v] += offsets[v]; + } + else + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + ++i) + { + operation.process_dof_gather(dof_indices, + *src[0], + 0, + values_dofs[comp][i], + vector_selector); + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + dof_indices[v] += offsets[v]; + } + } } else for (unsigned int comp = 0; comp < n_components; ++comp) { for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) operation.process_empty(values_dofs[comp][i]); - if (n_components == 1 || n_fe_components == 1) - for (unsigned int v = 0; v < vectorization_populated; ++v) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) - operation.process_dof(dof_indices[v] + i, - *src[comp], - values_dofs[comp][i][v]); + if (dof_info->index_storage_variants[ind][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + contiguous) + { + if (n_components == 1 || n_fe_components == 1) + for (unsigned int v = 0; v < vectorization_populated; ++v) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof(dof_indices[v] + i, + *src[comp], + values_dofs[comp][i][v]); + else + for (unsigned int v = 0; v < vectorization_populated; ++v) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof(dof_indices[v] + i + + comp * + data->dofs_per_component_on_cell, + *src[0], + values_dofs[comp][i][v]); + } else - for (unsigned int v = 0; v < vectorization_populated; ++v) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) - operation.process_dof(dof_indices[v] + i + - comp * data->dofs_per_component_on_cell, - *src[0], - values_dofs[comp][i][v]); + { + const unsigned int *offsets = + &dof_info->dof_indices_interleave_strides + [ind][VectorizedArray::n_array_elements * cell]; + for (unsigned int v = 0; v < vectorization_populated; ++v) + AssertIndexRange(offsets[v], + VectorizedArray::n_array_elements + 1); + if (n_components == 1 || n_fe_components == 1) + for (unsigned int v = 0; v < vectorization_populated; ++v) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof(dof_indices[v] + i * offsets[v], + *src[comp], + values_dofs[comp][i][v]); + else + for (unsigned int v = 0; v < vectorization_populated; ++v) + for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof( + dof_indices[v] + + (i + comp * data->dofs_per_component_on_cell) * + offsets[v], + *src[0], + values_dofs[comp][i][v]); + } } } @@ -6646,11 +6885,42 @@ FEEvaluation:: const bool evaluate_gradients, const bool evaluate_hessians) { - this->read_dof_values(input_vector); - evaluate(this->begin_dof_values(), - evaluate_values, - evaluate_gradients, - evaluate_hessians); + if (this->dof_info->index_storage_variants + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell] == internal::MatrixFreeFunctions::DoFInfo:: + IndexStorageVariants::interleaved_contiguous && + reinterpret_cast( + input_vector.begin() + + this->dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell * VectorizedArray::n_array_elements]) % + sizeof(VectorizedArray) == + 0) + { + const VectorizedArray *vec_values = + reinterpret_cast *>( + input_vector.begin() + + this->dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell * VectorizedArray::n_array_elements] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements); + + evaluate(vec_values, + evaluate_values, + evaluate_gradients, + evaluate_hessians); + } + else + { + this->read_dof_values(input_vector); + evaluate(this->begin_dof_values(), + evaluate_values, + evaluate_gradients, + evaluate_hessians); + } } @@ -6706,7 +6976,8 @@ FEEvaluation::integrate( ->gradients_quad[0][0], this->scratch_data, integrate_values, - integrate_gradients); + integrate_gradients, + false); # ifdef DEBUG this->dof_values_initialized = true; @@ -6727,8 +6998,49 @@ FEEvaluation:: const bool integrate_gradients, VectorType &destination) { - integrate(integrate_values, integrate_gradients, this->begin_dof_values()); - this->distribute_local_to_global(destination); + if (this->dof_info->index_storage_variants + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell] == internal::MatrixFreeFunctions::DoFInfo:: + IndexStorageVariants::interleaved_contiguous && + reinterpret_cast( + destination.begin() + + this->dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell * VectorizedArray::n_array_elements]) % + sizeof(VectorizedArray) == + 0) + { + VectorizedArray *vec_values = + reinterpret_cast *>( + destination.begin() + + this->dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [this->cell * VectorizedArray::n_array_elements] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements); + SelectEvaluator< + dim, + fe_degree, + n_q_points_1d, + n_components, + VectorizedArray>::integrate(*this->data, + vec_values, + this->values_quad[0], + this->gradients_quad[0][0], + this->scratch_data, + integrate_values, + integrate_gradients, + true); + } + else + { + integrate(integrate_values, + integrate_gradients, + this->begin_dof_values()); + this->distribute_local_to_global(destination); + } } @@ -7167,19 +7479,357 @@ FEFaceEvaluation:: temp1 = this->scratch_data; internal::VectorReader reader; + std::integral_constant< + bool, + std::is_same::value> + vector_selector; - 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->dof_access_index][this->cell] == - VectorizedArray::n_array_elements && - ((evaluate_gradients == false && + // case 1: contiguous and interleaved indices + if (((evaluate_gradients == false && this->data->nodal_at_cell_boundaries == true) || (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1))) + fe_degree > 1)) && + this->dof_info + ->index_storage_variants[this->dof_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous) + { + AssertDimension( + this->dof_info + ->n_vectorization_lanes_filled[this->dof_access_index][this->cell], + VectorizedArray::n_array_elements); + const unsigned int dof_index = + this->dof_info + ->dof_indices_contiguous[this->dof_access_index] + [this->cell * + VectorizedArray::n_array_elements] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements; + + if (fe_degree > 1 && evaluate_gradients == true) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 1 + side]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_hermite(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = index_array[2 * i]; + const unsigned int ind2 = index_array[2 * i + 1]; + AssertIndexRange(ind1, dofs_per_cell); + AssertIndexRange(ind2, dofs_per_cell); + for (unsigned int comp = 0; comp < n_components_; ++comp) + { + reader.process_dofs_vectorized( + 1, + dof_index + (ind1 + comp * static_dofs_per_component) * + VectorizedArray::n_array_elements, + input_vector, + temp1 + i + 2 * comp * dofs_per_face, + vector_selector); + reader.process_dofs_vectorized( + 1, + dof_index + (ind2 + comp * static_dofs_per_component) * + VectorizedArray::n_array_elements, + input_vector, + temp1 + dofs_per_face + i + 2 * comp * dofs_per_face, + vector_selector); + temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = + grad_weight * + (temp1[i + 2 * comp * dofs_per_face] - + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]); + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_nodal(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = index_array[i]; + for (unsigned int comp = 0; comp < n_components_; ++comp) + reader.process_dofs_vectorized( + 1, + dof_index + (ind + comp * static_dofs_per_component) * + VectorizedArray::n_array_elements, + input_vector, + temp1 + i + 2 * comp * dofs_per_face, + vector_selector); + } + } + } + + // case 2: contiguous and interleaved indices with fixed stride + else if (((evaluate_gradients == false && + this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + this->dof_info + ->index_storage_variants[this->dof_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_strided) + { + AssertDimension( + this->dof_info + ->n_vectorization_lanes_filled[this->dof_access_index][this->cell], + VectorizedArray::n_array_elements); + const unsigned int *indices = + &this->dof_info + ->dof_indices_contiguous[this->dof_access_index] + [this->cell * + VectorizedArray::n_array_elements]; + if (fe_degree > 1 && evaluate_gradients == true) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 1 + side]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &this->data->face_to_cell_index_hermite(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + index_array[2 * i] * VectorizedArray::n_array_elements; + const unsigned int ind2 = + index_array[2 * i + 1] * + VectorizedArray::n_array_elements; + for (unsigned int comp = 0; comp < n_components_; ++comp) + { + reader.process_dof_gather( + indices, + input_vector, + ind1 + + comp * static_dofs_per_component * + VectorizedArray::n_array_elements + + this->dof_info->component_dof_indices_offset + [this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements, + temp1[i + 2 * comp * dofs_per_face], + vector_selector); + VectorizedArray grad; + reader.process_dof_gather( + indices, + input_vector, + ind2 + + comp * static_dofs_per_component * + VectorizedArray::n_array_elements + + this->dof_info->component_dof_indices_offset + [this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements, + grad, + vector_selector); + temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = + grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad); + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_nodal(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = + index_array[i] * VectorizedArray::n_array_elements; + for (unsigned int comp = 0; comp < n_components_; ++comp) + reader.process_dof_gather( + indices, + input_vector, + ind + + comp * static_dofs_per_component * + VectorizedArray::n_array_elements + + this->dof_info->component_dof_indices_offset + [this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements, + temp1[i + 2 * comp * dofs_per_face], + vector_selector); + } + } + } + + // case 3: contiguous and interleaved indices with mixed stride + else if (((evaluate_gradients == false && + this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + this->dof_info + ->index_storage_variants[this->dof_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_mixed_strides) + { + const unsigned int *strides = + &this->dof_info->dof_indices_interleave_strides + [this->dof_access_index] + [this->cell * VectorizedArray::n_array_elements]; + unsigned int indices[VectorizedArray::n_array_elements]; + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; + ++v) + indices[v] = + this->dof_info->dof_indices_contiguous + [this->dof_access_index] + [this->cell * VectorizedArray::n_array_elements + v] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + strides[v]; + const unsigned int nvec = + this->dof_info + ->n_vectorization_lanes_filled[this->dof_access_index][this->cell]; + + if (fe_degree > 1 && evaluate_gradients == true) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 1 + side]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &this->data->face_to_cell_index_hermite(this->face_no, 0); + if (nvec == VectorizedArray::n_array_elements) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + unsigned int ind1[VectorizedArray::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + ind1[v] = indices[v] + (comp * static_dofs_per_component + + index_array[2 * i]) * + strides[v]; + unsigned int ind2[VectorizedArray::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + ind2[v] = indices[v] + (comp * static_dofs_per_component + + index_array[2 * i + 1]) * + strides[v]; + reader.process_dof_gather(ind1, + input_vector, + 0, + temp1[i + 2 * comp * dofs_per_face], + vector_selector); + VectorizedArray grad; + reader.process_dof_gather( + ind2, input_vector, 0, grad, vector_selector); + temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = + grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad); + } + else + { + for (unsigned int i = 0; i < n_components_ * 2 * dofs_per_face; + ++i) + temp1[i] = VectorizedArray(); + for (unsigned int v = 0; v < nvec; ++v) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i]) * + strides[v]; + const unsigned int ind2 = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i + 1]) * + strides[v]; + reader.process_dof( + ind1, + const_cast(input_vector), + temp1[i + 2 * comp * dofs_per_face][v]); + Number grad; + reader.process_dof(ind2, + const_cast(input_vector), + grad); + temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v] = + grad_weight[0] * + (temp1[i + 2 * comp * dofs_per_face][v] - grad); + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_nodal(this->face_no, 0); + if (nvec == VectorizedArray::n_array_elements) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + unsigned int ind[VectorizedArray::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + ind[v] = indices[v] + (comp * static_dofs_per_component + + index_array[i]) * + strides[v]; + reader.process_dof_gather(ind, + input_vector, + 0, + temp1[i + 2 * comp * dofs_per_face], + vector_selector); + } + else + { + for (unsigned int i = 0; i < n_components_ * dofs_per_face; ++i) + temp1[i] = VectorizedArray(); + for (unsigned int v = 0; v < nvec; ++v) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + indices[v] + + (comp * static_dofs_per_component + index_array[i]) * + strides[v]; + reader.process_dof( + ind1, + const_cast(input_vector), + temp1[i + 2 * comp * dofs_per_face][v]); + } + } + } + } + + // case 4: contiguous indices without interleaving + else if (((evaluate_gradients == false && + this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + 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->dof_access_index] + [this->cell] == + VectorizedArray::n_array_elements) { const unsigned int *indices = &this->dof_info @@ -7193,12 +7843,8 @@ FEFaceEvaluation:: // we know that the gradient weights for the Hermite case on the // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. - const VectorizedArray grad_weight0 = - (side ? -1. : 1.) * - this->data->shape_data_on_face[0][fe_degree + 1]; - const VectorizedArray grad_weight1 = - (side ? -1. : 1.) * - this->data->shape_data_on_face[0][fe_degree + 2]; + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 1 + side]; AssertDimension(this->data->face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -7217,10 +7863,7 @@ FEFaceEvaluation:: this->dof_info->component_dof_indices_offset [this->active_fe_index][this->first_selected_component], temp1[i + 2 * comp * dofs_per_face], - std::integral_constant< - bool, - std::is_same::value>()); + vector_selector); VectorizedArray grad; reader.process_dof_gather( indices, @@ -7229,13 +7872,9 @@ FEFaceEvaluation:: this->dof_info->component_dof_indices_offset [this->active_fe_index][this->first_selected_component], grad, - std::integral_constant< - bool, - std::is_same::value>()); + vector_selector); temp1[i + dofs_per_face + 2 * comp * dofs_per_face] = - grad_weight0 * temp1[i + 2 * comp * dofs_per_face] + - grad_weight1 * grad; + grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad); } } } @@ -7256,13 +7895,12 @@ FEFaceEvaluation:: this->dof_info->component_dof_indices_offset [this->active_fe_index][this->first_selected_component], temp1[i + comp * 2 * dofs_per_face], - std::integral_constant< - bool, - std::is_same::value>()); + vector_selector); } } } + + // case 5: default vector access else { this->read_dof_values(input_vector); @@ -7401,19 +8039,358 @@ FEFaceEvaluation:: # endif internal::VectorDistributorLocalToGlobal writer; + std::integral_constant< + bool, + std::is_same::value> + vector_selector; - 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->dof_access_index][this->cell] == - VectorizedArray::n_array_elements && - ((integrate_gradients == false && + // case 1: contiguous and interleaved indices + if (((integrate_gradients == false && this->data->nodal_at_cell_boundaries == true) || (this->data->element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && - fe_degree > 1))) + fe_degree > 1)) && + this->dof_info + ->index_storage_variants[this->dof_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous) + { + AssertDimension( + this->dof_info + ->n_vectorization_lanes_filled[this->dof_access_index][this->cell], + VectorizedArray::n_array_elements); + const unsigned int dof_index = + this->dof_info + ->dof_indices_contiguous[this->dof_access_index] + [this->cell * + VectorizedArray::n_array_elements] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements; + + if (fe_degree > 1 && integrate_gradients == true) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 2 - side]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_hermite(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = index_array[2 * i]; + const unsigned int ind2 = index_array[2 * i + 1]; + AssertIndexRange(ind1, dofs_per_cell); + AssertIndexRange(ind2, dofs_per_cell); + for (unsigned int comp = 0; comp < n_components_; ++comp) + { + VectorizedArray val = + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + VectorizedArray grad = + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + writer.process_dofs_vectorized( + 1, + dof_index + (ind1 + comp * static_dofs_per_component) * + VectorizedArray::n_array_elements, + destination, + &val, + vector_selector); + writer.process_dofs_vectorized( + 1, + dof_index + (ind2 + comp * static_dofs_per_component) * + VectorizedArray::n_array_elements, + destination, + &grad, + vector_selector); + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_nodal(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = index_array[i]; + for (unsigned int comp = 0; comp < n_components_; ++comp) + writer.process_dofs_vectorized( + 1, + dof_index + (ind + comp * static_dofs_per_component) * + VectorizedArray::n_array_elements, + destination, + temp1 + i + 2 * comp * dofs_per_face, + vector_selector); + } + } + } + + // case 2: contiguous and interleaved indices with fixed stride + else if (((integrate_gradients == false && + this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + this->dof_info + ->index_storage_variants[this->dof_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_strided) + { + AssertDimension( + this->dof_info + ->n_vectorization_lanes_filled[this->dof_access_index][this->cell], + VectorizedArray::n_array_elements); + const unsigned int *indices = + &this->dof_info + ->dof_indices_contiguous[this->dof_access_index] + [this->cell * + VectorizedArray::n_array_elements]; + if (fe_degree > 1 && integrate_gradients == true) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 2 - side]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &this->data->face_to_cell_index_hermite(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + index_array[2 * i] * VectorizedArray::n_array_elements; + const unsigned int ind2 = + index_array[2 * i + 1] * + VectorizedArray::n_array_elements; + for (unsigned int comp = 0; comp < n_components_; ++comp) + { + VectorizedArray val = + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + VectorizedArray grad = + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + writer.process_dof_gather( + indices, + destination, + ind1 + + comp * static_dofs_per_component * + VectorizedArray::n_array_elements + + this->dof_info->component_dof_indices_offset + [this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements, + val, + vector_selector); + writer.process_dof_gather( + indices, + destination, + ind2 + + comp * static_dofs_per_component * + VectorizedArray::n_array_elements + + this->dof_info->component_dof_indices_offset + [this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements, + grad, + vector_selector); + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_nodal(this->face_no, 0); + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind = + index_array[i] * VectorizedArray::n_array_elements; + for (unsigned int comp = 0; comp < n_components_; ++comp) + writer.process_dof_gather( + indices, + destination, + ind + + comp * static_dofs_per_component * + VectorizedArray::n_array_elements + + this->dof_info->component_dof_indices_offset + [this->active_fe_index] + [this->first_selected_component] * + VectorizedArray::n_array_elements, + temp1[i + 2 * comp * dofs_per_face], + vector_selector); + } + } + } + + // case 3: contiguous and interleaved indices with mixed stride + else if (((integrate_gradients == false && + this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + this->dof_info + ->index_storage_variants[this->dof_access_index][this->cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous_mixed_strides) + { + const unsigned int *strides = + &this->dof_info->dof_indices_interleave_strides + [this->dof_access_index] + [this->cell * VectorizedArray::n_array_elements]; + unsigned int indices[VectorizedArray::n_array_elements]; + for (unsigned int v = 0; v < VectorizedArray::n_array_elements; + ++v) + indices[v] = + this->dof_info->dof_indices_contiguous + [this->dof_access_index] + [this->cell * VectorizedArray::n_array_elements + v] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + strides[v]; + const unsigned int nvec = + this->dof_info + ->n_vectorization_lanes_filled[this->dof_access_index][this->cell]; + + if (fe_degree > 1 && integrate_gradients == true) + { + // we know that the gradient weights for the Hermite case on the + // right (side==1) are the negative from the value at the left + // (side==0), so we only read out one of them. + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 2 - side]; + AssertDimension(this->data->face_to_cell_index_hermite.size(1), + 2 * dofs_per_face); + + const unsigned int *index_array = + &this->data->face_to_cell_index_hermite(this->face_no, 0); + if (nvec == VectorizedArray::n_array_elements) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + unsigned int ind1[VectorizedArray::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + ind1[v] = indices[v] + (comp * static_dofs_per_component + + index_array[2 * i]) * + strides[v]; + unsigned int ind2[VectorizedArray::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + ind2[v] = indices[v] + (comp * static_dofs_per_component + + index_array[2 * i + 1]) * + strides[v]; + VectorizedArray val = + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + VectorizedArray grad = + grad_weight * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; + writer.process_dof_gather( + ind1, destination, 0, val, vector_selector); + writer.process_dof_gather( + ind2, destination, 0, grad, vector_selector); + } + else + { + for (unsigned int v = 0; v < nvec; ++v) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i]) * + strides[v]; + const unsigned int ind2 = + indices[v] + (comp * static_dofs_per_component + + index_array[2 * i + 1]) * + strides[v]; + Number val = + temp1[i + 2 * comp * dofs_per_face][v] - + grad_weight[0] * temp1[i + dofs_per_face + + 2 * comp * dofs_per_face][v]; + Number grad = + grad_weight[0] * + temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v]; + writer.process_dof(ind1, destination, val); + writer.process_dof(ind2, destination, grad); + } + } + } + else + { + AssertDimension(this->data->face_to_cell_index_nodal.size(1), + dofs_per_face); + const unsigned int *index_array = + &this->data->face_to_cell_index_nodal(this->face_no, 0); + if (nvec == VectorizedArray::n_array_elements) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + unsigned int ind[VectorizedArray::n_array_elements]; + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int v = 0; + v < VectorizedArray::n_array_elements; + ++v) + ind[v] = indices[v] + (comp * static_dofs_per_component + + index_array[i]) * + strides[v]; + writer.process_dof_gather(ind, + destination, + 0, + temp1[i + 2 * comp * dofs_per_face], + vector_selector); + } + else + { + for (unsigned int v = 0; v < nvec; ++v) + for (unsigned int comp = 0; comp < n_components_; ++comp) + for (unsigned int i = 0; i < dofs_per_face; ++i) + { + const unsigned int ind1 = + indices[v] + + (comp * static_dofs_per_component + index_array[i]) * + strides[v]; + writer.process_dof( + ind1, + destination, + temp1[i + 2 * comp * dofs_per_face][v]); + } + } + } + } + + // case 4: contiguous indices without interleaving + else if (((integrate_gradients == false && + this->data->nodal_at_cell_boundaries == true) || + (this->data->element_type == + internal::MatrixFreeFunctions::tensor_symmetric_hermite && + fe_degree > 1)) && + 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->dof_access_index] + [this->cell] == + VectorizedArray::n_array_elements) { const unsigned int *indices = &this->dof_info @@ -7428,12 +8405,8 @@ FEFaceEvaluation:: // we know that the gradient weights for the Hermite case on the // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. - const VectorizedArray grad_weight0 = - (side ? -1. : 1.) * - this->data->shape_data_on_face[0][fe_degree + 1]; - const VectorizedArray grad_weight1 = - (side ? -1. : 1.) * - this->data->shape_data_on_face[0][fe_degree + 2]; + const VectorizedArray grad_weight = + this->data->shape_data_on_face[0][fe_degree + 2 - side]; AssertDimension(this->data->face_to_cell_index_hermite.size(1), 2 * dofs_per_face); const unsigned int *index_array = @@ -7445,11 +8418,11 @@ FEFaceEvaluation:: for (unsigned int comp = 0; comp < n_components_; ++comp) { VectorizedArray val = - temp1[i + 2 * comp * dofs_per_face] + - grad_weight0 * + temp1[i + 2 * comp * dofs_per_face] - + grad_weight * temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; VectorizedArray grad = - grad_weight1 * + grad_weight * temp1[i + dofs_per_face + 2 * comp * dofs_per_face]; writer.process_dof_gather( indices, @@ -7458,10 +8431,7 @@ FEFaceEvaluation:: this->dof_info->component_dof_indices_offset [this->active_fe_index][this->first_selected_component], val, - std::integral_constant< - bool, - std::is_same::value>()); + vector_selector); writer.process_dof_gather( indices, destination, @@ -7469,10 +8439,7 @@ FEFaceEvaluation:: this->dof_info->component_dof_indices_offset [this->active_fe_index][this->first_selected_component], grad, - std::integral_constant< - bool, - std::is_same::value>()); + vector_selector); } } } @@ -7493,13 +8460,12 @@ FEFaceEvaluation:: this->dof_info->component_dof_indices_offset [this->active_fe_index][this->first_selected_component], temp1[i + 2 * comp * dofs_per_face], - std::integral_constant< - bool, - std::is_same::value>()); + vector_selector); } } } + + // case 5: default vector access else { internal::FEFaceNormalEvaluationImpl::initialize_indices( { const unsigned int p = face_info.faces[f].cells_exterior[v]; - const unsigned int stride = 1; - unsigned int i = 0; + const unsigned int stride = + dof_info[no].dof_indices_interleave_strides[2][p]; + unsigned int i = 0; for (unsigned int e = 0; e < dof_info[no].n_base_elements; ++e) @@ -1657,8 +1658,9 @@ MatrixFree::initialize_indices( { const unsigned int p = face_info.faces[f].cells_exterior[v]; - const unsigned int stride = 1; - unsigned int i = 0; + const unsigned int stride = + dof_info[no].dof_indices_interleave_strides[2][p]; + unsigned int i = 0; for (unsigned int e = 0; e < dof_info[no].n_base_elements; ++e) diff --git a/source/matrix_free/evaluation_selector.inst.in b/source/matrix_free/evaluation_selector.inst.in index 5be002385d..e456056bb4 100644 --- a/source/matrix_free/evaluation_selector.inst.in +++ b/source/matrix_free/evaluation_selector.inst.in @@ -29,6 +29,7 @@ for (deal_II_dimension : DIMENSIONS; components : SPACE_DIMENSIONS; VectorizedArray *, VectorizedArray *, const bool, + const bool, const bool); template void SelectEvaluator