From b40ea0239d713ca5423163ce213881902ddedc0e Mon Sep 17 00:00:00 2001 From: Martin Kronbichler 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(begin, 0, index_space_size)); + Assert(end <= index_space_size, + ExcIndexRangeType(end, 0, index_space_size + 1)); + Assert(begin <= end, ExcIndexRangeType(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 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> 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(begin, 0, index_space_size)); - Assert(end <= index_space_size, - ExcIndexRangeType(end, 0, index_space_size + 1)); - Assert(begin <= end, ExcIndexRangeType(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