From 65cb5e67fcefa1b99a9172e4edd56e98013f6b5f Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 17 Jun 2016 14:58:10 -0700 Subject: [PATCH] fix DynamicSparsityPattern::begin(row) performance, add IndexSet::at begin(row) and end(row) would search through all empty rows leading to quadratic complexity when iterating row by row as it is done in the construction of PETSc's SparseMatrix. Fix this by finding the closest locally owned row using a binary search and then only iterate over the locally owned set. --- doc/news/changes.h | 13 +++ include/deal.II/base/index_set.h | 65 ++++++++++++++ .../deal.II/lac/dynamic_sparsity_pattern.h | 83 ++++++++++++------ source/lac/dynamic_sparsity_pattern.cc | 39 ++++++--- tests/base/index_set_iterate_02.cc | 73 ++++++++++++++++ tests/base/index_set_iterate_02.output | 11 +++ .../dynamic_sparsity_pattern_iterator_05.cc | 86 +++++++++++++++++++ ...ynamic_sparsity_pattern_iterator_05.output | 7 ++ .../dynamic_sparsity_pattern_iterator_06.cc | 83 ++++++++++++++++++ ...ynamic_sparsity_pattern_iterator_06.output | 3 + 10 files changed, 427 insertions(+), 36 deletions(-) create mode 100644 tests/base/index_set_iterate_02.cc create mode 100644 tests/base/index_set_iterate_02.output create mode 100644 tests/bits/dynamic_sparsity_pattern_iterator_05.cc create mode 100644 tests/bits/dynamic_sparsity_pattern_iterator_05.output create mode 100644 tests/bits/dynamic_sparsity_pattern_iterator_06.cc create mode 100644 tests/bits/dynamic_sparsity_pattern_iterator_06.output diff --git a/doc/news/changes.h b/doc/news/changes.h index dc50fed535..13a304311f 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -288,6 +288,19 @@ inconvenience this causes.
    +
  1. New: IndexSet::at(idx) returns an iterator pointing to the given index + or the next larger element in the set if idx is not contained. +
    + (Timo Heister, 2016/06/21) +
    2. + +
  2. Fixed: Performance of DynamicSparsityPattern::begin(r) and + DynamicSparsityPattern::end(r) has been improved dramatically in parallel + computations and if the pattern is empty. +
    + (Timo Heister, 2016/06/21) +
    3. +
  3. Fixed: FEFieldFunction now works correctly in distributed computations, where before exceptions of type ExcPointNotAvailableHere could occur for evaluation points on or close to a boundary to a ghost cell. diff --git a/include/deal.II/base/index_set.h b/include/deal.II/base/index_set.h index 01dc388ce7..6812433b0e 100644 --- a/include/deal.II/base/index_set.h +++ b/include/deal.II/base/index_set.h @@ -665,6 +665,22 @@ public: */ ElementIterator begin() const; + /** + * Return an element iterator pointing to the element with global index + * @p global_index or the next larger element if the index is not in the + * set. This is equivalent to + * @code + * auto p = begin(); + * while (*p (global_index, 0, size())); + + if (ranges.empty()) + return end(); + + std::vector::const_iterator main_range=ranges.begin()+largest_range; + + Range r (global_index, global_index+1); + // This optimization makes the bounds for lower_bound smaller by checking + // the largest range first. + std::vector::const_iterator range_begin, range_end; + if (global_indexbegin) + { + range_begin = ranges.begin(); + range_end = main_range; + } + else + { + range_begin = main_range; + range_end = ranges.end(); + } + + // This will give us the first range p=[a,b[ with b>=global_index using + // a binary search + const std::vector::const_iterator + p = Utilities::lower_bound(range_begin, range_end, r, Range::end_compare); + + // We couldn't find a range, which means we have no range that contains + // global_index and also no range behind it, meaning we need to return end(). + if (p == ranges.end()) + return end(); + + // Finally, we can have two cases: Either global_index is not in [a,b[, + // which means we need to return an iterator to a because global_index, ..., + // a-1 is not in the IndexSet (if branch). Alternatively, global_index is in + // [a,b[ and we will return an iterator pointing directly at global_index + // (else branch). + if (global_index < p->begin) + return IndexSet::ElementIterator(this, p-ranges.begin(), p->begin); + else + return IndexSet::ElementIterator(this, p-ranges.begin(), global_index); +} + + inline IndexSet::ElementIterator IndexSet::end() const { diff --git a/include/deal.II/lac/dynamic_sparsity_pattern.h b/include/deal.II/lac/dynamic_sparsity_pattern.h index 86882de5b7..c2e35819af 100644 --- a/include/deal.II/lac/dynamic_sparsity_pattern.h +++ b/include/deal.II/lac/dynamic_sparsity_pattern.h @@ -559,6 +559,11 @@ public: size_type memory_consumption () const; private: + /** + * A flag that stores whether any entries have been added so far. + */ + bool have_entries; + /** * Number of rows that this sparsity structure shall represent. */ @@ -947,6 +952,8 @@ DynamicSparsityPattern::add (const size_type i, if (rowset.size() > 0 && !rowset.is_element(i)) return; + have_entries = true; + const size_type rowindex = rowset.size()==0 ? i : rowset.index_within_set(i); lines[rowindex].add (j); @@ -967,6 +974,9 @@ DynamicSparsityPattern::add_entries (const size_type row, if (rowset.size() > 0 && !rowset.is_element(row)) return; + if (!have_entries && begin (row, 0, n_rows())); + + if (!have_entries) + return 0; + if (rowset.size() > 0 && !rowset.is_element(row)) return 0; @@ -1034,20 +1048,50 @@ DynamicSparsityPattern::begin (const size_type r) const { Assert (r(r,0,n_rows())); - // find the first row starting at r that has entries and return the - // begin iterator to it. also skip rows for which we do not have - // store anything based on the IndexSet given to the sparsity - // pattern - // - // note: row_length(row) returns zero if the row is not locally stored - // - // TODO: this is way too slow when used in parallel, so do not use it on - // non-owned rows + if (!have_entries) + return iterator(this); + + if (rowset.size() > 0) + { + // We have an IndexSet that describes the locally owned set. For + // performance reasons we need to make sure that we don't do a + // linear search over 0..n_rows(). Instead, find the first entry + // >= row r in the locally owned set (this is done in log + // n_ranges time inside at()). From there, we move forward until + // we find a non-empty row. By iterating over the IndexSet instead + // of incrementing the row index, we potentially skip over entries + // not in the rowset. + IndexSet::ElementIterator it = rowset.at(r); + if (it == rowset.end()) + return end(); // we don't own any row between r and the end + + // Instead of using row_length(*it)==0 in the while loop below, + // which involves an expensive index_within_set() call, we + // look at the lines vector directly. This works, because we are + // walking over this vector entry by entry anyways. + size_type rowindex = rowset.index_within_set(*it); + + while (it!=rowset.end() + && lines[rowindex].entries.size()==0) + { + ++it; + ++rowindex; + } + + if (it == rowset.end()) + return end(); + else + return iterator(this, *it, 0); + } + + // Without an index set we have to do a linear search starting at + // row r until we find a non-empty one size_type row = r; - while ((row(r,0,n_rows())); - // find the first row after r that has entries and return the begin - // iterator to it. also skip rows for which we do not have - // store anything based on the IndexSet given to the sparsity - // pattern - // - // note: row_length(row) returns zero if the row is not locally stored unsigned int row = r+1; - while ((row +#include +#include +#include + +#include + +void test (IndexSet &index_set, unsigned int n) +{ + deallog << "n=" << n; + + IndexSet::ElementIterator it = index_set.at(n); + + deallog << " end?" << (it == index_set.end()); + if (it!=index_set.end()) + { + deallog << " value=" << *it; + } + + deallog << std::endl; +} + +void test() +{ + IndexSet index_set (20); + index_set.add_range (2,5); + index_set.add_index (9); + + index_set.print(deallog); + + test(index_set, 0); + test(index_set, 2); + test(index_set, 3); + test(index_set, 4); + test(index_set, 7); + test(index_set, 9); + test(index_set, 15); + + IndexSet empty (42); + empty.print(deallog); + test(empty, 6); +} + + + + +int main() +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test (); +} diff --git a/tests/base/index_set_iterate_02.output b/tests/base/index_set_iterate_02.output new file mode 100644 index 0000000000..8c2ec743e4 --- /dev/null +++ b/tests/base/index_set_iterate_02.output @@ -0,0 +1,11 @@ + +DEAL::{[2,4], 9} +DEAL::n=0 end?0 value=2 +DEAL::n=2 end?0 value=2 +DEAL::n=3 end?0 value=3 +DEAL::n=4 end?0 value=4 +DEAL::n=7 end?0 value=9 +DEAL::n=9 end?0 value=9 +DEAL::n=15 end?1 +DEAL::{} +DEAL::n=6 end?1 diff --git a/tests/bits/dynamic_sparsity_pattern_iterator_05.cc b/tests/bits/dynamic_sparsity_pattern_iterator_05.cc new file mode 100644 index 0000000000..407b2db73b --- /dev/null +++ b/tests/bits/dynamic_sparsity_pattern_iterator_05.cc @@ -0,0 +1,86 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// investigate performance issues in DynamicSparsityPattern::begin(r) and +// end(r) for large sets with many empty rows. + +#include "../tests.h" +#include +#include +#include + + +void test (bool empty, bool large_gap) +{ + const int size = 100000000; + const int my_start = size/3; + IndexSet owned(size); + owned.add_range(my_start, my_start+5); + if (large_gap) + owned.add_range(size-1, size); + DynamicSparsityPattern sp (size, 5, owned); + if (!empty) + sp.add(my_start+1, 1); + + for (unsigned int i=my_start-10; irow() << ' ' << p->column() << std::endl; + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + try + { + test (false, false); + test (true, false); + test (false, true); + test (true, true); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/dynamic_sparsity_pattern_iterator_05.output b/tests/bits/dynamic_sparsity_pattern_iterator_05.output new file mode 100644 index 0000000000..0f16dbbfca --- /dev/null +++ b/tests/bits/dynamic_sparsity_pattern_iterator_05.output @@ -0,0 +1,7 @@ + +DEAL::33333334 1 +DEAL::OK +DEAL::OK +DEAL::33333334 1 +DEAL::OK +DEAL::OK diff --git a/tests/bits/dynamic_sparsity_pattern_iterator_06.cc b/tests/bits/dynamic_sparsity_pattern_iterator_06.cc new file mode 100644 index 0000000000..98d9af2769 --- /dev/null +++ b/tests/bits/dynamic_sparsity_pattern_iterator_06.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// investigate performance for DynamicSparsityPattern::begin(r) if the +// SP is empty + +#include "../tests.h" +#include +#include +#include + + +void test (bool have_set) +{ + const int size = 100000000; + const int my_start = size/3; + + IndexSet empty_set; + IndexSet owned(size); + owned.add_range(my_start, my_start+5); + + DynamicSparsityPattern sp (size, 5, have_set?owned:empty_set); + + for (unsigned int i=my_start; irow() << ' ' << p->column() << std::endl; + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + try + { + test (false); + test (true); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/dynamic_sparsity_pattern_iterator_06.output b/tests/bits/dynamic_sparsity_pattern_iterator_06.output new file mode 100644 index 0000000000..8b3b075900 --- /dev/null +++ b/tests/bits/dynamic_sparsity_pattern_iterator_06.output @@ -0,0 +1,3 @@ + +DEAL::OK +DEAL::OK -- 2.39.5