17#ifdef DEAL_II_TRILINOS_WITH_TPETRA
34 namespace TpetraWrappers
38 template <
typename MemorySpace>
44 if (
static_cast<std::size_t
>(this->a_row) == sparsity_pattern->n_rows())
51 if (!sparsity_pattern->is_compressed())
52 sparsity_pattern->compress();
55 std::make_shared<std::vector<::types::signed_global_dof_index>>(
56 sparsity_pattern->row_length(this->a_row));
58 if (colnum_cache->size() > 0)
64 column_indices_view(colnum_cache->data(), colnum_cache->size());
66 sparsity_pattern->graph->getGlobalRowCopy(this->a_row,
80 template <
typename MemorySpace>
92 graph->fillComplete();
97 template <
typename MemorySpace>
103 reinit(m, n, n_entries_per_row);
108 template <
typename MemorySpace>
112 const std::vector<size_type> &n_entries_per_row)
114 reinit(m, n, n_entries_per_row);
119 template <
typename MemorySpace>
123 , column_space_map(std::move(other.column_space_map))
124 , graph(std::move(other.graph))
125 , nonlocal_graph(std::move(other.nonlocal_graph))
131 template <
typename MemorySpace>
139 Utilities::Trilinos::tpetra_comm_self()))
141 TpetraTypes::GraphType<
MemorySpace>>(column_space_map,
145 (void)input_sparsity;
148 "Copy constructor only works for empty sparsity patterns."));
153 template <
typename MemorySpace>
155 const IndexSet ¶llel_partitioning,
159 reinit(parallel_partitioning,
160 parallel_partitioning,
167 template <
typename MemorySpace>
169 const IndexSet ¶llel_partitioning,
171 const std::vector<size_type> &n_entries_per_row)
173 reinit(parallel_partitioning,
174 parallel_partitioning,
181 template <
typename MemorySpace>
183 const IndexSet &row_parallel_partitioning,
184 const IndexSet &col_parallel_partitioning,
188 reinit(row_parallel_partitioning,
189 col_parallel_partitioning,
196 template <
typename MemorySpace>
198 const IndexSet &row_parallel_partitioning,
199 const IndexSet &col_parallel_partitioning,
201 const std::vector<size_type> &n_entries_per_row)
203 reinit(row_parallel_partitioning,
204 col_parallel_partitioning,
211 template <
typename MemorySpace>
213 const IndexSet &row_parallel_partitioning,
214 const IndexSet &col_parallel_partitioning,
219 reinit(row_parallel_partitioning,
220 col_parallel_partitioning,
223 n_max_entries_per_row);
228 template <
typename MemorySpace>
242 template <
typename MemorySpace>
247 const std::vector<size_type> &n_entries_per_row)
257 namespace SparsityPatternImpl
259 template <
typename MemorySpace>
262 template <
typename MemorySpace>
272 Assert(row_map->isOneToOne(),
273 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
274 "the maps of different processors."));
275 Assert(col_map->isOneToOne(),
276 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
277 "the maps of different processors."));
279 nonlocal_graph.reset();
281 column_space_map = col_map;
299 template <
typename MemorySpace>
309 Assert(row_map->isOneToOne(),
310 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
311 "the maps of different processors."));
312 Assert(col_map->isOneToOne(),
313 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
314 "the maps of different processors."));
317 nonlocal_graph.reset();
320 row_map->getGlobalNumElements());
322 column_space_map = col_map;
326 Kokkos::DualView<size_t *, typename MemorySpace::kokkos_space>
327 local_entries_per_row(
"local_entries_per_row",
328 row_map->getMaxGlobalIndex() -
329 row_map->getMinGlobalIndex());
331 auto local_entries_per_row_host =
332 local_entries_per_row
333 .template view<Kokkos::DefaultHostExecutionSpace>();
335 std::uint64_t total_size = 0;
336 for (
unsigned int i = 0; i < local_entries_per_row.extent(0); ++i)
338 local_entries_per_row_host(i) =
339 n_entries_per_row[row_map->getMinGlobalIndex() + i];
340 total_size += local_entries_per_row_host[i];
342 local_entries_per_row
343 .template modify<Kokkos::DefaultHostExecutionSpace>();
344 local_entries_per_row
345 .template sync<typename MemorySpace::kokkos_space>();
348 total_size <
static_cast<std::uint64_t
>(
352 "You are requesting to store more elements than global ordinal type allows."));
357 local_entries_per_row);
362 template <
typename SparsityPatternType,
typename MemorySpace>
367 const SparsityPatternType &sp,
368 [[maybe_unused]]
const bool exchange_data,
373 nonlocal_graph.reset();
382 Assert(row_map->isContiguous() ==
true,
384 "This function only works if the row map is contiguous."));
388 row_map->getMaxGlobalIndex() + 1;
389 Teuchos::Array<size_t> n_entries_per_row(last_row - first_row);
392 n_entries_per_row[row - first_row] = sp.row_length(row);
395 std::accumulate(n_entries_per_row.begin(),
396 n_entries_per_row.end(),
398 static_cast<std::uint64_t
>(std::numeric_limits<int>::max()),
400 "The TrilinosWrappers use Tpetra internally, and "
401 "Trilinos/Tpetra was compiled with 'local ordinate = int'. "
402 "Therefore, 'signed int' is used to represent local indices, "
403 "and only 2,147,483,647 nonzero matrix entries can be stored "
404 "on a single process, but you are requesting more than "
405 "that. Either use more MPI processes or recompile Trilinos "
406 "with 'local ordinate = long long' "));
408 if (row_map->getComm()->getSize() > 1)
420 std::vector<TrilinosWrappers::types::int_type> row_indices;
429 row_indices.resize(row_length, -1);
431 typename SparsityPatternType::iterator p = sp.begin(row);
434 for (
int col = 0; col < row_length;)
436 row_indices[col++] = p->column();
437 if (col < row_length)
441 graph->insertGlobalIndices(row, row_length, row_indices.data());
444 graph->globalAssemble();
449 template <
typename MemorySpace>
456 parallel_partitioning.
size());
457 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
458 parallel_partitioning
459 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
460 communicator,
false);
461 SparsityPatternImpl::reinit_sp<MemorySpace>(
462 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
467 template <
typename MemorySpace>
470 const IndexSet ¶llel_partitioning,
472 const std::vector<size_type> &n_entries_per_row)
475 parallel_partitioning.
size());
476 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
477 parallel_partitioning
478 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
479 communicator,
false);
480 SparsityPatternImpl::reinit_sp<MemorySpace>(
481 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
486 template <
typename MemorySpace>
489 const IndexSet &row_parallel_partitioning,
490 const IndexSet &col_parallel_partitioning,
495 col_parallel_partitioning.
size());
496 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
497 row_parallel_partitioning
498 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
499 communicator,
false);
500 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
501 col_parallel_partitioning
502 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
503 communicator,
false);
504 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
514 template <
typename MemorySpace>
517 const IndexSet &row_parallel_partitioning,
518 const IndexSet &col_parallel_partitioning,
520 const std::vector<size_type> &n_entries_per_row)
523 col_parallel_partitioning.
size());
524 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
525 row_parallel_partitioning
526 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
527 communicator,
false);
528 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
529 col_parallel_partitioning
530 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
531 communicator,
false);
532 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
542 template <
typename MemorySpace>
545 const IndexSet &row_parallel_partitioning,
546 const IndexSet &col_parallel_partitioning,
552 col_parallel_partitioning.
size());
553 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
554 row_parallel_partitioning
555 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
556 communicator,
false);
557 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
558 col_parallel_partitioning
559 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
560 communicator,
false);
561 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
568 IndexSet nonlocal_partitioner = writable_rows;
570 row_parallel_partitioning.size());
574 IndexSet tmp = writable_rows & row_parallel_partitioning;
575 Assert(tmp == row_parallel_partitioning,
577 "The set of writable rows passed to this method does not "
578 "contain the locally owned rows, which is not allowed."));
581 nonlocal_partitioner.
subtract_set(row_parallel_partitioning);
584 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> nonlocal_map =
586 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
597 template <
typename MemorySpace>
598 template <
typename SparsityPatternType>
601 const IndexSet &row_parallel_partitioning,
602 const IndexSet &col_parallel_partitioning,
603 const SparsityPatternType &nontrilinos_sparsity_pattern,
605 const bool exchange_data)
608 col_parallel_partitioning.
size());
609 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
610 row_parallel_partitioning
611 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
612 communicator,
false);
613 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
614 col_parallel_partitioning
615 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
616 communicator,
false);
617 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
620 nontrilinos_sparsity_pattern,
629 template <
typename MemorySpace>
630 template <
typename SparsityPatternType>
633 const IndexSet ¶llel_partitioning,
634 const SparsityPatternType &nontrilinos_sparsity_pattern,
636 const bool exchange_data)
639 parallel_partitioning.
size());
641 parallel_partitioning.
size());
643 parallel_partitioning.
size());
644 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
645 parallel_partitioning
646 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>(
647 communicator,
false);
648 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
651 nontrilinos_sparsity_pattern,
660 template <
typename MemorySpace>
671 template <
typename MemorySpace>
686 nonlocal_graph.reset();
691 template <
typename MemorySpace>
692 template <
typename SparsityPatternType>
697 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> rows =
703 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> columns =
710 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
711 rows, columns, sp,
false, column_space_map, graph, nonlocal_graph);
716 template <
typename MemorySpace>
733 graph->fillComplete();
735 nonlocal_graph.reset();
740 template <
typename MemorySpace>
745 if (nonlocal_graph.get() !=
nullptr)
747# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
748 if (nonlocal_graph->getRowMap()->getLocalNumElements() > 0 &&
749 column_space_map->getGlobalNumElements() > 0)
751 if (nonlocal_graph->getRowMap()->getNodeNumElements() > 0 &&
752 column_space_map->getGlobalNumElements() > 0)
758 nonlocal_graph->getRowMap()->getGlobalElement(0);
763 if (column_space_map->getGlobalNumElements() ==
764 graph->getRangeMap()->getGlobalNumElements())
770# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
771 else if (column_space_map->getLocalNumElements() > 0)
773 else if (column_space_map->getNodeNumElements() > 0)
775 column = column_space_map->getGlobalElement(0);
776 nonlocal_graph->insertGlobalIndices(row, 1, &column);
778# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
779 Assert(nonlocal_graph->getRowMap()->getLocalNumElements() == 0 ||
780 column_space_map->getGlobalNumElements() == 0,
783 Assert(nonlocal_graph->getRowMap()->getNodeNumElements() == 0 ||
784 column_space_map->getGlobalNumElements() == 0,
788 nonlocal_graph->fillComplete(column_space_map, graph->getRowMap());
790 graph->fillComplete(column_space_map, graph->getRowMap());
801 template <
typename MemorySpace>
805 return graph->getRowMap()->getLocalElement(i) !=
806 Teuchos::OrdinalTraits<int>::invalid();
811 template <
typename MemorySpace>
816 if (!row_is_stored_locally(i))
820 const auto trilinos_i = graph->getRowMap()->getLocalElement(i);
821 const auto trilinos_j = graph->getColMap()->getLocalElement(j);
827 graph->getLocalRowView(trilinos_i, col_indices);
831 std::find(col_indices.data(),
832 col_indices.data() + col_indices.size(),
836 return static_cast<std::size_t
>(local_col_index) != col_indices.size();
841 template <
typename MemorySpace>
846 for (
int i = 0; i < static_cast<int>(local_size()); ++i)
848# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
852 Teuchos::ArrayView<const int> indices;
855 graph->getLocalRowView(i, indices);
856 const auto num_entries = indices.size();
857 for (
unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
874 template <
typename MemorySpace>
878# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
879 return graph->getLocalNumRows();
881 return graph->getNodeNumRows();
887 template <
typename MemorySpace>
888 std::pair<typename SparsityPattern<MemorySpace>::size_type,
892 const size_type begin = graph->getRowMap()->getMinGlobalIndex();
893 const size_type end = graph->getRowMap()->getMaxGlobalIndex() + 1;
900 template <
typename MemorySpace>
904 return graph->getGlobalNumEntries();
909 template <
typename MemorySpace>
913# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
914 return graph->getLocalMaxNumRowEntries();
916 return graph->getNodeMaxNumRowEntries();
922 template <
typename MemorySpace>
931 graph->getRowMap()->getLocalElement(row);
936 return graph->getNumEntriesInLocalRow(local_row);
943 template <
typename MemorySpace>
946 const ::types::global_dof_index &row,
948 const bool indices_are_sorted)
950 add_entries(row, columns.
begin(), columns.
end(), indices_are_sorted);
955 template <
typename MemorySpace>
956 Teuchos::RCP<const typename TpetraTypes::MapType<MemorySpace>>
959 return graph->getDomainMap();
964 template <
typename MemorySpace>
965 Teuchos::RCP<const typename TpetraTypes::MapType<MemorySpace>>
968 return graph->getRangeMap();
973 template <
typename MemorySpace>
978 graph->getRangeMap()->getComm());
983 template <
typename MemorySpace>
984 Teuchos::RCP<const Teuchos::Comm<int>>
987 return graph->getRangeMap()->getComm();
995 template <
typename MemorySpace>
999 const bool write_extended_trilinos_info)
const
1001 if (write_extended_trilinos_info)
1005# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1006 for (
unsigned int i = 0; i < graph->getLocalNumRows(); ++i)
1008 for (
unsigned int i = 0; i < graph->getNodeNumRows(); ++i)
1011# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1015 Teuchos::ArrayView<const int> indices;
1017 graph->getLocalRowView(i, indices);
1018 int num_entries = indices.size();
1019 for (
int j = 0; j < num_entries; ++j)
1020 out <<
"(" << graph->getRowMap()->getGlobalElement(i) <<
","
1021 << graph->getColMap()->getGlobalElement(indices[j]) <<
") "
1031 template <
typename MemorySpace>
1037 for (
unsigned int row = 0; row < local_size(); ++row)
1039# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1043 Teuchos::ArrayView<const int> indices;
1045 graph->getLocalRowView(row, indices);
1046 int num_entries = indices.size();
1050 const ::types::signed_global_dof_index num_entries_ =
1059 out <<
static_cast<int>(
1060 graph->getColMap()->getGlobalElement(indices[j]))
1062 << -
static_cast<int>(graph->getRowMap()->getGlobalElement(row))
936 return graph->getNumEntriesInLocalRow(local_row); {
…}
1070 template <
typename MemorySpace>
1085 const ::SparsityPattern &);
1088 const ::DynamicSparsityPattern &);
1093 const ::SparsityPattern &,
1099 const ::DynamicSparsityPattern &,
1108 const ::SparsityPattern &,
1115 const ::DynamicSparsityPattern &,
1124 const ::SparsityPattern &);
1127 const ::DynamicSparsityPattern &);
1132 const ::SparsityPattern &,
1138 const ::DynamicSparsityPattern &,
1147 const ::SparsityPattern &,
1154 const ::DynamicSparsityPattern &,
832 col_indices.data() + col_indices.size(), {
…}
790 graph->fillComplete(column_space_map, graph->getRowMap()); {
…}
758 nonlocal_graph->getRowMap()->getGlobalElement(0); {
…}
751 if (nonlocal_graph->getRowMap()->getNodeNumElements() > 0 && {
…}
660 template <
typename MemorySpace> {
…}
644 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map = {
…}
617 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>( {
…}
586 .template make_tpetra_map_rcp<TpetraTypes::NodeType<MemorySpace>>( {
…}
528 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map = {
…}
497 row_parallel_partitioning {
…}
481 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph); {
…}
449 template <
typename MemorySpace> {
…}
408 if (row_map->getComm()->getSize() > 1) {
…}
368 [[maybe_unused]]
const bool exchange_data, {
…}
348 total_size <
static_cast<std::uint64_t
>( {
…}
342 local_entries_per_row {
…}
322 column_space_map = col_map; {
…}
257 namespace SparsityPatternImpl {
…}
size_type n_elements() const
void subtract_set(const IndexSet &other)
Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > graph
Teuchos::RCP< const Teuchos::Comm< int > > get_teuchos_mpi_communicator() const
std::pair< size_type, size_type > local_range() const
Teuchos::RCP< const TpetraTypes::MapType< MemorySpace > > range_partitioner() const
void copy_from(const SparsityPattern< MemorySpace > &input_sparsity_pattern)
SparsityPattern< MemorySpace > & operator=(const SparsityPattern< MemorySpace > &input_sparsity_pattern)
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
unsigned int max_entries_per_row() const
Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > nonlocal_graph
void print_gnuplot(std::ostream &out) const
MPI_Comm get_mpi_communicator() const
Teuchos::RCP< const TpetraTypes::MapType< MemorySpace > > domain_partitioner() const
std::uint64_t n_nonzero_elements() const
size_type bandwidth() const
bool exists(const size_type i, const size_type j) const
Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > column_space_map
unsigned int local_size() const
virtual void add_row_entries(const ::types::global_dof_index &row, const ArrayView< const ::types::global_dof_index > &columns, const bool indices_are_sorted=false) override
size_type row_length(const size_type row) const
std::size_t memory_consumption() const
bool row_is_stored_locally(const size_type i) const
void reinit(const size_type m, const size_type n, const size_type n_entries_per_row=0)
virtual void resize(const size_type rows, const size_type cols)
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
IndexSet complete_index_set(const IndexSet::size_type N)
void reinit_sp(const Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > &row_map, const Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > &col_map, const size_type< MemorySpace > n_entries_per_row, Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > &column_space_map, Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > &graph, Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > &nonlocal_graph)
typename SparsityPattern< MemorySpace >::size_type size_type
Tpetra::CrsGraph< LO, GO, NodeType< MemorySpace > > GraphType
Tpetra::Map< LO, GO, NodeType< MemorySpace > > MapType
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
Teuchos::RCP< T > make_rcp(Args &&...args)
MPI_Comm teuchos_comm_to_mpi_comm(const Teuchos::RCP< const Teuchos::Comm< int > > &teuchos_comm)
const Teuchos::RCP< const Teuchos::Comm< int > > & tpetra_comm_self()
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)