From 90c9ee642e4d14904227f01460f24d47a02e1ff0 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@uni-a.de>
Date: Fri, 8 Jul 2022 08:35:38 +0200
Subject: [PATCH] MatrixFree DoFInfo: Break out from some expensive loops

---
 source/matrix_free/dof_info.cc | 62 ++++++++++++++++++++++++----------
 1 file changed, 44 insertions(+), 18 deletions(-)

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;
+              }
           }
     }
 
-- 
2.39.5