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