Reference documentation for deal.II version GIT relicensing-489-g2d48aca8cc 2024-04-28 17:30:02+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_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#ifdef DEAL_II_TRILINOS_WITH_TPETRA
21
25
29
30# include <Tpetra_CrsGraph.hpp>
31
32# include <cmath>
33# include <memory>
34# include <vector>
35
36
38
39// forward declarations
40# ifndef DOXYGEN
42
43namespace LinearAlgebra
44{
45 namespace TpetraWrappers
46 {
47 template <typename MemorySpace>
48 class SparsityPattern;
49
50 template <typename Number, typename MemorySpace>
51 class SparseMatrix;
52
54 {
55 template <typename MemorySpace>
56 class Iterator;
57 }
58 } // namespace TpetraWrappers
59} // namespace LinearAlgebra
60# endif
61
62namespace LinearAlgebra
63{
64 namespace TpetraWrappers
65 {
67 {
79 template <typename MemorySpace = ::MemorySpace::Host>
81 {
82 public:
87
91# if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0)
92 using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode<
93 typename MemorySpace::kokkos_space::execution_space,
94 typename MemorySpace::kokkos_space>;
95# else
96 using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode<
97 typename MemorySpace::kokkos_space::execution_space,
98 typename MemorySpace::kokkos_space>;
99# endif
100
105 const size_type row,
106 const size_type index);
107
112 row() const;
113
118 index() const;
119
124 column() const;
125
130
135 size_type,
136 size_type,
137 size_type,
138 << "You tried to access row " << arg1
139 << " of a distributed sparsity pattern, "
140 << " but only rows " << arg2 << " through " << arg3
141 << " are stored locally and can be accessed.");
142
143 private:
148
153
158
171 std::shared_ptr<std::vector<::types::signed_global_dof_index>>
173
179 void
181
182 // Make enclosing class a friend.
183 friend class Iterator<MemorySpace>;
184 };
185
191 template <typename MemorySpace = ::MemorySpace::Host>
193 {
194 public:
198 using size_type = size_t;
199
203# if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0)
204 using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode<
205 typename MemorySpace::kokkos_space::execution_space,
206 typename MemorySpace::kokkos_space>;
207# else
208 using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode<
209 typename MemorySpace::kokkos_space::execution_space,
210 typename MemorySpace::kokkos_space>;
211# endif
212
217 Iterator(const SparsityPattern<MemorySpace> *sparsity_pattern,
218 const size_type row,
219 const size_type index);
220
225
231
237
242 operator*() const;
243
248 operator->() const;
249
254 bool
256
260 bool
262
268 bool
269 operator<(const Iterator<MemorySpace> &) const;
270
275 size_type,
276 size_type,
277 << "Attempt to access element " << arg2 << " of row "
278 << arg1 << " which doesn't have that many elements.");
279
280 private:
285
287 };
288
289 } // namespace SparsityPatternIterators
290
291
309 template <typename MemorySpace = ::MemorySpace::Host>
311 {
312 public:
317
321# if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0)
322 using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode<
323 typename MemorySpace::kokkos_space::execution_space,
324 typename MemorySpace::kokkos_space>;
325# else
326 using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode<
327 typename MemorySpace::kokkos_space::execution_space,
328 typename MemorySpace::kokkos_space>;
329# endif
330
335
339 using MapType =
341
345 using GraphType =
347
348
357
373 const size_type n,
374 const size_type n_entries_per_row = 0);
375
385 const size_type n,
386 const std::vector<size_type> &n_entries_per_row);
387
393
400
404 virtual ~SparsityPattern() override = default;
405
417 void
418 reinit(const size_type m,
419 const size_type n,
420 const size_type n_entries_per_row = 0);
421
430 void
431 reinit(const size_type m,
432 const size_type n,
433 const std::vector<size_type> &n_entries_per_row);
434
439 void
441
447 template <typename SparsityPatternType>
448 void
450
459
467 void
468 clear();
469
479 void
480 compress();
500 const MPI_Comm communicator = MPI_COMM_WORLD,
501 const size_type n_entries_per_row = 0);
502
514 const MPI_Comm communicator,
515 const std::vector<size_type> &n_entries_per_row);
516
533 const MPI_Comm communicator = MPI_COMM_WORLD,
534 const size_type n_entries_per_row = 0);
535
549 const MPI_Comm communicator,
550 const std::vector<size_type> &n_entries_per_row);
551
580 const IndexSet &writable_rows,
581 const MPI_Comm communicator = MPI_COMM_WORLD,
582 const size_type n_entries_per_row = 0);
583
599 void
601 const MPI_Comm communicator = MPI_COMM_WORLD,
602 const size_type n_entries_per_row = 0);
603
614 void
616 const MPI_Comm communicator,
617 const std::vector<size_type> &n_entries_per_row);
618
635 void
638 const MPI_Comm communicator = MPI_COMM_WORLD,
639 const size_type n_entries_per_row = 0);
640
666 void
670 const MPI_Comm communicator = MPI_COMM_WORLD,
671 const size_type n_entries_per_row = 0);
672
677 void
680 const MPI_Comm communicator,
681 const std::vector<size_type> &n_entries_per_row);
682
692 template <typename SparsityPatternType>
693 void
697 const MPI_Comm communicator = MPI_COMM_WORLD,
698 const bool exchange_data = false);
699
708 template <typename SparsityPatternType>
709 void
712 const MPI_Comm communicator = MPI_COMM_WORLD,
713 const bool exchange_data = false);
724 bool
726
730 unsigned int
731 max_entries_per_row() const;
732
742 unsigned int
743 local_size() const;
744
753 std::pair<size_type, size_type>
754 local_range() const;
755
760 bool
761 in_local_range(const size_type index) const;
762
766 std::uint64_t
767 n_nonzero_elements() const;
768
778 row_length(const size_type row) const;
779
787 bandwidth() const;
788
793 bool
794 empty() const;
795
800 bool
801 exists(const size_type i, const size_type j) const;
802
807 bool
808 row_is_stored_locally(const size_type i) const;
809
814 std::size_t
815 memory_consumption() const;
816
825 void
826 add(const size_type i, const size_type j);
827
828
832 template <typename ForwardIterator>
833 void
837 const bool indices_are_sorted = false);
838
839 virtual void
841 const ::types::global_dof_index &row,
843 const bool indices_are_sorted = false) override;
844
846
857 Teuchos::RCP<GraphType>
859
867 domain_partitioner() const;
868
876 range_partitioner() const;
877
882 get_mpi_communicator() const;
883
889
904
912
924 begin() const;
925
930 end() const;
931
941 begin(const size_type r) const;
942
952 end(const size_type r) const;
953
967 void
968 print(std::ostream &out,
969 const bool write_extended_trilinos_info = false) const;
970
985 void
986 print_gnuplot(std::ostream &out) const;
987
997 int,
998 << "An error with error number " << arg1
999 << " occurred while calling a Trilinos function");
1000
1005 size_type,
1006 size_type,
1007 << "The entry with index <" << arg1 << ',' << arg2
1008 << "> does not exist.");
1009
1014 size_type,
1015 size_type,
1016 size_type,
1017 size_type,
1018 << "You tried to access element (" << arg1 << '/' << arg2
1019 << ')'
1020 << " of a distributed matrix, but only rows in range ["
1021 << arg3 << ',' << arg4
1022 << "] are stored locally and can be accessed.");
1023
1028 size_type,
1029 size_type,
1030 << "You tried to access element (" << arg1 << '/' << arg2
1031 << ')' << " of a sparse matrix, but it appears to not"
1032 << " exist in the Trilinos sparsity pattern.");
1034 private:
1039 Teuchos::RCP<MapType> column_space_map;
1040
1046 Teuchos::RCP<GraphType> graph;
1047
1054 Teuchos::RCP<GraphType> nonlocal_graph;
1055
1056 // TODO: currently only for double
1057 friend class SparseMatrix<double, MemorySpace>;
1060 };
1061
1062
1063
1064 // ---------------- inline and template functions -----------------
1065
1066
1067# ifndef DOXYGEN
1068
1069 namespace SparsityPatternIterators
1070 {
1071 template <typename MemorySpace>
1074 const size_type row,
1075 const size_type index)
1076 : sparsity_pattern(const_cast<SparsityPattern<MemorySpace> *>(sp))
1077 , a_row(row)
1078 , a_index(index)
1079 {
1080 visit_present_row();
1081 }
1082
1083
1084
1085 template <typename MemorySpace>
1086 inline typename Accessor<MemorySpace>::size_type
1087 Accessor<MemorySpace>::row() const
1088 {
1090 ExcBeyondEndOfSparsityPattern());
1091 return a_row;
1092 }
1093
1094
1095
1096 template <typename MemorySpace>
1097 inline typename Accessor<MemorySpace>::size_type
1098 Accessor<MemorySpace>::column() const
1099 {
1101 ExcBeyondEndOfSparsityPattern());
1102 return (*colnum_cache)[a_index];
1103 }
1104
1105
1106
1107 template <typename MemorySpace>
1108 inline typename Accessor<MemorySpace>::size_type
1109 Accessor<MemorySpace>::index() const
1110 {
1112 ExcBeyondEndOfSparsityPattern());
1113 return a_index;
1114 }
1115
1116
1117
1118 template <typename MemorySpace>
1119 inline Iterator<MemorySpace>::Iterator(
1121 const size_type row,
1122 const size_type index)
1123 : accessor(sp, row, index)
1124 {}
1125
1126
1127
1128 template <typename MemorySpace>
1129 inline Iterator<MemorySpace>::Iterator(const Iterator<MemorySpace> &) =
1130 default;
1131
1132
1133
1134 template <typename MemorySpace>
1135 inline Iterator<MemorySpace> &
1136 Iterator<MemorySpace>::operator++()
1137 {
1138 Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1140
1141 ++accessor.a_index;
1142
1143 // If at end of line: do one step, then cycle until we find a row with a
1144 // nonzero number of entries that is stored locally.
1145 if (accessor.a_index >=
1146 static_cast<::types::signed_global_dof_index>(
1147 accessor.colnum_cache->size()))
1148 {
1149 accessor.a_index = 0;
1150 ++accessor.a_row;
1151
1152 while (accessor.a_row <
1153 static_cast<::types::signed_global_dof_index>(
1154 accessor.sparsity_pattern->n_rows()))
1155 {
1156 const auto row_length =
1157 accessor.sparsity_pattern->row_length(accessor.a_row);
1158 if (row_length == 0 ||
1159 !accessor.sparsity_pattern->row_is_stored_locally(
1160 accessor.a_row))
1161 ++accessor.a_row;
1162 else
1163 break;
1164 }
1165
1166 accessor.visit_present_row();
1167 }
1168 return *this;
1169 }
1170
1171
1172
1173 template <typename MemorySpace>
1175 Iterator<MemorySpace>::operator++(int)
1176 {
1177 const Iterator<MemorySpace> old_state = *this;
1178 ++(*this);
1179 return old_state;
1180 }
1181
1182
1183
1184 template <typename MemorySpace>
1185 inline const Accessor<MemorySpace> &
1186 Iterator<MemorySpace>::operator*() const
1187 {
1188 return accessor;
1189 }
1190
1191
1192
1193 template <typename MemorySpace>
1194 inline const Accessor<MemorySpace> *
1195 Iterator<MemorySpace>::operator->() const
1196 {
1197 return &accessor;
1198 }
1199
1200
1201
1202 template <typename MemorySpace>
1203 inline bool
1204 Iterator<MemorySpace>::operator==(
1205 const Iterator<MemorySpace> &other) const
1206 {
1207 return (accessor.a_row == other.accessor.a_row &&
1208 accessor.a_index == other.accessor.a_index);
1209 }
1210
1211
1212
1213 template <typename MemorySpace>
1214 inline bool
1215 Iterator<MemorySpace>::operator!=(
1216 const Iterator<MemorySpace> &other) const
1217 {
1218 return !(*this == other);
1219 }
1220
1221
1222
1223 template <typename MemorySpace>
1224 inline bool
1225 Iterator<MemorySpace>::operator<(const Iterator<MemorySpace> &other) const
1226 {
1227 return (accessor.row() < other.accessor.row() ||
1228 (accessor.row() == other.accessor.row() &&
1229 accessor.index() < other.accessor.index()));
1230 }
1231
1232 } // namespace SparsityPatternIterators
1233
1234
1235
1236 template <typename MemorySpace>
1239 {
1240 const size_type first_valid_row = this->local_range().first;
1241 return const_iterator(this, first_valid_row, 0);
1242 }
1243
1244
1245
1246 template <typename MemorySpace>
1249 {
1250 return const_iterator(this, n_rows(), 0);
1251 }
1252
1253
1254
1255 template <typename MemorySpace>
1257 SparsityPattern<MemorySpace>::begin(const size_type r) const
1258 {
1259 AssertIndexRange(r, n_rows());
1260 if (row_length(r) > 0)
1261 return const_iterator(this, r, 0);
1262 else
1263 return end(r);
1264 }
1265
1266
1267
1268 template <typename MemorySpace>
1270 SparsityPattern<MemorySpace>::end(const size_type r) const
1271 {
1272 AssertIndexRange(r, n_rows());
1273
1274 // place the iterator on the first entry
1275 // past this line, or at the end of the
1276 // matrix
1277 for (size_type i = r + 1;
1279 ++i)
1280 if (row_length(i) > 0)
1281 return const_iterator(this, i, 0);
1282
1283 // if there is no such line, then take the
1284 // end iterator of the matrix
1285 return end();
1286 }
1287
1288
1289
1290 template <typename MemorySpace>
1291 inline bool
1292 SparsityPattern<MemorySpace>::in_local_range(const size_type index) const
1293 {
1295 graph->getRowMap()->getMinGlobalIndex();
1297 graph->getRowMap()->getMaxGlobalIndex() + 1;
1298
1299 return ((index >= static_cast<size_type>(begin)) &&
1300 (index < static_cast<size_type>(end)));
1301 }
1302
1303
1304
1305 template <typename MemorySpace>
1306 inline bool
1308 {
1309 return graph->isFillComplete();
1310 }
1311
1312
1313
1314 template <typename MemorySpace>
1315 inline bool
1317 {
1318 return ((n_rows() == 0) && (n_cols() == 0));
1319 }
1320
1321
1322
1323 template <typename MemorySpace>
1324 inline void
1325 SparsityPattern<MemorySpace>::add(const size_type i, const size_type j)
1326 {
1327 add_entries(i, &j, &j + 1);
1328 }
1329
1330
1331
1332 template <typename MemorySpace>
1333 template <typename ForwardIterator>
1334 inline void
1336 ForwardIterator begin,
1337 ForwardIterator end,
1338 const bool /*indices_are_sorted*/)
1339 {
1340 if (begin == end)
1341 return;
1342
1343 // verify that the size of the data type Trilinos expects matches that the
1344 // iterator points to. we allow for some slippage between signed and
1345 // unsigned and only compare that they are both either 32 or 64 bit. to
1346 // write this test properly, not that we cannot compare the size of
1347 // '*begin' because 'begin' may be an iterator and '*begin' may be an
1348 // accessor class. consequently, we need to somehow get an actual value
1349 // from it which we can by evaluating an expression such as when
1350 // multiplying the value produced by 2
1351 Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1353
1355 reinterpret_cast<TrilinosWrappers::types::int_type *>(
1356 const_cast<typename std::decay<decltype(*begin)>::type *>(&*begin));
1357
1359 reinterpret_cast<TrilinosWrappers::types::int_type *>(
1360 const_cast<typename std::decay<decltype(*end)>::type *>(&*end));
1361
1362 // Check at least for the first index that the conversion actually works
1366
1367 // TODO: The following line creates an array by copying the entries.
1368 // Perhaps there is a way to only create a 'view' of these arrays
1369 // and pass that to Tpetra?
1372
1373 if (row_is_stored_locally(row))
1374 graph->insertGlobalIndices(trilinos_row_index, array());
1375 else if (nonlocal_graph.get() != nullptr)
1376 {
1377 // this is the case when we have explicitly set the off-processor rows
1378 // and want to create a separate matrix object for them (to retain
1379 // thread-safety)
1380 Assert(nonlocal_graph->getRowMap()->getLocalElement(row) !=
1383 ExcMessage("Attempted to write into off-processor matrix row "
1384 "that has not be specified as being writable upon "
1385 "initialization"));
1386 nonlocal_graph->insertGlobalIndices(trilinos_row_index, array);
1387 }
1388 else
1389 graph->insertGlobalIndices(trilinos_row_index, array);
1390 }
1391
1392
1393
1394 template <typename MemorySpace>
1395 inline Teuchos::RCP<
1400 {
1401 return graph;
1402 }
1403
1404
1405
1406 template <typename MemorySpace>
1407 inline IndexSet
1409 {
1410 return IndexSet(graph->getDomainMap().getConst());
1411 }
1412
1413
1414
1415 template <typename MemorySpace>
1416 inline IndexSet
1418 {
1419 return IndexSet(graph->getRangeMap().getConst());
1420 }
1421
1422# endif // DOXYGEN
1423 } // namespace TpetraWrappers
1424
1425} // namespace LinearAlgebra
1426
1427
1429
1430
1431#endif // DEAL_II_TRILINOS_WITH_TPETRA
1432
1433#endif
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
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
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
Teuchos::RCP< const Teuchos::Comm< int > > get_teuchos_mpi_communicator() 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
Tpetra::CrsGraph< int, ::types::signed_global_dof_index, NodeType > GraphType
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)
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)
Teuchos::RCP< GraphType > trilinos_sparsity_pattern() const
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
const_iterator begin(const size_type r) const
Tpetra::Map< int, ::types::signed_global_dof_index, NodeType > MapType
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
friend class Tensor
Definition tensor.h:882
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:585
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)
Definition exceptions.h:539
#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)
Definition exceptions.h:562
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
types::global_dof_index size_type
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const Iterator & begin
Definition parallel.h:609
int signed_global_dof_index
Definition types.h:92