From: Martin Kronbichler Date: Mon, 18 Sep 2023 10:28:45 +0000 (+0200) Subject: Matrix-free evaluation: Avoid pointer aliasing X-Git-Tag: relicensing~492^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c399ab5dac1a12735054cfd506be559322d3606f;p=dealii.git Matrix-free evaluation: Avoid pointer aliasing --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 9c3bcfdb33..4bbcf8de5d 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3403,6 +3403,16 @@ FEEvaluationBase:: ->component_dof_indices_offset[this->active_fe_index] [this->first_selected_component] * n_lanes; + + std::array src_ptrs; + if (n_components == 1 || this->n_fe_components == 1) + for (unsigned int comp = 0; comp < n_components; ++comp) + src_ptrs[comp] = + const_cast(src[comp]->begin()); + else + src_ptrs[0] = + const_cast(src[0]->begin()); + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int i = 0; i < dofs_per_component; ++i, dof_indices += n_lanes) @@ -3410,14 +3420,19 @@ FEEvaluationBase:: operation.process_dof_gather(dof_indices, *src[comp], 0, + src_ptrs[comp], values_dofs[comp][i], vector_selector); else for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_component; ++i, dof_indices += n_lanes) - operation.process_dof_gather( - dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector); + operation.process_dof_gather(dof_indices, + *src[0], + 0, + src_ptrs[0], + values_dofs[comp][i], + vector_selector); return; } @@ -3962,6 +3977,15 @@ FEEvaluationBase:: internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous_strided) { + std::array src_ptrs; + if (n_components == 1 || this->n_fe_components == 1) + for (unsigned int comp = 0; comp < n_components; ++comp) + src_ptrs[comp] = const_cast( + src[comp]->begin()); + else + src_ptrs[0] = + const_cast(src[0]->begin()); + if (n_components == 1 || this->n_fe_components == 1) for (unsigned int i = 0; i < dofs_per_component; ++i) { @@ -3969,6 +3993,7 @@ FEEvaluationBase:: operation.process_dof_gather(dof_indices.data(), *src[comp], i * n_lanes, + src_ptrs[comp] + i * n_lanes, values_dofs[comp][i], vector_selector); } @@ -3976,12 +4001,13 @@ FEEvaluationBase:: for (unsigned int comp = 0; comp < n_components; ++comp) for (unsigned int i = 0; i < dofs_per_component; ++i) { - operation.process_dof_gather(dof_indices.data(), - *src[0], - (comp * dofs_per_component + i) * - n_lanes, - values_dofs[comp][i], - vector_selector); + operation.process_dof_gather( + dof_indices.data(), + *src[0], + (comp * dofs_per_component + i) * n_lanes, + src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes, + values_dofs[comp][i], + vector_selector); } } else @@ -3990,6 +4016,15 @@ FEEvaluationBase:: internal::MatrixFreeFunctions::DoFInfo:: IndexStorageVariants::interleaved_contiguous_mixed_strides, ExcNotImplemented()); + std::array src_ptrs; + if (n_components == 1 || this->n_fe_components == 1) + for (unsigned int comp = 0; comp < n_components; ++comp) + src_ptrs[comp] = const_cast( + src[comp]->begin()); + else + src_ptrs[0] = + const_cast(src[0]->begin()); + const unsigned int *offsets = &dof_info.dof_indices_interleave_strides[ind][n_lanes * this->cell]; if (n_components == 1 || this->n_fe_components == 1) @@ -3999,6 +4034,7 @@ FEEvaluationBase:: operation.process_dof_gather(dof_indices.data(), *src[comp], 0, + src_ptrs[comp], values_dofs[comp][i], vector_selector); DEAL_II_OPENMP_SIMD_PRAGMA @@ -4012,6 +4048,7 @@ FEEvaluationBase:: operation.process_dof_gather(dof_indices.data(), *src[0], 0, + src_ptrs[0], values_dofs[comp][i], vector_selector); DEAL_II_OPENMP_SIMD_PRAGMA diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h index 9b5535d209..1504b5b3d8 100644 --- a/include/deal.II/matrix_free/vector_access_internal.h +++ b/include/deal.II/matrix_free/vector_access_internal.h @@ -439,22 +439,28 @@ namespace internal // gather template void - process_dof_gather(const unsigned int *indices, - VectorType &vec, - const unsigned int constant_offset, - VectorizedArrayType &res, + process_dof_gather(const unsigned int *indices, + VectorType &vec, + const unsigned int constant_offset, + typename VectorType::value_type *vec_ptr, + VectorizedArrayType &res, std::integral_constant) const { + (void)constant_offset; + (void)vec; + #ifdef DEBUG // in debug mode, run non-vectorized version because this path // has additional checks (e.g., regarding ghosting) + Assert(vec_ptr == vec.begin() + constant_offset, ExcInternalError()); process_dof_gather(indices, vec, constant_offset, + vec_ptr, res, std::integral_constant()); #else - res.gather(vec.begin() + constant_offset, indices); + res.gather(vec_ptr, indices); #endif } @@ -464,9 +470,10 @@ namespace internal // manually load the data template void - process_dof_gather(const unsigned int *indices, - const VectorType &vec, - const unsigned int constant_offset, + process_dof_gather(const unsigned int *indices, + const VectorType &vec, + const unsigned int constant_offset, + typename VectorType::value_type *, VectorizedArrayType &res, std::integral_constant) const { @@ -683,21 +690,26 @@ namespace internal // scatter template void - process_dof_gather(const unsigned int *indices, - VectorType &vec, - const unsigned int constant_offset, - VectorizedArrayType &res, + process_dof_gather(const unsigned int *indices, + VectorType &vec, + const unsigned int constant_offset, + typename VectorType::value_type *vec_ptr, + VectorizedArrayType &res, std::integral_constant) const { + (void)constant_offset; + (void)vec_ptr; + (void)vec; + #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS < 512 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) vector_access(vec, indices[v] + constant_offset) += res[v]; #else // only use gather in case there is also scatter. VectorizedArrayType tmp; - tmp.gather(vec.begin() + constant_offset, indices); + tmp.gather(vec_ptr, indices); tmp += res; - tmp.scatter(indices, vec.begin() + constant_offset); + tmp.scatter(indices, vec_ptr); #endif } @@ -707,9 +719,10 @@ namespace internal // manually append all data template void - process_dof_gather(const unsigned int *indices, - VectorType &vec, - const unsigned int constant_offset, + process_dof_gather(const unsigned int *indices, + VectorType &vec, + const unsigned int constant_offset, + typename VectorType::value_type *, VectorizedArrayType &res, std::integral_constant) const { @@ -913,22 +926,28 @@ namespace internal template void - process_dof_gather(const unsigned int *indices, - VectorType &vec, - const unsigned int constant_offset, - VectorizedArrayType &res, + process_dof_gather(const unsigned int *indices, + VectorType &vec, + const unsigned int constant_offset, + typename VectorType::value_type *vec_ptr, + VectorizedArrayType &res, std::integral_constant) const { - res.scatter(indices, vec.begin() + constant_offset); + (void)vec_ptr; + (void)constant_offset; + (void)vec; + Assert(vec_ptr == vec.begin() + constant_offset, ExcInternalError()); + res.scatter(indices, vec_ptr); } template void - process_dof_gather(const unsigned int *indices, - VectorType &vec, - const unsigned int constant_offset, + process_dof_gather(const unsigned int *indices, + VectorType &vec, + const unsigned int constant_offset, + typename VectorType::value_type *, VectorizedArrayType &res, std::integral_constant) const {