Reference documentation for deal.II version 9.4.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 - 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
26
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
43class SparsityPattern;
45
46namespace TrilinosWrappers
47{
48 class SparsityPattern;
49 class SparseMatrix;
50
52 {
53 class Iterator;
54 }
55} // namespace TrilinosWrappers
56# endif
57
58namespace TrilinosWrappers
59{
61 {
74 {
75 public:
80
85 const size_type row,
86 const size_type index);
87
92 row() const;
93
98 index() const;
99
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
171 {
172 public:
177
182 Iterator(const SparsityPattern *sparsity_pattern,
183 const size_type row,
184 const size_type index);
185
190
194 Iterator &
196
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
295
311 const size_type n,
312 const size_type n_entries_per_row = 0);
313
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
663
667 unsigned int
668 max_entries_per_row() const;
669
674 n_rows() const;
675
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
727 row_length(const size_type row) const;
728
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
822 get_mpi_communicator() const;
824
829
837
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(),
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
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 {
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 {
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
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),
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
const Epetra_FECrsGraph & trilinos_sparsity_pattern() 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
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, typenameProductType< 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 & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcBeyondEndOfSparsityPattern()
std::unique_ptr< Epetra_FECrsGraph > graph
#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
static ::ExceptionBase & ExcNotImplemented()
void print_gnuplot(std::ostream &out) const
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcIteratorPastEnd()
std::unique_ptr< Epetra_CrsGraph > nonlocal_graph
static ::ExceptionBase & ExcAccessToNonlocalRow(size_type arg1, size_type arg2, size_type arg3)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcSourceEqualsDestination()
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
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
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 * end(VectorType &V)
VectorType::value_type * begin(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)