Reference documentation for deal.II version GIT 4e92747666 2022-12-02 08:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_sparsity_pattern.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_sparsity_pattern_h
17 #define dealii_trilinos_sparsity_pattern_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_WITH_TRILINOS
22 
23 # include <deal.II/base/index_set.h>
24 # include <deal.II/base/mpi_stub.h>
26 
27 # include <deal.II/lac/exceptions.h>
29 
30 # include <Epetra_FECrsGraph.h>
31 # include <Epetra_Map.h>
32 # include <Epetra_MpiComm.h>
33 
34 # include <cmath>
35 # include <memory>
36 # include <vector>
37 
38 
40 
41 // forward declarations
42 # ifndef DOXYGEN
43 class SparsityPattern;
45 
46 namespace TrilinosWrappers
47 {
48  class SparsityPattern;
49  class SparseMatrix;
50 
51  namespace SparsityPatternIterators
52  {
53  class Iterator;
54  }
55 } // namespace TrilinosWrappers
56 # endif
57 
58 namespace TrilinosWrappers
59 {
61  {
73  class Accessor
74  {
75  public:
80 
85  const size_type row,
86  const size_type index);
87 
91  size_type
92  row() const;
93 
97  size_type
98  index() const;
99 
103  size_type
104  column() const;
105 
110 
115  size_type,
116  size_type,
117  size_type,
118  << "You tried to access row " << arg1
119  << " of a distributed sparsity pattern, "
120  << " but only rows " << arg2 << " through " << arg3
121  << " are stored locally and can be accessed.");
122 
123  private:
128 
133 
138 
151  std::shared_ptr<const std::vector<size_type>> colnum_cache;
152 
158  void
160 
161  // Make enclosing class a friend.
162  friend class Iterator;
163  };
164 
170  class Iterator
171  {
172  public:
177 
182  Iterator(const SparsityPattern *sparsity_pattern,
183  const size_type row,
184  const size_type index);
185 
189  Iterator(const Iterator &i);
190 
194  Iterator &
196 
200  Iterator
202 
206  const Accessor &
207  operator*() const;
208 
212  const Accessor *
213  operator->() const;
214 
219  bool
220  operator==(const Iterator &) const;
221 
225  bool
226  operator!=(const Iterator &) const;
227 
233  bool
234  operator<(const Iterator &) const;
235 
240  size_type,
241  size_type,
242  << "Attempt to access element " << arg2 << " of row "
243  << arg1 << " which doesn't have that many elements.");
244 
245  private:
250 
252  };
253 
254  } // namespace SparsityPatternIterators
255 
256 
275  {
276  public:
281 
286 
294  SparsityPattern();
295 
310  SparsityPattern(const size_type m,
311  const size_type n,
312  const size_type n_entries_per_row = 0);
313 
322  SparsityPattern(const size_type m,
323  const size_type n,
324  const std::vector<size_type> &n_entries_per_row);
325 
330  SparsityPattern(SparsityPattern &&other) noexcept;
331 
336  SparsityPattern(const SparsityPattern &input_sparsity_pattern);
337 
341  virtual ~SparsityPattern() override = default;
342 
354  void
355  reinit(const size_type m,
356  const size_type n,
357  const size_type n_entries_per_row = 0);
358 
367  void
368  reinit(const size_type m,
369  const size_type n,
370  const std::vector<size_type> &n_entries_per_row);
371 
376  void
377  copy_from(const SparsityPattern &input_sparsity_pattern);
378 
384  template <typename SparsityPatternType>
385  void
386  copy_from(const SparsityPatternType &nontrilinos_sparsity_pattern);
387 
395  operator=(const SparsityPattern &input_sparsity_pattern);
396 
404  void
405  clear();
406 
416  void
417  compress();
436  SparsityPattern(const IndexSet &parallel_partitioning,
437  const MPI_Comm &communicator = MPI_COMM_WORLD,
438  const size_type n_entries_per_row = 0);
439 
450  SparsityPattern(const IndexSet & parallel_partitioning,
451  const MPI_Comm & communicator,
452  const std::vector<size_type> &n_entries_per_row);
453 
468  SparsityPattern(const IndexSet &row_parallel_partitioning,
469  const IndexSet &col_parallel_partitioning,
470  const MPI_Comm &communicator = MPI_COMM_WORLD,
471  const size_type n_entries_per_row = 0);
472 
484  SparsityPattern(const IndexSet & row_parallel_partitioning,
485  const IndexSet & col_parallel_partitioning,
486  const MPI_Comm & communicator,
487  const std::vector<size_type> &n_entries_per_row);
488 
515  SparsityPattern(const IndexSet &row_parallel_partitioning,
516  const IndexSet &col_parallel_partitioning,
517  const IndexSet &writable_rows,
518  const MPI_Comm &communicator = MPI_COMM_WORLD,
519  const size_type n_entries_per_row = 0);
520 
536  void
537  reinit(const IndexSet &parallel_partitioning,
538  const MPI_Comm &communicator = MPI_COMM_WORLD,
539  const size_type n_entries_per_row = 0);
540 
551  void
552  reinit(const IndexSet & parallel_partitioning,
553  const MPI_Comm & communicator,
554  const std::vector<size_type> &n_entries_per_row);
555 
572  void
573  reinit(const IndexSet &row_parallel_partitioning,
574  const IndexSet &col_parallel_partitioning,
575  const MPI_Comm &communicator = MPI_COMM_WORLD,
576  const size_type n_entries_per_row = 0);
577 
603  void
604  reinit(const IndexSet &row_parallel_partitioning,
605  const IndexSet &col_parallel_partitioning,
606  const IndexSet &writeable_rows,
607  const MPI_Comm &communicator = MPI_COMM_WORLD,
608  const size_type n_entries_per_row = 0);
609 
614  void
615  reinit(const IndexSet & row_parallel_partitioning,
616  const IndexSet & col_parallel_partitioning,
617  const MPI_Comm & communicator,
618  const std::vector<size_type> &n_entries_per_row);
619 
629  template <typename SparsityPatternType>
630  void
631  reinit(const IndexSet & row_parallel_partitioning,
632  const IndexSet & col_parallel_partitioning,
633  const SparsityPatternType &nontrilinos_sparsity_pattern,
634  const MPI_Comm & communicator = MPI_COMM_WORLD,
635  const bool exchange_data = false);
636 
645  template <typename SparsityPatternType>
646  void
647  reinit(const IndexSet & parallel_partitioning,
648  const SparsityPatternType &nontrilinos_sparsity_pattern,
649  const MPI_Comm & communicator = MPI_COMM_WORLD,
650  const bool exchange_data = false);
661  bool
662  is_compressed() const;
663 
667  unsigned int
668  max_entries_per_row() const;
669 
679  unsigned int
680  local_size() const;
681 
690  std::pair<size_type, size_type>
691  local_range() const;
692 
697  bool
698  in_local_range(const size_type index) const;
699 
703  std::uint64_t
704  n_nonzero_elements() const;
705 
714  size_type
715  row_length(const size_type row) const;
716 
723  size_type
724  bandwidth() const;
725 
730  bool
731  empty() const;
732 
737  bool
738  exists(const size_type i, const size_type j) const;
739 
744  bool
745  row_is_stored_locally(const size_type i) const;
746 
751  std::size_t
752  memory_consumption() const;
753 
762  void
763  add(const size_type i, const size_type j);
764 
765 
769  template <typename ForwardIterator>
770  void
772  ForwardIterator begin,
773  ForwardIterator end,
774  const bool indices_are_sorted = false);
775 
776  virtual void
777  add_row_entries(const size_type & row,
778  const ArrayView<const size_type> &columns,
779  const bool indices_are_sorted = false) override;
780 
782 
793  const Epetra_FECrsGraph &
795 
802  const Epetra_Map &
803  domain_partitioner() const;
804 
811  const Epetra_Map &
812  range_partitioner() const;
813 
817  MPI_Comm
818  get_mpi_communicator() const;
831  IndexSet
833 
839  IndexSet
841 
853  begin() const;
854 
859  end() const;
860 
870  begin(const size_type r) const;
871 
881  end(const size_type r) const;
882 
894  void
895  write_ascii();
896 
904  void
905  print(std::ostream &out,
906  const bool write_extended_trilinos_info = false) const;
907 
922  void
923  print_gnuplot(std::ostream &out) const;
924 
934  int,
935  << "An error with error number " << arg1
936  << " occurred while calling a Trilinos function");
937 
942  size_type,
943  size_type,
944  << "The entry with index <" << arg1 << ',' << arg2
945  << "> does not exist.");
946 
951  size_type,
952  size_type,
953  size_type,
954  size_type,
955  << "You tried to access element (" << arg1 << '/' << arg2
956  << ')'
957  << " of a distributed matrix, but only rows in range ["
958  << arg3 << ',' << arg4
959  << "] are stored locally and can be accessed.");
960 
965  size_type,
966  size_type,
967  << "You tried to access element (" << arg1 << '/' << arg2
968  << ')' << " of a sparse matrix, but it appears to not"
969  << " exist in the Trilinos sparsity pattern.");
971  private:
976  std::unique_ptr<Epetra_Map> column_space_map;
977 
983  std::unique_ptr<Epetra_FECrsGraph> graph;
984 
991  std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;
992 
996  };
997 
998 
999 
1000  // ----------------------- inline and template functions --------------------
1001 
1002 
1003 # ifndef DOXYGEN
1004 
1005  namespace SparsityPatternIterators
1006  {
1007  inline Accessor::Accessor(const SparsityPattern *sp,
1008  const size_type row,
1009  const size_type index)
1010  : sparsity_pattern(const_cast<SparsityPattern *>(sp))
1011  , a_row(row)
1012  , a_index(index)
1013  {
1014  visit_present_row();
1015  }
1016 
1017 
1018 
1019  inline Accessor::size_type
1020  Accessor::row() const
1021  {
1022  Assert(a_row < sparsity_pattern->n_rows(),
1023  ExcBeyondEndOfSparsityPattern());
1024  return a_row;
1025  }
1026 
1027 
1028 
1029  inline Accessor::size_type
1030  Accessor::column() const
1031  {
1032  Assert(a_row < sparsity_pattern->n_rows(),
1033  ExcBeyondEndOfSparsityPattern());
1034  return (*colnum_cache)[a_index];
1035  }
1036 
1037 
1038 
1039  inline Accessor::size_type
1040  Accessor::index() const
1041  {
1042  Assert(a_row < sparsity_pattern->n_rows(),
1043  ExcBeyondEndOfSparsityPattern());
1044  return a_index;
1045  }
1046 
1047 
1048 
1049  inline Iterator::Iterator(const SparsityPattern *sp,
1050  const size_type row,
1051  const size_type index)
1052  : accessor(sp, row, index)
1053  {}
1054 
1055 
1056 
1057  inline Iterator::Iterator(const Iterator &) = default;
1058 
1059 
1060 
1061  inline Iterator &
1063  {
1064  Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1065  ExcIteratorPastEnd());
1066 
1067  ++accessor.a_index;
1068 
1069  // If at end of line: do one step, then cycle until we find a row with a
1070  // nonzero number of entries that is stored locally.
1071  if (accessor.a_index >= accessor.colnum_cache->size())
1072  {
1073  accessor.a_index = 0;
1074  ++accessor.a_row;
1075 
1076  while (accessor.a_row < accessor.sparsity_pattern->n_rows())
1077  {
1078  const auto row_length =
1079  accessor.sparsity_pattern->row_length(accessor.a_row);
1080  if (row_length == 0 ||
1081  !accessor.sparsity_pattern->row_is_stored_locally(
1082  accessor.a_row))
1083  ++accessor.a_row;
1084  else
1085  break;
1086  }
1087 
1088  accessor.visit_present_row();
1089  }
1090  return *this;
1091  }
1092 
1093 
1094 
1095  inline Iterator
1097  {
1098  const Iterator old_state = *this;
1099  ++(*this);
1100  return old_state;
1101  }
1102 
1103 
1104 
1105  inline const Accessor &
1106  Iterator::operator*() const
1107  {
1108  return accessor;
1109  }
1110 
1111 
1112 
1113  inline const Accessor *
1114  Iterator::operator->() const
1115  {
1116  return &accessor;
1117  }
1118 
1119 
1120 
1121  inline bool
1122  Iterator::operator==(const Iterator &other) const
1123  {
1124  return (accessor.a_row == other.accessor.a_row &&
1125  accessor.a_index == other.accessor.a_index);
1126  }
1127 
1128 
1129 
1130  inline bool
1131  Iterator::operator!=(const Iterator &other) const
1132  {
1133  return !(*this == other);
1134  }
1135 
1136 
1137 
1138  inline bool
1139  Iterator::operator<(const Iterator &other) const
1140  {
1141  return (accessor.row() < other.accessor.row() ||
1142  (accessor.row() == other.accessor.row() &&
1143  accessor.index() < other.accessor.index()));
1144  }
1145 
1146  } // namespace SparsityPatternIterators
1147 
1148 
1149 
1151  SparsityPattern::begin() const
1152  {
1153  const size_type first_valid_row = this->local_range().first;
1154  return const_iterator(this, first_valid_row, 0);
1155  }
1156 
1157 
1158 
1160  SparsityPattern::end() const
1161  {
1162  return const_iterator(this, n_rows(), 0);
1163  }
1164 
1165 
1166 
1168  SparsityPattern::begin(const size_type r) const
1169  {
1170  AssertIndexRange(r, n_rows());
1171  if (row_length(r) > 0)
1172  return const_iterator(this, r, 0);
1173  else
1174  return end(r);
1175  }
1176 
1177 
1178 
1180  SparsityPattern::end(const size_type r) const
1181  {
1182  AssertIndexRange(r, n_rows());
1183 
1184  // place the iterator on the first entry
1185  // past this line, or at the end of the
1186  // matrix
1187  for (size_type i = r + 1; i < n_rows(); ++i)
1188  if (row_length(i) > 0)
1189  return const_iterator(this, i, 0);
1190 
1191  // if there is no such line, then take the
1192  // end iterator of the matrix
1193  return end();
1194  }
1195 
1196 
1197 
1198  inline bool
1199  SparsityPattern::in_local_range(const size_type index) const
1200  {
1202 # ifndef DEAL_II_WITH_64BIT_INDICES
1203  begin = graph->RowMap().MinMyGID();
1204  end = graph->RowMap().MaxMyGID() + 1;
1205 # else
1206  begin = graph->RowMap().MinMyGID64();
1207  end = graph->RowMap().MaxMyGID64() + 1;
1208 # endif
1209 
1210  return ((index >= static_cast<size_type>(begin)) &&
1211  (index < static_cast<size_type>(end)));
1212  }
1213 
1214 
1215 
1216  inline bool
1218  {
1219  return graph->Filled();
1220  }
1221 
1222 
1223 
1224  inline bool
1225  SparsityPattern::empty() const
1226  {
1227  return ((n_rows() == 0) && (n_cols() == 0));
1228  }
1229 
1230 
1231 
1232  inline void
1233  SparsityPattern::add(const size_type i, const size_type j)
1234  {
1235  add_entries(i, &j, &j + 1);
1236  }
1237 
1238 
1239 
1240  template <typename ForwardIterator>
1241  inline void
1243  ForwardIterator begin,
1244  ForwardIterator end,
1245  const bool /*indices_are_sorted*/)
1246  {
1247  if (begin == end)
1248  return;
1249 
1250  // verify that the size of the data type Trilinos expects matches that the
1251  // iterator points to. we allow for some slippage between signed and
1252  // unsigned and only compare that they are both either 32 or 64 bit. to
1253  // write this test properly, not that we cannot compare the size of
1254  // '*begin' because 'begin' may be an iterator and '*begin' may be an
1255  // accessor class. consequently, we need to somehow get an actual value
1256  // from it which we can by evaluating an expression such as when
1257  // multiplying the value produced by 2
1258  Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1259  ExcNotImplemented());
1260 
1261  TrilinosWrappers::types::int_type *col_index_ptr =
1262  reinterpret_cast<TrilinosWrappers::types::int_type *>(
1263  const_cast<typename std::decay<decltype(*begin)>::type *>(&*begin));
1264  // Check at least for the first index that the conversion actually works
1265  AssertDimension(*col_index_ptr, *begin);
1266  TrilinosWrappers::types::int_type trilinos_row_index = row;
1267  const int n_cols = static_cast<int>(end - begin);
1268 
1269  int ierr;
1270  if (row_is_stored_locally(row))
1271  ierr =
1272  graph->InsertGlobalIndices(trilinos_row_index, n_cols, col_index_ptr);
1273  else if (nonlocal_graph.get() != nullptr)
1274  {
1275  // this is the case when we have explicitly set the off-processor rows
1276  // and want to create a separate matrix object for them (to retain
1277  // thread-safety)
1278  Assert(nonlocal_graph->RowMap().LID(
1279  static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1280  ExcMessage("Attempted to write into off-processor matrix row "
1281  "that has not be specified as being writable upon "
1282  "initialization"));
1283  ierr = nonlocal_graph->InsertGlobalIndices(trilinos_row_index,
1284  n_cols,
1285  col_index_ptr);
1286  }
1287  else
1288  ierr = graph->InsertGlobalIndices(1,
1289  &trilinos_row_index,
1290  n_cols,
1291  col_index_ptr);
1292 
1293  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1294  }
1295 
1296 
1297 
1298  inline const Epetra_FECrsGraph &
1299  SparsityPattern::trilinos_sparsity_pattern() const
1300  {
1301  return *graph;
1302  }
1303 
1304 
1305 
1306  inline IndexSet
1307  SparsityPattern::locally_owned_domain_indices() const
1308  {
1309  return IndexSet(graph->DomainMap());
1310  }
1311 
1312 
1313 
1314  inline IndexSet
1315  SparsityPattern::locally_owned_range_indices() const
1316  {
1317  return IndexSet(graph->RangeMap());
1318  }
1319 
1320 # endif // DOXYGEN
1321 } // namespace TrilinosWrappers
1322 
1323 
1325 
1326 
1327 #endif // DEAL_II_WITH_TRILINOS
1328 
1329 #endif
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
size_type n_cols() const
bool is_compressed() const
iterator begin() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
void add(const size_type i, const size_type j)
SparsityPatternIterators::Iterator const_iterator
bool empty() const
iterator end() const
types::global_dof_index size_type
unsigned int row_length(const size_type row) const
std::shared_ptr< const std::vector< size_type > > colnum_cache
Accessor(const SparsityPattern *sparsity_pattern, const size_type row, const size_type index)
Iterator(const SparsityPattern *sparsity_pattern, const size_type row, const size_type index)
IndexSet locally_owned_domain_indices() const
size_type row_length(const size_type row) const
std::unique_ptr< Epetra_FECrsGraph > graph
void add(const size_type i, const size_type j)
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void print_gnuplot(std::ostream &out) const
const_iterator end() const
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
const_iterator begin(const size_type r) const
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
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
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
void copy_from(const SparsityPattern &input_sparsity_pattern)
IndexSet locally_owned_range_indices() const
std::unique_ptr< Epetra_Map > column_space_map
const_iterator begin() const
const_iterator end(const size_type r) const
bool in_local_range(const size_type index) const
const Epetra_FECrsGraph & trilinos_sparsity_pattern() const
virtual ~SparsityPattern() override=default
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)
std::enable_if_t< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
static ::ExceptionBase & ExcAccessToNonlocalRow(size_type arg1, size_type arg2, size_type arg3)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcBeyondEndOfSparsityPattern()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
types::global_dof_index size_type
Definition: cuda_kernels.h:45
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
unsigned int global_dof_index
Definition: types.h:81
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)