24#ifdef DEAL_II_WITH_TRILINOS
25# ifdef DEAL_II_WITH_MPI
26# include <Epetra_MpiComm.h>
28# include <Epetra_Map.h>
29# include <Epetra_SerialComm.h>
30# ifdef DEAL_II_TRILINOS_WITH_TPETRA
31# include <Tpetra_Map.hpp>
39#ifdef DEAL_II_WITH_TRILINOS
41# ifdef DEAL_II_TRILINOS_WITH_TPETRA
43template <
typename NodeType>
46 const Tpetra::Map<int, types::signed_global_dof_index, NodeType>> &map)
48 , index_space_size(1 + map->getMaxAllGlobalIndex())
49 , largest_range(
numbers::invalid_unsigned_int)
51 Assert(map->getMinAllGlobalIndex() == 0,
53 "The Tpetra::Map does not contain the global index 0, "
54 "which means some entries are not present on any processor."));
57 if (map->isContiguous())
62# if DEAL_II_TRILINOS_VERSION_GTE(13, 4, 0)
63 const size_type n_indices = map->getLocalNumElements();
65 const size_type n_indices = map->getNodeNumElements();
68 map->getMyGlobalIndices().data();
80# ifdef DEAL_II_WITH_64BIT_INDICES
84 , index_space_size(1 + map.MaxAllGID64())
85 , largest_range(
numbers::invalid_unsigned_int)
87 Assert(map.MinAllGID64() == 0,
89 "The Epetra_BlockMap does not contain the global index 0, "
90 "which means some entries are not present on any processor."));
97 const size_type n_indices = map.NumMyElements();
99 reinterpret_cast<size_type *
>(map.MyGlobalElements64());
110 : is_compressed(true)
111 , index_space_size(1 + map.MaxAllGID())
112 , largest_range(
numbers::invalid_unsigned_int)
114 Assert(map.MinAllGID() == 0,
116 "The Epetra_BlockMap does not contain the global index 0, "
117 "which means some entries are not present on any processor."));
124 const size_type n_indices = map.NumMyElements();
125 unsigned int *indices =
126 reinterpret_cast<unsigned int *
>(map.MyGlobalElements());
153 std::vector<Range>::iterator store =
ranges.begin();
154 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();)
156 std::vector<Range>::iterator next = i;
163 while (next !=
ranges.end() && (next->begin <= last_index))
165 last_index =
std::max(last_index, next->end);
171 *store =
Range(first_index, last_index);
175 if (store !=
ranges.end())
177 std::vector<Range> new_ranges(
ranges.begin(), store);
182 size_type next_index = 0, largest_range_size = 0;
183 for (std::vector<Range>::iterator i =
ranges.begin(); i !=
ranges.end();
188 i->nth_index_in_set = next_index;
189 next_index += (i->end - i->begin);
190 if (i->end - i->begin > largest_range_size)
192 largest_range_size = i->end - i->begin;
208 for (
const auto &range :
ranges)
210 ExcMessage(
"In the process of creating the current IndexSet "
211 "object, you added indices beyond the size of the index "
212 "space. Specifically, you added elements that form the "
214 std::to_string(range.begin) +
"," +
215 std::to_string(range.end) +
216 "), but the size of the index space is only " +
232 std::vector<Range>::const_iterator r1 =
ranges.begin(),
240 if (r1->end <= r2->begin)
242 else if (r2->end <= r1->begin)
247 Assert(((r1->begin <= r2->begin) && (r1->end > r2->begin)) ||
248 ((r2->begin <= r1->begin) && (r2->end > r1->begin)),
252 result.add_range(
std::max(r1->begin, r2->begin),
258 if (r1->end <= r2->end)
276 ExcMessage(
"End index needs to be larger or equal to begin index!"));
278 ExcMessage(
"You are asking for a view into an IndexSet object "
279 "that would cover the sub-range [" +
280 std::to_string(
begin) +
',' + std::to_string(
end) +
281 "). But this is not a subset of the range "
282 "of the current object, which is [0," +
283 std::to_string(
size()) +
")."));
286 std::vector<Range>::const_iterator r1 =
ranges.begin();
288 while (r1 !=
ranges.end())
290 if ((r1->end >
begin) && (r1->begin <
end))
295 else if (r1->begin >=
end)
311 ExcMessage(
"The mask must have the same size index space "
312 "as the index set it is applied to."));
324 if (mask.ranges.size() == 1)
325 return get_view(mask.ranges[0].begin, mask.ranges[0].end);
334 std::vector<Range> new_ranges;
336 std::vector<Range>::iterator own_it =
ranges.begin();
337 std::vector<Range>::iterator mask_it = mask.ranges.begin();
339 while ((own_it !=
ranges.end()) && (mask_it != mask.ranges.end()))
345 if (own_it->end <= mask_it->begin)
353 if (mask_it->end <= own_it->begin)
372 if ((own_it->begin <= mask_it->begin) && (own_it->end <= mask_it->end))
374 new_ranges.emplace_back(mask_it->begin - mask_it->nth_index_in_set,
375 own_it->end - mask_it->nth_index_in_set);
379 if ((mask_it->begin <= own_it->begin) && (mask_it->end <= own_it->end))
381 const size_type offset_within_mask_interval =
382 own_it->begin - mask_it->begin;
383 new_ranges.emplace_back(mask_it->nth_index_in_set +
384 offset_within_mask_interval,
385 mask_it->nth_index_in_set +
386 (mask_it->end - mask_it->begin));
390 if ((own_it->begin <= mask_it->begin) &&
391 (own_it->end >= mask_it->end))
393 new_ranges.emplace_back(mask_it->begin -
394 mask_it->nth_index_in_set,
395 mask_it->end - mask_it->nth_index_in_set);
399 if ((mask_it->begin <= own_it->begin) &&
400 (mask_it->end >= own_it->end))
402 const size_type offset_within_mask_interval =
403 own_it->begin - mask_it->begin;
404 new_ranges.emplace_back(mask_it->nth_index_in_set +
405 offset_within_mask_interval,
406 mask_it->nth_index_in_set +
407 offset_within_mask_interval +
408 (own_it->end - own_it->begin));
417 if (own_it->end < mask_it->end)
419 else if (mask_it->end < own_it->end)
434 for (
const auto &range : new_ranges)
435 result.
add_range(range.begin, range.end);
445 const std::vector<types::global_dof_index> &n_indices_per_block)
const
447 std::vector<IndexSet> partitioned;
448 const unsigned int n_blocks = n_indices_per_block.size();
450 partitioned.reserve(n_blocks);
452 for (
const auto n_block_indices : n_indices_per_block)
454 partitioned.push_back(this->
get_view(start, start + n_block_indices));
455 start += n_block_indices;
460 for (
const auto &partition : partitioned)
462 sum += partition.size();
482 std::vector<Range> new_ranges;
484 std::vector<Range>::iterator own_it =
ranges.begin();
485 std::vector<Range>::iterator other_it = other.
ranges.begin();
487 while (own_it !=
ranges.end() && other_it != other.
ranges.end())
490 if (own_it->end <= other_it->begin)
492 new_ranges.push_back(*own_it);
497 if (own_it->begin >= other_it->end)
505 if (own_it->begin < other_it->begin)
507 Range r(own_it->begin, other_it->begin);
509 new_ranges.push_back(r);
514 own_it->begin = other_it->end;
515 if (own_it->begin > own_it->end)
517 own_it->begin = own_it->end;
526 for (; own_it !=
ranges.end(); ++own_it)
527 new_ranges.push_back(*own_it);
532 const std::vector<Range>::iterator
end = new_ranges.end();
533 for (std::vector<Range>::iterator it = new_ranges.begin(); it !=
end; ++it)
545 for (
const auto el : *
this)
558 "pop_back() failed, because this IndexSet contains no entries."));
576 "pop_front() failed, because this IndexSet contains no entries."));
599 const auto insert_position =
601 if (insert_position ==
ranges.end() ||
602 insert_position->begin > new_range.
begin ||
603 insert_position->end < new_range.
end)
604 ranges.insert(insert_position, new_range);
611 boost::container::small_vector<std::pair<size_type, size_type>, 200>
613 const bool ranges_are_sorted)
615 if (!ranges_are_sorted)
616 std::sort(tmp_ranges.begin(), tmp_ranges.end());
626 if (tmp_ranges.size() > 9)
629 tmp_set.
ranges.reserve(tmp_ranges.size());
630 for (
const auto &i : tmp_ranges)
635 if (this->
ranges.size() <= 1)
637 if (this->
ranges.size() == 1)
639 std::swap(*
this, tmp_set);
645 for (
const auto &i : tmp_ranges)
654 if ((
this == &other) && (offset == 0))
657 if (other.
ranges.size() != 0)
665 std::vector<Range>::const_iterator r1 =
ranges.begin(),
666 r2 = other.
ranges.begin();
668 std::vector<Range> new_ranges;
675 if (r2 == other.
ranges.end() ||
676 (r1 !=
ranges.end() && r1->end < (r2->begin + offset)))
678 new_ranges.push_back(*r1);
681 else if (r1 ==
ranges.end() || (r2->end + offset) < r1->begin)
683 new_ranges.emplace_back(r2->begin + offset, r2->end + offset);
691 std::max(r1->end, r2->end + offset));
692 new_ranges.push_back(next);
709 ExcMessage(
"One index set can only be a subset of another if they "
710 "describe index spaces of the same size. The ones in "
711 "question here have sizes " +
712 std::to_string(
size()) +
" and " +
713 std::to_string(other.
size()) +
"."));
733 out <<
size() <<
" ";
734 out <<
ranges.size() << std::endl;
735 std::vector<Range>::const_iterator r =
ranges.begin();
736 for (; r !=
ranges.end(); ++r)
738 out << r->begin <<
" " << r->end << std::endl;
750 unsigned int n_ranges;
755 for (
unsigned int i = 0; i < n_ranges; ++i)
772 std::size_t n_ranges =
ranges.size();
773 out.write(
reinterpret_cast<const char *
>(&n_ranges),
sizeof(n_ranges));
774 if (
ranges.empty() ==
false)
775 out.write(
reinterpret_cast<const char *
>(&*
ranges.begin()),
784 std::size_t n_ranges;
785 in.read(
reinterpret_cast<char *
>(&
size),
sizeof(
size));
786 in.read(
reinterpret_cast<char *
>(&n_ranges),
sizeof(n_ranges));
792 in.read(
reinterpret_cast<char *
>(&*
ranges.begin()),
815 std::vector<Range>::const_iterator p = std::upper_bound(
823 return ((index >= p->begin) && (index < p->end));
831 return (p->end > index);
848 return p->begin + (n - p->nth_index_in_set);
860 std::vector<Range>::const_iterator p =
864 if (p ==
ranges.end() || p->end == n || p->begin > n)
870 return (n - p->begin) + p->nth_index_in_set;
884 std::vector<Range>::const_iterator main_range =
887 Range r(global_index, global_index + 1);
890 std::vector<Range>::const_iterator range_begin, range_end;
891 if (global_index < main_range->
begin)
893 range_begin =
ranges.begin();
894 range_end = main_range;
898 range_begin = main_range;
904 const std::vector<Range>::const_iterator p =
917 if (global_index < p->
begin)
920 return {
this,
static_cast<size_type>(p -
ranges.begin()), global_index};
925std::vector<IndexSet::size_type>
930 std::vector<size_type> indices;
933 for (
const auto &range :
ranges)
934 for (
size_type entry = range.begin; entry < range.end; ++entry)
935 indices.push_back(entry);
952#ifdef DEAL_II_WITH_TRILINOS
953# ifdef DEAL_II_TRILINOS_WITH_TPETRA
955template <
typename NodeType>
956Tpetra::Map<int, types::signed_global_dof_index, NodeType>
960 return *make_tpetra_map_rcp<NodeType>(communicator,
overlapping);
965template <
typename NodeType>
966Teuchos::RCP<Tpetra::Map<int, types::signed_global_dof_index, NodeType>>
979 ExcMessage(
"You are trying to create an Tpetra::Map object "
980 "that partitions elements of an index set "
981 "between processors. However, the union of the "
982 "index sets on different processors does not "
983 "contain all indices exactly once: the sum of "
984 "the number of entries the various processors "
985 "want to store locally is " +
986 std::to_string(n_global_elements) +
987 " whereas the total size of the object to be "
989 std::to_string(
size()) +
990 ". In other words, there are "
991 "either indices that are not spoken for "
992 "by any processor, or there are indices that are "
993 "claimed by multiple processors."));
1003 Tpetra::Map<int, types::signed_global_dof_index, NodeType>>(
1007# ifdef DEAL_II_WITH_MPI
1008 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1017 std::vector<types::signed_global_dof_index> int_indices(indices.size());
1018 std::copy(indices.begin(), indices.end(), int_indices.begin());
1019 const Teuchos::ArrayView<types::signed_global_dof_index> arr_view(
1023 Tpetra::Map<int, types::signed_global_dof_index, NodeType>>(
1027# ifdef DEAL_II_WITH_MPI
1028 Utilities::Trilinos::internal::make_rcp<Teuchos::MpiComm<int>>(
1053 ExcMessage(
"You are trying to create an Epetra_Map object "
1054 "that partitions elements of an index set "
1055 "between processors. However, the union of the "
1056 "index sets on different processors does not "
1057 "contain all indices exactly once: the sum of "
1058 "the number of entries the various processors "
1059 "want to store locally is " +
1060 std::to_string(n_global_elements) +
1061 " whereas the total size of the object to be "
1063 std::to_string(
size()) +
1064 ". In other words, there are "
1065 "either indices that are not spoken for "
1066 "by any processor, or there are indices that are "
1067 "claimed by multiple processors."));
1080# ifdef DEAL_II_WITH_MPI
1081 Epetra_MpiComm(communicator)
1097# ifdef DEAL_II_WITH_MPI
1098 Epetra_MpiComm(communicator)
1108#ifdef DEAL_II_WITH_PETSC
1112 std::vector<size_type> indices;
1118 const auto local_size =
static_cast<PetscInt
>(
n_elements());
1122 std::vector<PetscInt> petsc_indices(
n_elements());
1123 for (
const auto &index : *
this)
1124 petsc_indices[i] =
static_cast<PetscInt
>(index);
1127 PetscErrorCode ierr = ISCreateGeneral(
1128 communicator, local_size, petsc_indices.data(), PETSC_COPY_VALUES, &is);
1144 if (n_global_elements !=
size())
1147 if (n_global_elements == 0)
1150#ifdef DEAL_II_WITH_MPI
1152 const bool all_contiguous =
1154 if (!all_contiguous)
1157 bool is_globally_ascending =
true;
1164 const std::vector<types::global_dof_index> global_dofs =
1174 for (; index < global_dofs.size(); ++index)
1179 if (new_dof <= old_dof)
1181 is_globally_ascending =
false;
1191 int is_ascending = is_globally_ascending ? 1 : 0;
1192 int ierr = MPI_Bcast(&is_ascending, 1, MPI_INT, 0, communicator);
1195 return (is_ascending == 1);
1215# ifdef DEAL_II_WITH_TRILINOS
1216# ifdef DEAL_II_TRILINOS_WITH_TPETRA
1219 const Teuchos::RCP<
const Tpetra::Map<
1225# if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || \
1226 defined(KOKKOS_ENABLE_SYCL)
1228 const Teuchos::RCP<
const Tpetra::Map<
1241# if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || \
1242 defined(KOKKOS_ENABLE_SYCL)
1251template Teuchos::RCP<
1258# if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) || \
1259 defined(KOKKOS_ENABLE_SYCL)
1260template 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
#define DEAL_II_NAMESPACE_CLOSE
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)
#define AssertIntegerConversion(index1, index2)
#define DEAL_II_ASSERT_UNREACHABLE()
#define AssertThrowIntegerConversion(index1, index2)
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)
const 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