index_set.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2009 - 2018 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.
//
// ---------------------------------------------------------------------

#ifndef dealii_index_set_h
#define dealii_index_set_h

#include <deal.II/base/config.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/thread_management.h>
#include <boost/serialization/vector.hpp>
#include <vector>
#include <algorithm>
#include <iterator>


#ifdef DEAL_II_WITH_TRILINOS
#  include <Epetra_Map.h>
#endif

#if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
#include <mpi.h>
#else
typedef int MPI_Comm;
#  ifndef MPI_COMM_WORLD
#    define MPI_COMM_WORLD 0
#  endif
#endif

DEAL_II_NAMESPACE_OPEN

class IndexSet
{
public:
  // forward declarations:
  class ElementIterator;
  class IntervalIterator;

  typedef types::global_dof_index size_type;

  typedef signed int value_type;


  IndexSet ();

  explicit IndexSet (const size_type size);

  IndexSet (const IndexSet &) = default;

  IndexSet &operator= (const IndexSet &) = default;

  IndexSet (IndexSet &&is) noexcept;

  IndexSet &operator= (IndexSet &&is) noexcept;

#ifdef DEAL_II_WITH_TRILINOS

  explicit IndexSet(const Epetra_Map &map);
#endif

  void clear ();

  void set_size (const size_type size);

  size_type size () const;

  void add_range (const size_type begin,
                  const size_type end);

  void add_index (const size_type index);

  template <typename ForwardIterator>
  void add_indices (const ForwardIterator &begin,
                    const ForwardIterator &end);

  void add_indices(const IndexSet &other,
                   const unsigned int offset = 0);

  bool is_element (const size_type index) const;

  bool is_contiguous () const;

  bool is_empty () const;

  bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const;

  size_type n_elements () const;

  size_type nth_index_in_set (const unsigned int local_index) const;

  size_type index_within_set (const size_type global_index) const;

  unsigned int n_intervals () const;

  unsigned int largest_range_starting_index() const;

  void compress () const;

  bool operator == (const IndexSet &is) const;

  bool operator != (const IndexSet &is) const;

  IndexSet operator & (const IndexSet &is) const;

  IndexSet get_view (const size_type begin,
                     const size_type end) const;

  void subtract_set (const IndexSet &other);

  size_type pop_back ();

  size_type pop_front ();

  void fill_index_vector(std::vector<size_type> &indices) const;

  template <typename VectorType>
  void fill_binary_vector (VectorType &vector) const;

  template <class StreamType>
  void print(StreamType &out) const;

  void write(std::ostream &out) const;

  void read(std::istream &in);

  void block_write(std::ostream &out) const;

  void block_read(std::istream &in);

#ifdef DEAL_II_WITH_TRILINOS

  Epetra_Map make_trilinos_map (const MPI_Comm &communicator = MPI_COMM_WORLD,
                                const bool      overlapping  = false) const;
#endif


  std::size_t memory_consumption () const;

  DeclException1 (ExcIndexNotPresent, size_type,
                  << "The global index " << arg1
                  << " is not an element of this set.");

  template <class Archive>
  void serialize (Archive &ar, const unsigned int version);


  class IntervalAccessor
  {
  public:
    IntervalAccessor(const IndexSet *idxset, const size_type range_idx);

    explicit IntervalAccessor(const IndexSet *idxset);

    size_type n_elements() const;

    bool is_valid() const;

    ElementIterator begin() const;

    ElementIterator end() const;

    size_type last() const;

  private:
    IntervalAccessor(const IntervalAccessor &other);
    IntervalAccessor &operator = (const IntervalAccessor &other);

    bool operator == (const IntervalAccessor &other) const;
    bool operator < (const IntervalAccessor &other) const;
    void advance ();
    const IndexSet *index_set;

    size_type range_idx;

    friend class IntervalIterator;
  };

  class IntervalIterator
  {
  public:
    IntervalIterator(const IndexSet *idxset, const size_type range_idx);

    explicit IntervalIterator(const IndexSet *idxset);

    IntervalIterator();

    IntervalIterator(const IntervalIterator &other);

    IntervalIterator &operator = (const IntervalIterator &other);

    IntervalIterator &operator++ ();

    IntervalIterator operator++ (int);

    const IntervalAccessor &operator* () const;

    const IntervalAccessor *operator-> () const;

    bool operator == (const IntervalIterator &) const;

    bool operator != (const IntervalIterator &) const;

    bool operator < (const IntervalIterator &) const;

    int operator - (const IntervalIterator &p) const;

    typedef std::forward_iterator_tag iterator_category;
    typedef IntervalAccessor          value_type;
    typedef std::ptrdiff_t            difference_type;
    typedef IntervalAccessor         *pointer;
    typedef IntervalAccessor         &reference;

  private:
    IntervalAccessor accessor;
  };

  class ElementIterator
  {
  public:
    ElementIterator(const IndexSet *idxset,
                    const size_type range_idx,
                    const size_type index);

    explicit ElementIterator(const IndexSet *idxset);

    size_type operator* () const;

    bool is_valid () const;

    ElementIterator &operator++ ();

    ElementIterator operator++ (int);

    bool operator == (const ElementIterator &) const;

    bool operator != (const ElementIterator &) const;

    bool operator < (const ElementIterator &) const;

    std::ptrdiff_t operator - (const ElementIterator &p) const;

    typedef std::forward_iterator_tag iterator_category;
    typedef size_type                 value_type;
    typedef std::ptrdiff_t            difference_type;
    typedef size_type                *pointer;
    typedef size_type                &reference;

  private:
    void advance ();

    const IndexSet *index_set;
    size_type range_idx;
    size_type idx;
  };

  ElementIterator begin() const;

  ElementIterator at(const size_type global_index) const;

  ElementIterator end() const;

  IntervalIterator begin_intervals() const;

  IntervalIterator end_intervals() const;

private:
  struct Range
  {
    size_type begin;
    size_type end;

    size_type nth_index_in_set;

    Range ();

    Range (const size_type i1,
           const size_type i2);

    friend
    inline bool operator< (const Range &range_1,
                           const Range &range_2)
    {
      return ((range_1.begin < range_2.begin)
              ||
              ((range_1.begin == range_2.begin)
               &&
               (range_1.end < range_2.end)));
    }

    static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
    {
      return x.end < y.end;
    }

    static bool nth_index_compare (const IndexSet::Range &x,
                                   const IndexSet::Range &y)
    {
      return (x.nth_index_in_set+(x.end-x.begin) <
              y.nth_index_in_set+(y.end-y.begin));
    }

    friend
    inline bool operator== (const Range &range_1,
                            const Range &range_2)
    {
      return ((range_1.begin == range_2.begin)
              &&
              (range_1.end == range_2.end));
    }

    static std::size_t memory_consumption ()
    {
      return sizeof(Range);
    }

    template <class Archive>
    void serialize (Archive &ar, const unsigned int version);
  };

  mutable std::vector<Range> ranges;

  mutable bool is_compressed;

  size_type index_space_size;

  mutable size_type largest_range;

  mutable Threads::Mutex compress_mutex;

  void do_compress() const;
};


inline
IndexSet complete_index_set (const unsigned int N)
{
  IndexSet is (N);
  is.add_range(0, N);
  is.compress();
  return is;
}

/* ------------------ inline functions ------------------ */


/* IntervalAccessor */

inline
IndexSet::IntervalAccessor::
IntervalAccessor(const IndexSet *idxset,
                 const IndexSet::size_type range_idx)
  :
  index_set(idxset),
  range_idx(range_idx)
{
  Assert(range_idx < idxset->n_intervals(), ExcInternalError("Invalid range index"));
}



inline
IndexSet::IntervalAccessor::
IntervalAccessor(const IndexSet *idxset)
  :
  index_set(idxset),
  range_idx(numbers::invalid_dof_index)
{}



inline
IndexSet::IntervalAccessor::
IntervalAccessor(const IndexSet::IntervalAccessor &other)
  :
  index_set (other.index_set),
  range_idx(other.range_idx)
{
  Assert( range_idx == numbers::invalid_dof_index || is_valid(),
          ExcMessage("invalid iterator"));
}



inline
IndexSet::size_type
IndexSet::IntervalAccessor::n_elements() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));
  return index_set->ranges[range_idx].end - index_set->ranges[range_idx].begin;
}



inline
bool
IndexSet::IntervalAccessor::is_valid() const
{
  return index_set != nullptr && range_idx < index_set->n_intervals();
}



inline
IndexSet::ElementIterator
IndexSet::IntervalAccessor::begin() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));
  return IndexSet::ElementIterator(index_set, range_idx,
                                   index_set->ranges[range_idx].begin);
}



inline
IndexSet::ElementIterator
IndexSet::IntervalAccessor::end() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));

  // point to first index in next interval unless we are the last interval.
  if (range_idx < index_set->ranges.size()-1)
    return IndexSet::ElementIterator(index_set, range_idx+1, index_set->ranges[range_idx+1].begin);
  else
    return index_set->end();
}



inline
IndexSet::size_type
IndexSet::IntervalAccessor::last() const
{
  Assert(is_valid(), ExcMessage("invalid iterator"));

  return index_set->ranges[range_idx].end-1;
}



inline
IndexSet::IntervalAccessor &
IndexSet::IntervalAccessor::operator = (const IndexSet::IntervalAccessor &other)
{
  index_set = other.index_set;
  range_idx = other.range_idx;
  Assert( range_idx == numbers::invalid_dof_index || is_valid(),
          ExcMessage("invalid iterator"));
  return *this;
}



inline
bool
IndexSet::IntervalAccessor::operator == (const IndexSet::IntervalAccessor &other) const
{
  Assert (index_set == other.index_set,
          ExcMessage("Can not compare accessors pointing to different IndexSets"));
  return range_idx == other.range_idx;
}



inline
bool
IndexSet::IntervalAccessor::operator < (const IndexSet::IntervalAccessor &other) const
{
  Assert (index_set == other.index_set,
          ExcMessage("Can not compare accessors pointing to different IndexSets"));
  return range_idx < other.range_idx;
}



inline
void IndexSet::IntervalAccessor::advance ()
{
  Assert (is_valid(),
          ExcMessage("Impossible to advance an IndexSet::IntervalIterator that is invalid"));
  ++range_idx;

  // set ourselves to invalid if we walk off the end
  if (range_idx>=index_set->ranges.size())
    range_idx = numbers::invalid_dof_index;
}


/* IntervalIterator */

inline
IndexSet::IntervalIterator::IntervalIterator(const IndexSet *idxset,
                                             const IndexSet::size_type range_idx)
  :
  accessor(idxset, range_idx)
{}



inline
IndexSet::IntervalIterator::IntervalIterator()
  :
  accessor(nullptr)
{}



inline
IndexSet::IntervalIterator::IntervalIterator(const IndexSet *idxset)
  :
  accessor(idxset)
{}



inline
IndexSet::IntervalIterator::IntervalIterator(const IndexSet::IntervalIterator &other)
  :
  accessor(other.accessor)
{}



inline
IndexSet::IntervalIterator &
IndexSet::IntervalIterator::operator = (const IntervalIterator &other)
{
  accessor = other.accessor;
  return *this;
}



inline
IndexSet::IntervalIterator &
IndexSet::IntervalIterator::operator++ ()
{
  accessor.advance();
  return *this;
}



inline
IndexSet::IntervalIterator
IndexSet::IntervalIterator::operator++ (int)
{
  const IndexSet::IntervalIterator iter = *this;
  accessor.advance ();
  return iter;
}



inline
const IndexSet::IntervalAccessor &
IndexSet::IntervalIterator::operator* () const
{
  return accessor;
}



inline
const IndexSet::IntervalAccessor *
IndexSet::IntervalIterator::operator-> () const
{
  return &accessor;
}



inline
bool
IndexSet::IntervalIterator::operator == (const IndexSet::IntervalIterator &other) const
{
  return accessor == other.accessor;
}



inline
bool
IndexSet::IntervalIterator::operator != (const IndexSet::IntervalIterator &other) const
{
  return !(*this == other);
}



inline
bool
IndexSet::IntervalIterator::operator < (const IndexSet::IntervalIterator &other) const
{
  return accessor < other.accessor;
}



inline
int
IndexSet::IntervalIterator::operator - (const IndexSet::IntervalIterator &other) const
{
  Assert (accessor.index_set == other.accessor.index_set,
          ExcMessage("Can not compare iterators belonging to different IndexSets"));

  const size_type lhs = (accessor.range_idx == numbers::invalid_dof_index)
                        ?
                        accessor.index_set->ranges.size()
                        :
                        accessor.range_idx;
  const size_type rhs = (other.accessor.range_idx == numbers::invalid_dof_index)
                        ?
                        accessor.index_set->ranges.size()
                        :
                        other.accessor.range_idx;

  if (lhs > rhs)
    return static_cast<int>(lhs - rhs);
  else
    return -static_cast<int>(rhs - lhs);
}



/* ElementIterator */

inline
IndexSet::ElementIterator::
ElementIterator(const IndexSet *idxset,
                const IndexSet::size_type range_idx,
                const IndexSet::size_type index)
  :
  index_set(idxset),
  range_idx(range_idx),
  idx(index)
{
  Assert(range_idx < index_set->ranges.size(),
         ExcMessage("Invalid range index for IndexSet::ElementIterator constructor."));
  Assert(idx >= index_set->ranges[range_idx].begin
         &&
         idx < index_set->ranges[range_idx].end,
         ExcInternalError("Invalid index argument for IndexSet::ElementIterator constructor."));
}



inline
IndexSet::ElementIterator::ElementIterator(const IndexSet *idxset)
  :
  index_set(idxset),
  range_idx(numbers::invalid_dof_index),
  idx(numbers::invalid_dof_index)
{}



inline
bool
IndexSet::ElementIterator::is_valid() const
{
  Assert((range_idx == numbers::invalid_dof_index && idx == numbers::invalid_dof_index)
         ||
         (range_idx < index_set->ranges.size() && idx<index_set->ranges[range_idx].end),
         ExcInternalError("Invalid ElementIterator state."));

  return (range_idx < index_set->ranges.size() &&
          idx<index_set->ranges[range_idx].end);
}




inline
IndexSet::size_type
IndexSet::ElementIterator::operator* () const
{
  Assert (is_valid(),
          ExcMessage("Impossible to dereference an IndexSet::ElementIterator that is invalid"));
  return idx;
}



inline
bool
IndexSet::ElementIterator::operator == (const IndexSet::ElementIterator &other) const
{
  Assert (index_set == other.index_set,
          ExcMessage("Can not compare iterators belonging to different IndexSets"));
  return range_idx == other.range_idx && idx==other.idx;
}



inline
void
IndexSet::ElementIterator::advance ()
{
  Assert (is_valid(),
          ExcMessage("Impossible to advance an IndexSet::ElementIterator that is invalid"));
  if (idx < index_set->ranges[range_idx].end)
    ++idx;
  // end of this range?
  if (idx == index_set->ranges[range_idx].end)
    {
      // point to first element in next interval if possible
      if (range_idx < index_set->ranges.size()-1)
        {
          ++range_idx;
          idx = index_set->ranges[range_idx].begin;
        }
      else
        {
          // we just fell off the end, set to invalid:
          range_idx = numbers::invalid_dof_index;
          idx = numbers::invalid_dof_index;
        }
    }
}



inline
IndexSet::ElementIterator &
IndexSet::ElementIterator::operator++ ()
{
  advance();
  return *this;
}



inline
IndexSet::ElementIterator
IndexSet::ElementIterator::operator++ (int)
{
  const IndexSet::ElementIterator it = *this;
  advance();
  return it;
}



inline
bool IndexSet::ElementIterator::operator != (const IndexSet::ElementIterator &other) const
{
  return !(*this == other);
}



inline
bool IndexSet::ElementIterator::operator < (const IndexSet::ElementIterator &other) const
{
  Assert (index_set == other.index_set,
          ExcMessage("Can not compare iterators belonging to different IndexSets"));
  return range_idx < other.range_idx || (range_idx == other.range_idx && idx<other.idx);
}



inline
std::ptrdiff_t
IndexSet::ElementIterator::operator - (const IndexSet::ElementIterator &other) const
{
  Assert (index_set == other.index_set,
          ExcMessage("Can not compare iterators belonging to different IndexSets"));
  if (*this == other)
    return 0;
  if (!(*this < other))
    return -(other-*this);

  // only other can be equal to end() because of the checks above.
  Assert (is_valid(), ExcInternalError());

  // Note: we now compute how far advance *this in "*this < other" to get other, so we need to return -c at the end.

  // first finish the current range:
  std::ptrdiff_t c = index_set->ranges[range_idx].end-idx;

  // now walk in steps of ranges (need to start one behind our current one):
  for (size_type range=range_idx+1; range<index_set->ranges.size() && range<=other.range_idx; ++range)
    c += index_set->ranges[range].end-index_set->ranges[range].begin;

  Assert(other.range_idx < index_set->ranges.size() || other.range_idx == numbers::invalid_dof_index,
         ExcMessage("Inconsistent iterator state. Did you invalidate iterators by modifying the IndexSet?"));

  // We might have walked too far because we went until the end of other.range_idx, so walk backwards to other.idx:
  if (other.range_idx != numbers::invalid_dof_index)
    c -= index_set->ranges[other.range_idx].end - other.idx;

  return -c;
}


/* Range */

inline
IndexSet::Range::Range ()
  :
  begin(numbers::invalid_dof_index),
  end(numbers::invalid_dof_index),
  nth_index_in_set(numbers::invalid_dof_index)
{}



inline
IndexSet::Range::Range (const size_type i1,
                        const size_type i2)
  :
  begin(i1),
  end(i2),
  nth_index_in_set(numbers::invalid_dof_index)
{}



/* IndexSet itself */

inline
IndexSet::IndexSet ()
  :
  is_compressed (true),
  index_space_size (0),
  largest_range (numbers::invalid_unsigned_int)
{}



inline
IndexSet::IndexSet (const size_type size)
  :
  is_compressed (true),
  index_space_size (size),
  largest_range (numbers::invalid_unsigned_int)
{}



inline
IndexSet::IndexSet (IndexSet &&is) noexcept
:
ranges (std::move(is.ranges)),
       is_compressed (is.is_compressed),
       index_space_size (is.index_space_size),
       largest_range (is.largest_range)
{
  is.ranges.clear ();
  is.is_compressed = true;
  is.index_space_size = 0;
  is.largest_range = numbers::invalid_unsigned_int;

  compress ();
}



inline
IndexSet &IndexSet::operator= (IndexSet &&is) noexcept
{
  ranges = std::move (is.ranges);
  is_compressed = is.is_compressed;
  index_space_size = is.index_space_size;
  largest_range = is.largest_range;

  is.ranges.clear ();
  is.is_compressed = true;
  is.index_space_size = 0;
  is.largest_range = numbers::invalid_unsigned_int;

  compress ();

  return *this;
}



inline
IndexSet::ElementIterator IndexSet::begin() const
{
  compress();
  if (ranges.size()>0)
    return IndexSet::ElementIterator(this, 0, ranges[0].begin);
  else
    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
{
  compress();
  return IndexSet::ElementIterator(this);
}



inline
IndexSet::IntervalIterator
IndexSet::begin_intervals() const
{
  compress();
  if (ranges.size()>0)
    return IndexSet::IntervalIterator(this, 0);
  else
    return end_intervals();
}



inline
IndexSet::IntervalIterator
IndexSet::end_intervals() const
{
  compress();
  return IndexSet::IntervalIterator(this);
}



inline
void
IndexSet::clear ()
{
  // reset so that there are no indices in the set any more; however,
  // as documented, the index set retains its size
  ranges.clear ();
  is_compressed = true;
  largest_range = numbers::invalid_unsigned_int;
}



inline
void
IndexSet::set_size (const size_type sz)
{
  Assert (ranges.empty(),
          ExcMessage ("This function can only be called if the current "
                      "object does not yet contain any elements."));
  index_space_size = sz;
  is_compressed = true;
}



inline
IndexSet::size_type
IndexSet::size () const
{
  return index_space_size;
}



inline
void
IndexSet::compress () const
{
  if (is_compressed == true)
    return;

  do_compress();
}



inline
void
IndexSet::add_index (const size_type index)
{
  Assert (index < index_space_size,
          ExcIndexRangeType<size_type> (index, 0, index_space_size));

  const Range new_range(index, index+1);
  if (ranges.size() == 0 || index > ranges.back().end)
    ranges.push_back(new_range);
  else if (index == ranges.back().end)
    ranges.back().end++;
  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
  for (ForwardIterator p=begin; p!=end;)
    {
      const size_type begin_index = *p;
      size_type       end_index   = begin_index + 1;
      ForwardIterator q = p;
      ++q;
      while ((q != end) && (*q == end_index))
        {
          ++end_index;
          ++q;
        }

      add_range (begin_index, end_index);
      p = q;
    }
}



inline
bool
IndexSet::is_element (const size_type index) const
{
  if (ranges.empty() == false)
    {
      compress ();

      // fast check whether the index is in the largest range
      Assert (largest_range < ranges.size(), ExcInternalError());
      if (index >= ranges[largest_range].begin &&
          index < ranges[largest_range].end)
        return true;

      // get the element after which we would have to insert a range that
      // consists of all elements from this element to the end of the index
      // range plus one. after this call we know that if p!=end() then
      // p->begin<=index unless there is no such range at all
      //
      // if the searched for element is an element of this range, then we're
      // done. otherwise, the element can't be in one of the following ranges
      // because otherwise p would be a different iterator
      //
      // since we already know the position relative to the largest range (we
      // called compress!), we can perform the binary search on ranges with
      // lower/higher number compared to the largest range
      std::vector<Range>::const_iterator
      p = std::upper_bound (ranges.begin() + (index<ranges[largest_range].begin?
                                              0 : largest_range+1),
                            index<ranges[largest_range].begin ?
                            ranges.begin() + largest_range:
                            ranges.end(),
                            Range (index, size()+1));

      if (p == ranges.begin())
        return ((index >= p->begin) && (index < p->end));

      Assert ((p == ranges.end()) || (p->begin > index),
              ExcInternalError());

      // now move to that previous range
      --p;
      Assert (p->begin <= index, ExcInternalError());

      return (p->end > index);
    }

  // didn't find this index, so it's not in the set
  return false;
}



inline
bool
IndexSet::is_contiguous () const
{
  compress ();
  return (ranges.size() <= 1);
}



inline
bool
IndexSet::is_empty () const
{
  return ranges.empty();
}



inline
IndexSet::size_type
IndexSet::n_elements () const
{
  // make sure we have non-overlapping ranges
  compress ();

  size_type v = 0;
  if (!ranges.empty())
    {
      Range &r = ranges.back();
      v = r.nth_index_in_set + r.end - r.begin;
    }

#ifdef DEBUG
  size_type s = 0;
  for (std::vector<Range>::iterator range = ranges.begin();
       range != ranges.end();
       ++range)
    s += (range->end - range->begin);
  Assert(s==v, ExcInternalError());
#endif

  return v;
}



inline
unsigned int
IndexSet::n_intervals () const
{
  compress ();
  return ranges.size();
}



inline
unsigned int
IndexSet::largest_range_starting_index() const
{
  Assert(ranges.empty()==false, ExcMessage("IndexSet cannot be empty."));

  compress();
  const std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;

  return main_range->nth_index_in_set;
}



inline
IndexSet::size_type
IndexSet::nth_index_in_set (const unsigned int n) const
{
  Assert (n < n_elements(), ExcIndexRangeType<size_type> (n, 0, n_elements()));

  compress ();

  // first check whether the index is in the largest range
  Assert (largest_range < ranges.size(), ExcInternalError());
  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
  if (n>=main_range->nth_index_in_set &&
      n<main_range->nth_index_in_set+(main_range->end-main_range->begin))
    return main_range->begin + (n-main_range->nth_index_in_set);

  // find out which chunk the local index n belongs to by using a binary
  // search. the comparator is based on the end of the ranges. Use the
  // position relative to main_range to subdivide the ranges
  Range r (n,n+1);
  r.nth_index_in_set = n;
  std::vector<Range>::const_iterator range_begin, range_end;
  if (n<main_range->nth_index_in_set)
    {
      range_begin = ranges.begin();
      range_end   = main_range;
    }
  else
    {
      range_begin = main_range + 1;
      range_end   = ranges.end();
    }

  const std::vector<Range>::const_iterator
  p = Utilities::lower_bound(range_begin, range_end, r,
                             Range::nth_index_compare);

  Assert (p != ranges.end(), ExcInternalError());
  return p->begin + (n-p->nth_index_in_set);
}



inline
IndexSet::size_type
IndexSet::index_within_set (const size_type n) const
{
  // to make this call thread-safe, compress() must not be called through this
  // function
  Assert (is_compressed == true, ExcMessage ("IndexSet must be compressed."));
  Assert (n < size(), ExcIndexRangeType<size_type> (n, 0, size()));

  // return immediately if the index set is empty
  if (is_empty())
    return numbers::invalid_dof_index;

  // check whether the index is in the largest range. use the result to
  // perform a one-sided binary search afterward
  Assert (largest_range < ranges.size(), ExcInternalError());
  std::vector<Range>::const_iterator main_range=ranges.begin()+largest_range;
  if (n >= main_range->begin && n < main_range->end)
    return (n-main_range->begin) + main_range->nth_index_in_set;

  Range r(n, n);
  std::vector<Range>::const_iterator range_begin, range_end;
  if (n<main_range->begin)
    {
      range_begin = ranges.begin();
      range_end   = main_range;
    }
  else
    {
      range_begin = main_range + 1;
      range_end   = ranges.end();
    }

  std::vector<Range>::const_iterator
  p = Utilities::lower_bound(range_begin, range_end, r,
                             Range::end_compare);

  // if n is not in this set
  if (p==range_end || p->end == n || p->begin > n)
    return numbers::invalid_dof_index;

  Assert(p!=ranges.end(), ExcInternalError());
  Assert(p->begin<=n, ExcInternalError());
  Assert(n<p->end, ExcInternalError());
  return (n-p->begin) + p->nth_index_in_set;
}



inline
bool
IndexSet::operator == (const IndexSet &is) const
{
  Assert (size() == is.size(),
          ExcDimensionMismatch (size(), is.size()));

  compress ();
  is.compress ();

  return ranges == is.ranges;
}



inline
bool
IndexSet::operator != (const IndexSet &is) const
{
  Assert (size() == is.size(),
          ExcDimensionMismatch (size(), is.size()));

  compress ();
  is.compress ();

  return ranges != is.ranges;
}



template <typename Vector>
void
IndexSet::fill_binary_vector (Vector &vector) const
{
  Assert (vector.size() == size(),
          ExcDimensionMismatch (vector.size(), size()));

  compress();
  // first fill all elements of the vector with zeroes.
  std::fill (vector.begin(), vector.end(), 0);

  // then write ones into the elements whose indices are contained in the
  // index set
  for (std::vector<Range>::iterator it = ranges.begin();
       it != ranges.end();
       ++it)
    for (size_type i=it->begin; i<it->end; ++i)
      vector[i] = 1;
}



template <class StreamType>
inline
void
IndexSet::print (StreamType &out) const
{
  compress();
  out << "{";
  std::vector<Range>::const_iterator p;
  for (p = ranges.begin(); p != ranges.end(); ++p)
    {
      if (p->end-p->begin==1)
        out << p->begin;
      else
        out << "[" << p->begin << "," << p->end-1 << "]";

      if (p != --ranges.end())
        out << ", ";
    }
  out << "}" << std::endl;
}



template <class Archive>
inline
void
IndexSet::Range::serialize (Archive &ar, const unsigned int)
{
  ar &begin &end &nth_index_in_set;
}



template <class Archive>
inline
void
IndexSet::serialize (Archive &ar, const unsigned int)
{
  ar &ranges &is_compressed &index_space_size &largest_range;
}

DEAL_II_NAMESPACE_CLOSE

#endif