Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2740-g4da21c751f 2025-03-02 17:30:00+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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
trilinos_tpetra_sparsity_pattern.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 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_tpetra_sparsity_pattern_h
16#define dealii_trilinos_tpetra_sparsity_pattern_h
17
18#include <deal.II/base/config.h>
19
20#include <deal.II/base/types.h>
21
23
24#ifdef DEAL_II_TRILINOS_WITH_TPETRA
25
29
32
33# include <Tpetra_CrsGraph.hpp>
34
35# include <cmath>
36# include <memory>
37# include <vector>
38
39
41
42// forward declarations
43# ifndef DOXYGEN
45
46namespace LinearAlgebra
47{
48 namespace TpetraWrappers
49 {
50 template <typename MemorySpace>
51 class SparsityPattern;
52
53 template <typename Number, typename MemorySpace>
54 class SparseMatrix;
55
57 {
58 template <typename MemorySpace>
59 class Iterator;
60 }
61 } // namespace TpetraWrappers
62} // namespace LinearAlgebra
63# endif
64
65namespace LinearAlgebra
66{
67 namespace TpetraWrappers
68 {
70 {
82 template <typename MemorySpace = ::MemorySpace::Host>
84 {
85 public:
90
95 const size_type row,
96 const size_type index);
97
102 row() const;
103
108 index() const;
109
114 column() const;
115
120
125 size_type,
126 size_type,
127 size_type,
128 << "You tried to access row " << arg1
129 << " of a distributed sparsity pattern, "
130 << " but only rows " << arg2 << " through " << arg3
131 << " are stored locally and can be accessed.");
132
133 private:
138
143
148
161 std::shared_ptr<std::vector<::types::signed_global_dof_index>>
163
169 void
171
172 // Make enclosing class a friend.
173 friend class Iterator<MemorySpace>;
174 };
175
181 template <typename MemorySpace = ::MemorySpace::Host>
183 {
184 public:
189
194 Iterator(const SparsityPattern<MemorySpace> *sparsity_pattern,
195 const size_type row,
196 const size_type index);
197
202
208
214
219 operator*() const;
220
225 operator->() const;
226
231 bool
233
237 bool
239
245 bool
246 operator<(const Iterator<MemorySpace> &) const;
247
252 size_type,
253 size_type,
254 << "Attempt to access element " << arg2 << " of row "
255 << arg1 << " which doesn't have that many elements.");
256
257 private:
262
264 };
265
266 } // namespace SparsityPatternIterators
267
268
286 template <typename MemorySpace = ::MemorySpace::Host>
288 {
289 public:
298
307
323 const size_type n,
324 const size_type n_entries_per_row = 0);
325
335 const size_type n,
336 const std::vector<size_type> &n_entries_per_row);
337
343
349 const SparsityPattern<MemorySpace> &input_sparsity_pattern);
350
354 virtual ~SparsityPattern() override = default;
355
367 void
368 reinit(const size_type m,
369 const size_type n,
370 const size_type n_entries_per_row = 0);
371
380 void
381 reinit(const size_type m,
382 const size_type n,
383 const std::vector<size_type> &n_entries_per_row);
384
389 void
390 copy_from(const SparsityPattern<MemorySpace> &input_sparsity_pattern);
391
397 template <typename SparsityPatternType>
398 void
399 copy_from(const SparsityPatternType &nontrilinos_sparsity_pattern);
400
408 operator=(const SparsityPattern<MemorySpace> &input_sparsity_pattern);
409
417 void
418 clear();
419
429 void
430 compress();
449 SparsityPattern(const IndexSet &parallel_partitioning,
450 const MPI_Comm communicator = MPI_COMM_WORLD,
451 const size_type n_entries_per_row = 0);
452
463 SparsityPattern(const IndexSet &parallel_partitioning,
464 const MPI_Comm communicator,
465 const std::vector<size_type> &n_entries_per_row);
466
481 SparsityPattern(const IndexSet &row_parallel_partitioning,
482 const IndexSet &col_parallel_partitioning,
483 const MPI_Comm communicator = MPI_COMM_WORLD,
484 const size_type n_entries_per_row = 0);
485
497 SparsityPattern(const IndexSet &row_parallel_partitioning,
498 const IndexSet &col_parallel_partitioning,
499 const MPI_Comm communicator,
500 const std::vector<size_type> &n_entries_per_row);
501
528 SparsityPattern(const IndexSet &row_parallel_partitioning,
529 const IndexSet &col_parallel_partitioning,
530 const IndexSet &writable_rows,
531 const MPI_Comm communicator = MPI_COMM_WORLD,
532 const size_type n_entries_per_row = 0);
533
549 void
550 reinit(const IndexSet &parallel_partitioning,
551 const MPI_Comm communicator = MPI_COMM_WORLD,
552 const size_type n_entries_per_row = 0);
553
564 void
565 reinit(const IndexSet &parallel_partitioning,
566 const MPI_Comm communicator,
567 const std::vector<size_type> &n_entries_per_row);
568
585 void
586 reinit(const IndexSet &row_parallel_partitioning,
587 const IndexSet &col_parallel_partitioning,
588 const MPI_Comm communicator = MPI_COMM_WORLD,
589 const size_type n_entries_per_row = 0);
590
616 void
617 reinit(const IndexSet &row_parallel_partitioning,
618 const IndexSet &col_parallel_partitioning,
619 const IndexSet &writeable_rows,
620 const MPI_Comm communicator = MPI_COMM_WORLD,
621 const size_type n_entries_per_row = 0);
622
627 void
628 reinit(const IndexSet &row_parallel_partitioning,
629 const IndexSet &col_parallel_partitioning,
630 const MPI_Comm communicator,
631 const std::vector<size_type> &n_entries_per_row);
632
642 template <typename SparsityPatternType>
643 void
644 reinit(const IndexSet &row_parallel_partitioning,
645 const IndexSet &col_parallel_partitioning,
646 const SparsityPatternType &nontrilinos_sparsity_pattern,
647 const MPI_Comm communicator = MPI_COMM_WORLD,
648 const bool exchange_data = false);
649
658 template <typename SparsityPatternType>
659 void
660 reinit(const IndexSet &parallel_partitioning,
661 const SparsityPatternType &nontrilinos_sparsity_pattern,
662 const MPI_Comm communicator = MPI_COMM_WORLD,
663 const bool exchange_data = false);
674 bool
676
680 unsigned int
681 max_entries_per_row() const;
682
692 unsigned int
693 local_size() const;
694
703 std::pair<size_type, size_type>
704 local_range() const;
705
710 bool
711 in_local_range(const size_type index) const;
712
716 std::uint64_t
717 n_nonzero_elements() const;
718
728 row_length(const size_type row) const;
729
737 bandwidth() const;
738
743 bool
744 empty() const;
745
750 bool
751 exists(const size_type i, const size_type j) const;
752
757 bool
758 row_is_stored_locally(const size_type i) const;
759
764 std::size_t
765 memory_consumption() const;
766
775 void
776 add(const size_type i, const size_type j);
777
778
782 template <typename ForwardIterator>
783 void
785 ForwardIterator begin,
786 ForwardIterator end,
787 const bool indices_are_sorted = false);
788
789 virtual void
791 const ::types::global_dof_index &row,
793 const bool indices_are_sorted = false) override;
794
796
807 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>>
809
816 Teuchos::RCP<const TpetraTypes::MapType<MemorySpace>>
817 domain_partitioner() const;
818
825 Teuchos::RCP<const TpetraTypes::MapType<MemorySpace>>
826 range_partitioner() const;
827
832 get_mpi_communicator() const;
833
837 Teuchos::RCP<const Teuchos::Comm<int>>
839
854
862
874 begin() const;
875
880 end() const;
881
891 begin(const size_type r) const;
892
902 end(const size_type r) const;
903
917 void
918 print(std::ostream &out,
919 const bool write_extended_trilinos_info = false) const;
920
935 void
936 print_gnuplot(std::ostream &out) const;
937
947 int,
948 << "An error with error number " << arg1
949 << " occurred while calling a Trilinos function");
950
955 size_type,
956 size_type,
957 << "The entry with index <" << arg1 << ',' << arg2
958 << "> does not exist.");
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 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> column_space_map;
990
996 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> graph;
997
1004 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> nonlocal_graph;
1005
1006 // TODO: currently only for double
1007 friend class SparseMatrix<double, MemorySpace>;
1010 };
1011
1012
1013
1014 // ---------------- inline and template functions -----------------
1015
1016
1017# ifndef DOXYGEN
1018
1019 namespace SparsityPatternIterators
1020 {
1021 template <typename MemorySpace>
1024 const size_type row,
1025 const size_type index)
1026 : sparsity_pattern(const_cast<SparsityPattern<MemorySpace> *>(sp))
1027 , a_row(row)
1028 , a_index(index)
1029 {
1030 visit_present_row();
1031 }
1032
1033
1034
1035 template <typename MemorySpace>
1036 inline typename Accessor<MemorySpace>::size_type
1037 Accessor<MemorySpace>::row() const
1038 {
1039 Assert(a_row < sparsity_pattern->n_rows(),
1040 ExcBeyondEndOfSparsityPattern());
1041 return a_row;
1042 }
1043
1044
1045
1046 template <typename MemorySpace>
1047 inline typename Accessor<MemorySpace>::size_type
1048 Accessor<MemorySpace>::column() const
1049 {
1050 Assert(a_row < sparsity_pattern->n_rows(),
1051 ExcBeyondEndOfSparsityPattern());
1052 return (*colnum_cache)[a_index];
1053 }
1054
1055
1056
1057 template <typename MemorySpace>
1058 inline typename Accessor<MemorySpace>::size_type
1059 Accessor<MemorySpace>::index() const
1060 {
1061 Assert(a_row < sparsity_pattern->n_rows(),
1062 ExcBeyondEndOfSparsityPattern());
1063 return a_index;
1064 }
1065
1066
1067
1068 template <typename MemorySpace>
1069 inline Iterator<MemorySpace>::Iterator(
1071 const size_type row,
1072 const size_type index)
1073 : accessor(sp, row, index)
1074 {}
1075
1076
1077
1078 template <typename MemorySpace>
1079 inline Iterator<MemorySpace>::Iterator(const Iterator<MemorySpace> &) =
1080 default;
1081
1082
1083
1084 template <typename MemorySpace>
1085 inline Iterator<MemorySpace> &
1086 Iterator<MemorySpace>::operator++()
1087 {
1088 Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1090
1091 ++accessor.a_index;
1092
1093 // If at end of line: do one step, then cycle until we find a row with a
1094 // nonzero number of entries that is stored locally.
1095 if (accessor.a_index >=
1096 static_cast<::types::signed_global_dof_index>(
1097 accessor.colnum_cache->size()))
1098 {
1099 accessor.a_index = 0;
1100 ++accessor.a_row;
1101
1102 while (accessor.a_row <
1103 static_cast<::types::signed_global_dof_index>(
1104 accessor.sparsity_pattern->n_rows()))
1105 {
1106 const auto row_length =
1107 accessor.sparsity_pattern->row_length(accessor.a_row);
1108 if (row_length == 0 ||
1109 !accessor.sparsity_pattern->row_is_stored_locally(
1110 accessor.a_row))
1111 ++accessor.a_row;
1112 else
1113 break;
1114 }
1115
1116 accessor.visit_present_row();
1117 }
1118 return *this;
1119 }
1120
1121
1122
1123 template <typename MemorySpace>
1124 inline Iterator<MemorySpace>
1125 Iterator<MemorySpace>::operator++(int)
1126 {
1127 const Iterator<MemorySpace> old_state = *this;
1128 ++(*this);
1129 return old_state;
1130 }
1131
1132
1133
1134 template <typename MemorySpace>
1135 inline const Accessor<MemorySpace> &
1136 Iterator<MemorySpace>::operator*() const
1137 {
1138 return accessor;
1139 }
1140
1141
1142
1143 template <typename MemorySpace>
1144 inline const Accessor<MemorySpace> *
1145 Iterator<MemorySpace>::operator->() const
1146 {
1147 return &accessor;
1148 }
1149
1150
1151
1152 template <typename MemorySpace>
1153 inline bool
1154 Iterator<MemorySpace>::operator==(
1155 const Iterator<MemorySpace> &other) const
1156 {
1157 return (accessor.a_row == other.accessor.a_row &&
1158 accessor.a_index == other.accessor.a_index);
1159 }
1160
1161
1162
1163 template <typename MemorySpace>
1164 inline bool
1165 Iterator<MemorySpace>::operator!=(
1166 const Iterator<MemorySpace> &other) const
1167 {
1168 return !(*this == other);
1169 }
1170
1171
1172
1173 template <typename MemorySpace>
1174 inline bool
1175 Iterator<MemorySpace>::operator<(const Iterator<MemorySpace> &other) const
1176 {
1177 return (accessor.row() < other.accessor.row() ||
1178 (accessor.row() == other.accessor.row() &&
1179 accessor.index() < other.accessor.index()));
1180 }
1181
1182 } // namespace SparsityPatternIterators
1183
1184
1185
1186 template <typename MemorySpace>
1189 {
1190 const size_type first_valid_row = this->local_range().first;
1191 return const_iterator(this, first_valid_row, 0);
1192 }
1193
1194
1195
1196 template <typename MemorySpace>
1199 {
1200 return const_iterator(this, n_rows(), 0);
1201 }
1202
1203
1204
1205 template <typename MemorySpace>
1207 SparsityPattern<MemorySpace>::begin(const size_type r) const
1208 {
1209 AssertIndexRange(r, n_rows());
1210 if (row_length(r) > 0)
1211 return const_iterator(this, r, 0);
1212 else
1213 return end(r);
1214 }
1215
1216
1217
1218 template <typename MemorySpace>
1220 SparsityPattern<MemorySpace>::end(const size_type r) const
1221 {
1222 AssertIndexRange(r, n_rows());
1223
1224 // place the iterator on the first entry
1225 // past this line, or at the end of the
1226 // matrix
1227 for (size_type i = r + 1; i < n_rows(); ++i)
1228 if (row_length(i) > 0)
1229 return const_iterator(this, i, 0);
1230
1231 // if there is no such line, then take the
1232 // end iterator of the matrix
1233 return end();
1234 }
1235
1236
1237
1238 template <typename MemorySpace>
1239 inline bool
1240 SparsityPattern<MemorySpace>::in_local_range(const size_type index) const
1241 {
1243 graph->getRowMap()->getMinGlobalIndex();
1245 graph->getRowMap()->getMaxGlobalIndex() + 1;
1246
1247 return ((index >= static_cast<size_type>(begin)) &&
1248 (index < static_cast<size_type>(end)));
1249 }
1250
1251
1252
1253 template <typename MemorySpace>
1254 inline bool
1256 {
1257 return graph->isFillComplete();
1258 }
1259
1260
1261
1262 template <typename MemorySpace>
1263 inline bool
1265 {
1266 return ((n_rows() == 0) && (n_cols() == 0));
1267 }
1268
1269
1270
1271 template <typename MemorySpace>
1272 inline void
1273 SparsityPattern<MemorySpace>::add(const size_type i, const size_type j)
1274 {
1275 add_entries(i, &j, &j + 1);
1276 }
1277
1278
1279
1280 template <typename MemorySpace>
1281 template <typename ForwardIterator>
1282 inline void
1284 ForwardIterator begin,
1285 ForwardIterator end,
1286 const bool /*indices_are_sorted*/)
1287 {
1288 if (begin == end)
1289 return;
1290
1291 // verify that the size of the data type Trilinos expects matches that the
1292 // iterator points to. we allow for some slippage between signed and
1293 // unsigned and only compare that they are both either 32 or 64 bit. to
1294 // write this test properly, not that we cannot compare the size of
1295 // '*begin' because 'begin' may be an iterator and '*begin' may be an
1296 // accessor class. consequently, we need to somehow get an actual value
1297 // from it which we can by evaluating an expression such as when
1298 // multiplying the value produced by 2
1299 Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1301
1302 const TrilinosWrappers::types::int_type *col_index_ptr_begin =
1303 reinterpret_cast<TrilinosWrappers::types::int_type *>(
1304 const_cast<std::decay_t<decltype(*begin)> *>(&*begin));
1305
1306 const TrilinosWrappers::types::int_type *col_index_ptr_end =
1307 reinterpret_cast<TrilinosWrappers::types::int_type *>(
1308 const_cast<std::decay_t<decltype(*end)> *>(&*end));
1309
1310 // Check at least for the first index that the conversion actually works
1311 AssertDimension(*col_index_ptr_begin, *begin);
1312 AssertDimension(*col_index_ptr_end, *end);
1313 TrilinosWrappers::types::int_type trilinos_row_index = row;
1314
1315 // TODO: The following line creates an array by copying the entries.
1316 // Perhaps there is a way to only create a 'view' of these arrays
1317 // and pass that to Tpetra?
1318 Teuchos::Array<TrilinosWrappers::types::int_type> array(
1319 col_index_ptr_begin, col_index_ptr_end);
1320
1321 if (row_is_stored_locally(row))
1322 graph->insertGlobalIndices(trilinos_row_index, array());
1323 else if (nonlocal_graph.get() != nullptr)
1324 {
1325 // this is the case when we have explicitly set the off-processor rows
1326 // and want to create a separate matrix object for them (to retain
1327 // thread-safety)
1328 Assert(nonlocal_graph->getRowMap()->getLocalElement(row) !=
1329 Teuchos::OrdinalTraits<
1331 ExcMessage("Attempted to write into off-processor matrix row "
1332 "that has not be specified as being writable upon "
1333 "initialization"));
1334 nonlocal_graph->insertGlobalIndices(trilinos_row_index, array);
1335 }
1336 else
1337 graph->insertGlobalIndices(trilinos_row_index, array);
1338 }
1339
1340
1341
1342 template <typename MemorySpace>
1343 inline Teuchos::RCP<Tpetra::CrsGraph<int,
1345 TpetraTypes::NodeType<MemorySpace>>>
1347 {
1348 return graph;
1349 }
1350
1351
1352
1353 template <typename MemorySpace>
1354 inline IndexSet
1356 {
1357 return IndexSet(graph->getDomainMap().getConst());
1358 }
1359
1360
1361
1362 template <typename MemorySpace>
1363 inline IndexSet
1365 {
1366 return IndexSet(graph->getRangeMap().getConst());
1367 }
1368
1369# endif // DOXYGEN
1370 } // namespace TpetraWrappers
1371
1372} // namespace LinearAlgebra
1373
1374
1376
1377
1378#endif // DEAL_II_TRILINOS_WITH_TPETRA
1379
1380#endif
Accessor(const SparsityPattern< MemorySpace > *sparsity_pattern, const size_type row, const size_type index)
std::shared_ptr< std::vector<::types::signed_global_dof_index > > colnum_cache
Iterator(const SparsityPattern< MemorySpace > *sparsity_pattern, const size_type row, const size_type index)
bool operator!=(const Iterator< MemorySpace > &) const
bool operator<(const Iterator< MemorySpace > &) const
bool operator==(const Iterator< MemorySpace > &) const
Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > graph
Teuchos::RCP< const Teuchos::Comm< int > > get_teuchos_mpi_communicator() const
Teuchos::RCP< const TpetraTypes::MapType< MemorySpace > > range_partitioner() const
Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > trilinos_sparsity_pattern() const
void copy_from(const SparsityPattern< MemorySpace > &input_sparsity_pattern)
SparsityPattern< MemorySpace > & operator=(const SparsityPattern< MemorySpace > &input_sparsity_pattern)
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
Teuchos::RCP< TpetraTypes::GraphType< MemorySpace > > nonlocal_graph
Teuchos::RCP< const TpetraTypes::MapType< MemorySpace > > domain_partitioner() const
const_iterator end(const size_type r) const
bool exists(const size_type i, const size_type j) const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
Teuchos::RCP< TpetraTypes::MapType< MemorySpace > > column_space_map
virtual void add_row_entries(const ::types::global_dof_index &row, const ArrayView< const ::types::global_dof_index > &columns, const bool indices_are_sorted=false) override
void add(const size_type i, const size_type j)
const_iterator begin(const size_type r) const
virtual ~SparsityPattern() override=default
bool in_local_range(const size_type index) const
void reinit(const size_type m, const size_type n, const size_type n_entries_per_row=0)
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
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)
iterator end() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException0(Exception0)
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcIteratorPastEnd()
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcAccessToNonlocalRow(size_type arg1, size_type arg2, size_type arg3)
#define DeclException3(Exception3, type1, type2, type3, outsequence)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::pair< types::global_dof_index, types::global_dof_index > local_range
Definition mpi.cc:824
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
int signed_global_dof_index
Definition types.h:101
unsigned int global_dof_index
Definition types.h:90