Reference documentation for deal.II version Git 8ba89c874e 2021-08-04 22:06:33 -0600
\(\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 - 2021 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 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_TRILINOS
23 
24 # include <deal.II/base/index_set.h>
26 
27 # include <deal.II/lac/exceptions.h>
28 
29 # include <Epetra_FECrsGraph.h>
30 # include <Epetra_Map.h>
31 # include <Epetra_MpiComm.h>
32 # include <mpi.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 
84  Accessor(const SparsityPattern *sparsity_pattern,
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 
109  DeclException0(ExcBeyondEndOfSparsityPattern);
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
159  visit_present_row();
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 &
195  operator++();
196 
200  Iterator
201  operator++(int);
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 
239  DeclException2(ExcInvalidIndexWithinRow,
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();
419 
424 
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);
652 
656 
661  bool
662  is_compressed() const;
663 
667  unsigned int
668  max_entries_per_row() const;
669 
673  size_type
674  n_rows() const;
675 
679  size_type
680  n_cols() const;
681 
691  unsigned int
692  local_size() const;
693 
702  std::pair<size_type, size_type>
703  local_range() const;
704 
709  bool
710  in_local_range(const size_type index) const;
711 
715  size_type
716  n_nonzero_elements() const;
717 
726  size_type
727  row_length(const size_type row) const;
728 
735  size_type
736  bandwidth() const;
737 
742  bool
743  empty() const;
744 
749  bool
750  exists(const size_type i, const size_type j) const;
751 
756  bool
757  row_is_stored_locally(const size_type i) const;
758 
763  std::size_t
764  memory_consumption() const;
765 
767 
774  void
775  add(const size_type i, const size_type j);
776 
777 
781  template <typename ForwardIterator>
782  void
783  add_entries(const size_type row,
784  ForwardIterator begin,
785  ForwardIterator end,
786  const bool indices_are_sorted = false);
788 
792 
797  const Epetra_FECrsGraph &
798  trilinos_sparsity_pattern() const;
799 
806  const Epetra_Map &
807  domain_partitioner() const;
808 
815  const Epetra_Map &
816  range_partitioner() const;
817 
821  MPI_Comm
822  get_mpi_communicator() const;
824 
829 
835  IndexSet
836  locally_owned_domain_indices() const;
837 
843  IndexSet
844  locally_owned_range_indices() const;
845 
847 
852 
857  begin() const;
858 
863  end() const;
864 
874  begin(const size_type r) const;
875 
885  end(const size_type r) const;
886 
888 
892 
898  void
899  write_ascii();
900 
908  void
909  print(std::ostream &out,
910  const bool write_extended_trilinos_info = false) const;
911 
926  void
927  print_gnuplot(std::ostream &out) const;
928 
930 
938  int,
939  << "An error with error number " << arg1
940  << " occurred while calling a Trilinos function");
941 
945  DeclException2(ExcInvalidIndex,
946  size_type,
947  size_type,
948  << "The entry with index <" << arg1 << ',' << arg2
949  << "> does not exist.");
950 
955  ExcSourceEqualsDestination,
956  "You are attempting an operation on two sparsity patterns that "
957  "are the same object, but the operation requires that the "
958  "two objects are in fact different.");
959 
963  DeclException4(ExcAccessToNonLocalElement,
964  size_type,
965  size_type,
966  size_type,
967  size_type,
968  << "You tried to access element (" << arg1 << "/" << arg2
969  << ")"
970  << " of a distributed matrix, but only rows in range ["
971  << arg3 << "," << arg4
972  << "] are stored locally and can be accessed.");
973 
977  DeclException2(ExcAccessToNonPresentElement,
978  size_type,
979  size_type,
980  << "You tried to access element (" << arg1 << "/" << arg2
981  << ")"
982  << " of a sparse matrix, but it appears to not"
983  << " exist in the Trilinos sparsity pattern.");
985  private:
990  std::unique_ptr<Epetra_Map> column_space_map;
991 
997  std::unique_ptr<Epetra_FECrsGraph> graph;
998 
1005  std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;
1006 
1010  };
1011 
1012 
1013 
1014  // ----------------------- inline and template functions --------------------
1015 
1016 
1017 # ifndef DOXYGEN
1018 
1019  namespace SparsityPatternIterators
1020  {
1021  inline Accessor::Accessor(const SparsityPattern *sp,
1022  const size_type row,
1023  const size_type index)
1024  : sparsity_pattern(const_cast<SparsityPattern *>(sp))
1025  , a_row(row)
1026  , a_index(index)
1027  {
1028  visit_present_row();
1029  }
1030 
1031 
1032 
1033  inline Accessor::size_type
1034  Accessor::row() const
1035  {
1036  Assert(a_row < sparsity_pattern->n_rows(),
1037  ExcBeyondEndOfSparsityPattern());
1038  return a_row;
1039  }
1040 
1041 
1042 
1043  inline Accessor::size_type
1044  Accessor::column() const
1045  {
1046  Assert(a_row < sparsity_pattern->n_rows(),
1047  ExcBeyondEndOfSparsityPattern());
1048  return (*colnum_cache)[a_index];
1049  }
1050 
1051 
1052 
1053  inline Accessor::size_type
1054  Accessor::index() const
1055  {
1056  Assert(a_row < sparsity_pattern->n_rows(),
1057  ExcBeyondEndOfSparsityPattern());
1058  return a_index;
1059  }
1060 
1061 
1062 
1063  inline Iterator::Iterator(const SparsityPattern *sp,
1064  const size_type row,
1065  const size_type index)
1066  : accessor(sp, row, index)
1067  {}
1068 
1069 
1070 
1071  inline Iterator::Iterator(const Iterator &) = default;
1072 
1073 
1074 
1075  inline Iterator &
1077  {
1078  Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1079  ExcIteratorPastEnd());
1080 
1081  ++accessor.a_index;
1082 
1083  // If at end of line: do one step, then cycle until we find a row with a
1084  // nonzero number of entries that is stored locally.
1085  if (accessor.a_index >= accessor.colnum_cache->size())
1086  {
1087  accessor.a_index = 0;
1088  ++accessor.a_row;
1089 
1090  while (accessor.a_row < accessor.sparsity_pattern->n_rows())
1091  {
1092  const auto row_length =
1093  accessor.sparsity_pattern->row_length(accessor.a_row);
1094  if (row_length == 0 ||
1095  !accessor.sparsity_pattern->row_is_stored_locally(
1096  accessor.a_row))
1097  ++accessor.a_row;
1098  else
1099  break;
1100  }
1101 
1102  accessor.visit_present_row();
1103  }
1104  return *this;
1105  }
1106 
1107 
1108 
1109  inline Iterator
1111  {
1112  const Iterator old_state = *this;
1113  ++(*this);
1114  return old_state;
1115  }
1116 
1117 
1118 
1119  inline const Accessor &
1120  Iterator::operator*() const
1121  {
1122  return accessor;
1123  }
1124 
1125 
1126 
1127  inline const Accessor *
1128  Iterator::operator->() const
1129  {
1130  return &accessor;
1131  }
1132 
1133 
1134 
1135  inline bool
1136  Iterator::operator==(const Iterator &other) const
1137  {
1138  return (accessor.a_row == other.accessor.a_row &&
1139  accessor.a_index == other.accessor.a_index);
1140  }
1141 
1142 
1143 
1144  inline bool
1145  Iterator::operator!=(const Iterator &other) const
1146  {
1147  return !(*this == other);
1148  }
1149 
1150 
1151 
1152  inline bool
1153  Iterator::operator<(const Iterator &other) const
1154  {
1155  return (accessor.row() < other.accessor.row() ||
1156  (accessor.row() == other.accessor.row() &&
1157  accessor.index() < other.accessor.index()));
1158  }
1159 
1160  } // namespace SparsityPatternIterators
1161 
1162 
1163 
1165  SparsityPattern::begin() const
1166  {
1167  const size_type first_valid_row = this->local_range().first;
1168  return const_iterator(this, first_valid_row, 0);
1169  }
1170 
1171 
1172 
1174  SparsityPattern::end() const
1175  {
1176  return const_iterator(this, n_rows(), 0);
1177  }
1178 
1179 
1180 
1182  SparsityPattern::begin(const size_type r) const
1183  {
1184  AssertIndexRange(r, n_rows());
1185  if (row_length(r) > 0)
1186  return const_iterator(this, r, 0);
1187  else
1188  return end(r);
1189  }
1190 
1191 
1192 
1194  SparsityPattern::end(const size_type r) const
1195  {
1196  AssertIndexRange(r, n_rows());
1197 
1198  // place the iterator on the first entry
1199  // past this line, or at the end of the
1200  // matrix
1201  for (size_type i = r + 1; i < n_rows(); ++i)
1202  if (row_length(i) > 0)
1203  return const_iterator(this, i, 0);
1204 
1205  // if there is no such line, then take the
1206  // end iterator of the matrix
1207  return end();
1208  }
1209 
1210 
1211 
1212  inline bool
1213  SparsityPattern::in_local_range(const size_type index) const
1214  {
1216 # ifndef DEAL_II_WITH_64BIT_INDICES
1217  begin = graph->RowMap().MinMyGID();
1218  end = graph->RowMap().MaxMyGID() + 1;
1219 # else
1220  begin = graph->RowMap().MinMyGID64();
1221  end = graph->RowMap().MaxMyGID64() + 1;
1222 # endif
1223 
1224  return ((index >= static_cast<size_type>(begin)) &&
1225  (index < static_cast<size_type>(end)));
1226  }
1227 
1228 
1229 
1230  inline bool
1232  {
1233  return graph->Filled();
1234  }
1235 
1236 
1237 
1238  inline bool
1239  SparsityPattern::empty() const
1240  {
1241  return ((n_rows() == 0) && (n_cols() == 0));
1242  }
1243 
1244 
1245 
1246  inline void
1247  SparsityPattern::add(const size_type i, const size_type j)
1248  {
1249  add_entries(i, &j, &j + 1);
1250  }
1251 
1252 
1253 
1254  template <typename ForwardIterator>
1255  inline void
1257  ForwardIterator begin,
1258  ForwardIterator end,
1259  const bool /*indices_are_sorted*/)
1260  {
1261  if (begin == end)
1262  return;
1263 
1264  // verify that the size of the data type Trilinos expects matches that the
1265  // iterator points to. we allow for some slippage between signed and
1266  // unsigned and only compare that they are both either 32 or 64 bit. to
1267  // write this test properly, not that we cannot compare the size of
1268  // '*begin' because 'begin' may be an iterator and '*begin' may be an
1269  // accessor class. consequently, we need to somehow get an actual value
1270  // from it which we can by evaluating an expression such as when
1271  // multiplying the value produced by 2
1272  Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1273  ExcNotImplemented());
1274 
1275  TrilinosWrappers::types::int_type *col_index_ptr =
1276  reinterpret_cast<TrilinosWrappers::types::int_type *>(
1277  const_cast<typename std::decay<decltype(*begin)>::type *>(&*begin));
1278  // Check at least for the first index that the conversion actually works
1279  AssertDimension(*col_index_ptr, *begin);
1280  TrilinosWrappers::types::int_type trilinos_row_index = row;
1281  const int n_cols = static_cast<int>(end - begin);
1282 
1283  int ierr;
1284  if (row_is_stored_locally(row))
1285  ierr =
1286  graph->InsertGlobalIndices(trilinos_row_index, n_cols, col_index_ptr);
1287  else if (nonlocal_graph.get() != nullptr)
1288  {
1289  // this is the case when we have explicitly set the off-processor rows
1290  // and want to create a separate matrix object for them (to retain
1291  // thread-safety)
1292  Assert(nonlocal_graph->RowMap().LID(
1293  static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1294  ExcMessage("Attempted to write into off-processor matrix row "
1295  "that has not be specified as being writable upon "
1296  "initialization"));
1297  ierr = nonlocal_graph->InsertGlobalIndices(trilinos_row_index,
1298  n_cols,
1299  col_index_ptr);
1300  }
1301  else
1302  ierr = graph->InsertGlobalIndices(1,
1303  &trilinos_row_index,
1304  n_cols,
1305  col_index_ptr);
1306 
1307  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1308  }
1309 
1310 
1311 
1312  inline const Epetra_FECrsGraph &
1314  {
1315  return *graph;
1316  }
1317 
1318 
1319 
1320  inline IndexSet
1322  {
1323  return IndexSet(graph->DomainMap());
1324  }
1325 
1326 
1327 
1328  inline IndexSet
1330  {
1331  return IndexSet(graph->RangeMap());
1332  }
1333 
1334 # endif // DOXYGEN
1335 } // namespace TrilinosWrappers
1336 
1337 
1339 
1340 
1341 # endif // DEAL_II_WITH_TRILINOS
1342 
1343 
1344 /*-------------------- trilinos_sparsity_pattern.h --------------------*/
1345 
1346 #endif
1347 /*-------------------- trilinos_sparsity_pattern.h --------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1659
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::shared_ptr< const std::vector< size_type > > colnum_cache
size_type row_length(const size_type row) const
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
IndexSet locally_owned_range_indices() const
const Epetra_FECrsGraph & trilinos_sparsity_pattern() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1724
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1577
void add(const size_type i, const size_type j)
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
SparsityPatternIterators::Iterator const_iterator
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
const_iterator begin() const
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::string compress(const std::string &input)
Definition: utilities.cc:392
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
std::unique_ptr< Epetra_Map > column_space_map
#define Assert(cond, exc)
Definition: exceptions.h:1467
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
VectorType::value_type * end(VectorType &V)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
const_iterator end() const
static ::ExceptionBase & ExcIteratorPastEnd()
unsigned int global_dof_index
Definition: types.h:76
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
static ::ExceptionBase & ExcNotImplemented()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:561
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
bool in_local_range(const size_type index) const
IndexSet locally_owned_domain_indices() const
std::pair< size_type, size_type > local_range() const
std::unique_ptr< Epetra_FECrsGraph > graph
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Accessor(const SparsityPattern *sparsity_pattern, const size_type row, const size_type index)
bool row_is_stored_locally(const size_type i) const