22#include <boost/container/small_vector.hpp>
26#ifdef DEAL_II_WITH_TRILINOS
27# ifdef DEAL_II_WITH_MPI
28# include <Epetra_MpiComm.h>
30# include <Epetra_Map.h>
31# include <Epetra_SerialComm.h>
32# ifdef DEAL_II_TRILINOS_WITH_TPETRA
33# include <Tpetra_Map.hpp>
41#ifdef DEAL_II_WITH_TRILINOS
43# ifdef DEAL_II_TRILINOS_WITH_TPETRA
45template <
typename NodeType>
48 const Tpetra::Map<int, types::signed_global_dof_index, NodeType>> &map)
50 , index_space_size(1 + map->getMaxAllGlobalIndex())
51 , largest_range(
numbers::invalid_unsigned_int)
53 Assert(map->getMinAllGlobalIndex() == 0,
55 "The Tpetra::Map does not contain the global index 0, "
56 "which means some entries are not present on any processor."));
59 if (map->isContiguous())
64# if DEAL_II_TRILINOS_VERSION_GTE(13, 4, 0)
65 const size_type n_indices = map->getLocalNumElements();
67 const size_type n_indices = map->getNodeNumElements();
70 map->getMyGlobalIndices().data();
82# ifdef DEAL_II_WITH_64BIT_INDICES
86 , index_space_size(1 + map.MaxAllGID64())
87 , largest_range(
numbers::invalid_unsigned_int)
89 Assert(map.MinAllGID64() == 0,
91 "The Epetra_BlockMap does not contain the global index 0, "
92 "which means some entries are not present on any processor."));
99 const size_type n_indices = map.NumMyElements();
101 reinterpret_cast<size_type *
>(map.MyGlobalElements64());
112 : is_compressed(true)
113 , index_space_size(1 + map.MaxAllGID())
114 , largest_range(
numbers::invalid_unsigned_int)
116 Assert(map.MinAllGID() == 0,
118 "The Epetra_BlockMap does not contain the global index 0, "
119 "which means some entries are not present on any processor."));
126 const size_type n_indices = map.NumMyElements();
127 unsigned int *indices =
128 reinterpret_cast<unsigned int *
>(map.MyGlobalElements());
155 std::vector<Range>::iterator store =
ranges.begin();
156 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();)
158 std::vector<Range>::iterator next = i;
165 while (next !=
ranges.end() && (next->begin <= last_index))
167 last_index =
std::max(last_index, next->end);
173 *store =
Range(first_index, last_index);
177 if (store !=
ranges.end())
179 std::vector<Range> new_ranges(
ranges.begin(), store);
184 size_type next_index = 0, largest_range_size = 0;
185 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();
190 i->nth_index_in_set = next_index;
191 next_index += (i->end - i->begin);
192 if (i->end - i->begin > largest_range_size)
194 largest_range_size = i->end - i->begin;
211 for (
const auto &range :
ranges)
215 "In the process of creating the current IndexSet "
216 "object, you added indices beyond the size of the index "
217 "space. Specifically, you added elements that form the "
219 std::to_string(range.begin) +
"," + std::to_string(range.end) +
220 "), but the size of the index space is only " +
236 std::vector<Range>::const_iterator r1 =
ranges.begin(),
244 if (r1->end <= r2->begin)
246 else if (r2->end <= r1->begin)
251 Assert(((r1->begin <= r2->begin) && (r1->end > r2->begin)) ||
252 ((r2->begin <= r1->begin) && (r2->end > r1->begin)),
256 result.add_range(
std::max(r1->begin, r2->begin),
262 if (r1->end <= r2->end)
280 ExcMessage(
"End index needs to be larger or equal to begin index!"));
282 ExcMessage(
"You are asking for a view into an IndexSet object "
283 "that would cover the sub-range [" +
284 std::to_string(
begin) +
',' + std::to_string(
end) +
285 "). But this is not a subset of the range "
286 "of the current object, which is [0," +
287 std::to_string(
size()) +
")."));
290 std::vector<Range>::const_iterator r1 =
ranges.begin();
292 while (r1 !=
ranges.end())
294 if ((r1->end >
begin) && (r1->begin <
end))
299 else if (r1->begin >=
end)
315 ExcMessage(
"The mask must have the same size index space "
316 "as the index set it is applied to."));
328 if (mask.ranges.size() == 1)
329 return get_view(mask.ranges[0].begin, mask.ranges[0].end);
338 std::vector<Range> new_ranges;
340 std::vector<Range>::iterator own_it =
ranges.begin();
341 std::vector<Range>::iterator mask_it = mask.ranges.begin();
343 while ((own_it !=
ranges.end()) && (mask_it != mask.ranges.end()))
349 if (own_it->end <= mask_it->begin)
357 if (mask_it->end <= own_it->begin)
376 if ((own_it->begin <= mask_it->begin) && (own_it->end <= mask_it->end))
378 new_ranges.emplace_back(mask_it->begin - mask_it->nth_index_in_set,
379 own_it->end - mask_it->nth_index_in_set);
383 if ((mask_it->begin <= own_it->begin) && (mask_it->end <= own_it->end))
385 const size_type offset_within_mask_interval =
386 own_it->begin - mask_it->begin;
387 new_ranges.emplace_back(mask_it->nth_index_in_set +
388 offset_within_mask_interval,
389 mask_it->nth_index_in_set +
390 (mask_it->end - mask_it->begin));
394 if ((own_it->begin <= mask_it->begin) &&
395 (own_it->end >= mask_it->end))
397 new_ranges.emplace_back(mask_it->begin -
398 mask_it->nth_index_in_set,
399 mask_it->end - mask_it->nth_index_in_set);
403 if ((mask_it->begin <= own_it->begin) &&
404 (mask_it->end >= own_it->end))
406 const size_type offset_within_mask_interval =
407 own_it->begin - mask_it->begin;
408 new_ranges.emplace_back(mask_it->nth_index_in_set +
409 offset_within_mask_interval,
410 mask_it->nth_index_in_set +
411 offset_within_mask_interval +
412 (own_it->end - own_it->begin));
421 if (own_it->end < mask_it->end)
423 else if (mask_it->end < own_it->end)
438 for (
const auto &range : new_ranges)
439 result.
add_range(range.begin, range.end);
449 const std::vector<types::global_dof_index> &n_indices_per_block)
const
451 std::vector<IndexSet> partitioned;
452 const unsigned int n_blocks = n_indices_per_block.size();
454 partitioned.reserve(n_blocks);
456 for (
const auto n_block_indices : n_indices_per_block)
458 partitioned.push_back(this->
get_view(start, start + n_block_indices));
459 start += n_block_indices;
465 for (
const auto &partition : partitioned)
467 sum += partition.size();
487 std::vector<Range> new_ranges;
489 std::vector<Range>::iterator own_it =
ranges.begin();
490 std::vector<Range>::iterator other_it = other.
ranges.begin();
492 while (own_it !=
ranges.end() && other_it != other.
ranges.end())
495 if (own_it->end <= other_it->begin)
497 new_ranges.push_back(*own_it);
502 if (own_it->begin >= other_it->end)
510 if (own_it->begin < other_it->begin)
512 Range r(own_it->begin, other_it->begin);
514 new_ranges.push_back(r);
519 own_it->begin = other_it->end;
520 if (own_it->begin > own_it->end)
522 own_it->begin = own_it->end;
531 for (; own_it !=
ranges.end(); ++own_it)
532 new_ranges.push_back(*own_it);
537 const std::vector<Range>::iterator
end = new_ranges.end();
538 for (std::vector<Range>::iterator it = new_ranges.begin(); it !=
end; ++it)
550 for (
const auto el : *
this)
563 "pop_back() failed, because this IndexSet contains no entries."));
581 "pop_front() failed, because this IndexSet contains no entries."));
604 const auto insert_position =
606 if (insert_position ==
ranges.end() ||
607 insert_position->begin > new_range.
begin ||
608 insert_position->end < new_range.
end)
609 ranges.insert(insert_position, new_range);
616 boost::container::small_vector<std::pair<size_type, size_type>, 200>
618 const bool ranges_are_sorted)
620 if (!ranges_are_sorted)
621 std::sort(tmp_ranges.begin(), tmp_ranges.end());
631 if (tmp_ranges.size() > 9)
634 tmp_set.
ranges.reserve(tmp_ranges.size());
635 for (
const auto &i : tmp_ranges)
640 if (this->
ranges.size() <= 1)
642 if (this->
ranges.size() == 1)
644 std::swap(*
this, tmp_set);
650 for (
const auto &i : tmp_ranges)
659 if ((
this == &other) && (offset == 0))
662 if (other.
ranges.size() != 0)
670 std::vector<Range>::const_iterator r1 =
ranges.begin(),
671 r2 = other.
ranges.begin();
673 std::vector<Range> new_ranges;
680 if (r2 == other.
ranges.end() ||
681 (r1 !=
ranges.end() && r1->end < (r2->begin + offset)))
683 new_ranges.push_back(*r1);
686 else if (r1 ==
ranges.end() || (r2->end + offset) < r1->begin)
688 new_ranges.emplace_back(r2->begin + offset, r2->end + offset);
696 std::max(r1->end, r2->end + offset));
697 new_ranges.push_back(next);
714 ExcMessage(
"One index set can only be a subset of another if they "
715 "describe index spaces of the same size. The ones in "
716 "question here have sizes " +
717 std::to_string(
size()) +
" and " +
718 std::to_string(other.
size()) +
"."));
738 out <<
size() <<
" ";
739 out <<
ranges.size() << std::endl;
740 std::vector<Range>::const_iterator r =
ranges.begin();
741 for (; r !=
ranges.end(); ++r)
743 out << r->begin <<
" " << r->end << std::endl;
755 unsigned int n_ranges;
760 for (
unsigned int i = 0; i < n_ranges; ++i)
777 std::size_t n_ranges =
ranges.size();
778 out.write(
reinterpret_cast<const char *
>(&n_ranges),
sizeof(n_ranges));
779 if (
ranges.empty() ==
false)
780 out.write(
reinterpret_cast<const char *
>(&*
ranges.begin()),
789 std::size_t n_ranges;
790 in.read(
reinterpret_cast<char *
>(&
size),
sizeof(
size));
791 in.read(
reinterpret_cast<char *
>(&n_ranges),
sizeof(n_ranges));
797 in.read(
reinterpret_cast<char *
>(&*
ranges.begin()),
820 std::vector<Range>::const_iterator p = std::upper_bound(
828 return ((index >= p->begin) && (index < p->end));
836 return (p->end > index);
853 return p->begin + (n - p->nth_index_in_set);
865 std::vector<Range>::const_iterator p =
869 if (p ==
ranges.end() || p->end == n || p->begin > n)
875 return (n - p->begin) + p->nth_index_in_set;
889 std::vector<Range>::const_iterator main_range =
892 Range r(global_index, global_index + 1);
895 std::vector<Range>::const_iterator range_begin, range_end;
896 if (global_index < main_range->
begin)
898 range_begin =
ranges.begin();
899 range_end = main_range;
903 range_begin = main_range;
909 const std::vector<Range>::const_iterator p =
922 if (global_index < p->
begin)
925 return {
this,
static_cast<size_type>(p -
ranges.begin()), global_index};
930std::vector<IndexSet::size_type>
935 std::vector<size_type> indices;
938 for (
const auto &range :
ranges)
939 for (
size_type entry = range.begin; entry < range.end; ++entry)
940 indices.push_back(entry);
957#ifdef DEAL_II_WITH_TRILINOS
958# ifdef DEAL_II_TRILINOS_WITH_TPETRA
960template <
typename NodeType>
961Tpetra::Map<int, types::signed_global_dof_index, NodeType>
965 return *make_tpetra_map_rcp<NodeType>(communicator,
overlapping);
970template <
typename NodeType>
971Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index, NodeType>>
985 ExcMessage(
"You are trying to create an Tpetra::Map object "
986 "that partitions elements of an index set "
987 "between processors. However, the union of the "
988 "index sets on different processors does not "
989 "contain all indices exactly once: the sum of "
990 "the number of entries the various processors "
991 "want to store locally is " +
992 std::to_string(n_global_elements) +
993 " whereas the total size of the object to be "
995 std::to_string(
size()) +
996 ". In other words, there are "
997 "either indices that are not spoken for "
998 "by any processor, or there are indices that are "
999 "claimed by multiple processors."));
1009 Tpetra::Map<int, types::signed_global_dof_index, NodeType>>(
1013# ifdef DEAL_II_WITH_MPI
1014 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1023 std::vector<types::signed_global_dof_index> int_indices(indices.size());
1024 std::copy(indices.begin(), indices.end(), int_indices.begin());
1025 const Teuchos::ArrayView<types::signed_global_dof_index> arr_view(
1029 Tpetra::Map<int, types::signed_global_dof_index, NodeType>>(
1033# ifdef DEAL_II_WITH_MPI
1034 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1060 ExcMessage(
"You are trying to create an Epetra_Map object "
1061 "that partitions elements of an index set "
1062 "between processors. However, the union of the "
1063 "index sets on different processors does not "
1064 "contain all indices exactly once: the sum of "
1065 "the number of entries the various processors "
1066 "want to store locally is " +
1067 std::to_string(n_global_elements) +
1068 " whereas the total size of the object to be "
1070 std::to_string(
size()) +
1071 ". In other words, there are "
1072 "either indices that are not spoken for "
1073 "by any processor, or there are indices that are "
1074 "claimed by multiple processors."));
1087# ifdef DEAL_II_WITH_MPI
1088 Epetra_MpiComm(communicator)
1104# ifdef DEAL_II_WITH_MPI
1105 Epetra_MpiComm(communicator)
1115#ifdef DEAL_II_WITH_PETSC
1119 std::vector<size_type> indices;
1125 const auto local_size =
static_cast<PetscInt
>(
n_elements());
1129 std::vector<PetscInt> petsc_indices(
n_elements());
1130 for (
const auto &index : *
this)
1131 petsc_indices[i] =
static_cast<PetscInt
>(index);
1134 PetscErrorCode ierr = ISCreateGeneral(
1135 communicator, local_size, petsc_indices.data(), PETSC_COPY_VALUES, &is);
1151 if (n_global_elements !=
size())
1154 if (n_global_elements == 0)
1157#ifdef DEAL_II_WITH_MPI
1159 const bool all_contiguous =
1161 if (!all_contiguous)
1164 bool is_globally_ascending =
true;
1171 const std::vector<types::global_dof_index> global_dofs =
1181 for (; index < global_dofs.size(); ++index)
1186 if (new_dof <= old_dof)
1188 is_globally_ascending =
false;
1198 int is_ascending = is_globally_ascending ? 1 : 0;
1199 int ierr = MPI_Bcast(&is_ascending, 1, MPI_INT, 0, communicator);
1202 return (is_ascending == 1);
1222# ifdef DEAL_II_WITH_TRILINOS
1223# ifdef DEAL_II_TRILINOS_WITH_TPETRA
1226 const Teuchos::RCP<
const Tpetra::Map<
1232# if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || \
1233 defined(KOKKOS_ENABLE_SYCL)
1235 const Teuchos::RCP<
const Tpetra::Map<
1248# if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || \
1249 defined(KOKKOS_ENABLE_SYCL)
1258template Teuchos::RCP<
1265# if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || \
1266 defined(KOKKOS_ENABLE_SYCL)
1267template Teuchos::RCP<
ElementIterator begin() const
bool is_subset_of(const IndexSet &other) const
bool is_element_binary_search(const size_type local_index) const
size_type index_within_set_binary_search(const size_type global_index) const
IS make_petsc_is(const MPI_Comm communicator=MPI_COMM_WORLD) const
bool is_ascending_and_one_to_one(const MPI_Comm communicator) const
bool is_contiguous() const
Tpetra::Map< int, types::signed_global_dof_index, NodeType > make_tpetra_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
ElementIterator at(const size_type global_index) const
void fill_index_vector(std::vector< size_type > &indices) const
std::vector< IndexSet > split_by_block(const std::vector< types::global_dof_index > &n_indices_per_block) const
size_type n_elements() const
void add_range_lower_bound(const Range &range)
ElementIterator begin() const
void set_size(const size_type size)
void read(std::istream &in)
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
IndexSet tensor_product(const IndexSet &other) const
void write(std::ostream &out) const
void block_read(std::istream &in)
void add_ranges_internal(boost::container::small_vector< std::pair< size_type, size_type >, 200 > &tmp_ranges, const bool ranges_are_sorted)
IntervalIterator begin_intervals() const
std::vector< Range > ranges
void subtract_set(const IndexSet &other)
ElementIterator end() const
Threads::Mutex compress_mutex
size_type index_space_size
void block_write(std::ostream &out) const
IndexSet get_view(const size_type begin, const size_type end) const
void add_range(const size_type begin, const size_type end)
std::size_t memory_consumption() const
size_type nth_index_in_set_binary_search(const size_type local_index) const
std::vector< size_type > get_index_vector() const
Teuchos::RCP< Tpetra::Map< int, types::signed_global_dof_index, NodeType > > make_tpetra_map_rcp(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
types::global_dof_index size_type
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
IndexSet operator&(const IndexSet &is) const
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define AssertIntegerConversion(index1, index2)
#define DEAL_II_ASSERT_UNREACHABLE()
#define AssertThrowIntegerConversion(index1, index2)
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrowMPI(error_code)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
const unsigned int my_rank
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
Tpetra::Map< LO, GO, NodeType< MemorySpace > > MapType
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
T sum(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::vector< T > gather(const MPI_Comm comm, const T &object_to_send, const unsigned int root_process=0)
Teuchos::RCP< T > make_rcp(Args &&...args)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
constexpr types::global_dof_index invalid_dof_index
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
static bool end_compare(const IndexSet::Range &x, const IndexSet::Range &y)
static bool nth_index_compare(const IndexSet::Range &x, const IndexSet::Range &y)
size_type nth_index_in_set