19#ifdef DEAL_II_WITH_TRILINOS
27# include <Epetra_Export.h>
64 AssertThrow(
static_cast<std::vector<size_type>::size_type
>(ncols) ==
89 graph = std::make_unique<Epetra_FECrsGraph>(View,
93 graph->FillComplete();
102 reinit(m, n, n_entries_per_row);
110 const std::vector<size_type> &n_entries_per_row)
112 reinit(m, n, n_entries_per_row);
119 , column_space_map(std::move(other.column_space_map))
120 , graph(std::move(other.graph))
121 , nonlocal_graph(std::move(other.nonlocal_graph))
134 new Epetra_FECrsGraph(View, *column_space_map, *column_space_map, 0))
136 (void)input_sparsity;
139 "Copy constructor only works for empty sparsity patterns."));
148 reinit(parallel_partitioning,
149 parallel_partitioning,
157 const IndexSet & parallel_partitioning,
159 const std::vector<size_type> &n_entries_per_row)
161 reinit(parallel_partitioning,
162 parallel_partitioning,
170 const IndexSet &col_parallel_partitioning,
174 reinit(row_parallel_partitioning,
175 col_parallel_partitioning,
183 const IndexSet & row_parallel_partitioning,
184 const IndexSet & col_parallel_partitioning,
186 const std::vector<size_type> &n_entries_per_row)
188 reinit(row_parallel_partitioning,
189 col_parallel_partitioning,
197 const IndexSet &col_parallel_partitioning,
202 reinit(row_parallel_partitioning,
203 col_parallel_partitioning,
206 n_max_entries_per_row);
227 const std::vector<size_type> &n_entries_per_row)
242 reinit_sp(
const Epetra_Map & row_map,
243 const Epetra_Map & col_map,
244 const size_type n_entries_per_row,
245 std::unique_ptr<Epetra_Map> & column_space_map,
246 std::unique_ptr<Epetra_FECrsGraph> &graph,
247 std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
249 Assert(row_map.IsOneToOne(),
250 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
251 "the maps of different processors."));
252 Assert(col_map.IsOneToOne(),
253 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
254 "the maps of different processors."));
256 nonlocal_graph.reset();
258 column_space_map = std::make_unique<Epetra_Map>(col_map);
262 std::uint64_t(n_entries_per_row) <
263 static_cast<std::uint64_t
>(std::numeric_limits<int>::max()),
264 ExcMessage(
"The TrilinosWrappers use Epetra internally which "
265 "uses 'signed int' to represent local indices. "
266 "Therefore, only 2,147,483,647 nonzero matrix "
267 "entries can be stored on a single process, "
268 "but you are requesting more than that. "
269 "If possible, use more MPI processes."));
280 if (row_map.Comm().NumProc() > 1)
281 graph = std::make_unique<Epetra_FECrsGraph>(
282 Copy, row_map, n_entries_per_row,
false
290 graph = std::make_unique<Epetra_FECrsGraph>(
291 Copy, row_map, col_map, n_entries_per_row,
false);
297 reinit_sp(
const Epetra_Map & row_map,
298 const Epetra_Map & col_map,
299 const std::vector<size_type> & n_entries_per_row,
300 std::unique_ptr<Epetra_Map> & column_space_map,
301 std::unique_ptr<Epetra_FECrsGraph> &graph,
302 std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
304 Assert(row_map.IsOneToOne(),
305 ExcMessage(
"Row map must be 1-to-1, i.e., no overlap between "
306 "the maps of different processors."));
307 Assert(col_map.IsOneToOne(),
308 ExcMessage(
"Column map must be 1-to-1, i.e., no overlap between "
309 "the maps of different processors."));
312 nonlocal_graph.reset();
317 column_space_map = std::make_unique<Epetra_Map>(col_map);
321 std::vector<int> local_entries_per_row(
324 for (
unsigned int i = 0; i < local_entries_per_row.size(); ++i)
325 local_entries_per_row[i] =
328 AssertThrow(std::accumulate(local_entries_per_row.begin(),
329 local_entries_per_row.end(),
331 static_cast<std::uint64_t
>(std::numeric_limits<int>::max()),
332 ExcMessage(
"The TrilinosWrappers use Epetra internally which "
333 "uses 'signed int' to represent local indices. "
334 "Therefore, only 2,147,483,647 nonzero matrix "
335 "entries can be stored on a single process, "
336 "but you are requesting more than that. "
337 "If possible, use more MPI processes."));
339 if (row_map.Comm().NumProc() > 1)
340 graph = std::make_unique<Epetra_FECrsGraph>(
341 Copy, row_map, local_entries_per_row.data(),
false
349 graph = std::make_unique<Epetra_FECrsGraph>(
350 Copy, row_map, col_map, local_entries_per_row.data(),
false);
355 template <
typename SparsityPatternType>
357 reinit_sp(
const Epetra_Map & row_map,
358 const Epetra_Map & col_map,
359 const SparsityPatternType & sp,
360 const bool exchange_data,
361 std::unique_ptr<Epetra_Map> & column_space_map,
362 std::unique_ptr<Epetra_FECrsGraph> &graph,
363 std::unique_ptr<Epetra_CrsGraph> & nonlocal_graph)
365 nonlocal_graph.reset();
373 column_space_map = std::make_unique<Epetra_Map>(col_map);
375 Assert(row_map.LinearMap() ==
true,
377 "This function only works if the row map is contiguous."));
381 std::vector<int> n_entries_per_row(last_row - first_row);
385 for (size_type row = first_row; row < last_row; ++row)
386 n_entries_per_row[row - first_row] =
387 static_cast<int>(sp.row_length(row));
389 AssertThrow(std::accumulate(n_entries_per_row.begin(),
390 n_entries_per_row.end(),
392 static_cast<std::uint64_t
>(std::numeric_limits<int>::max()),
393 ExcMessage(
"The TrilinosWrappers use Epetra internally which "
394 "uses 'signed int' to represent local indices. "
395 "Therefore, only 2,147,483,647 nonzero matrix "
396 "entries can be stored on a single process, "
397 "but you are requesting more than that. "
398 "If possible, use more MPI processes."));
400 if (row_map.Comm().NumProc() > 1)
401 graph = std::make_unique<Epetra_FECrsGraph>(Copy,
403 n_entries_per_row.data(),
406 graph = std::make_unique<Epetra_FECrsGraph>(
407 Copy, row_map, col_map, n_entries_per_row.data(),
false);
411 std::vector<TrilinosWrappers::types::int_type> row_indices;
415 if (exchange_data ==
false)
416 for (size_type row = first_row; row < last_row; ++row)
423 row_indices.resize(row_length, -1);
425 typename SparsityPatternType::iterator p = sp.begin(row);
428 for (
int col = 0; col < row_length;)
430 row_indices[col++] = p->column();
431 if (col < row_length)
435 graph->Epetra_CrsGraph::InsertGlobalIndices(row,
440 for (size_type row = 0; row < sp.n_rows(); ++row)
447 row_indices.resize(row_length, -1);
449 typename SparsityPatternType::iterator p = sp.begin(row);
452 for (
int col = 0; col < row_length;)
454 row_indices[col++] = p->column();
455 if (col < row_length)
460 graph->InsertGlobalIndices(1,
467 const auto &range_map =
468 static_cast<const Epetra_Map &
>(graph->RangeMap());
469 int ierr = graph->GlobalAssemble(*column_space_map, range_map,
true);
472 ierr = graph->OptimizeStorage();
495 const std::vector<size_type> &n_entries_per_row)
507 const IndexSet &col_parallel_partitioning,
527 const IndexSet &col_parallel_partitioning,
529 const std::vector<size_type> &n_entries_per_row)
547 const IndexSet &col_parallel_partitioning,
563 IndexSet nonlocal_partitioner = writable_rows;
565 row_parallel_partitioning.
size());
568 IndexSet tmp = writable_rows & row_parallel_partitioning;
569 Assert(tmp == row_parallel_partitioning,
571 "The set of writable rows passed to this method does not "
572 "contain the locally owned rows, which is not allowed."));
575 nonlocal_partitioner.
subtract_set(row_parallel_partitioning);
578 Epetra_Map nonlocal_map =
581 std::make_unique<Epetra_CrsGraph>(Copy, nonlocal_map, 0);
589 template <
typename SparsityPatternType>
592 const IndexSet & row_parallel_partitioning,
593 const IndexSet & col_parallel_partitioning,
594 const SparsityPatternType &nontrilinos_sparsity_pattern,
596 const bool exchange_data)
604 nontrilinos_sparsity_pattern,
613 template <
typename SparsityPatternType>
616 const IndexSet & parallel_partitioning,
617 const SparsityPatternType &nontrilinos_sparsity_pattern,
619 const bool exchange_data)
625 nontrilinos_sparsity_pattern,
647 graph = std::make_unique<Epetra_FECrsGraph>(*sp.
graph);
657 template <
typename SparsityPatternType>
684 graph = std::make_unique<Epetra_FECrsGraph>(View,
688 graph->FillComplete();
744 const auto &range_map =
745 static_cast<const Epetra_Map &
>(
graph->RangeMap());
752 ierr =
graph->OptimizeStorage();
754 catch (
const int error_code)
759 "The Epetra_CrsGraph::OptimizeStorage() function "
760 "has thrown an error with code " +
761 std::to_string(error_code) +
762 ". You will have to look up the exact meaning of this error "
763 "in the Trilinos source code, but oftentimes, this function "
764 "throwing an error indicates that you are trying to allocate "
765 "more than 2,147,483,647 nonzero entries in the sparsity "
766 "pattern on the local process; this will not work because "
767 "Epetra indexes entries with a simple 'signed int'."));
777 return graph->RowMap().LID(
802 if (
graph->Filled() ==
false)
804 int nnz_present =
graph->NumGlobalIndices(i);
814 int ierr =
graph->ExtractGlobalRowView(trilinos_i,
819 Assert(nnz_present == nnz_extracted,
823 const std::ptrdiff_t local_col_index =
824 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
827 if (local_col_index == nnz_present)
834 int nnz_present =
graph->NumGlobalIndices(i);
842 graph->ExtractMyRowView(trilinos_i, nnz_extracted, col_indices);
846 Assert(nnz_present == nnz_extracted,
850 const std::ptrdiff_t local_col_index =
851 std::find(col_indices, col_indices + nnz_present, trilinos_j) -
854 if (local_col_index == nnz_present)
869 for (
int i = 0; i < static_cast<int>(
local_size()); ++i)
873 graph->ExtractMyRowView(i, num_entries, indices);
874 for (
unsigned int j = 0; j < static_cast<unsigned int>(num_entries);
903 if (
graph->Filled() ==
true)
923 std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
940 return static_cast<std::uint64_t
>(nnz);
948 int nnz =
graph->MaxNumIndices();
950 return static_cast<unsigned int>(nnz);
968 return graph->NumMyIndices(local_row);
979 const auto &domain_map =
980 static_cast<const Epetra_Map &
>(
graph->DomainMap());
990 const auto &range_map =
991 static_cast<const Epetra_Map &
>(
graph->RangeMap());
1000 const Epetra_MpiComm *mpi_comm =
1001 dynamic_cast<const Epetra_MpiComm *
>(&
graph->RangeMap().Comm());
1003 return mpi_comm->Comm();
1021 const bool write_extended_trilinos_info)
const
1023 if (write_extended_trilinos_info)
1030 for (
int i = 0; i <
graph->NumMyRows(); ++i)
1032 graph->ExtractMyRowView(i, num_entries, indices);
1033 for (
int j = 0; j < num_entries; ++j)
1037 <<
") " << std::endl;
1054 graph->ExtractMyRowView(row, num_entries, indices);
1058 const ::types::global_dof_index num_entries_ = num_entries;
1065 out <<
static_cast<int>(
1068 << -
static_cast<int>(
1095 const ::SparsityPattern &,
1100 const ::DynamicSparsityPattern &,
1108 const ::SparsityPattern &,
1114 const ::DynamicSparsityPattern &,
size_type n_elements() const
void subtract_set(const IndexSet &other)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
std::shared_ptr< const std::vector< size_type > > colnum_cache
SparsityPattern * sparsity_pattern
size_type row_length(const size_type row) const
unsigned int max_entries_per_row() const
size_type bandwidth() const
const_iterator end() const
MPI_Comm get_mpi_communicator() const
std::uint64_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
const Epetra_Map & domain_partitioner() const
std::pair< size_type, size_type > local_range() const
::types::global_dof_index size_type
bool is_compressed() const
unsigned int local_size() const
void copy_from(const SparsityPattern &input_sparsity_pattern)
const Epetra_Map & range_partitioner() const
const_iterator begin() const
std::size_t memory_consumption() const
SparsityPattern & operator=(const SparsityPattern &input_sparsity_pattern)
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)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcIO()
std::unique_ptr< Epetra_FECrsGraph > graph
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
static ::ExceptionBase & ExcNotImplemented()
void print_gnuplot(std::ostream &out) const
#define Assert(cond, exc)
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
#define AssertDimension(dim1, dim2)
std::unique_ptr< Epetra_Map > column_space_map
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
IndexSet complete_index_set(const IndexSet::size_type N)
types::global_dof_index size_type
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
TrilinosWrappers::types::int_type n_global_rows(const Epetra_CrsGraph &graph)
TrilinosWrappers::types::int_type min_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int64_type n_global_entries(const Epetra_CrsGraph &graph)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type max_my_gid(const Epetra_BlockMap &map)
TrilinosWrappers::types::int_type n_global_cols(const Epetra_CrsGraph &graph)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
const Epetra_Comm & comm_self()
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)