<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.
*/
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.
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
{
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.
*/
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);
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);
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;
{
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);
{
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);
}
DynamicSparsityPattern::DynamicSparsityPattern ()
:
- rows(0),
- cols(0),
- rowset(0)
+ have_entries (false),
+ rows (0),
+ cols (0),
+ rowset (0)
{}
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());
const IndexSet &rowset_
)
:
- rows(0),
- cols(0),
- rowset(0)
+ have_entries (false),
+ rows (0),
+ cols (0),
+ rowset (0)
{
reinit (m,n, rowset_);
}
DynamicSparsityPattern::DynamicSparsityPattern (const IndexSet &rowset_)
:
+ have_entries (false),
rows(0),
cols(0),
rowset(0)
DynamicSparsityPattern::DynamicSparsityPattern (const size_type n)
:
+ have_entries (false),
rows(0),
cols(0),
rowset(0)
const size_type n,
const IndexSet &rowset_)
{
+ have_entries = false;
rows = m;
cols = n;
rowset=rowset_;
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)
{
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);
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)
{
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]);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 ();
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::33333334 1
+DEAL::OK
+DEAL::OK
+DEAL::33333334 1
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
projects / dealii.git / commitdiff