From: Martin Kronbichler Date: Fri, 8 Jul 2022 06:35:38 +0000 (+0200) Subject: MatrixFree DoFInfo: Break out from some expensive loops X-Git-Tag: v9.5.0-rc1~1075^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=90c9ee642e4d14904227f01460f24d47a02e1ff0;p=dealii.git MatrixFree DoFInfo: Break out from some expensive loops --- diff --git a/source/matrix_free/dof_info.cc b/source/matrix_free/dof_info.cc index 482614408e..8cf256b08a 100644 --- a/source/matrix_free/dof_info.cc +++ b/source/matrix_free/dof_info.cc @@ -612,7 +612,9 @@ namespace internal const unsigned int *dof_indices = this->dof_indices.data() + row_starts[i * vectorization_length * n_components].first; - for (unsigned int k = 0; k < ndofs; ++k) + for (unsigned int k = 0; + k < ndofs && indices_are_interleaved_and_contiguous; + ++k) for (unsigned int j = 0; j < n_comp; ++j) if (dof_indices[j * ndofs + k] != dof_indices[0] + k * n_comp + j) @@ -663,7 +665,9 @@ namespace internal 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 k = 0; + k < ndofs && indices_are_interleaved_and_mixed != 0; + ++k) for (unsigned int j = 0; j < n_comp; ++j) // the first if case is to avoid negative offsets // (invalid) @@ -712,19 +716,36 @@ namespace internal 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 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; - } + // do not use interleaved storage if two entries within + // vectorized array point to the same index (scatter not + // possible due to race condition) + for (const unsigned int *indices = dof_indices; + indices != dof_indices + ndofs; + ++indices) + { + bool is_sorted = true; + unsigned int previous = indices[0]; + for (unsigned int l = 1; l < n_comp; ++l) + { + const unsigned int current = indices[l * ndofs]; + if (current <= previous) + is_sorted = false; + + // the simple check failed, must compare all + // indices manually - due to short sizes this + // O(n^2) algorithm is better than sorting + if (!is_sorted) + for (unsigned int j = 0; j < l; ++j) + if (indices[j * ndofs] == current) + { + index_storage_variants + [dof_access_cell][i] = + IndexStorageVariants::full; + break; + } + previous = current; + } + } } } } @@ -848,9 +869,14 @@ namespace internal ndofs * vectorization_length, this->dof_indices_interleaved.size() + 1); 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]; + { + const unsigned int *my_dof_indices = dof_indices + k; + const unsigned int *end = + interleaved_dof_indices + vectorization_length; + for (; interleaved_dof_indices != end; + ++interleaved_dof_indices, my_dof_indices += ndofs) + *interleaved_dof_indices = *my_dof_indices; + } } }