Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18:55: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 
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 
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();
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  std::uint64_t
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
784  ForwardIterator begin,
785  ForwardIterator end,
786  const bool indices_are_sorted = false);
788 
792 
797  const Epetra_FECrsGraph &
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
837 
843  IndexSet
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 
946  size_type,
947  size_type,
948  << "The entry with index <" << arg1 << ',' << arg2
949  << "> does not exist.");
950 
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 
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 
978  size_type,
979  size_type,
980  << "You tried to access element (" << arg1 << '/' << arg2
981  << ')' << " of a sparse matrix, but it appears to not"
982  << " exist in the Trilinos sparsity pattern.");
984  private:
989  std::unique_ptr<Epetra_Map> column_space_map;
990 
996  std::unique_ptr<Epetra_FECrsGraph> graph;
997 
1004  std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;
1005 
1009  };
1010 
1011 
1012 
1013  // ----------------------- inline and template functions --------------------
1014 
1015 
1016 # ifndef DOXYGEN
1017 
1018  namespace SparsityPatternIterators
1019  {
1020  inline Accessor::Accessor(const SparsityPattern *sp,
1021  const size_type row,
1022  const size_type index)
1023  : sparsity_pattern(const_cast<SparsityPattern *>(sp))
1024  , a_row(row)
1025  , a_index(index)
1026  {
1027  visit_present_row();
1028  }
1029 
1030 
1031 
1032  inline Accessor::size_type
1033  Accessor::row() const
1034  {
1035  Assert(a_row < sparsity_pattern->n_rows(),
1036  ExcBeyondEndOfSparsityPattern());
1037  return a_row;
1038  }
1039 
1040 
1041 
1042  inline Accessor::size_type
1043  Accessor::column() const
1044  {
1045  Assert(a_row < sparsity_pattern->n_rows(),
1046  ExcBeyondEndOfSparsityPattern());
1047  return (*colnum_cache)[a_index];
1048  }
1049 
1050 
1051 
1052  inline Accessor::size_type
1053  Accessor::index() const
1054  {
1055  Assert(a_row < sparsity_pattern->n_rows(),
1056  ExcBeyondEndOfSparsityPattern());
1057  return a_index;
1058  }
1059 
1060 
1061 
1062  inline Iterator::Iterator(const SparsityPattern *sp,
1063  const size_type row,
1064  const size_type index)
1065  : accessor(sp, row, index)
1066  {}
1067 
1068 
1069 
1070  inline Iterator::Iterator(const Iterator &) = default;
1071 
1072 
1073 
1074  inline Iterator &
1076  {
1077  Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1078  ExcIteratorPastEnd());
1079 
1080  ++accessor.a_index;
1081 
1082  // If at end of line: do one step, then cycle until we find a row with a
1083  // nonzero number of entries that is stored locally.
1084  if (accessor.a_index >= accessor.colnum_cache->size())
1085  {
1086  accessor.a_index = 0;
1087  ++accessor.a_row;
1088 
1089  while (accessor.a_row < accessor.sparsity_pattern->n_rows())
1090  {
1091  const auto row_length =
1092  accessor.sparsity_pattern->row_length(accessor.a_row);
1093  if (row_length == 0 ||
1094  !accessor.sparsity_pattern->row_is_stored_locally(
1095  accessor.a_row))
1096  ++accessor.a_row;
1097  else
1098  break;
1099  }
1100 
1101  accessor.visit_present_row();
1102  }
1103  return *this;
1104  }
1105 
1106 
1107 
1108  inline Iterator
1110  {
1111  const Iterator old_state = *this;
1112  ++(*this);
1113  return old_state;
1114  }
1115 
1116 
1117 
1118  inline const Accessor &
1119  Iterator::operator*() const
1120  {
1121  return accessor;
1122  }
1123 
1124 
1125 
1126  inline const Accessor *
1127  Iterator::operator->() const
1128  {
1129  return &accessor;
1130  }
1131 
1132 
1133 
1134  inline bool
1135  Iterator::operator==(const Iterator &other) const
1136  {
1137  return (accessor.a_row == other.accessor.a_row &&
1138  accessor.a_index == other.accessor.a_index);
1139  }
1140 
1141 
1142 
1143  inline bool
1144  Iterator::operator!=(const Iterator &other) const
1145  {
1146  return !(*this == other);
1147  }
1148 
1149 
1150 
1151  inline bool
1152  Iterator::operator<(const Iterator &other) const
1153  {
1154  return (accessor.row() < other.accessor.row() ||
1155  (accessor.row() == other.accessor.row() &&
1156  accessor.index() < other.accessor.index()));
1157  }
1158 
1159  } // namespace SparsityPatternIterators
1160 
1161 
1162 
1164  SparsityPattern::begin() const
1165  {
1166  const size_type first_valid_row = this->local_range().first;
1167  return const_iterator(this, first_valid_row, 0);
1168  }
1169 
1170 
1171 
1173  SparsityPattern::end() const
1174  {
1175  return const_iterator(this, n_rows(), 0);
1176  }
1177 
1178 
1179 
1181  SparsityPattern::begin(const size_type r) const
1182  {
1183  AssertIndexRange(r, n_rows());
1184  if (row_length(r) > 0)
1185  return const_iterator(this, r, 0);
1186  else
1187  return end(r);
1188  }
1189 
1190 
1191 
1193  SparsityPattern::end(const size_type r) const
1194  {
1195  AssertIndexRange(r, n_rows());
1196 
1197  // place the iterator on the first entry
1198  // past this line, or at the end of the
1199  // matrix
1200  for (size_type i = r + 1; i < n_rows(); ++i)
1201  if (row_length(i) > 0)
1202  return const_iterator(this, i, 0);
1203 
1204  // if there is no such line, then take the
1205  // end iterator of the matrix
1206  return end();
1207  }
1208 
1209 
1210 
1211  inline bool
1212  SparsityPattern::in_local_range(const size_type index) const
1213  {
1215 # ifndef DEAL_II_WITH_64BIT_INDICES
1216  begin = graph->RowMap().MinMyGID();
1217  end = graph->RowMap().MaxMyGID() + 1;
1218 # else
1219  begin = graph->RowMap().MinMyGID64();
1220  end = graph->RowMap().MaxMyGID64() + 1;
1221 # endif
1222 
1223  return ((index >= static_cast<size_type>(begin)) &&
1224  (index < static_cast<size_type>(end)));
1225  }
1226 
1227 
1228 
1229  inline bool
1231  {
1232  return graph->Filled();
1233  }
1234 
1235 
1236 
1237  inline bool
1238  SparsityPattern::empty() const
1239  {
1240  return ((n_rows() == 0) && (n_cols() == 0));
1241  }
1242 
1243 
1244 
1245  inline void
1246  SparsityPattern::add(const size_type i, const size_type j)
1247  {
1248  add_entries(i, &j, &j + 1);
1249  }
1250 
1251 
1252 
1253  template <typename ForwardIterator>
1254  inline void
1256  ForwardIterator begin,
1257  ForwardIterator end,
1258  const bool /*indices_are_sorted*/)
1259  {
1260  if (begin == end)
1261  return;
1262 
1263  // verify that the size of the data type Trilinos expects matches that the
1264  // iterator points to. we allow for some slippage between signed and
1265  // unsigned and only compare that they are both either 32 or 64 bit. to
1266  // write this test properly, not that we cannot compare the size of
1267  // '*begin' because 'begin' may be an iterator and '*begin' may be an
1268  // accessor class. consequently, we need to somehow get an actual value
1269  // from it which we can by evaluating an expression such as when
1270  // multiplying the value produced by 2
1271  Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1272  ExcNotImplemented());
1273 
1274  TrilinosWrappers::types::int_type *col_index_ptr =
1275  reinterpret_cast<TrilinosWrappers::types::int_type *>(
1276  const_cast<typename std::decay<decltype(*begin)>::type *>(&*begin));
1277  // Check at least for the first index that the conversion actually works
1278  AssertDimension(*col_index_ptr, *begin);
1279  TrilinosWrappers::types::int_type trilinos_row_index = row;
1280  const int n_cols = static_cast<int>(end - begin);
1281 
1282  int ierr;
1283  if (row_is_stored_locally(row))
1284  ierr =
1285  graph->InsertGlobalIndices(trilinos_row_index, n_cols, col_index_ptr);
1286  else if (nonlocal_graph.get() != nullptr)
1287  {
1288  // this is the case when we have explicitly set the off-processor rows
1289  // and want to create a separate matrix object for them (to retain
1290  // thread-safety)
1291  Assert(nonlocal_graph->RowMap().LID(
1292  static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1293  ExcMessage("Attempted to write into off-processor matrix row "
1294  "that has not be specified as being writable upon "
1295  "initialization"));
1296  ierr = nonlocal_graph->InsertGlobalIndices(trilinos_row_index,
1297  n_cols,
1298  col_index_ptr);
1299  }
1300  else
1301  ierr = graph->InsertGlobalIndices(1,
1302  &trilinos_row_index,
1303  n_cols,
1304  col_index_ptr);
1305 
1306  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1307  }
1308 
1309 
1310 
1311  inline const Epetra_FECrsGraph &
1312  SparsityPattern::trilinos_sparsity_pattern() const
1313  {
1314  return *graph;
1315  }
1316 
1317 
1318 
1319  inline IndexSet
1320  SparsityPattern::locally_owned_domain_indices() const
1321  {
1322  return IndexSet(graph->DomainMap());
1323  }
1324 
1325 
1326 
1327  inline IndexSet
1328  SparsityPattern::locally_owned_range_indices() const
1329  {
1330  return IndexSet(graph->RangeMap());
1331  }
1332 
1333 # endif // DOXYGEN
1334 } // namespace TrilinosWrappers
1335 
1336 
1338 
1339 
1340 # endif // DEAL_II_WITH_TRILINOS
1341 
1342 
1343 /*-------------------- trilinos_sparsity_pattern.h --------------------*/
1344 
1345 #endif
1346 /*-------------------- trilinos_sparsity_pattern.h --------------------*/
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void add(const size_type i, const size_type j)
unsigned int row_length(const size_type row) const
size_type n_rows() const
iterator end() const
bool is_compressed() const
iterator begin() const
size_type n_cols() const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
SparsityPatternBase::const_iterator const_iterator
SparsityPatternBase::size_type size_type
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
void add(const size_type i, const size_type j)
const_iterator end() const
const_iterator begin(const size_type r) 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
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
const Epetra_Map & range_partitioner() const
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< 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)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
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)
std::unique_ptr< Epetra_FECrsGraph > graph
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void print_gnuplot(std::ostream &out) const
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
std::unique_ptr< Epetra_Map > column_space_map
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
#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:1583
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:76
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)