Reference documentation for deal.II version GIT relicensing-756-gc37f912904 2024-05-25 08:20: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\}}\)
Loading...
Searching...
No Matches
trilinos_sparsity_pattern.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_trilinos_sparsity_pattern_h
16#define dealii_trilinos_sparsity_pattern_h
17
18#include <deal.II/base/config.h>
19
20#ifdef DEAL_II_WITH_TRILINOS
21
25
28
29# include <Epetra_FECrsGraph.h>
30# include <Epetra_Map.h>
31# include <Epetra_MpiComm.h>
32
33# include <cmath>
34# include <memory>
35# include <vector>
36
37
39
40// forward declarations
41# ifndef DOXYGEN
42class SparsityPattern;
44
45namespace TrilinosWrappers
46{
47 class SparsityPattern;
48 class SparseMatrix;
49
51 {
52 class Iterator;
53 }
54} // namespace TrilinosWrappers
55# endif
56
57namespace TrilinosWrappers
58{
60 {
73 {
74 public:
79
84 const size_type row,
85 const size_type index);
86
91 row() const;
92
97 index() const;
98
103 column() const;
104
109
114 size_type,
115 size_type,
116 size_type,
117 << "You tried to access row " << arg1
118 << " of a distributed sparsity pattern, "
119 << " but only rows " << arg2 << " through " << arg3
120 << " are stored locally and can be accessed.");
121
122 private:
127
132
137
150 std::shared_ptr<const std::vector<size_type>> colnum_cache;
151
157 void
159
160 // Make enclosing class a friend.
161 friend class Iterator;
162 };
163
170 {
171 public:
176
181 Iterator(const SparsityPattern *sparsity_pattern,
182 const size_type row,
183 const size_type index);
184
189
193 Iterator &
195
201
205 const Accessor &
206 operator*() const;
207
211 const Accessor *
212 operator->() const;
213
218 bool
219 operator==(const Iterator &) const;
220
224 bool
225 operator!=(const Iterator &) const;
226
232 bool
233 operator<(const Iterator &) const;
234
239 size_type,
240 size_type,
241 << "Attempt to access element " << arg2 << " of row "
242 << arg1 << " which doesn't have that many elements.");
243
244 private:
249
251 };
252
253 } // namespace SparsityPatternIterators
254
255
274 {
275 public:
280
285
294
310 const size_type n,
311 const size_type n_entries_per_row = 0);
312
322 const size_type n,
323 const std::vector<size_type> &n_entries_per_row);
324
329 SparsityPattern(SparsityPattern &&other) noexcept;
330
335 SparsityPattern(const SparsityPattern &input_sparsity_pattern);
336
340 virtual ~SparsityPattern() override = default;
341
353 void
354 reinit(const size_type m,
355 const size_type n,
356 const size_type n_entries_per_row = 0);
357
366 void
367 reinit(const size_type m,
368 const size_type n,
369 const std::vector<size_type> &n_entries_per_row);
370
375 void
376 copy_from(const SparsityPattern &input_sparsity_pattern);
377
383 template <typename SparsityPatternType>
384 void
385 copy_from(const SparsityPatternType &nontrilinos_sparsity_pattern);
386
394 operator=(const SparsityPattern &input_sparsity_pattern);
395
403 void
404 clear();
405
415 void
416 compress();
435 SparsityPattern(const IndexSet &parallel_partitioning,
436 const MPI_Comm communicator = MPI_COMM_WORLD,
437 const size_type n_entries_per_row = 0);
438
449 SparsityPattern(const IndexSet &parallel_partitioning,
450 const MPI_Comm communicator,
451 const std::vector<size_type> &n_entries_per_row);
452
467 SparsityPattern(const IndexSet &row_parallel_partitioning,
468 const IndexSet &col_parallel_partitioning,
469 const MPI_Comm communicator = MPI_COMM_WORLD,
470 const size_type n_entries_per_row = 0);
471
483 SparsityPattern(const IndexSet &row_parallel_partitioning,
484 const IndexSet &col_parallel_partitioning,
485 const MPI_Comm communicator,
486 const std::vector<size_type> &n_entries_per_row);
487
514 SparsityPattern(const IndexSet &row_parallel_partitioning,
515 const IndexSet &col_parallel_partitioning,
516 const IndexSet &writable_rows,
517 const MPI_Comm communicator = MPI_COMM_WORLD,
518 const size_type n_entries_per_row = 0);
519
535 void
536 reinit(const IndexSet &parallel_partitioning,
537 const MPI_Comm communicator = MPI_COMM_WORLD,
538 const size_type n_entries_per_row = 0);
539
550 void
551 reinit(const IndexSet &parallel_partitioning,
552 const MPI_Comm communicator,
553 const std::vector<size_type> &n_entries_per_row);
554
571 void
572 reinit(const IndexSet &row_parallel_partitioning,
573 const IndexSet &col_parallel_partitioning,
574 const MPI_Comm communicator = MPI_COMM_WORLD,
575 const size_type n_entries_per_row = 0);
576
602 void
603 reinit(const IndexSet &row_parallel_partitioning,
604 const IndexSet &col_parallel_partitioning,
605 const IndexSet &writeable_rows,
606 const MPI_Comm communicator = MPI_COMM_WORLD,
607 const size_type n_entries_per_row = 0);
608
613 void
614 reinit(const IndexSet &row_parallel_partitioning,
615 const IndexSet &col_parallel_partitioning,
616 const MPI_Comm communicator,
617 const std::vector<size_type> &n_entries_per_row);
618
628 template <typename SparsityPatternType>
629 void
630 reinit(const IndexSet &row_parallel_partitioning,
631 const IndexSet &col_parallel_partitioning,
632 const SparsityPatternType &nontrilinos_sparsity_pattern,
633 const MPI_Comm communicator = MPI_COMM_WORLD,
634 const bool exchange_data = false);
635
644 template <typename SparsityPatternType>
645 void
646 reinit(const IndexSet &parallel_partitioning,
647 const SparsityPatternType &nontrilinos_sparsity_pattern,
648 const MPI_Comm communicator = MPI_COMM_WORLD,
649 const bool exchange_data = false);
660 bool
662
666 unsigned int
667 max_entries_per_row() const;
668
678 unsigned int
679 local_size() const;
680
689 std::pair<size_type, size_type>
690 local_range() const;
691
696 bool
697 in_local_range(const size_type index) const;
698
702 std::uint64_t
703 n_nonzero_elements() const;
704
714 row_length(const size_type row) const;
715
723 bandwidth() const;
724
729 bool
730 empty() const;
731
736 bool
737 exists(const size_type i, const size_type j) const;
738
743 bool
744 row_is_stored_locally(const size_type i) const;
745
750 std::size_t
751 memory_consumption() const;
752
761 void
762 add(const size_type i, const size_type j);
763
764
768 template <typename ForwardIterator>
769 void
771 ForwardIterator begin,
772 ForwardIterator end,
773 const bool indices_are_sorted = false);
774
775 virtual void
776 add_row_entries(const size_type &row,
777 const ArrayView<const size_type> &columns,
778 const bool indices_are_sorted = false) override;
779
781
792 const Epetra_FECrsGraph &
794
801 const Epetra_Map &
802 domain_partitioner() const;
803
810 const Epetra_Map &
811 range_partitioner() const;
812
817 get_mpi_communicator() const;
818
833
841
853 begin() const;
854
859 end() const;
860
870 begin(const size_type r) const;
871
881 end(const size_type r) const;
882
894 void
895 write_ascii();
896
904 void
905 print(std::ostream &out,
906 const bool write_extended_trilinos_info = false) const;
907
922 void
923 print_gnuplot(std::ostream &out) const;
924
934 int,
935 << "An error with error number " << arg1
936 << " occurred while calling a Trilinos function");
937
942 size_type,
943 size_type,
944 << "The entry with index <" << arg1 << ',' << arg2
945 << "> does not exist.");
946
951 size_type,
952 size_type,
953 size_type,
954 size_type,
955 << "You tried to access element (" << arg1 << '/' << arg2
956 << ')'
957 << " of a distributed matrix, but only rows in range ["
958 << arg3 << ',' << arg4
959 << "] are stored locally and can be accessed.");
960
965 size_type,
966 size_type,
967 << "You tried to access element (" << arg1 << '/' << arg2
968 << ')' << " of a sparse matrix, but it appears to not"
969 << " exist in the Trilinos sparsity pattern.");
971 private:
976 std::unique_ptr<Epetra_Map> column_space_map;
977
983 std::unique_ptr<Epetra_FECrsGraph> graph;
984
991 std::unique_ptr<Epetra_CrsGraph> nonlocal_graph;
992
996 };
997
998
999
1000 // ----------------------- inline and template functions --------------------
1001
1002
1003# ifndef DOXYGEN
1004
1005 namespace SparsityPatternIterators
1006 {
1007 inline Accessor::Accessor(const SparsityPattern *sp,
1008 const size_type row,
1009 const size_type index)
1010 : sparsity_pattern(const_cast<SparsityPattern *>(sp))
1011 , a_row(row)
1012 , a_index(index)
1013 {
1014 visit_present_row();
1015 }
1016
1017
1018
1019 inline Accessor::size_type
1020 Accessor::row() const
1021 {
1022 Assert(a_row < sparsity_pattern->n_rows(),
1023 ExcBeyondEndOfSparsityPattern());
1024 return a_row;
1025 }
1026
1027
1028
1029 inline Accessor::size_type
1030 Accessor::column() const
1031 {
1032 Assert(a_row < sparsity_pattern->n_rows(),
1033 ExcBeyondEndOfSparsityPattern());
1034 return (*colnum_cache)[a_index];
1035 }
1036
1037
1038
1039 inline Accessor::size_type
1040 Accessor::index() const
1041 {
1042 Assert(a_row < sparsity_pattern->n_rows(),
1043 ExcBeyondEndOfSparsityPattern());
1044 return a_index;
1045 }
1046
1047
1048
1049 inline Iterator::Iterator(const SparsityPattern *sp,
1050 const size_type row,
1051 const size_type index)
1052 : accessor(sp, row, index)
1053 {}
1054
1055
1056
1057 inline Iterator::Iterator(const Iterator &) = default;
1058
1059
1060
1061 inline Iterator &
1062 Iterator::operator++()
1063 {
1064 Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1066
1067 ++accessor.a_index;
1068
1069 // If at end of line: do one step, then cycle until we find a row with a
1070 // nonzero number of entries that is stored locally.
1071 if (accessor.a_index >= accessor.colnum_cache->size())
1072 {
1073 accessor.a_index = 0;
1074 ++accessor.a_row;
1075
1076 while (accessor.a_row < accessor.sparsity_pattern->n_rows())
1077 {
1078 const auto row_length =
1079 accessor.sparsity_pattern->row_length(accessor.a_row);
1080 if (row_length == 0 ||
1081 !accessor.sparsity_pattern->row_is_stored_locally(
1082 accessor.a_row))
1083 ++accessor.a_row;
1084 else
1085 break;
1086 }
1087
1088 accessor.visit_present_row();
1089 }
1090 return *this;
1091 }
1092
1093
1094
1095 inline Iterator
1096 Iterator::operator++(int)
1097 {
1098 const Iterator old_state = *this;
1099 ++(*this);
1100 return old_state;
1101 }
1102
1103
1104
1105 inline const Accessor &
1106 Iterator::operator*() const
1107 {
1108 return accessor;
1109 }
1110
1111
1112
1113 inline const Accessor *
1114 Iterator::operator->() const
1115 {
1116 return &accessor;
1117 }
1118
1119
1120
1121 inline bool
1122 Iterator::operator==(const Iterator &other) const
1123 {
1124 return (accessor.a_row == other.accessor.a_row &&
1125 accessor.a_index == other.accessor.a_index);
1126 }
1127
1128
1129
1130 inline bool
1131 Iterator::operator!=(const Iterator &other) const
1132 {
1133 return !(*this == other);
1134 }
1135
1136
1137
1138 inline bool
1139 Iterator::operator<(const Iterator &other) const
1140 {
1141 return (accessor.row() < other.accessor.row() ||
1142 (accessor.row() == other.accessor.row() &&
1143 accessor.index() < other.accessor.index()));
1144 }
1145
1146 } // namespace SparsityPatternIterators
1147
1148
1149
1152 {
1153 const size_type first_valid_row = this->local_range().first;
1154 return const_iterator(this, first_valid_row, 0);
1155 }
1156
1157
1158
1160 SparsityPattern::end() const
1161 {
1162 return const_iterator(this, n_rows(), 0);
1163 }
1164
1165
1166
1168 SparsityPattern::begin(const size_type r) const
1169 {
1171 if (row_length(r) > 0)
1172 return const_iterator(this, r, 0);
1173 else
1174 return end(r);
1175 }
1176
1177
1178
1180 SparsityPattern::end(const size_type r) const
1181 {
1183
1184 // place the iterator on the first entry
1185 // past this line, or at the end of the
1186 // matrix
1187 for (size_type i = r + 1; i < n_rows(); ++i)
1188 if (row_length(i) > 0)
1189 return const_iterator(this, i, 0);
1190
1191 // if there is no such line, then take the
1192 // end iterator of the matrix
1193 return end();
1194 }
1195
1196
1197
1198 inline bool
1199 SparsityPattern::in_local_range(const size_type index) const
1200 {
1202# ifndef DEAL_II_WITH_64BIT_INDICES
1203 begin = graph->RowMap().MinMyGID();
1204 end = graph->RowMap().MaxMyGID() + 1;
1205# else
1206 begin = graph->RowMap().MinMyGID64();
1207 end = graph->RowMap().MaxMyGID64() + 1;
1208# endif
1209
1210 return ((index >= static_cast<size_type>(begin)) &&
1211 (index < static_cast<size_type>(end)));
1212 }
1213
1214
1215
1216 inline bool
1218 {
1219 return graph->Filled();
1220 }
1221
1222
1223
1224 inline bool
1226 {
1227 return ((n_rows() == 0) && (n_cols() == 0));
1228 }
1229
1230
1231
1232 inline void
1233 SparsityPattern::add(const size_type i, const size_type j)
1234 {
1235 add_entries(i, &j, &j + 1);
1236 }
1237
1238
1239
1240 template <typename ForwardIterator>
1241 inline void
1242 SparsityPattern::add_entries(const size_type row,
1243 ForwardIterator begin,
1244 ForwardIterator end,
1245 const bool /*indices_are_sorted*/)
1246 {
1247 if (begin == end)
1248 return;
1249
1250 // verify that the size of the data type Trilinos expects matches that the
1251 // iterator points to. we allow for some slippage between signed and
1252 // unsigned and only compare that they are both either 32 or 64 bit. to
1253 // write this test properly, not that we cannot compare the size of
1254 // '*begin' because 'begin' may be an iterator and '*begin' may be an
1255 // accessor class. consequently, we need to somehow get an actual value
1256 // from it which we can by evaluating an expression such as when
1257 // multiplying the value produced by 2
1258 Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1260
1261 TrilinosWrappers::types::int_type *col_index_ptr =
1262 reinterpret_cast<TrilinosWrappers::types::int_type *>(
1263 const_cast<std::decay_t<decltype(*begin)> *>(&*begin));
1264 // Check at least for the first index that the conversion actually works
1265 AssertDimension(*col_index_ptr, *begin);
1266 TrilinosWrappers::types::int_type trilinos_row_index = row;
1267 const int n_cols = static_cast<int>(end - begin);
1268
1269 int ierr;
1270 if (row_is_stored_locally(row))
1271 ierr =
1272 graph->InsertGlobalIndices(trilinos_row_index, n_cols, col_index_ptr);
1273 else if (nonlocal_graph.get() != nullptr)
1274 {
1275 // this is the case when we have explicitly set the off-processor rows
1276 // and want to create a separate matrix object for them (to retain
1277 // thread-safety)
1278 Assert(nonlocal_graph->RowMap().LID(
1279 static_cast<TrilinosWrappers::types::int_type>(row)) != -1,
1280 ExcMessage("Attempted to write into off-processor matrix row "
1281 "that has not be specified as being writable upon "
1282 "initialization"));
1283 ierr = nonlocal_graph->InsertGlobalIndices(trilinos_row_index,
1284 n_cols,
1285 col_index_ptr);
1286 }
1287 else
1288 ierr = graph->InsertGlobalIndices(1,
1289 &trilinos_row_index,
1290 n_cols,
1291 col_index_ptr);
1292
1293 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1294 }
1295
1296
1297
1298 inline const Epetra_FECrsGraph &
1299 SparsityPattern::trilinos_sparsity_pattern() const
1300 {
1301 return *graph;
1302 }
1303
1304
1305
1306 inline IndexSet
1307 SparsityPattern::locally_owned_domain_indices() const
1308 {
1309 return IndexSet(graph->DomainMap());
1310 }
1311
1312
1313
1314 inline IndexSet
1315 SparsityPattern::locally_owned_range_indices() const
1316 {
1317 return IndexSet(graph->RangeMap());
1318 }
1319
1320# endif // DOXYGEN
1321} // namespace TrilinosWrappers
1322
1323
1325
1326
1327#endif // DEAL_II_WITH_TRILINOS
1328
1329#endif
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
size_type n_cols() const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
bool is_compressed() const
iterator begin() const
void add(const size_type i, const size_type j)
SparsityPatternIterators::Iterator const_iterator
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
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
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
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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcBeyondEndOfSparsityPattern()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:585
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcAccessToNonlocalRow(size_type arg1, size_type arg2, size_type arg3)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define AssertIndexRange(index, range)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition exceptions.h:562
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
unsigned int global_dof_index
Definition types.h:81