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)
65 column_indices_view(colnum_cache->data(), colnum_cache->size());
67 Teuchos::ArrayView<::types::signed_global_dof_index>
68 column_indices_view(colnum_cache->data(), colnum_cache->size());
70 sparsity_pattern->graph->getGlobalRowCopy(this->a_row,
84 template <
typename MemorySpace>
96 graph->fillComplete();
101 template <
typename MemorySpace>
107 reinit(m, n, n_entries_per_row);
112 template <
typename MemorySpace>
116 const std::vector<size_type> &n_entries_per_row)
118 reinit(m, n, n_entries_per_row);
123 template <
typename MemorySpace>
127 , column_space_map(std::move(other.column_space_map))
128 , graph(std::move(other.graph))
129 , nonlocal_graph(std::move(other.nonlocal_graph))
135 template <
typename MemorySpace>
143 Utilities::Trilinos::tpetra_comm_self()))
145 TpetraTypes::GraphType<
MemorySpace>>(column_space_map,
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>
276 Assert(row_map->isOneToOne(),
277 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
278 "the maps of different processors."));
279 Assert(col_map->isOneToOne(),
280 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
281 "the maps of different processors."));
283 nonlocal_graph.reset();
285 column_space_map = col_map;
303 template <
typename MemorySpace>
313 Assert(row_map->isOneToOne(),
314 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
315 "the maps of different processors."));
316 Assert(col_map->isOneToOne(),
317 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
318 "the maps of different processors."));
321 nonlocal_graph.reset();
324 row_map->getGlobalNumElements());
326 column_space_map = col_map;
330 Kokkos::DualView<size_t *, typename MemorySpace::kokkos_space>
331 local_entries_per_row(
"local_entries_per_row",
332 row_map->getMaxGlobalIndex() -
333 row_map->getMinGlobalIndex());
335 auto local_entries_per_row_host =
336 local_entries_per_row
337 .template view<Kokkos::DefaultHostExecutionSpace>();
339 std::uint64_t total_size = 0;
340 for (
unsigned int i = 0; i < local_entries_per_row.extent(0); ++i)
342 local_entries_per_row_host(i) =
343 n_entries_per_row[row_map->getMinGlobalIndex() + i];
344 total_size += local_entries_per_row_host[i];
346 local_entries_per_row
347 .template modify<Kokkos::DefaultHostExecutionSpace>();
348 local_entries_per_row
349 .template sync<typename MemorySpace::kokkos_space>();
352 total_size <
static_cast<std::uint64_t
>(
356 "You are requesting to store more elements than global ordinal type allows."));
361 local_entries_per_row);
366 template <
typename SparsityPatternType,
typename MemorySpace>
371 const SparsityPatternType &sp,
372 [[maybe_unused]]
const bool exchange_data,
377 nonlocal_graph.reset();
386 Assert(row_map->isContiguous() ==
true,
388 "This function only works if the row map is contiguous."));
392 row_map->getMaxGlobalIndex() + 1;
393 Teuchos::Array<size_t> n_entries_per_row(last_row - first_row);
396 n_entries_per_row[row - first_row] = sp.row_length(row);
399 std::accumulate(n_entries_per_row.begin(),
400 n_entries_per_row.end(),
402 static_cast<std::uint64_t
>(std::numeric_limits<int>::max()),
404 "The TrilinosWrappers use Tpetra internally, and "
405 "Trilinos/Tpetra was compiled with 'local ordinate = int'. "
406 "Therefore, 'signed int' is used to represent local indices, "
407 "and only 2,147,483,647 nonzero matrix entries can be stored "
408 "on a single process, but you are requesting more than "
409 "that. Either use more MPI processes or recompile Trilinos "
410 "with 'local ordinate = long long' "));
412# if DEAL_II_TRILINOS_VERSION_GTE(12, 16, 0)
413 if (row_map->getComm()->getSize() > 1)
422 if (row_map->getComm()->getSize() > 1)
425 row_map, Teuchos::arcpFromArray(n_entries_per_row));
429 row_map, col_map, Teuchos::arcpFromArray(n_entries_per_row));
436 std::vector<TrilinosWrappers::types::int_type> row_indices;
445 row_indices.resize(row_length, -1);
447 typename SparsityPatternType::iterator p = sp.begin(row);
450 for (
int col = 0; col < row_length;)
452 row_indices[col++] = p->column();
453 if (col < row_length)
457 graph->insertGlobalIndices(row, row_length, row_indices.data());
460 graph->globalAssemble();
465 template <
typename MemorySpace>
472 parallel_partitioning.
size());
473 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
475 SparsityPatternImpl::reinit_sp<MemorySpace>(
476 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
481 template <
typename MemorySpace>
484 const IndexSet ¶llel_partitioning,
486 const std::vector<size_type> &n_entries_per_row)
489 parallel_partitioning.
size());
490 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
492 SparsityPatternImpl::reinit_sp<MemorySpace>(
493 map, map, n_entries_per_row, column_space_map, graph, nonlocal_graph);
498 template <
typename MemorySpace>
501 const IndexSet &row_parallel_partitioning,
502 const IndexSet &col_parallel_partitioning,
507 col_parallel_partitioning.
size());
508 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
510 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
512 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
522 template <
typename MemorySpace>
525 const IndexSet &row_parallel_partitioning,
526 const IndexSet &col_parallel_partitioning,
528 const std::vector<size_type> &n_entries_per_row)
531 col_parallel_partitioning.
size());
532 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
534 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
536 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
546 template <
typename MemorySpace>
549 const IndexSet &row_parallel_partitioning,
550 const IndexSet &col_parallel_partitioning,
556 col_parallel_partitioning.
size());
557 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
559 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
561 SparsityPatternImpl::reinit_sp<MemorySpace>(row_map,
568 IndexSet nonlocal_partitioner = writable_rows;
570 row_parallel_partitioning.size());
573 IndexSet tmp = writable_rows & row_parallel_partitioning;
574 Assert(tmp == row_parallel_partitioning,
576 "The set of writable rows passed to this method does not "
577 "contain the locally owned rows, which is not allowed."));
580 nonlocal_partitioner.
subtract_set(row_parallel_partitioning);
583 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> nonlocal_map =
594 template <
typename MemorySpace>
595 template <
typename SparsityPatternType>
598 const IndexSet &row_parallel_partitioning,
599 const IndexSet &col_parallel_partitioning,
600 const SparsityPatternType &nontrilinos_sparsity_pattern,
602 const bool exchange_data)
605 col_parallel_partitioning.
size());
606 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> row_map =
608 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> col_map =
610 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
613 nontrilinos_sparsity_pattern,
622 template <
typename MemorySpace>
623 template <
typename SparsityPatternType>
626 const IndexSet ¶llel_partitioning,
627 const SparsityPatternType &nontrilinos_sparsity_pattern,
629 const bool exchange_data)
632 parallel_partitioning.
size());
634 parallel_partitioning.
size());
636 parallel_partitioning.
size());
637 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> map =
639 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
642 nontrilinos_sparsity_pattern,
651 template <
typename MemorySpace>
662 template <
typename MemorySpace>
677 nonlocal_graph.reset();
682 template <
typename MemorySpace>
683 template <
typename SparsityPatternType>
688 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> rows =
694 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> columns =
701 SparsityPatternImpl::reinit_sp<SparsityPatternType, MemorySpace>(
702 rows, columns, sp,
false, column_space_map, graph, nonlocal_graph);
707 template <
typename MemorySpace>
724 graph->fillComplete();
726 nonlocal_graph.reset();
731 template <
typename MemorySpace>
736 if (nonlocal_graph.get() !=
nullptr)
738# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
739 if (nonlocal_graph->getRowMap()->getLocalNumElements() > 0 &&
740 column_space_map->getGlobalNumElements() > 0)
742 if (nonlocal_graph->getRowMap()->getNodeNumElements() > 0 &&
743 column_space_map->getGlobalNumElements() > 0)
749 nonlocal_graph->getRowMap()->getGlobalElement(0);
754 if (column_space_map->getGlobalNumElements() ==
755 graph->getRangeMap()->getGlobalNumElements())
761# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
762 else if (column_space_map->getLocalNumElements() > 0)
764 else if (column_space_map->getNodeNumElements() > 0)
766 column = column_space_map->getGlobalElement(0);
767 nonlocal_graph->insertGlobalIndices(row, 1, &column);
769# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
770 Assert(nonlocal_graph->getRowMap()->getLocalNumElements() == 0 ||
771 column_space_map->getGlobalNumElements() == 0,
774 Assert(nonlocal_graph->getRowMap()->getNodeNumElements() == 0 ||
775 column_space_map->getGlobalNumElements() == 0,
779 nonlocal_graph->fillComplete(column_space_map, graph->getRowMap());
781 graph->fillComplete(column_space_map, graph->getRowMap());
792 template <
typename MemorySpace>
796 return graph->getRowMap()->getLocalElement(i) !=
797 Teuchos::OrdinalTraits<int>::invalid();
802 template <
typename MemorySpace>
807# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
808 if (!row_is_stored_locally(i))
812 const auto trilinos_i = graph->getRowMap()->getLocalElement(i);
813 const auto trilinos_j = graph->getColMap()->getLocalElement(j);
819 graph->getLocalRowView(trilinos_i, col_indices);
823 std::find(col_indices.data(),
824 col_indices.data() + col_indices.size(),
828 return static_cast<size_t>(local_col_index) != col_indices.size();
839 template <
typename MemorySpace>
844 for (
int i = 0; i < static_cast<int>(local_size()); ++i)
846# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
850 Teuchos::ArrayView<const int> indices;
853 graph->getLocalRowView(i, indices);
854 const auto num_entries = indices.size();
855 for (
unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
872 template <
typename MemorySpace>
876# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
877 return graph->getLocalNumRows();
879 return graph->getNodeNumRows();
885 template <
typename MemorySpace>
886 std::pair<typename SparsityPattern<MemorySpace>::size_type,
890 const size_type begin = graph->getRowMap()->getMinGlobalIndex();
891 const size_type end = graph->getRowMap()->getMaxGlobalIndex() + 1;
898 template <
typename MemorySpace>
902 return graph->getGlobalNumEntries();
907 template <
typename MemorySpace>
911# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
912 return graph->getLocalMaxNumRowEntries();
914 return graph->getNodeMaxNumRowEntries();
920 template <
typename MemorySpace>
929 graph->getRowMap()->getLocalElement(row);
934 return graph->getNumEntriesInLocalRow(local_row);
941 template <
typename MemorySpace>
944 const ::types::global_dof_index &row,
946 const bool indices_are_sorted)
948 add_entries(row, columns.
begin(), columns.
end(), indices_are_sorted);
953 template <
typename MemorySpace>
954 Teuchos::RCP<const typename TpetraTypes::MapType<MemorySpace>>
957 return graph->getDomainMap();
962 template <
typename MemorySpace>
963 Teuchos::RCP<const typename TpetraTypes::MapType<MemorySpace>>
966 return graph->getRangeMap();
971 template <
typename MemorySpace>
976 graph->getRangeMap()->getComm());
981 template <
typename MemorySpace>
982 Teuchos::RCP<const Teuchos::Comm<int>>
985 return graph->getRangeMap()->getComm();
993 template <
typename MemorySpace>
997 const bool write_extended_trilinos_info)
const
999 if (write_extended_trilinos_info)
1003# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1004 for (
unsigned int i = 0; i < graph->getLocalNumRows(); ++i)
1006 for (
unsigned int i = 0; i < graph->getNodeNumRows(); ++i)
1009# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1013 Teuchos::ArrayView<const int> indices;
1015 graph->getLocalRowView(i, indices);
1016 int num_entries = indices.size();
1017 for (
int j = 0; j < num_entries; ++j)
1018 out <<
"(" << graph->getRowMap()->getGlobalElement(i) <<
","
1019 << graph->getColMap()->getGlobalElement(indices[j]) <<
") "
1029 template <
typename MemorySpace>
1035 for (
unsigned int row = 0; row < local_size(); ++row)
1037# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
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
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
#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)
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 > &)