Reference documentation for deal.II version 9.2.0
\(\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  {
79  class Accessor
80  {
81  public:
86 
90  Accessor(const SparsityPattern *sparsity_pattern,
91  const size_type row,
92  const size_type index);
93 
97  size_type
98  row() const;
99 
103  size_type
104  index() const;
105 
109  size_type
110  column() const;
111 
115  DeclException0(ExcBeyondEndOfSparsityPattern);
116 
121  size_type,
122  size_type,
123  size_type,
124  << "You tried to access row " << arg1
125  << " of a distributed sparsity pattern, "
126  << " but only rows " << arg2 << " through " << arg3
127  << " are stored locally and can be accessed.");
128 
129  private:
134 
139 
144 
157  std::shared_ptr<const std::vector<size_type>> colnum_cache;
158 
164  void
165  visit_present_row();
166 
167  // Make enclosing class a friend.
168  friend class Iterator;
169  };
170 
176  class Iterator
177  {
178  public:
183 
188  Iterator(const SparsityPattern *sparsity_pattern,
189  const size_type row,
190  const size_type index);
191 
195  Iterator(const Iterator &i);
196 
200  Iterator &
201  operator++();
202 
206  Iterator
207  operator++(int);
208 
212  const Accessor &operator*() const;
213 
217  const Accessor *operator->() const;
218 
223  bool
224  operator==(const Iterator &) const;
225 
229  bool
230  operator!=(const Iterator &) const;
231 
237  bool
238  operator<(const Iterator &) const;
239 
243  DeclException2(ExcInvalidIndexWithinRow,
244  size_type,
245  size_type,
246  << "Attempt to access element " << arg2 << " of row "
247  << arg1 << " which doesn't have that many elements.");
248 
249  private:
254 
256  };
257 
258  } // namespace SparsityPatternIterators
259 
260 
280  {
281  public:
286 
291 
299  SparsityPattern();
300 
315  SparsityPattern(const size_type m,
316  const size_type n,
317  const size_type n_entries_per_row = 0);
318 
327  SparsityPattern(const size_type m,
328  const size_type n,
329  const std::vector<size_type> &n_entries_per_row);
330 
335  SparsityPattern(SparsityPattern &&other) noexcept;
336 
341  SparsityPattern(const SparsityPattern &input_sparsity_pattern);
342 
346  virtual ~SparsityPattern() override = default;
347 
359  void
360  reinit(const size_type m,
361  const size_type n,
362  const size_type n_entries_per_row = 0);
363 
372  void
373  reinit(const size_type m,
374  const size_type n,
375  const std::vector<size_type> &n_entries_per_row);
376 
381  void
382  copy_from(const SparsityPattern &input_sparsity_pattern);
383 
389  template <typename SparsityPatternType>
390  void
391  copy_from(const SparsityPatternType &nontrilinos_sparsity_pattern);
392 
400  operator=(const SparsityPattern &input_sparsity_pattern);
401 
409  void
410  clear();
411 
421  void
422  compress();
424 
429 
441  SparsityPattern(const IndexSet &parallel_partitioning,
442  const MPI_Comm &communicator = MPI_COMM_WORLD,
443  const size_type n_entries_per_row = 0);
444 
455  SparsityPattern(const IndexSet & parallel_partitioning,
456  const MPI_Comm & communicator,
457  const std::vector<size_type> &n_entries_per_row);
458 
473  SparsityPattern(const IndexSet &row_parallel_partitioning,
474  const IndexSet &col_parallel_partitioning,
475  const MPI_Comm &communicator = MPI_COMM_WORLD,
476  const size_type n_entries_per_row = 0);
477 
489  SparsityPattern(const IndexSet & row_parallel_partitioning,
490  const IndexSet & col_parallel_partitioning,
491  const MPI_Comm & communicator,
492  const std::vector<size_type> &n_entries_per_row);
493 
520  SparsityPattern(const IndexSet &row_parallel_partitioning,
521  const IndexSet &col_parallel_partitioning,
522  const IndexSet &writable_rows,
523  const MPI_Comm &communicator = MPI_COMM_WORLD,
524  const size_type n_entries_per_row = 0);
525 
541  void
542  reinit(const IndexSet &parallel_partitioning,
543  const MPI_Comm &communicator = MPI_COMM_WORLD,
544  const size_type n_entries_per_row = 0);
545 
556  void
557  reinit(const IndexSet & parallel_partitioning,
558  const MPI_Comm & communicator,
559  const std::vector<size_type> &n_entries_per_row);
560 
577  void
578  reinit(const IndexSet &row_parallel_partitioning,
579  const IndexSet &col_parallel_partitioning,
580  const MPI_Comm &communicator = MPI_COMM_WORLD,
581  const size_type n_entries_per_row = 0);
582 
608  void
609  reinit(const IndexSet &row_parallel_partitioning,
610  const IndexSet &col_parallel_partitioning,
611  const IndexSet &writeable_rows,
612  const MPI_Comm &communicator = MPI_COMM_WORLD,
613  const size_type n_entries_per_row = 0);
614 
619  void
620  reinit(const IndexSet & row_parallel_partitioning,
621  const IndexSet & col_parallel_partitioning,
622  const MPI_Comm & communicator,
623  const std::vector<size_type> &n_entries_per_row);
624 
634  template <typename SparsityPatternType>
635  void
636  reinit(const IndexSet & row_parallel_partitioning,
637  const IndexSet & col_parallel_partitioning,
638  const SparsityPatternType &nontrilinos_sparsity_pattern,
639  const MPI_Comm & communicator = MPI_COMM_WORLD,
640  const bool exchange_data = false);
641 
650  template <typename SparsityPatternType>
651  void
652  reinit(const IndexSet & parallel_partitioning,
653  const SparsityPatternType &nontrilinos_sparsity_pattern,
654  const MPI_Comm & communicator = MPI_COMM_WORLD,
655  const bool exchange_data = false);
657 
661 
666  bool
667  is_compressed() const;
668 
672  unsigned int
673  max_entries_per_row() const;
674 
678  size_type
679  n_rows() const;
680 
684  size_type
685  n_cols() const;
686 
696  unsigned int
697  local_size() const;
698 
707  std::pair<size_type, size_type>
708  local_range() const;
709 
714  bool
715  in_local_range(const size_type index) const;
716 
720  size_type
721  n_nonzero_elements() const;
722 
731  size_type
732  row_length(const size_type row) const;
733 
740  size_type
741  bandwidth() const;
742 
747  bool
748  empty() const;
749 
754  bool
755  exists(const size_type i, const size_type j) const;
756 
761  bool
762  row_is_stored_locally(const size_type i) const;
763 
768  std::size_t
769  memory_consumption() const;
770 
772 
779  void
780  add(const size_type i, const size_type j);
781 
782 
786  template <typename ForwardIterator>
787  void
788  add_entries(const size_type row,
789  ForwardIterator begin,
790  ForwardIterator end,
791  const bool indices_are_sorted = false);
793 
797 
802  const Epetra_FECrsGraph &
803  trilinos_sparsity_pattern() const;
804 
811  const Epetra_Map &
812  domain_partitioner() const;
813 
820  const Epetra_Map &
821  range_partitioner() const;
822 
826  MPI_Comm
827  get_mpi_communicator() const;
829 
834 
840  IndexSet
841  locally_owned_domain_indices() const;
842 
848  IndexSet
849  locally_owned_range_indices() const;
850 
852 
857 
862  begin() const;
863 
868  end() const;
869 
879  begin(const size_type r) const;
880 
890  end(const size_type r) const;
891 
893 
897 
903  void
904  write_ascii();
905 
913  void
914  print(std::ostream &out,
915  const bool write_extended_trilinos_info = false) const;
916 
931  void
932  print_gnuplot(std::ostream &out) const;
933 
935 
943  int,
944  << "An error with error number " << arg1
945  << " occurred while calling a Trilinos function");
946 
950  DeclException2(ExcInvalidIndex,
951  size_type,
952  size_type,
953  << "The entry with index <" << arg1 << ',' << arg2
954  << "> does not exist.");
955 
960  ExcSourceEqualsDestination,
961  "You are attempting an operation on two sparsity patterns that "
962  "are the same object, but the operation requires that the "
963  "two objects are in fact different.");
964 
968  DeclException4(ExcAccessToNonLocalElement,
969  size_type,
970  size_type,
971  size_type,
972  size_type,
973  << "You tried to access element (" << arg1 << "/" << arg2
974  << ")"
975  << " of a distributed matrix, but only rows " << arg3
976  << " through " << arg4
977  << " are stored locally and can be accessed.");
978 
982  DeclException2(ExcAccessToNonPresentElement,
983  size_type,
984  size_type,
985  << "You tried to access element (" << arg1 << "/" << arg2
986  << ")"
987  << " of a sparse matrix, but it appears to not"
988  << " exist in the Trilinos sparsity pattern.");
990  private:
995  std::unique_ptr<Epetra_Map> column_space_map;
996 
1002  std::unique_ptr<Epetra_FECrsGraph> graph;
1003 
1010  std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;
1011 
1015  };
1016 
1017 
1018 
1019  // ----------------------- inline and template functions --------------------
1020 
1021 
1022 # ifndef DOXYGEN
1023 
1024  namespace SparsityPatternIterators
1025  {
1026  inline Accessor::Accessor(const SparsityPattern *sp,
1027  const size_type row,
1028  const size_type index)
1029  : sparsity_pattern(const_cast<SparsityPattern *>(sp))
1030  , a_row(row)
1031  , a_index(index)
1032  {
1033  visit_present_row();
1034  }
1035 
1036 
1037 
1038  inline Accessor::size_type
1039  Accessor::row() const
1040  {
1041  Assert(a_row < sparsity_pattern->n_rows(),
1042  ExcBeyondEndOfSparsityPattern());
1043  return a_row;
1044  }
1045 
1046 
1047 
1048  inline Accessor::size_type
1049  Accessor::column() const
1050  {
1051  Assert(a_row < sparsity_pattern->n_rows(),
1052  ExcBeyondEndOfSparsityPattern());
1053  return (*colnum_cache)[a_index];
1054  }
1055 
1056 
1057 
1058  inline Accessor::size_type
1059  Accessor::index() const
1060  {
1061  Assert(a_row < sparsity_pattern->n_rows(),
1062  ExcBeyondEndOfSparsityPattern());
1063  return a_index;
1064  }
1065 
1066 
1067 
1068  inline Iterator::Iterator(const SparsityPattern *sp,
1069  const size_type row,
1070  const size_type index)
1071  : accessor(sp, row, index)
1072  {}
1073 
1074 
1075 
1076  inline Iterator::Iterator(const Iterator &) = default;
1077 
1078 
1079 
1080  inline Iterator &
1082  {
1083  Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1084  ExcIteratorPastEnd());
1085 
1086  ++accessor.a_index;
1087 
1088  // If at end of line: do one step, then cycle until we find a row with a
1089  // nonzero number of entries that is stored locally.
1090  if (accessor.a_index >= accessor.colnum_cache->size())
1091  {
1092  accessor.a_index = 0;
1093  ++accessor.a_row;
1094 
1095  while (accessor.a_row < accessor.sparsity_pattern->n_rows())
1096  {
1097  const auto row_length =
1098  accessor.sparsity_pattern->row_length(accessor.a_row);
1099  if (row_length == 0 ||
1100  !accessor.sparsity_pattern->row_is_stored_locally(
1101  accessor.a_row))
1102  ++accessor.a_row;
1103  else
1104  break;
1105  }
1106 
1107  accessor.visit_present_row();
1108  }
1109  return *this;
1110  }
1111 
1112 
1113 
1114  inline Iterator
1116  {
1117  const Iterator old_state = *this;
1118  ++(*this);
1119  return old_state;
1120  }
1121 
1122 
1123 
1124  inline const Accessor &Iterator::operator*() const
1125  {
1126  return accessor;
1127  }
1128 
1129 
1130 
1131  inline const Accessor *Iterator::operator->() const
1132  {
1133  return &accessor;
1134  }
1135 
1136 
1137 
1138  inline bool
1139  Iterator::operator==(const Iterator &other) const
1140  {
1141  return (accessor.a_row == other.accessor.a_row &&
1142  accessor.a_index == other.accessor.a_index);
1143  }
1144 
1145 
1146 
1147  inline bool
1148  Iterator::operator!=(const Iterator &other) const
1149  {
1150  return !(*this == other);
1151  }
1152 
1153 
1154 
1155  inline bool
1156  Iterator::operator<(const Iterator &other) const
1157  {
1158  return (accessor.row() < other.accessor.row() ||
1159  (accessor.row() == other.accessor.row() &&
1160  accessor.index() < other.accessor.index()));
1161  }
1162 
1163  } // namespace SparsityPatternIterators
1164 
1165 
1166 
1168  SparsityPattern::begin() const
1169  {
1170  const size_type first_valid_row = this->local_range().first;
1171  return const_iterator(this, first_valid_row, 0);
1172  }
1173 
1174 
1175 
1177  SparsityPattern::end() const
1178  {
1179  return const_iterator(this, n_rows(), 0);
1180  }
1181 
1182 
1183 
1185  SparsityPattern::begin(const size_type r) const
1186  {
1187  AssertIndexRange(r, n_rows());
1188  if (row_length(r) > 0)
1189  return const_iterator(this, r, 0);
1190  else
1191  return end(r);
1192  }
1193 
1194 
1195 
1197  SparsityPattern::end(const size_type r) const
1198  {
1199  AssertIndexRange(r, n_rows());
1200 
1201  // place the iterator on the first entry
1202  // past this line, or at the end of the
1203  // matrix
1204  for (size_type i = r + 1; i < n_rows(); ++i)
1205  if (row_length(i) > 0)
1206  return const_iterator(this, i, 0);
1207 
1208  // if there is no such line, then take the
1209  // end iterator of the matrix
1210  return end();
1211  }
1212 
1213 
1214 
1215  inline bool
1216  SparsityPattern::in_local_range(const size_type index) const
1217  {
1219 # ifndef DEAL_II_WITH_64BIT_INDICES
1220  begin = graph->RowMap().MinMyGID();
1221  end = graph->RowMap().MaxMyGID() + 1;
1222 # else
1223  begin = graph->RowMap().MinMyGID64();
1224  end = graph->RowMap().MaxMyGID64() + 1;
1225 # endif
1226 
1227  return ((index >= static_cast<size_type>(begin)) &&
1228  (index < static_cast<size_type>(end)));
1229  }
1230 
1231 
1232 
1233  inline bool
1235  {
1236  return graph->Filled();
1237  }
1238 
1239 
1240 
1241  inline bool
1242  SparsityPattern::empty() const
1243  {
1244  return ((n_rows() == 0) && (n_cols() == 0));
1245  }
1246 
1247 
1248 
1249  inline void
1250  SparsityPattern::add(const size_type i, const size_type j)
1251  {
1252  add_entries(i, &j, &j + 1);
1253  }
1254 
1255 
1256 
1257  template <typename ForwardIterator>
1258  inline void
1260  ForwardIterator begin,
1261  ForwardIterator end,
1262  const bool /*indices_are_sorted*/)
1263  {
1264  if (begin == end)
1265  return;
1266 
1267  // verify that the size of the data type Trilinos expects matches that the
1268  // iterator points to. we allow for some slippage between signed and
1269  // unsigned and only compare that they are both either 32 or 64 bit. to
1270  // write this test properly, not that we cannot compare the size of
1271  // '*begin' because 'begin' may be an iterator and '*begin' may be an
1272  // accessor class. consequently, we need to somehow get an actual value
1273  // from it which we can by evaluating an expression such as when
1274  // multiplying the value produced by 2
1275  Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1276  ExcNotImplemented());
1277 
1278  TrilinosWrappers::types::int_type *col_index_ptr =
1279  reinterpret_cast<TrilinosWrappers::types::int_type *>(
1280  const_cast<typename std::decay<decltype(*begin)>::type *>(&*begin));
1281  // Check at least for the first index that the conversion actually works
1282  AssertDimension(*col_index_ptr, *begin);
1283  TrilinosWrappers::types::int_type trilinos_row_index = row;
1284  const int n_cols = static_cast<int>(end - begin);
1285 
1286  int ierr;
1287  if (row_is_stored_locally(row))
1288  ierr =
1289  graph->InsertGlobalIndices(trilinos_row_index, n_cols, col_index_ptr);
1290  else if (nonlocal_graph.get() != nullptr)
1291  {
1292  // this is the case when we have explicitly set the off-processor rows
1293  // and want to create a separate matrix object for them (to retain
1294  // thread-safety)
1295  Assert(nonlocal_graph->RowMap().LID(
1296  static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1297  ExcMessage("Attempted to write into off-processor matrix row "
1298  "that has not be specified as being writable upon "
1299  "initialization"));
1300  ierr = nonlocal_graph->InsertGlobalIndices(trilinos_row_index,
1301  n_cols,
1302  col_index_ptr);
1303  }
1304  else
1305  ierr = graph->InsertGlobalIndices(1,
1306  &trilinos_row_index,
1307  n_cols,
1308  col_index_ptr);
1309 
1310  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1311  }
1312 
1313 
1314 
1315  inline const Epetra_FECrsGraph &
1317  {
1318  return *graph;
1319  }
1320 
1321 
1322 
1323  inline IndexSet
1325  {
1326  return IndexSet(graph->DomainMap());
1327  }
1328 
1329 
1330 
1331  inline IndexSet
1333  {
1334  return IndexSet(graph->RangeMap());
1335  }
1336 
1337 # endif // DOXYGEN
1338 } // namespace TrilinosWrappers
1339 
1340 
1342 
1343 
1344 # endif // DEAL_II_WITH_TRILINOS
1345 
1346 
1347 /*-------------------- trilinos_sparsity_pattern.h --------------------*/
1348 
1349 #endif
1350 /*-------------------- trilinos_sparsity_pattern.h --------------------*/
static ::ExceptionBase & ExcTrilinosError(int arg1)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
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:1649
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:1531
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:518
std::unique_ptr< Epetra_Map > column_space_map
#define Assert(cond, exc)
Definition: exceptions.h:1419
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
#define DeclException0(Exception0)
Definition: exceptions.h:473
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType::value_type * end(VectorType &V)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
const_iterator end() const
MatrixTableIterators::Accessor< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Accessor
Definition: table.h:1907
static ::ExceptionBase & ExcIteratorPastEnd()
unsigned int global_dof_index
Definition: types.h:76
MatrixTableIterators::Iterator< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Iterator
Definition: table.h:1913
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
static ::ExceptionBase & ExcNotImplemented()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:564
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