]> https://gitweb.dealii.org/ - dealii.git/commitdiff
fix DynamicSparsityPattern::begin(row) performance, add IndexSet::at
authorTimo Heister <timo.heister@gmail.com>
Fri, 17 Jun 2016 21:58:10 +0000 (14:58 -0700)
committerTimo Heister <timo.heister@gmail.com>
Tue, 21 Jun 2016 18:30:27 +0000 (14:30 -0400)
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
include/deal.II/base/index_set.h
include/deal.II/lac/dynamic_sparsity_pattern.h
source/lac/dynamic_sparsity_pattern.cc
tests/base/index_set_iterate_02.cc [new file with mode: 0644]
tests/base/index_set_iterate_02.output [new file with mode: 0644]
tests/bits/dynamic_sparsity_pattern_iterator_05.cc [new file with mode: 0644]
tests/bits/dynamic_sparsity_pattern_iterator_05.output [new file with mode: 0644]
tests/bits/dynamic_sparsity_pattern_iterator_06.cc [new file with mode: 0644]
tests/bits/dynamic_sparsity_pattern_iterator_06.output [new file with mode: 0644]

index dc50fed535d70556bb9ca095f1d290adb57e878f..13a304311fe22c54cf924441a389101eaf91cf9b 100644 (file)
@@ -288,6 +288,19 @@ inconvenience this causes.
 
 <ol>
 
+ <li> 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.
+ <br>
+ (Timo Heister, 2016/06/21)
+ </li>
+
+ <li> Fixed: Performance of DynamicSparsityPattern::begin(r) and
+ DynamicSparsityPattern::end(r) has been improved dramatically in parallel
+ computations and if the pattern is empty.
+ <br>
+ (Timo Heister, 2016/06/21)
+ </li>
+
  <li> 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.
index 01dc388ce75db80c2338cebef9ddf07139898999..6812433b0e2d0619c756435ab6f917cd97b478c4 100644 (file)
@@ -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)
+   *   ++p;
+   * return p;
+   * @endcode
+   *
+   * If there is no element in this IndexSet at or behind @p global_index,
+   * this method will return end().
+   */
+  ElementIterator at(const size_type global_index) const;
+
   /**
    * Return an iterator that points one after the last index that is contained
    * in this IndexSet.
@@ -1194,6 +1210,55 @@ IndexSet::ElementIterator IndexSet::begin() const
     return end();
 }
 
+inline
+IndexSet::ElementIterator IndexSet::at(const size_type global_index) const
+{
+  compress();
+  Assert (global_index < size(),
+          ExcIndexRangeType<size_type> (global_index, 0, size()));
+
+  if (ranges.empty())
+    return end();
+
+  std::vector<Range>::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<Range>::const_iterator range_begin, range_end;
+  if (global_index<main_range->begin)
+    {
+      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<Range>::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
 {
index 86882de5b723c3ea60820b17d92ab3bd3bde9ef9..c2e35819af685cbdd81b40a9cfb017cd55c8498f 100644 (file)
@@ -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<end)
+    have_entries = true;
+
   const size_type rowindex =
     rowset.size()==0 ? row : rowset.index_within_set(row);
   lines[rowindex].add_entries (begin, end, indices_are_sorted);
@@ -985,6 +995,10 @@ types::global_dof_index
 DynamicSparsityPattern::row_length (const size_type row) const
 {
   Assert (row < n_rows(), ExcIndexRangeType<size_type> (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<n_rows(), ExcIndexRangeType<size_type>(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<n_rows())
-         &&
-         (row_length(row)==0))
-    ++row;
+
+  while (row<n_rows() && row_length(row)==0)
+    {
+      ++row;
+    }
 
   if (row == n_rows())
     return iterator(this);
@@ -1063,22 +1107,11 @@ DynamicSparsityPattern::end (const size_type r) const
 {
   Assert (r<n_rows(), ExcIndexRangeType<size_type>(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<n_rows())
-         &&
-         (row_length(row)==0))
-    ++row;
-
   if (row == n_rows())
     return iterator(this);
   else
-    return iterator(this, row, 0);
+    return begin(row);
 }
 
 
index 94c3be446f0eaf50dedf480e5609eb5a12d93b81..d097f52f07c5fae05fbffe79ff9b95001a57f2a0 100644 (file)
@@ -216,9 +216,10 @@ DynamicSparsityPattern::Line::memory_consumption () const
 
 DynamicSparsityPattern::DynamicSparsityPattern ()
   :
-  rows(0),
-  cols(0),
-  rowset(0)
+  have_entries (false),
+  rows (0),
+  cols (0),
+  rowset (0)
 {}
 
 
@@ -227,9 +228,10 @@ DynamicSparsityPattern::
 DynamicSparsityPattern (const DynamicSparsityPattern &s)
   :
   Subscriptor(),
-  rows(0),
-  cols(0),
-  rowset(0)
+  have_entries (false),
+  rows (0),
+  cols (0),
+  rowset (0)
 {
   (void)s;
   Assert (s.rows == 0, ExcInvalidConstructorCall());
@@ -243,9 +245,10 @@ DynamicSparsityPattern::DynamicSparsityPattern (const size_type m,
                                                 const IndexSet &rowset_
                                                )
   :
-  rows(0),
-  cols(0),
-  rowset(0)
+  have_entries (false),
+  rows (0),
+  cols (0),
+  rowset (0)
 {
   reinit (m,n, rowset_);
 }
@@ -253,6 +256,7 @@ DynamicSparsityPattern::DynamicSparsityPattern (const size_type m,
 
 DynamicSparsityPattern::DynamicSparsityPattern (const IndexSet &rowset_)
   :
+  have_entries (false),
   rows(0),
   cols(0),
   rowset(0)
@@ -263,6 +267,7 @@ DynamicSparsityPattern::DynamicSparsityPattern (const IndexSet &rowset_)
 
 DynamicSparsityPattern::DynamicSparsityPattern (const size_type n)
   :
+  have_entries (false),
   rows(0),
   cols(0),
   rowset(0)
@@ -292,6 +297,7 @@ DynamicSparsityPattern::reinit (const size_type m,
                                 const size_type n,
                                 const IndexSet &rowset_)
 {
+  have_entries = false;
   rows = m;
   cols = n;
   rowset=rowset_;
@@ -321,6 +327,9 @@ DynamicSparsityPattern::empty () const
 DynamicSparsityPattern::size_type
 DynamicSparsityPattern::max_entries_per_row () const
 {
+  if (!have_entries)
+    return 0;
+
   size_type m = 0;
   for (size_type i=0; i<lines.size(); ++i)
     {
@@ -340,6 +349,9 @@ DynamicSparsityPattern::exists (const size_type i,
   Assert (j<cols, ExcIndexRange(j, 0, cols));
   Assert( rowset.size()==0 || rowset.is_element(i), ExcInternalError());
 
+  if (!have_entries)
+    return false;
+
   const size_type rowindex =
     rowset.size()==0 ? i : rowset.index_within_set(i);
 
@@ -455,6 +467,9 @@ DynamicSparsityPattern::bandwidth () const
 DynamicSparsityPattern::size_type
 DynamicSparsityPattern::n_nonzero_elements () const
 {
+  if (!have_entries)
+    return 0;
+
   size_type n=0;
   for (size_type i=0; i<lines.size(); ++i)
     {
@@ -468,8 +483,10 @@ DynamicSparsityPattern::n_nonzero_elements () const
 DynamicSparsityPattern::size_type
 DynamicSparsityPattern::memory_consumption () const
 {
-  //TODO: IndexSet...
-  size_type mem = sizeof(DynamicSparsityPattern);
+  size_type mem = sizeof(DynamicSparsityPattern)
+                  + MemoryConsumption::memory_consumption(rowset)
+                  - sizeof(rowset);
+
   for (size_type i=0; i<lines.size(); ++i)
     mem += MemoryConsumption::memory_consumption (lines[i]);
 
diff --git a/tests/base/index_set_iterate_02.cc b/tests/base/index_set_iterate_02.cc
new file mode 100644 (file)
index 0000000..f8d6eb0
--- /dev/null
@@ -0,0 +1,73 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test IndexSet::at()
+
+#include "../tests.h"
+#include <iomanip>
+#include <iomanip>
+#include <fstream>
+#include <cmath>
+
+#include <deal.II/base/index_set.h>
+
+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 (file)
index 0000000..8c2ec74
--- /dev/null
@@ -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 (file)
index 0000000..407b2db
--- /dev/null
@@ -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 <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <fstream>
+#include <iomanip>
+
+
+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; i<my_start+10; ++i)
+    for (DynamicSparsityPattern::iterator p=sp.begin(i);
+         p != sp.end(i); ++p)
+      deallog << p->row() << ' ' << 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 (file)
index 0000000..0f16dbb
--- /dev/null
@@ -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 (file)
index 0000000..98d9af2
--- /dev/null
@@ -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 <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <fstream>
+#include <iomanip>
+
+
+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; i<my_start+5; ++i)
+    for (DynamicSparsityPattern::iterator p=sp.begin(i);
+         p != sp.end(i); ++p)
+      deallog << p->row() << ' ' << 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 (file)
index 0000000..8b3b075
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.