Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00:01+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 - 2023 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;
819 
832  IndexSet
834 
840  IndexSet
842 
854  begin() const;
855 
860  end() const;
861 
871  begin(const size_type r) const;
872 
882  end(const size_type r) const;
883 
895  void
896  write_ascii();
897 
905  void
906  print(std::ostream &out,
907  const bool write_extended_trilinos_info = false) const;
908 
923  void
924  print_gnuplot(std::ostream &out) const;
925 
935  int,
936  << "An error with error number " << arg1
937  << " occurred while calling a Trilinos function");
938 
943  size_type,
944  size_type,
945  << "The entry with index <" << arg1 << ',' << arg2
946  << "> does not exist.");
947 
952  size_type,
953  size_type,
954  size_type,
955  size_type,
956  << "You tried to access element (" << arg1 << '/' << arg2
957  << ')'
958  << " of a distributed matrix, but only rows in range ["
959  << arg3 << ',' << arg4
960  << "] are stored locally and can be accessed.");
961 
966  size_type,
967  size_type,
968  << "You tried to access element (" << arg1 << '/' << arg2
969  << ')' << " of a sparse matrix, but it appears to not"
970  << " exist in the Trilinos sparsity pattern.");
972  private:
977  std::unique_ptr<Epetra_Map> column_space_map;
978 
984  std::unique_ptr<Epetra_FECrsGraph> graph;
985 
992  std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;
993 
997  };
998 
999 
1000 
1001  // ----------------------- inline and template functions --------------------
1002 
1003 
1004 # ifndef DOXYGEN
1005 
1006  namespace SparsityPatternIterators
1007  {
1008  inline Accessor::Accessor(const SparsityPattern *sp,
1009  const size_type row,
1010  const size_type index)
1011  : sparsity_pattern(const_cast<SparsityPattern *>(sp))
1012  , a_row(row)
1013  , a_index(index)
1014  {
1015  visit_present_row();
1016  }
1017 
1018 
1019 
1020  inline Accessor::size_type
1021  Accessor::row() const
1022  {
1023  Assert(a_row < sparsity_pattern->n_rows(),
1024  ExcBeyondEndOfSparsityPattern());
1025  return a_row;
1026  }
1027 
1028 
1029 
1030  inline Accessor::size_type
1031  Accessor::column() const
1032  {
1033  Assert(a_row < sparsity_pattern->n_rows(),
1034  ExcBeyondEndOfSparsityPattern());
1035  return (*colnum_cache)[a_index];
1036  }
1037 
1038 
1039 
1040  inline Accessor::size_type
1041  Accessor::index() const
1042  {
1043  Assert(a_row < sparsity_pattern->n_rows(),
1044  ExcBeyondEndOfSparsityPattern());
1045  return a_index;
1046  }
1047 
1048 
1049 
1050  inline Iterator::Iterator(const SparsityPattern *sp,
1051  const size_type row,
1052  const size_type index)
1053  : accessor(sp, row, index)
1054  {}
1055 
1056 
1057 
1058  inline Iterator::Iterator(const Iterator &) = default;
1059 
1060 
1061 
1062  inline Iterator &
1064  {
1065  Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1066  ExcIteratorPastEnd());
1067 
1068  ++accessor.a_index;
1069 
1070  // If at end of line: do one step, then cycle until we find a row with a
1071  // nonzero number of entries that is stored locally.
1072  if (accessor.a_index >= accessor.colnum_cache->size())
1073  {
1074  accessor.a_index = 0;
1075  ++accessor.a_row;
1076 
1077  while (accessor.a_row < accessor.sparsity_pattern->n_rows())
1078  {
1079  const auto row_length =
1080  accessor.sparsity_pattern->row_length(accessor.a_row);
1081  if (row_length == 0 ||
1082  !accessor.sparsity_pattern->row_is_stored_locally(
1083  accessor.a_row))
1084  ++accessor.a_row;
1085  else
1086  break;
1087  }
1088 
1089  accessor.visit_present_row();
1090  }
1091  return *this;
1092  }
1093 
1094 
1095 
1096  inline Iterator
1098  {
1099  const Iterator old_state = *this;
1100  ++(*this);
1101  return old_state;
1102  }
1103 
1104 
1105 
1106  inline const Accessor &
1107  Iterator::operator*() const
1108  {
1109  return accessor;
1110  }
1111 
1112 
1113 
1114  inline const Accessor *
1115  Iterator::operator->() const
1116  {
1117  return &accessor;
1118  }
1119 
1120 
1121 
1122  inline bool
1123  Iterator::operator==(const Iterator &other) const
1124  {
1125  return (accessor.a_row == other.accessor.a_row &&
1126  accessor.a_index == other.accessor.a_index);
1127  }
1128 
1129 
1130 
1131  inline bool
1132  Iterator::operator!=(const Iterator &other) const
1133  {
1134  return !(*this == other);
1135  }
1136 
1137 
1138 
1139  inline bool
1140  Iterator::operator<(const Iterator &other) const
1141  {
1142  return (accessor.row() < other.accessor.row() ||
1143  (accessor.row() == other.accessor.row() &&
1144  accessor.index() < other.accessor.index()));
1145  }
1146 
1147  } // namespace SparsityPatternIterators
1148 
1149 
1150 
1152  SparsityPattern::begin() const
1153  {
1154  const size_type first_valid_row = this->local_range().first;
1155  return const_iterator(this, first_valid_row, 0);
1156  }
1157 
1158 
1159 
1161  SparsityPattern::end() const
1162  {
1163  return const_iterator(this, n_rows(), 0);
1164  }
1165 
1166 
1167 
1169  SparsityPattern::begin(const size_type r) const
1170  {
1171  AssertIndexRange(r, n_rows());
1172  if (row_length(r) > 0)
1173  return const_iterator(this, r, 0);
1174  else
1175  return end(r);
1176  }
1177 
1178 
1179 
1181  SparsityPattern::end(const size_type r) const
1182  {
1183  AssertIndexRange(r, n_rows());
1184 
1185  // place the iterator on the first entry
1186  // past this line, or at the end of the
1187  // matrix
1188  for (size_type i = r + 1; i < n_rows(); ++i)
1189  if (row_length(i) > 0)
1190  return const_iterator(this, i, 0);
1191 
1192  // if there is no such line, then take the
1193  // end iterator of the matrix
1194  return end();
1195  }
1196 
1197 
1198 
1199  inline bool
1200  SparsityPattern::in_local_range(const size_type index) const
1201  {
1203 # ifndef DEAL_II_WITH_64BIT_INDICES
1204  begin = graph->RowMap().MinMyGID();
1205  end = graph->RowMap().MaxMyGID() + 1;
1206 # else
1207  begin = graph->RowMap().MinMyGID64();
1208  end = graph->RowMap().MaxMyGID64() + 1;
1209 # endif
1210 
1211  return ((index >= static_cast<size_type>(begin)) &&
1212  (index < static_cast<size_type>(end)));
1213  }
1214 
1215 
1216 
1217  inline bool
1219  {
1220  return graph->Filled();
1221  }
1222 
1223 
1224 
1225  inline bool
1226  SparsityPattern::empty() const
1227  {
1228  return ((n_rows() == 0) && (n_cols() == 0));
1229  }
1230 
1231 
1232 
1233  inline void
1234  SparsityPattern::add(const size_type i, const size_type j)
1235  {
1236  add_entries(i, &j, &j + 1);
1237  }
1238 
1239 
1240 
1241  template <typename ForwardIterator>
1242  inline void
1244  ForwardIterator begin,
1245  ForwardIterator end,
1246  const bool /*indices_are_sorted*/)
1247  {
1248  if (begin == end)
1249  return;
1250 
1251  // verify that the size of the data type Trilinos expects matches that the
1252  // iterator points to. we allow for some slippage between signed and
1253  // unsigned and only compare that they are both either 32 or 64 bit. to
1254  // write this test properly, not that we cannot compare the size of
1255  // '*begin' because 'begin' may be an iterator and '*begin' may be an
1256  // accessor class. consequently, we need to somehow get an actual value
1257  // from it which we can by evaluating an expression such as when
1258  // multiplying the value produced by 2
1259  Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1260  ExcNotImplemented());
1261 
1262  TrilinosWrappers::types::int_type *col_index_ptr =
1263  reinterpret_cast<TrilinosWrappers::types::int_type *>(
1264  const_cast<std::decay_t<decltype(*begin)> *>(&*begin));
1265  // Check at least for the first index that the conversion actually works
1266  AssertDimension(*col_index_ptr, *begin);
1267  TrilinosWrappers::types::int_type trilinos_row_index = row;
1268  const int n_cols = static_cast<int>(end - begin);
1269 
1270  int ierr;
1271  if (row_is_stored_locally(row))
1272  ierr =
1273  graph->InsertGlobalIndices(trilinos_row_index, n_cols, col_index_ptr);
1274  else if (nonlocal_graph.get() != nullptr)
1275  {
1276  // this is the case when we have explicitly set the off-processor rows
1277  // and want to create a separate matrix object for them (to retain
1278  // thread-safety)
1279  Assert(nonlocal_graph->RowMap().LID(
1280  static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1281  ExcMessage("Attempted to write into off-processor matrix row "
1282  "that has not be specified as being writable upon "
1283  "initialization"));
1284  ierr = nonlocal_graph->InsertGlobalIndices(trilinos_row_index,
1285  n_cols,
1286  col_index_ptr);
1287  }
1288  else
1289  ierr = graph->InsertGlobalIndices(1,
1290  &trilinos_row_index,
1291  n_cols,
1292  col_index_ptr);
1293 
1294  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1295  }
1296 
1297 
1298 
1299  inline const Epetra_FECrsGraph &
1300  SparsityPattern::trilinos_sparsity_pattern() const
1301  {
1302  return *graph;
1303  }
1304 
1305 
1306 
1307  inline IndexSet
1308  SparsityPattern::locally_owned_domain_indices() const
1309  {
1310  return IndexSet(graph->DomainMap());
1311  }
1312 
1313 
1314 
1315  inline IndexSet
1316  SparsityPattern::locally_owned_range_indices() const
1317  {
1318  return IndexSet(graph->RangeMap());
1319  }
1320 
1321 # endif // DOXYGEN
1322 } // namespace TrilinosWrappers
1323 
1324 
1326 
1327 
1328 #endif // DEAL_II_WITH_TRILINOS
1329 
1330 #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_v< T > &&std::is_floating_point_v< U >, 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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
#define DeclException0(Exception0)
Definition: exceptions.h:472
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:586
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:540
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:563
static ::ExceptionBase & ExcAccessToNonlocalRow(size_type arg1, size_type arg2, size_type arg3)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:517
static ::ExceptionBase & ExcBeyondEndOfSparsityPattern()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
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:82
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)