Reference documentation for deal.II version 9.6.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\}}\)
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
26
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
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();
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
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
715 row_length(const size_type row) const;
716
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
818 get_mpi_communicator() const;
819
834
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 &
1063 Iterator::operator++()
1064 {
1065 Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
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
1097 Iterator::operator++(int)
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
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 {
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 {
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
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
1243 SparsityPattern::add_entries(const size_type row,
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),
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 " +
1283 " that has not been specified as being writable upon "
1284 "initialization."));
1285 ierr = nonlocal_graph->InsertGlobalIndices(trilinos_row_index,
1286 n_cols,
1287 col_index_ptr);
1288 }
1289 else
1290 ierr = graph->InsertGlobalIndices(1,
1291 &trilinos_row_index,
1292 n_cols,
1293 col_index_ptr);
1294
1295 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
1296 }
1297
1298
1299
1300 inline const Epetra_FECrsGraph &
1301 SparsityPattern::trilinos_sparsity_pattern() const
1302 {
1303 return *graph;
1304 }
1305
1306
1307
1308 inline IndexSet
1309 SparsityPattern::locally_owned_domain_indices() const
1310 {
1311 return IndexSet(graph->DomainMap());
1312 }
1313
1314
1315
1316 inline IndexSet
1317 SparsityPattern::locally_owned_range_indices() const
1318 {
1319 return IndexSet(graph->RangeMap());
1320 }
1321
1322# endif // DOXYGEN
1323} // namespace TrilinosWrappers
1324
1325
1327
1328
1329#endif // DEAL_II_WITH_TRILINOS
1330
1331#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
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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:479
unsigned int global_dof_index
Definition types.h:81