17#ifdef DEAL_II_TRILINOS_WITH_TPETRA
34 namespace TpetraWrappers
38 template <
typename MemorySpace>
44 if (
static_cast<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)
62# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
63 typename Tpetra::CrsGraph<
66 NodeType>::nonconst_global_inds_host_view_type
67 column_indices_view(colnum_cache->data(), colnum_cache->size());
69 Teuchos::ArrayView<::types::signed_global_dof_index>
70 column_indices_view(colnum_cache->data(), colnum_cache->size());
72 sparsity_pattern->graph->getGlobalRowCopy(this->a_row,
86 template <
typename MemorySpace>
89 column_space_map = Utilities::Trilinos::internal::make_rcp<MapType>(
94 Utilities::Trilinos::internal::make_rcp<GraphType>(column_space_map,
97 graph->fillComplete();
102 template <
typename MemorySpace>
108 reinit(m, n, n_entries_per_row);
113 template <
typename MemorySpace>
117 const std::vector<size_type> &n_entries_per_row)
119 reinit(m, n, n_entries_per_row);
124 template <
typename MemorySpace>
128 , column_space_map(std::move(other.column_space_map))
129 , graph(std::move(other.graph))
130 , nonlocal_graph(std::move(other.nonlocal_graph))
136 template <
typename MemorySpace>
143 Utilities::Trilinos::tpetra_comm_self()))
149 (void)input_sparsity;
152 "Copy constructor only works for empty sparsity patterns."));
157 template <
typename MemorySpace>
159 const IndexSet ¶llel_partitioning,
163 reinit(parallel_partitioning,
164 parallel_partitioning,
171 template <
typename MemorySpace>
173 const IndexSet ¶llel_partitioning,
175 const std::vector<size_type> &n_entries_per_row)
177 reinit(parallel_partitioning,
178 parallel_partitioning,
185 template <
typename MemorySpace>
187 const IndexSet &row_parallel_partitioning,
188 const IndexSet &col_parallel_partitioning,
192 reinit(row_parallel_partitioning,
193 col_parallel_partitioning,
200 template <
typename MemorySpace>
202 const IndexSet &row_parallel_partitioning,
203 const IndexSet &col_parallel_partitioning,
205 const std::vector<size_type> &n_entries_per_row)
207 reinit(row_parallel_partitioning,
208 col_parallel_partitioning,
215 template <
typename MemorySpace>
217 const IndexSet &row_parallel_partitioning,
218 const IndexSet &col_parallel_partitioning,
223 reinit(row_parallel_partitioning,
224 col_parallel_partitioning,
227 n_max_entries_per_row);
232 template <
typename MemorySpace>
246 template <
typename MemorySpace>
251 const std::vector<size_type> &n_entries_per_row)
261 namespace SparsityPatternImpl
263 template <
typename MemorySpace>
266 template <
typename MemorySpace>
272 template <
typename MemorySpace>
274 Tpetra::CrsGraph<
int,
278 template <
typename MemorySpace>
287 Assert(row_map->isOneToOne(),
288 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
289 "the maps of different processors."));
290 Assert(col_map->isOneToOne(),
291 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
292 "the maps of different processors."));
294 nonlocal_graph.reset();
296 column_space_map = col_map;
306 graph = Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
307 row_map, row_map, n_entries_per_row);
312 template <
typename MemorySpace>
321 Assert(row_map->isOneToOne(),
322 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
323 "the maps of different processors."));
324 Assert(col_map->isOneToOne(),
325 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
326 "the maps of different processors."));
329 nonlocal_graph.reset();
332 row_map->getGlobalNumElements());
334 column_space_map = col_map;
338 Kokkos::DualView<size_t *, typename MemorySpace::kokkos_space>
339 local_entries_per_row(
"local_entries_per_row",
340 row_map->getMaxGlobalIndex() -
341 row_map->getMinGlobalIndex());
343 auto local_entries_per_row_host =
344 local_entries_per_row
345 .template view<Kokkos::DefaultHostExecutionSpace>();
347 std::uint64_t total_size = 0;
348 for (
unsigned int i = 0; i < local_entries_per_row.extent(0); ++i)
350 local_entries_per_row_host(i) =
351 n_entries_per_row[row_map->getMinGlobalIndex() + i];
352 total_size += local_entries_per_row_host[i];
354 local_entries_per_row
355 .template modify<Kokkos::DefaultHostExecutionSpace>();
356 local_entries_per_row
357 .template sync<typename MemorySpace::kokkos_space>();
360 total_size <
static_cast<std::uint64_t
>(
364 "You are requesting to store more elements than global ordinal type allows."));
366 graph = Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
367 row_map, col_map, local_entries_per_row);
372 template <
typename SparsityPatternType,
typename MemorySpace>
376 const SparsityPatternType &sp,
377 [[maybe_unused]]
const bool exchange_data,
382 nonlocal_graph.reset();
389 Utilities::Trilinos::internal::make_rcp<MapType<MemorySpace>>(
394 "This function only works if the row map is contiguous."));
398 row_map->getMaxGlobalIndex() + 1;
399 Teuchos::Array<size_t> n_entries_per_row(last_row - first_row);
402 n_entries_per_row[row - first_row] = sp.row_length(row);
405 std::accumulate(n_entries_per_row.begin(),
406 n_entries_per_row.end(),
408 static_cast<std::uint64_t
>(std::numeric_limits<int>::max()),
410 "The TrilinosWrappers use Tpetra internally, and "
411 "Trilinos/Tpetra was compiled with 'local ordinate = int'. "
412 "Therefore, 'signed int' is used to represent local indices, "
413 "and only 2,147,483,647 nonzero matrix entries can be stored "
414 "on a single process, but you are requesting more than "
415 "that. Either use more MPI processes or recompile Trilinos "
416 "with 'local ordinate = long long' "));
418# if DEAL_II_TRILINOS_VERSION_GTE(12, 16, 0)
419 if (row_map->getComm()->getSize() > 1)
421 Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
422 row_map, n_entries_per_row);
425 Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
426 row_map, col_map, n_entries_per_row);
428 if (row_map->getComm()->getSize() > 1)
431 row_map, Teuchos::arcpFromArray(n_entries_per_row));
434 Utilities::Trilinos::internal::make_rcp<GraphType<MemorySpace>>(
435 row_map, col_map, Teuchos::arcpFromArray(n_entries_per_row));
442 std::vector<TrilinosWrappers::types::int_type> row_indices;
451 row_indices.resize(row_length, -1);
453 typename SparsityPatternType::iterator p = sp.begin(row);
456 for (
int col = 0; col < row_length;)
458 row_indices[col++] = p->column();
459 if (col < row_length)
463 graph->insertGlobalIndices(row, row_length, row_indices.data());
466 graph->globalAssemble();
471 template <
typename MemorySpace>
478 parallel_partitioning.
size());
479 Teuchos::RCP<MapType> map =
481 SparsityPatternImpl::reinit_sp<MemorySpace>(
482 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
487 template <
typename MemorySpace>
490 const IndexSet ¶llel_partitioning,
492 const std::vector<size_type> &n_entries_per_row)
495 parallel_partitioning.
size());
496 Teuchos::RCP<MapType> map =
498 SparsityPatternImpl::reinit_sp<MemorySpace>(
499 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
504 template <
typename MemorySpace>
507 const IndexSet &row_parallel_partitioning,
508 const IndexSet &col_parallel_partitioning,
513 col_parallel_partitioning.
size());
514 Teuchos::RCP<MapType> row_map =
515 row_parallel_partitioning.make_tpetra_map_rcp(communicator,
false);
516 Teuchos::RCP<MapType> col_map =
517 col_parallel_partitioning.make_tpetra_map_rcp(communicator,
false);
518 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
528 template <
typename MemorySpace>
532 const IndexSet &col_parallel_partitioning,
534 const std::vector<size_type> &n_entries_per_row)
537 col_parallel_partitioning.
size());
538 Teuchos::RCP<MapType> row_map =
540 Teuchos::RCP<MapType> col_map =
542 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
552 template <
typename MemorySpace>
555 const IndexSet &row_parallel_partitioning,
556 const IndexSet &col_parallel_partitioning,
562 col_parallel_partitioning.
size());
563 Teuchos::RCP<MapType> row_map =
565 Teuchos::RCP<MapType> col_map =
567 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
574 IndexSet nonlocal_partitioner = writable_rows;
576 row_parallel_partitioning.
size());
579 IndexSet tmp = writable_rows & row_parallel_partitioning;
580 Assert(tmp == row_parallel_partitioning,
582 "The set of writable rows passed to this method does not "
583 "contain the locally owned rows, which is not allowed."));
586 nonlocal_partitioner.
subtract_set(row_parallel_partitioning);
589 Teuchos::RCP<MapType> nonlocal_map =
592 Utilities::Trilinos::internal::make_rcp<GraphType>(nonlocal_map,
602 template <
typename MemorySpace>
603 template <
typename SparsityPatternType>
606 const IndexSet &row_parallel_partitioning,
607 const IndexSet &col_parallel_partitioning,
608 const SparsityPatternType &nontrilinos_sparsity_pattern,
610 const bool exchange_data)
613 col_parallel_partitioning.
size());
614 Teuchos::RCP<MapType> row_map =
616 Teuchos::RCP<MapType> col_map =
618 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
621 nontrilinos_sparsity_pattern,
630 template <
typename MemorySpace>
631 template <
typename SparsityPatternType>
634 const IndexSet ¶llel_partitioning,
635 const SparsityPatternType &nontrilinos_sparsity_pattern,
637 const bool exchange_data)
640 parallel_partitioning.
size());
642 parallel_partitioning.
size());
644 parallel_partitioning.
size());
645 Teuchos::RCP<MapType> map =
647 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
650 nontrilinos_sparsity_pattern,
659 template <
typename MemorySpace>
670 template <
typename MemorySpace>
678 graph = Utilities::Trilinos::internal::make_rcp<GraphType>(*sp.
graph);
681 nonlocal_graph = Utilities::Trilinos::internal::make_rcp<GraphType>(
684 nonlocal_graph.reset();
689 template <
typename MemorySpace>
690 template <
typename SparsityPatternType>
695 Teuchos::RCP<MapType> rows =
696 Utilities::Trilinos::internal::make_rcp<MapType>(
700 Teuchos::RCP<MapType> columns =
701 Utilities::Trilinos::internal::make_rcp<MapType>(
706 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
707 rows, columns, sp,
false, column_space_map, graph, nonlocal_graph);
712 template <
typename MemorySpace>
720 column_space_map = Utilities::Trilinos::internal::make_rcp<MapType>(
725 Utilities::Trilinos::internal::make_rcp<GraphType>(column_space_map,
728 graph->fillComplete();
730 nonlocal_graph.reset();
735 template <
typename MemorySpace>
740 if (nonlocal_graph.get() !=
nullptr)
742# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
743 if (nonlocal_graph->getRowMap()->getLocalNumElements() > 0 &&
744 column_space_map->getGlobalNumElements() > 0)
746 if (nonlocal_graph->getRowMap()->getNodeNumElements() > 0 &&
747 column_space_map->getGlobalNumElements() > 0)
753 nonlocal_graph->getRowMap()->getGlobalElement(0);
758 if (column_space_map->getGlobalNumElements() ==
759 graph->getRangeMap()->getGlobalNumElements())
765# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
766 else if (column_space_map->getLocalNumElements() > 0)
768 else if (column_space_map->getNodeNumElements() > 0)
770 column = column_space_map->getGlobalElement(0);
771 nonlocal_graph->insertGlobalIndices(row, 1, &column);
773# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
774 Assert(nonlocal_graph->getRowMap()->getLocalNumElements() == 0 ||
775 column_space_map->getGlobalNumElements() == 0,
778 Assert(nonlocal_graph->getRowMap()->getNodeNumElements() == 0 ||
779 column_space_map->getGlobalNumElements() == 0,
783 nonlocal_graph->fillComplete(column_space_map, graph->getRowMap());
785 graph->fillComplete(column_space_map, graph->getRowMap());
796 template <
typename MemorySpace>
800 return graph->getRowMap()->getLocalElement(i) !=
801 Teuchos::OrdinalTraits<int>::invalid();
806 template <
typename MemorySpace>
811# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
812 if (!row_is_stored_locally(i))
816 const auto trilinos_i = graph->getRowMap()->getLocalElement(i);
817 const auto trilinos_j = graph->getColMap()->getLocalElement(j);
819 typename GraphType::local_inds_host_view_type col_indices;
822 graph->getLocalRowView(trilinos_i, col_indices);
826 std::find(col_indices.data(),
827 col_indices.data() + col_indices.size(),
831 return static_cast<size_t>(local_col_index) != col_indices.size();
842 template <
typename MemorySpace>
847 for (
int i = 0; i < static_cast<int>(local_size()); ++i)
849# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
850 typename GraphType::local_inds_host_view_type indices;
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 SparsityPattern<MemorySpace>::MapType>
959 return graph->getDomainMap();
964 template <
typename MemorySpace>
965 Teuchos::RCP<const typename SparsityPattern<MemorySpace>::MapType>
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)
1012 typename GraphType::local_inds_host_view_type indices;
1014 Teuchos::ArrayView<const int> indices;
1016 graph->getLocalRowView(i, indices);
1017 int num_entries = indices.size();
1018 for (
int j = 0; j < num_entries; ++j)
1019 out <<
"(" << graph->getRowMap()->getGlobalElement(i) <<
","
1020 << graph->getColMap()->getGlobalElement(indices[j]) <<
") "
1030 template <
typename MemorySpace>
1036 for (
unsigned int row = 0; row < local_size(); ++row)
1038# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1039 typename GraphType::local_inds_host_view_type indices;
1041 Teuchos::ArrayView<const int> indices;
1043 graph->getLocalRowView(row, indices);
1044 int num_entries = indices.size();
1048 const ::types::signed_global_dof_index num_entries_ =
1057 out <<
static_cast<int>(
1058 graph->getColMap()->getGlobalElement(indices[j]))
1060 << -
static_cast<int>(graph->getRowMap()->getGlobalElement(row))
1068 template <
typename MemorySpace>
1083 const ::SparsityPattern &);
1086 const ::DynamicSparsityPattern &);
1091 const ::SparsityPattern &,
1097 const ::DynamicSparsityPattern &,
1106 const ::SparsityPattern &,
1113 const ::DynamicSparsityPattern &,
size_type n_elements() const
void subtract_set(const IndexSet &other)
Teuchos::RCP< Tpetra::Map< int, types::signed_global_dof_index > > make_tpetra_map_rcp(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
Teuchos::RCP< MapType > column_space_map
Teuchos::RCP< const Teuchos::Comm< int > > get_teuchos_mpi_communicator() const
Teuchos::RCP< const MapType > domain_partitioner() const
std::pair< size_type, size_type > local_range() const
Teuchos::RCP< GraphType > nonlocal_graph
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
void print_gnuplot(std::ostream &out) const
MPI_Comm get_mpi_communicator() const
Tpetra::CrsGraph< int, ::types::signed_global_dof_index, NodeType > GraphType
std::uint64_t n_nonzero_elements() const
size_type bandwidth() const
bool exists(const size_type i, const size_type j) const
Teuchos::RCP< const MapType > range_partitioner() const
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
Teuchos::RCP< GraphType > graph
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
Tpetra::Map< int, ::types::signed_global_dof_index, NodeType > MapType
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
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
IndexSet complete_index_set(const IndexSet::size_type N)
Tpetra::CrsGraph< int, ::types::signed_global_dof_index, typename SparsityPattern< MemorySpace >::NodeType > GraphType
Tpetra::Map< int, ::types::signed_global_dof_index, typename SparsityPattern< MemorySpace >::NodeType > MapType
void reinit_sp(const Teuchos::RCP< MapType< MemorySpace > > &row_map, const Teuchos::RCP< MapType< MemorySpace > > &col_map, const size_type< MemorySpace > n_entries_per_row, Teuchos::RCP< MapType< MemorySpace > > &column_space_map, Teuchos::RCP< GraphType< MemorySpace > > &graph, Teuchos::RCP< GraphType< MemorySpace > > &nonlocal_graph)
typename SparsityPattern< MemorySpace >::size_type size_type
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 > &)
int signed_global_dof_index