Reference documentation for deal.II version Git 2375602 2018-06-23 12:24:17 -0500
index_set.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/index_set.h>
17 #include <deal.II/base/memory_consumption.h>
18 #include <deal.II/base/mpi.h>
19 
20 #include <vector>
21 
22 #ifdef DEAL_II_WITH_TRILINOS
23 # ifdef DEAL_II_WITH_MPI
24 # include <Epetra_MpiComm.h>
25 # endif
26 # include <Epetra_Map.h>
27 # include <Epetra_SerialComm.h>
28 #endif
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 
33 
34 #ifdef DEAL_II_WITH_TRILINOS
35 
36 // the 64-bit path uses a few different names, so put that into a separate
37 // implementation
38 
39 # ifdef DEAL_II_WITH_64BIT_INDICES
40 
41 IndexSet::IndexSet(const Epetra_Map &map)
42  : is_compressed(true)
43  , index_space_size(1 + map.MaxAllGID64())
44  , largest_range(numbers::invalid_unsigned_int)
45 {
46  Assert(map.MinAllGID64() == 0,
47  ExcMessage("The Epetra_Map does not contain the global index 0, which "
48  "means some entries are not present on any processor."));
49 
50  // For a contiguous map, we do not need to go through the whole data...
51  if (map.LinearMap())
52  add_range(size_type(map.MinMyGID64()), size_type(map.MaxMyGID64() + 1));
53  else
54  {
55  const size_type n_indices = map.NumMyElements();
56  size_type * indices =
57  reinterpret_cast<size_type *>(map.MyGlobalElements64());
58  add_indices(indices, indices + n_indices);
59  }
60  compress();
61 }
62 
63 # else
64 
65 // this is the standard 32-bit implementation
66 
67 IndexSet::IndexSet(const Epetra_Map &map)
68  : is_compressed(true)
69  , index_space_size(1 + map.MaxAllGID())
70  , largest_range(numbers::invalid_unsigned_int)
71 {
72  Assert(map.MinAllGID() == 0,
73  ExcMessage("The Epetra_Map does not contain the global index 0, which "
74  "means some entries are not present on any processor."));
75 
76  // For a contiguous map, we do not need to go through the whole data...
77  if (map.LinearMap())
78  add_range(size_type(map.MinMyGID()), size_type(map.MaxMyGID() + 1));
79  else
80  {
81  const size_type n_indices = map.NumMyElements();
82  unsigned int * indices =
83  reinterpret_cast<unsigned int *>(map.MyGlobalElements());
84  add_indices(indices, indices + n_indices);
85  }
86  compress();
87 }
88 
89 # endif
90 
91 #endif // ifdef DEAL_II_WITH_TRILINOS
92 
93 
94 
95 void
96 IndexSet::add_range(const size_type begin, const size_type end)
97 {
98  Assert((begin < index_space_size) ||
99  ((begin == index_space_size) && (end == index_space_size)),
100  ExcIndexRangeType<size_type>(begin, 0, index_space_size));
101  Assert(end <= index_space_size,
102  ExcIndexRangeType<size_type>(end, 0, index_space_size + 1));
103  Assert(begin <= end, ExcIndexRangeType<size_type>(begin, 0, end));
104 
105  if (begin != end)
106  {
107  const Range new_range(begin, end);
108 
109  // the new index might be larger than the last index present in the
110  // ranges. Then we can skip the binary search
111  if (ranges.size() == 0 || begin > ranges.back().end)
112  ranges.push_back(new_range);
113  else
114  ranges.insert(Utilities::lower_bound(ranges.begin(),
115  ranges.end(),
116  new_range),
117  new_range);
118  is_compressed = false;
119  }
120 }
121 
122 
123 
124 void
126 {
127  // we will, in the following, modify mutable variables. this can only
128  // work in multithreaded applications if we lock the data structures
129  // via a mutex, so that users can call 'const' functions from threads
130  // in parallel (and these 'const' functions can then call compress()
131  // which itself calls the current function)
133 
134  // see if any of the contiguous ranges can be merged. do not use
135  // std::vector::erase in-place as it is quadratic in the number of
136  // ranges. since the ranges are sorted by their first index, determining
137  // overlap isn't all that hard
138  std::vector<Range>::iterator store = ranges.begin();
139  for (std::vector<Range>::iterator i = ranges.begin(); i != ranges.end();)
140  {
141  std::vector<Range>::iterator next = i;
142  ++next;
143 
144  size_type first_index = i->begin;
145  size_type last_index = i->end;
146 
147  // see if we can merge any of the following ranges
148  while (next != ranges.end() && (next->begin <= last_index))
149  {
150  last_index = std::max(last_index, next->end);
151  ++next;
152  }
153  i = next;
154 
155  // store the new range in the slot we last occupied
156  *store = Range(first_index, last_index);
157  ++store;
158  }
159  // use a compact array with exactly the right amount of storage
160  if (store != ranges.end())
161  {
162  std::vector<Range> new_ranges(ranges.begin(), store);
163  ranges.swap(new_ranges);
164  }
165 
166  // now compute indices within set and the range with most elements
167  size_type next_index = 0, largest_range_size = 0;
168  for (std::vector<Range>::iterator i = ranges.begin(); i != ranges.end(); ++i)
169  {
170  Assert(i->begin < i->end, ExcInternalError());
171 
172  i->nth_index_in_set = next_index;
173  next_index += (i->end - i->begin);
174  if (i->end - i->begin > largest_range_size)
175  {
176  largest_range_size = i->end - i->begin;
177  largest_range = i - ranges.begin();
178  }
179  }
180  is_compressed = true;
181 
182  // check that next_index is correct. needs to be after the previous
183  // statement because we otherwise will get into an endless loop
184  Assert(next_index == n_elements(), ExcInternalError());
185 }
186 
187 
188 
190 {
191  Assert(size() == is.size(), ExcDimensionMismatch(size(), is.size()));
192 
193  compress();
194  is.compress();
195 
196  std::vector<Range>::const_iterator r1 = ranges.begin(),
197  r2 = is.ranges.begin();
198  IndexSet result(size());
199 
200  while ((r1 != ranges.end()) && (r2 != is.ranges.end()))
201  {
202  // if r1 and r2 do not overlap at all, then move the pointer that sits
203  // to the left of the other up by one
204  if (r1->end <= r2->begin)
205  ++r1;
206  else if (r2->end <= r1->begin)
207  ++r2;
208  else
209  {
210  // the ranges must overlap somehow
211  Assert(((r1->begin <= r2->begin) && (r1->end > r2->begin)) ||
212  ((r2->begin <= r1->begin) && (r2->end > r1->begin)),
213  ExcInternalError());
214 
215  // add the overlapping range to the result
216  result.add_range(std::max(r1->begin, r2->begin),
217  std::min(r1->end, r2->end));
218 
219  // now move that iterator that ends earlier one up. note that it has
220  // to be this one because a subsequent range may still have a chance
221  // of overlapping with the range that ends later
222  if (r1->end <= r2->end)
223  ++r1;
224  else
225  ++r2;
226  }
227  }
228 
229  result.compress();
230  return result;
231 }
232 
233 
234 
235 IndexSet
236 IndexSet::get_view(const size_type begin, const size_type end) const
237 {
238  Assert(begin <= end,
239  ExcMessage("End index needs to be larger or equal to begin index!"));
240  Assert(end <= size(), ExcMessage("Given range exceeds index set dimension"));
241 
242  IndexSet result(end - begin);
243  std::vector<Range>::const_iterator r1 = ranges.begin();
244 
245  while (r1 != ranges.end())
246  {
247  if ((r1->end > begin) && (r1->begin < end))
248  {
249  result.add_range(std::max(r1->begin, begin) - begin,
250  std::min(r1->end, end) - begin);
251  }
252  else if (r1->begin >= end)
253  break;
254 
255  ++r1;
256  }
257 
258  result.compress();
259  return result;
260 }
261 
262 
263 
264 void
266 {
267  compress();
268  other.compress();
269  is_compressed = false;
270 
271 
272  // we save new ranges to be added to our IndexSet in an temporary vector and
273  // add all of them in one go at the end.
274  std::vector<Range> new_ranges;
275 
276  std::vector<Range>::iterator own_it = ranges.begin();
277  std::vector<Range>::iterator other_it = other.ranges.begin();
278 
279  while (own_it != ranges.end() && other_it != other.ranges.end())
280  {
281  // advance own iterator until we get an overlap
282  if (own_it->end <= other_it->begin)
283  {
284  ++own_it;
285  continue;
286  }
287  // we are done with other_it, so advance
288  if (own_it->begin >= other_it->end)
289  {
290  ++other_it;
291  continue;
292  }
293 
294  // Now own_it and other_it overlap. First save the part of own_it that
295  // is before other_it (if not empty).
296  if (own_it->begin < other_it->begin)
297  {
298  Range r(own_it->begin, other_it->begin);
299  r.nth_index_in_set = 0; // fix warning of unused variable
300  new_ranges.push_back(r);
301  }
302  // change own_it to the sub range behind other_it. Do not delete own_it
303  // in any case. As removal would invalidate iterators, we just shrink
304  // the range to an empty one.
305  own_it->begin = other_it->end;
306  if (own_it->begin > own_it->end)
307  {
308  own_it->begin = own_it->end;
309  ++own_it;
310  }
311 
312  // continue without advancing iterators, the right one will be advanced
313  // next.
314  }
315 
316  // Now delete all empty ranges we might
317  // have created.
318  for (std::vector<Range>::iterator it = ranges.begin(); it != ranges.end();)
319  {
320  if (it->begin >= it->end)
321  it = ranges.erase(it);
322  else
323  ++it;
324  }
325 
326  // done, now add the temporary ranges
327  const std::vector<Range>::iterator end = new_ranges.end();
328  for (std::vector<Range>::iterator it = new_ranges.begin(); it != end; ++it)
329  add_range(it->begin, it->end);
330 
331  compress();
332 }
333 
334 
335 
338 {
339  Assert(is_empty() == false,
340  ExcMessage(
341  "pop_back() failed, because this IndexSet contains no entries."));
342 
343  const size_type index = ranges.back().end - 1;
344  --ranges.back().end;
345 
346  if (ranges.back().begin == ranges.back().end)
347  ranges.pop_back();
348 
349  return index;
350 }
351 
352 
353 
356 {
357  Assert(is_empty() == false,
358  ExcMessage(
359  "pop_front() failed, because this IndexSet contains no entries."));
360 
361  const size_type index = ranges.front().begin;
362  ++ranges.front().begin;
363 
364  if (ranges.front().begin == ranges.front().end)
365  ranges.erase(ranges.begin());
366 
367  // We have to set this in any case, because nth_index_in_set is no longer
368  // up to date for all but the first range
369  is_compressed = false;
370 
371  return index;
372 }
373 
374 
375 
376 void
377 IndexSet::add_indices(const IndexSet &other, const unsigned int offset)
378 {
379  if ((this == &other) && (offset == 0))
380  return;
381 
382  Assert(other.ranges.size() == 0 ||
383  other.ranges.back().end - 1 < index_space_size,
384  ExcIndexRangeType<size_type>(other.ranges.back().end - 1,
385  0,
387 
388  compress();
389  other.compress();
390 
391  std::vector<Range>::const_iterator r1 = ranges.begin(),
392  r2 = other.ranges.begin();
393 
394  std::vector<Range> new_ranges;
395  // just get the start and end of the ranges right in this method, everything
396  // else will be done in compress()
397  while (r1 != ranges.end() || r2 != other.ranges.end())
398  {
399  // the two ranges do not overlap or we are at the end of one of the
400  // ranges
401  if (r2 == other.ranges.end() ||
402  (r1 != ranges.end() && r1->end < (r2->begin + offset)))
403  {
404  new_ranges.push_back(*r1);
405  ++r1;
406  }
407  else if (r1 == ranges.end() || (r2->end + offset) < r1->begin)
408  {
409  new_ranges.emplace_back(r2->begin + offset, r2->end + offset);
410  ++r2;
411  }
412  else
413  {
414  // ok, we do overlap, so just take the combination of the current
415  // range (do not bother to merge with subsequent ranges)
416  Range next(std::min(r1->begin, r2->begin + offset),
417  std::max(r1->end, r2->end + offset));
418  new_ranges.push_back(next);
419  ++r1;
420  ++r2;
421  }
422  }
423  ranges.swap(new_ranges);
424 
425  is_compressed = false;
426  compress();
427 }
428 
429 
430 
431 void
432 IndexSet::write(std::ostream &out) const
433 {
434  compress();
435  out << size() << " ";
436  out << ranges.size() << std::endl;
437  std::vector<Range>::const_iterator r = ranges.begin();
438  for (; r != ranges.end(); ++r)
439  {
440  out << r->begin << " " << r->end << std::endl;
441  }
442 }
443 
444 
445 
446 void
447 IndexSet::read(std::istream &in)
448 {
449  AssertThrow(in, ExcIO());
450 
451  size_type s;
452  unsigned int n_ranges;
453 
454  in >> s >> n_ranges;
455  ranges.clear();
456  set_size(s);
457  for (unsigned int i = 0; i < n_ranges; ++i)
458  {
459  AssertThrow(in, ExcIO());
460 
461  size_type b, e;
462  in >> b >> e;
463  add_range(b, e);
464  }
465 }
466 
467 
468 void
469 IndexSet::block_write(std::ostream &out) const
470 {
471  AssertThrow(out, ExcIO());
472  out.write(reinterpret_cast<const char *>(&index_space_size),
473  sizeof(index_space_size));
474  size_t n_ranges = ranges.size();
475  out.write(reinterpret_cast<const char *>(&n_ranges), sizeof(n_ranges));
476  if (ranges.empty() == false)
477  out.write(reinterpret_cast<const char *>(&*ranges.begin()),
478  ranges.size() * sizeof(Range));
479  AssertThrow(out, ExcIO());
480 }
481 
482 void
483 IndexSet::block_read(std::istream &in)
484 {
485  size_type size;
486  size_t n_ranges;
487  in.read(reinterpret_cast<char *>(&size), sizeof(size));
488  in.read(reinterpret_cast<char *>(&n_ranges), sizeof(n_ranges));
489  // we have to clear ranges first
490  ranges.clear();
491  set_size(size);
492  ranges.resize(n_ranges, Range(0, 0));
493  if (n_ranges)
494  in.read(reinterpret_cast<char *>(&*ranges.begin()),
495  ranges.size() * sizeof(Range));
496 
497  do_compress(); // needed so that largest_range can be recomputed
498 }
499 
500 
501 
502 void
503 IndexSet::fill_index_vector(std::vector<size_type> &indices) const
504 {
505  compress();
506 
507  indices.clear();
508  indices.reserve(n_elements());
509 
510  for (std::vector<Range>::iterator it = ranges.begin(); it != ranges.end();
511  ++it)
512  for (size_type i = it->begin; i < it->end; ++i)
513  indices.push_back(i);
514 
515  Assert(indices.size() == n_elements(), ExcInternalError());
516 }
517 
518 
519 
520 #ifdef DEAL_II_WITH_TRILINOS
521 
522 Epetra_Map
523 IndexSet::make_trilinos_map(const MPI_Comm &communicator,
524  const bool overlapping) const
525 {
526  compress();
527  (void)communicator;
528 
529 # ifdef DEBUG
530  if (!overlapping)
531  {
532  const size_type n_global_elements =
533  Utilities::MPI::sum(n_elements(), communicator);
534  Assert(n_global_elements == size(),
535  ExcMessage("You are trying to create an Epetra_Map object "
536  "that partitions elements of an index set "
537  "between processors. However, the union of the "
538  "index sets on different processors does not "
539  "contain all indices exactly once: the sum of "
540  "the number of entries the various processors "
541  "want to store locally is " +
542  Utilities::to_string(n_global_elements) +
543  " whereas the total size of the object to be "
544  "allocated is " +
546  ". In other words, there are "
547  "either indices that are not spoken for "
548  "by any processor, or there are indices that are "
549  "claimed by multiple processors."));
550  }
551 # endif
552 
553  // Find out if the IndexSet is ascending and 1:1. This corresponds to a
554  // linear EpetraMap. Overlapping IndexSets are never 1:1.
555  const bool linear =
556  overlapping ? false : is_ascending_and_one_to_one(communicator);
557 
558  if (linear)
559  return Epetra_Map(TrilinosWrappers::types::int_type(size()),
560  TrilinosWrappers::types::int_type(n_elements()),
561  0,
562 # ifdef DEAL_II_WITH_MPI
563  Epetra_MpiComm(communicator)
564 # else
565  Epetra_SerialComm()
566 # endif
567  );
568  else
569  {
570  std::vector<size_type> indices;
571  fill_index_vector(indices);
572  return Epetra_Map(
573  TrilinosWrappers::types::int_type(-1),
574  TrilinosWrappers::types::int_type(n_elements()),
575  (n_elements() > 0 ?
576  reinterpret_cast<TrilinosWrappers::types::int_type *>(
577  indices.data()) :
578  nullptr),
579  0,
580 # ifdef DEAL_II_WITH_MPI
581  Epetra_MpiComm(communicator)
582 # else
583  Epetra_SerialComm()
584 # endif
585  );
586  }
587 }
588 #endif
589 
590 
591 
592 bool
593 IndexSet::is_ascending_and_one_to_one(const MPI_Comm &communicator) const
594 {
595  // If the sum of local elements does not add up to the total size,
596  // the IndexSet can't be complete.
597  const size_type n_global_elements =
598  Utilities::MPI::sum(n_elements(), communicator);
599  if (n_global_elements != size())
600  return false;
601 
602 #ifdef DEAL_II_WITH_MPI
603  // Non-contiguous IndexSets can't be linear.
604  const bool all_contiguous =
605  (Utilities::MPI::min(is_contiguous() ? 1 : 0, communicator) == 1);
606  if (!all_contiguous)
607  return false;
608 
609  bool is_globally_ascending = true;
610  // we know that there is only one interval
611  types::global_dof_index first_local_dof = (n_elements() > 0) ?
612  *(begin_intervals()->begin()) :
614 
615  const unsigned int my_rank = Utilities::MPI::this_mpi_process(communicator);
616  const std::vector<types::global_dof_index> global_dofs =
617  Utilities::MPI::gather(communicator, first_local_dof, 0);
618 
619  if (my_rank == 0)
620  {
621  // find out if the received std::vector is ascending
622  types::global_dof_index index = 0;
623  while (global_dofs[index] == numbers::invalid_dof_index)
624  ++index;
625  types::global_dof_index old_dof = global_dofs[index++];
626  for (; index < global_dofs.size(); ++index)
627  {
628  const types::global_dof_index new_dof = global_dofs[index];
629  if (new_dof != numbers::invalid_dof_index)
630  {
631  if (new_dof <= old_dof)
632  {
633  is_globally_ascending = false;
634  break;
635  }
636  else
637  old_dof = new_dof;
638  }
639  }
640  }
641 
642  // now broadcast the result
643  int is_ascending = is_globally_ascending ? 1 : 0;
644  int ierr = MPI_Bcast(&is_ascending, 1, MPI_INT, 0, communicator);
645  AssertThrowMPI(ierr);
646 
647  return (is_ascending == 1);
648 #else
649  return true;
650 #endif // DEAL_II_WITH_MPI
651 }
652 
653 
654 
655 std::size_t
657 {
661  sizeof(compress_mutex));
662 }
663 
664 
665 
666 DEAL_II_NAMESPACE_CLOSE
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:912
IntervalIterator begin_intervals() const
Definition: index_set.h:1566
static::ExceptionBase & ExcIO()
void block_read(std::istream &in)
Definition: index_set.cc:483
bool is_ascending_and_one_to_one(const MPI_Comm &communicator) const
Definition: index_set.cc:593
size_type n_elements() const
Definition: index_set.h:1743
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1652
#define AssertThrow(cond, exc)
Definition: exceptions.h:1329
bool is_empty() const
Definition: index_set.h:1735
size_type size() const
Definition: index_set.h:1611
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:105
size_type index_space_size
Definition: index_set.h:924
ElementIterator begin() const
Definition: index_set.h:1030
void block_write(std::ostream &out) const
Definition: index_set.cc:469
bool is_contiguous() const
Definition: index_set.h:1726
static::ExceptionBase & ExcMessage(std::string arg1)
void write(std::ostream &out) const
Definition: index_set.cc:432
std::vector< Range > ranges
Definition: index_set.h:908
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:89
void subtract_set(const IndexSet &other)
Definition: index_set.cc:265
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:503
#define Assert(cond, exc)
Definition: exceptions.h:1227
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
IndexSet operator&(const IndexSet &is) const
Definition: index_set.cc:189
std::vector< T > gather(const MPI_Comm &comm, const T &object_to_send, const unsigned int root_process=0)
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:236
void compress() const
Definition: index_set.h:1619
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:523
size_type pop_back()
Definition: index_set.cc:337
void add_range(const size_type begin, const size_type end)
Definition: index_set.cc:96
void set_size(const size_type size)
Definition: index_set.h:1599
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1438
Threads::Mutex compress_mutex
Definition: index_set.h:941
void do_compress() const
Definition: index_set.cc:125
T min(const T &t, const MPI_Comm &mpi_communicator)
size_type pop_front()
Definition: index_set.cc:355
bool is_compressed
Definition: index_set.h:918
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:80
std::size_t memory_consumption() const
Definition: index_set.cc:656
ElementIterator end() const
Definition: index_set.h:1557
const types::global_dof_index invalid_dof_index
Definition: types.h:188
types::global_dof_index size_type
Definition: index_set.h:76
void read(std::istream &in)
Definition: index_set.cc:447
size_type largest_range
Definition: index_set.h:935
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static::ExceptionBase & ExcInternalError()