From b40ea0239d713ca5423163ce213881902ddedc0e Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 19 Sep 2019 12:58:28 +0200
Subject: [PATCH] Speed up IndexSet::add_indices()

---
 include/deal.II/base/index_set.h | 68 ++++++++++++++++++++++++++++++--
 source/base/index_set.cc         | 29 --------------
 2 files changed, 65 insertions(+), 32 deletions(-)

diff --git a/include/deal.II/base/index_set.h b/include/deal.II/base/index_set.h
index 883dd1a997..f4e6dfe24a 100644
--- a/include/deal.II/base/index_set.h
+++ b/include/deal.II/base/index_set.h
@@ -1636,12 +1636,49 @@ IndexSet::add_index(const size_type index)
 
 
 
+inline void
+IndexSet::add_range(const size_type begin, const size_type end)
+{
+  Assert((begin < index_space_size) ||
+           ((begin == index_space_size) && (end == index_space_size)),
+         ExcIndexRangeType<size_type>(begin, 0, index_space_size));
+  Assert(end <= index_space_size,
+         ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
+  Assert(begin <= end, ExcIndexRangeType<size_type>(begin, 0, end));
+
+  if (begin != end)
+    {
+      const Range new_range(begin, end);
+
+      // the new index might be larger than the last index present in the
+      // ranges. Then we can skip the binary search
+      if (ranges.size() == 0 || begin > ranges.back().end)
+        ranges.push_back(new_range);
+      else
+        ranges.insert(Utilities::lower_bound(ranges.begin(),
+                                             ranges.end(),
+                                             new_range),
+                      new_range);
+      is_compressed = false;
+    }
+}
+
+
+
 template <typename ForwardIterator>
 inline void
 IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
 {
-  // insert each element of the range. if some of them happen to be
-  // consecutive, merge them to a range
+  if (begin == end)
+    return;
+
+  // identify ranges in the given iterator range by checking whether some
+  // indices happen to be consecutive. to avoid quadratic complexity when
+  // calling add_range many times (as add_range() going into the middle of an
+  // already existing range must shift entries around), we first collect a
+  // vector of ranges.
+  std::vector<std::pair<size_type, size_type>> tmp_ranges;
+  bool                                         ranges_are_sorted = true;
   for (ForwardIterator p = begin; p != end;)
     {
       const size_type begin_index = *p;
@@ -1654,9 +1691,34 @@ IndexSet::add_indices(const ForwardIterator &begin, const ForwardIterator &end)
           ++q;
         }
 
-      add_range(begin_index, end_index);
+      tmp_ranges.emplace_back(begin_index, end_index);
       p = q;
+      if (p != end && *p < end_index)
+        ranges_are_sorted = false;
     }
+
+  if (!ranges_are_sorted)
+    std::sort(tmp_ranges.begin(), tmp_ranges.end());
+
+  // if we have many ranges, we first construct a temporary index set (where
+  // we add ranges in a consecutive way, so fast), otherwise, we work with
+  // add_range(). the number 9 is chosen heuristically given the fact that
+  // there are typically up to 8 independent ranges when adding the degrees of
+  // freedom on a 3D cell or 9 when adding degrees of freedom of faces. if
+  // doing cell-by-cell additions, we want to avoid repeated calls to
+  // IndexSet::compress() which gets called upon merging two index sets, so we
+  // want to be in the other branch then.
+  if (tmp_ranges.size() > 9)
+    {
+      IndexSet tmp_set(size());
+      tmp_set.ranges.reserve(tmp_ranges.size());
+      for (const auto &i : tmp_ranges)
+        tmp_set.add_range(i.first, i.second);
+      this->add_indices(tmp_set);
+    }
+  else
+    for (const auto &i : tmp_ranges)
+      add_range(i.first, i.second);
 }
 
 
diff --git a/source/base/index_set.cc b/source/base/index_set.cc
index 273a46c9b3..ef7ff15d95 100644
--- a/source/base/index_set.cc
+++ b/source/base/index_set.cc
@@ -94,35 +94,6 @@ IndexSet::IndexSet(const Epetra_BlockMap &map)
 
 
 
-void
-IndexSet::add_range(const size_type begin, const size_type end)
-{
-  Assert((begin < index_space_size) ||
-           ((begin == index_space_size) && (end == index_space_size)),
-         ExcIndexRangeType<size_type>(begin, 0, index_space_size));
-  Assert(end <= index_space_size,
-         ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
-  Assert(begin <= end, ExcIndexRangeType<size_type>(begin, 0, end));
-
-  if (begin != end)
-    {
-      const Range new_range(begin, end);
-
-      // the new index might be larger than the last index present in the
-      // ranges. Then we can skip the binary search
-      if (ranges.size() == 0 || begin > ranges.back().end)
-        ranges.push_back(new_range);
-      else
-        ranges.insert(Utilities::lower_bound(ranges.begin(),
-                                             ranges.end(),
-                                             new_range),
-                      new_range);
-      is_compressed = false;
-    }
-}
-
-
-
 void
 IndexSet::do_compress() const
 {
-- 
2.39.5