Reference documentation for deal.II version Git d435e92be9 2020-09-23 09:04:26 -0400
\(\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 - 2020 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 
32 # include <cmath>
33 # include <memory>
34 # include <vector>
35 # ifdef DEAL_II_WITH_MPI
36 # include <Epetra_MpiComm.h>
37 # include <mpi.h>
38 # else
39 # include <Epetra_SerialComm.h>
40 # endif
41 
42 
44 
45 // forward declarations
46 # ifndef DOXYGEN
47 class SparsityPattern;
49 
50 namespace TrilinosWrappers
51 {
52  class SparsityPattern;
53  class SparseMatrix;
54 
55  namespace SparsityPatternIterators
56  {
57  class Iterator;
58  }
59 } // namespace TrilinosWrappers
60 # endif
61 
62 namespace TrilinosWrappers
63 {
65  {
77  class Accessor
78  {
79  public:
84 
88  Accessor(const SparsityPattern *sparsity_pattern,
89  const size_type row,
90  const size_type index);
91 
95  size_type
96  row() const;
97 
101  size_type
102  index() const;
103 
107  size_type
108  column() const;
109 
113  DeclException0(ExcBeyondEndOfSparsityPattern);
114 
119  size_type,
120  size_type,
121  size_type,
122  << "You tried to access row " << arg1
123  << " of a distributed sparsity pattern, "
124  << " but only rows " << arg2 << " through " << arg3
125  << " are stored locally and can be accessed.");
126 
127  private:
132 
137 
142 
155  std::shared_ptr<const std::vector<size_type>> colnum_cache;
156 
162  void
163  visit_present_row();
164 
165  // Make enclosing class a friend.
166  friend class Iterator;
167  };
168 
174  class Iterator
175  {
176  public:
181 
186  Iterator(const SparsityPattern *sparsity_pattern,
187  const size_type row,
188  const size_type index);
189 
193  Iterator(const Iterator &i);
194 
198  Iterator &
199  operator++();
200 
204  Iterator
205  operator++(int);
206 
210  const Accessor &operator*() const;
211 
215  const Accessor *operator->() const;
216 
221  bool
222  operator==(const Iterator &) const;
223 
227  bool
228  operator!=(const Iterator &) const;
229 
235  bool
236  operator<(const Iterator &) const;
237 
241  DeclException2(ExcInvalidIndexWithinRow,
242  size_type,
243  size_type,
244  << "Attempt to access element " << arg2 << " of row "
245  << arg1 << " which doesn't have that many elements.");
246 
247  private:
252 
254  };
255 
256  } // namespace SparsityPatternIterators
257 
258 
277  {
278  public:
283 
288 
296  SparsityPattern();
297 
312  SparsityPattern(const size_type m,
313  const size_type n,
314  const size_type n_entries_per_row = 0);
315 
324  SparsityPattern(const size_type m,
325  const size_type n,
326  const std::vector<size_type> &n_entries_per_row);
327 
332  SparsityPattern(SparsityPattern &&other) noexcept;
333 
338  SparsityPattern(const SparsityPattern &input_sparsity_pattern);
339 
343  virtual ~SparsityPattern() override = default;
344 
356  void
357  reinit(const size_type m,
358  const size_type n,
359  const size_type n_entries_per_row = 0);
360 
369  void
370  reinit(const size_type m,
371  const size_type n,
372  const std::vector<size_type> &n_entries_per_row);
373 
378  void
379  copy_from(const SparsityPattern &input_sparsity_pattern);
380 
386  template <typename SparsityPatternType>
387  void
388  copy_from(const SparsityPatternType &nontrilinos_sparsity_pattern);
389 
397  operator=(const SparsityPattern &input_sparsity_pattern);
398 
406  void
407  clear();
408 
418  void
419  compress();
421 
426 
438  SparsityPattern(const IndexSet &parallel_partitioning,
439  const MPI_Comm &communicator = MPI_COMM_WORLD,
440  const size_type n_entries_per_row = 0);
441 
452  SparsityPattern(const IndexSet & parallel_partitioning,
453  const MPI_Comm & communicator,
454  const std::vector<size_type> &n_entries_per_row);
455 
470  SparsityPattern(const IndexSet &row_parallel_partitioning,
471  const IndexSet &col_parallel_partitioning,
472  const MPI_Comm &communicator = MPI_COMM_WORLD,
473  const size_type n_entries_per_row = 0);
474 
486  SparsityPattern(const IndexSet & row_parallel_partitioning,
487  const IndexSet & col_parallel_partitioning,
488  const MPI_Comm & communicator,
489  const std::vector<size_type> &n_entries_per_row);
490 
517  SparsityPattern(const IndexSet &row_parallel_partitioning,
518  const IndexSet &col_parallel_partitioning,
519  const IndexSet &writable_rows,
520  const MPI_Comm &communicator = MPI_COMM_WORLD,
521  const size_type n_entries_per_row = 0);
522 
538  void
539  reinit(const IndexSet &parallel_partitioning,
540  const MPI_Comm &communicator = MPI_COMM_WORLD,
541  const size_type n_entries_per_row = 0);
542 
553  void
554  reinit(const IndexSet & parallel_partitioning,
555  const MPI_Comm & communicator,
556  const std::vector<size_type> &n_entries_per_row);
557 
574  void
575  reinit(const IndexSet &row_parallel_partitioning,
576  const IndexSet &col_parallel_partitioning,
577  const MPI_Comm &communicator = MPI_COMM_WORLD,
578  const size_type n_entries_per_row = 0);
579 
605  void
606  reinit(const IndexSet &row_parallel_partitioning,
607  const IndexSet &col_parallel_partitioning,
608  const IndexSet &writeable_rows,
609  const MPI_Comm &communicator = MPI_COMM_WORLD,
610  const size_type n_entries_per_row = 0);
611 
616  void
617  reinit(const IndexSet & row_parallel_partitioning,
618  const IndexSet & col_parallel_partitioning,
619  const MPI_Comm & communicator,
620  const std::vector<size_type> &n_entries_per_row);
621 
631  template <typename SparsityPatternType>
632  void
633  reinit(const IndexSet & row_parallel_partitioning,
634  const IndexSet & col_parallel_partitioning,
635  const SparsityPatternType &nontrilinos_sparsity_pattern,
636  const MPI_Comm & communicator = MPI_COMM_WORLD,
637  const bool exchange_data = false);
638 
647  template <typename SparsityPatternType>
648  void
649  reinit(const IndexSet & parallel_partitioning,
650  const SparsityPatternType &nontrilinos_sparsity_pattern,
651  const MPI_Comm & communicator = MPI_COMM_WORLD,
652  const bool exchange_data = false);
654 
658 
663  bool
664  is_compressed() const;
665 
669  unsigned int
670  max_entries_per_row() const;
671 
675  size_type
676  n_rows() const;
677 
681  size_type
682  n_cols() const;
683 
693  unsigned int
694  local_size() const;
695 
704  std::pair<size_type, size_type>
705  local_range() const;
706 
711  bool
712  in_local_range(const size_type index) const;
713 
717  size_type
718  n_nonzero_elements() const;
719 
728  size_type
729  row_length(const size_type row) const;
730 
737  size_type
738  bandwidth() const;
739 
744  bool
745  empty() const;
746 
751  bool
752  exists(const size_type i, const size_type j) const;
753 
758  bool
759  row_is_stored_locally(const size_type i) const;
760 
765  std::size_t
766  memory_consumption() const;
767 
769 
776  void
777  add(const size_type i, const size_type j);
778 
779 
783  template <typename ForwardIterator>
784  void
785  add_entries(const size_type row,
786  ForwardIterator begin,
787  ForwardIterator end,
788  const bool indices_are_sorted = false);
790 
794 
799  const Epetra_FECrsGraph &
800  trilinos_sparsity_pattern() const;
801 
808  const Epetra_Map &
809  domain_partitioner() const;
810 
817  const Epetra_Map &
818  range_partitioner() const;
819 
823  MPI_Comm
824  get_mpi_communicator() const;
826 
831 
837  IndexSet
838  locally_owned_domain_indices() const;
839 
845  IndexSet
846  locally_owned_range_indices() const;
847 
849 
854 
859  begin() const;
860 
865  end() const;
866 
876  begin(const size_type r) const;
877 
887  end(const size_type r) const;
888 
890 
894 
900  void
901  write_ascii();
902 
910  void
911  print(std::ostream &out,
912  const bool write_extended_trilinos_info = false) const;
913 
928  void
929  print_gnuplot(std::ostream &out) const;
930 
932 
940  int,
941  << "An error with error number " << arg1
942  << " occurred while calling a Trilinos function");
943 
947  DeclException2(ExcInvalidIndex,
948  size_type,
949  size_type,
950  << "The entry with index <" << arg1 << ',' << arg2
951  << "> does not exist.");
952 
957  ExcSourceEqualsDestination,
958  "You are attempting an operation on two sparsity patterns that "
959  "are the same object, but the operation requires that the "
960  "two objects are in fact different.");
961 
965  DeclException4(ExcAccessToNonLocalElement,
966  size_type,
967  size_type,
968  size_type,
969  size_type,
970  << "You tried to access element (" << arg1 << "/" << arg2
971  << ")"
972  << " of a distributed matrix, but only rows " << arg3
973  << " through " << arg4
974  << " are stored locally and can be accessed.");
975 
979  DeclException2(ExcAccessToNonPresentElement,
980  size_type,
981  size_type,
982  << "You tried to access element (" << arg1 << "/" << arg2
983  << ")"
984  << " of a sparse matrix, but it appears to not"
985  << " exist in the Trilinos sparsity pattern.");
987  private:
992  std::unique_ptr<Epetra_Map> column_space_map;
993 
999  std::unique_ptr<Epetra_FECrsGraph> graph;
1000 
1007  std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;
1008 
1012  };
1013 
1014 
1015 
1016  // ----------------------- inline and template functions --------------------
1017 
1018 
1019 # ifndef DOXYGEN
1020 
1021  namespace SparsityPatternIterators
1022  {
1023  inline Accessor::Accessor(const SparsityPattern *sp,
1024  const size_type row,
1025  const size_type index)
1026  : sparsity_pattern(const_cast<SparsityPattern *>(sp))
1027  , a_row(row)
1028  , a_index(index)
1029  {
1030  visit_present_row();
1031  }
1032 
1033 
1034 
1035  inline Accessor::size_type
1036  Accessor::row() const
1037  {
1038  Assert(a_row < sparsity_pattern->n_rows(),
1039  ExcBeyondEndOfSparsityPattern());
1040  return a_row;
1041  }
1042 
1043 
1044 
1045  inline Accessor::size_type
1046  Accessor::column() const
1047  {
1048  Assert(a_row < sparsity_pattern->n_rows(),
1049  ExcBeyondEndOfSparsityPattern());
1050  return (*colnum_cache)[a_index];
1051  }
1052 
1053 
1054 
1055  inline Accessor::size_type
1056  Accessor::index() const
1057  {
1058  Assert(a_row < sparsity_pattern->n_rows(),
1059  ExcBeyondEndOfSparsityPattern());
1060  return a_index;
1061  }
1062 
1063 
1064 
1065  inline Iterator::Iterator(const SparsityPattern *sp,
1066  const size_type row,
1067  const size_type index)
1068  : accessor(sp, row, index)
1069  {}
1070 
1071 
1072 
1073  inline Iterator::Iterator(const Iterator &) = default;
1074 
1075 
1076 
1077  inline Iterator &
1079  {
1080  Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1081  ExcIteratorPastEnd());
1082 
1083  ++accessor.a_index;
1084 
1085  // If at end of line: do one step, then cycle until we find a row with a
1086  // nonzero number of entries that is stored locally.
1087  if (accessor.a_index >= accessor.colnum_cache->size())
1088  {
1089  accessor.a_index = 0;
1090  ++accessor.a_row;
1091 
1092  while (accessor.a_row < accessor.sparsity_pattern->n_rows())
1093  {
1094  const auto row_length =
1095  accessor.sparsity_pattern->row_length(accessor.a_row);
1096  if (row_length == 0 ||
1097  !accessor.sparsity_pattern->row_is_stored_locally(
1098  accessor.a_row))
1099  ++accessor.a_row;
1100  else
1101  break;
1102  }
1103 
1104  accessor.visit_present_row();
1105  }
1106  return *this;
1107  }
1108 
1109 
1110 
1111  inline Iterator
1113  {
1114  const Iterator old_state = *this;
1115  ++(*this);
1116  return old_state;
1117  }
1118 
1119 
1120 
1121  inline const Accessor &Iterator::operator*() const
1122  {
1123  return accessor;
1124  }
1125 
1126 
1127 
1128  inline const Accessor *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:1568
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:1636
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:1521
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:1411
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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:362
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