deal.II version GIT relicensing-2069-gd7d5dbb18d 2024-11-01 22:40: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\}}\)
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#include "deal.II/base/types.h"
21
23
24#ifdef DEAL_II_TRILINOS_WITH_TPETRA
25
29
33
34# include <Tpetra_CrsGraph.hpp>
35
36# include <cmath>
37# include <memory>
38# include <vector>
39
40
42
43// forward declarations
44# ifndef DOXYGEN
46
47namespace LinearAlgebra
48{
49 namespace TpetraWrappers
50 {
51 template <typename MemorySpace>
52 class SparsityPattern;
53
54 template <typename Number, typename MemorySpace>
55 class SparseMatrix;
56
58 {
59 template <typename MemorySpace>
60 class Iterator;
61 }
62 } // namespace TpetraWrappers
63} // namespace LinearAlgebra
64# endif
65
66namespace LinearAlgebra
67{
68 namespace TpetraWrappers
69 {
71 {
83 template <typename MemorySpace = ::MemorySpace::Host>
85 {
86 public:
91
96 const size_type row,
97 const size_type index);
98
103 row() const;
104
109 index() const;
110
115 column() const;
116
121
126 size_type,
127 size_type,
128 size_type,
129 << "You tried to access row " << arg1
130 << " of a distributed sparsity pattern, "
131 << " but only rows " << arg2 << " through " << arg3
132 << " are stored locally and can be accessed.");
133
134 private:
139
144
149
162 std::shared_ptr<std::vector<::types::signed_global_dof_index>>
164
170 void
172
173 // Make enclosing class a friend.
174 friend class Iterator<MemorySpace>;
175 };
176
182 template <typename MemorySpace = ::MemorySpace::Host>
184 {
185 public:
190
195 Iterator(const SparsityPattern<MemorySpace> *sparsity_pattern,
196 const size_type row,
197 const size_type index);
198
203
209
215
220 operator*() const;
221
226 operator->() const;
227
232 bool
234
238 bool
240
246 bool
247 operator<(const Iterator<MemorySpace> &) const;
248
253 size_type,
254 size_type,
255 << "Attempt to access element " << arg2 << " of row "
256 << arg1 << " which doesn't have that many elements.");
257
258 private:
263
265 };
266
267 } // namespace SparsityPatternIterators
268
269
287 template <typename MemorySpace = ::MemorySpace::Host>
289 {
290 public:
299
308
324 const size_type n,
325 const size_type n_entries_per_row = 0);
326
336 const size_type n,
337 const std::vector<size_type> &n_entries_per_row);
338
344
350 const SparsityPattern<MemorySpace> &input_sparsity_pattern);
351
355 virtual ~SparsityPattern() override = default;
356
368 void
369 reinit(const size_type m,
370 const size_type n,
371 const size_type n_entries_per_row = 0);
372
381 void
382 reinit(const size_type m,
383 const size_type n,
384 const std::vector<size_type> &n_entries_per_row);
385
390 void
391 copy_from(const SparsityPattern<MemorySpace> &input_sparsity_pattern);
392
398 template <typename SparsityPatternType>
399 void
400 copy_from(const SparsityPatternType &nontrilinos_sparsity_pattern);
401
409 operator=(const SparsityPattern<MemorySpace> &input_sparsity_pattern);
410
418 void
419 clear();
420
430 void
431 compress();
450 SparsityPattern(const IndexSet &parallel_partitioning,
451 const MPI_Comm communicator = MPI_COMM_WORLD,
452 const size_type n_entries_per_row = 0);
453
464 SparsityPattern(const IndexSet &parallel_partitioning,
465 const MPI_Comm communicator,
466 const std::vector<size_type> &n_entries_per_row);
467
482 SparsityPattern(const IndexSet &row_parallel_partitioning,
483 const IndexSet &col_parallel_partitioning,
484 const MPI_Comm communicator = MPI_COMM_WORLD,
485 const size_type n_entries_per_row = 0);
486
498 SparsityPattern(const IndexSet &row_parallel_partitioning,
499 const IndexSet &col_parallel_partitioning,
500 const MPI_Comm communicator,
501 const std::vector<size_type> &n_entries_per_row);
502
529 SparsityPattern(const IndexSet &row_parallel_partitioning,
530 const IndexSet &col_parallel_partitioning,
531 const IndexSet &writable_rows,
532 const MPI_Comm communicator = MPI_COMM_WORLD,
533 const size_type n_entries_per_row = 0);
534
550 void
551 reinit(const IndexSet &parallel_partitioning,
552 const MPI_Comm communicator = MPI_COMM_WORLD,
553 const size_type n_entries_per_row = 0);
554
565 void
566 reinit(const IndexSet &parallel_partitioning,
567 const MPI_Comm communicator,
568 const std::vector<size_type> &n_entries_per_row);
569
586 void
587 reinit(const IndexSet &row_parallel_partitioning,
588 const IndexSet &col_parallel_partitioning,
589 const MPI_Comm communicator = MPI_COMM_WORLD,
590 const size_type n_entries_per_row = 0);
591
617 void
618 reinit(const IndexSet &row_parallel_partitioning,
619 const IndexSet &col_parallel_partitioning,
620 const IndexSet &writeable_rows,
621 const MPI_Comm communicator = MPI_COMM_WORLD,
622 const size_type n_entries_per_row = 0);
623
628 void
629 reinit(const IndexSet &row_parallel_partitioning,
630 const IndexSet &col_parallel_partitioning,
631 const MPI_Comm communicator,
632 const std::vector<size_type> &n_entries_per_row);
633
643 template <typename SparsityPatternType>
644 void
645 reinit(const IndexSet &row_parallel_partitioning,
646 const IndexSet &col_parallel_partitioning,
647 const SparsityPatternType &nontrilinos_sparsity_pattern,
648 const MPI_Comm communicator = MPI_COMM_WORLD,
649 const bool exchange_data = false);
650
659 template <typename SparsityPatternType>
660 void
661 reinit(const IndexSet &parallel_partitioning,
662 const SparsityPatternType &nontrilinos_sparsity_pattern,
663 const MPI_Comm communicator = MPI_COMM_WORLD,
664 const bool exchange_data = false);
675 bool
677
681 unsigned int
682 max_entries_per_row() const;
683
693 unsigned int
694 local_size() const;
695
704 std::pair<size_type, size_type>
705 local_range() const;
706
711 bool
712 in_local_range(const size_type index) const;
713
717 std::uint64_t
718 n_nonzero_elements() const;
719
729 row_length(const size_type row) const;
730
738 bandwidth() const;
739
744 bool
745 empty() const;
746
751 bool
752 exists(const size_type i, const size_type j) const;
753
758 bool
759 row_is_stored_locally(const size_type i) const;
760
765 std::size_t
766 memory_consumption() const;
767
776 void
777 add(const size_type i, const size_type j);
778
779
783 template <typename ForwardIterator>
784 void
786 ForwardIterator begin,
787 ForwardIterator end,
788 const bool indices_are_sorted = false);
789
790 virtual void
792 const ::types::global_dof_index &row,
794 const bool indices_are_sorted = false) override;
795
797
808 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>>
810
817 Teuchos::RCP<const TpetraTypes::MapType<MemorySpace>>
818 domain_partitioner() const;
819
826 Teuchos::RCP<const TpetraTypes::MapType<MemorySpace>>
827 range_partitioner() const;
828
833 get_mpi_communicator() const;
834
838 Teuchos::RCP<const Teuchos::Comm<int>>
840
855
863
875 begin() const;
876
881 end() const;
882
892 begin(const size_type r) const;
893
903 end(const size_type r) const;
904
918 void
919 print(std::ostream &out,
920 const bool write_extended_trilinos_info = false) const;
921
936 void
937 print_gnuplot(std::ostream &out) const;
938
948 int,
949 << "An error with error number " << arg1
950 << " occurred while calling a Trilinos function");
951
956 size_type,
957 size_type,
958 << "The entry with index <" << arg1 << ',' << arg2
959 << "> does not exist.");
960
965 size_type,
966 size_type,
967 size_type,
968 size_type,
969 << "You tried to access element (" << arg1 << '/' << arg2
970 << ')'
971 << " of a distributed matrix, but only rows in range ["
972 << arg3 << ',' << arg4
973 << "] are stored locally and can be accessed.");
974
979 size_type,
980 size_type,
981 << "You tried to access element (" << arg1 << '/' << arg2
982 << ')' << " of a sparse matrix, but it appears to not"
983 << " exist in the Trilinos sparsity pattern.");
985 private:
990 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> column_space_map;
991
997 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> graph;
998
1005 Teuchos::RCP<TpetraTypes::GraphType<MemorySpace>> nonlocal_graph;
1006
1007 // TODO: currently only for double
1008 friend class SparseMatrix<double, MemorySpace>;
1011 };
1012
1013
1014
1015 // ---------------- inline and template functions -----------------
1016
1017
1018# ifndef DOXYGEN
1019
1020 namespace SparsityPatternIterators
1021 {
1022 template <typename MemorySpace>
1025 const size_type row,
1026 const size_type index)
1027 : sparsity_pattern(const_cast<SparsityPattern<MemorySpace> *>(sp))
1028 , a_row(row)
1029 , a_index(index)
1030 {
1031 visit_present_row();
1032 }
1033
1034
1035
1036 template <typename MemorySpace>
1037 inline typename Accessor<MemorySpace>::size_type
1038 Accessor<MemorySpace>::row() const
1039 {
1040 Assert(a_row < sparsity_pattern->n_rows(),
1041 ExcBeyondEndOfSparsityPattern());
1042 return a_row;
1043 }
1044
1045
1046
1047 template <typename MemorySpace>
1048 inline typename Accessor<MemorySpace>::size_type
1049 Accessor<MemorySpace>::column() const
1050 {
1051 Assert(a_row < sparsity_pattern->n_rows(),
1052 ExcBeyondEndOfSparsityPattern());
1053 return (*colnum_cache)[a_index];
1054 }
1055
1056
1057
1058 template <typename MemorySpace>
1059 inline typename Accessor<MemorySpace>::size_type
1060 Accessor<MemorySpace>::index() const
1061 {
1062 Assert(a_row < sparsity_pattern->n_rows(),
1063 ExcBeyondEndOfSparsityPattern());
1064 return a_index;
1065 }
1066
1067
1068
1069 template <typename MemorySpace>
1070 inline Iterator<MemorySpace>::Iterator(
1072 const size_type row,
1073 const size_type index)
1074 : accessor(sp, row, index)
1075 {}
1076
1077
1078
1079 template <typename MemorySpace>
1080 inline Iterator<MemorySpace>::Iterator(const Iterator<MemorySpace> &) =
1081 default;
1082
1083
1084
1085 template <typename MemorySpace>
1086 inline Iterator<MemorySpace> &
1087 Iterator<MemorySpace>::operator++()
1088 {
1089 Assert(accessor.a_row < accessor.sparsity_pattern->n_rows(),
1091
1092 ++accessor.a_index;
1093
1094 // If at end of line: do one step, then cycle until we find a row with a
1095 // nonzero number of entries that is stored locally.
1096 if (accessor.a_index >=
1097 static_cast<::types::signed_global_dof_index>(
1098 accessor.colnum_cache->size()))
1099 {
1100 accessor.a_index = 0;
1101 ++accessor.a_row;
1102
1103 while (accessor.a_row <
1104 static_cast<::types::signed_global_dof_index>(
1105 accessor.sparsity_pattern->n_rows()))
1106 {
1107 const auto row_length =
1108 accessor.sparsity_pattern->row_length(accessor.a_row);
1109 if (row_length == 0 ||
1110 !accessor.sparsity_pattern->row_is_stored_locally(
1111 accessor.a_row))
1112 ++accessor.a_row;
1113 else
1114 break;
1115 }
1116
1117 accessor.visit_present_row();
1118 }
1119 return *this;
1120 }
1121
1122
1123
1124 template <typename MemorySpace>
1125 inline Iterator<MemorySpace>
1126 Iterator<MemorySpace>::operator++(int)
1127 {
1128 const Iterator<MemorySpace> old_state = *this;
1129 ++(*this);
1130 return old_state;
1131 }
1132
1133
1134
1135 template <typename MemorySpace>
1136 inline const Accessor<MemorySpace> &
1137 Iterator<MemorySpace>::operator*() const
1138 {
1139 return accessor;
1140 }
1141
1142
1143
1144 template <typename MemorySpace>
1145 inline const Accessor<MemorySpace> *
1146 Iterator<MemorySpace>::operator->() const
1147 {
1148 return &accessor;
1149 }
1150
1151
1152
1153 template <typename MemorySpace>
1154 inline bool
1155 Iterator<MemorySpace>::operator==(
1156 const Iterator<MemorySpace> &other) const
1157 {
1158 return (accessor.a_row == other.accessor.a_row &&
1159 accessor.a_index == other.accessor.a_index);
1160 }
1161
1162
1163
1164 template <typename MemorySpace>
1165 inline bool
1166 Iterator<MemorySpace>::operator!=(
1167 const Iterator<MemorySpace> &other) const
1168 {
1169 return !(*this == other);
1170 }
1171
1172
1173
1174 template <typename MemorySpace>
1175 inline bool
1176 Iterator<MemorySpace>::operator<(const Iterator<MemorySpace> &other) const
1177 {
1178 return (accessor.row() < other.accessor.row() ||
1179 (accessor.row() == other.accessor.row() &&
1180 accessor.index() < other.accessor.index()));
1181 }
1182
1183 } // namespace SparsityPatternIterators
1184
1185
1186
1187 template <typename MemorySpace>
1190 {
1191 const size_type first_valid_row = this->local_range().first;
1192 return const_iterator(this, first_valid_row, 0);
1193 }
1194
1195
1196
1197 template <typename MemorySpace>
1200 {
1201 return const_iterator(this, n_rows(), 0);
1202 }
1203
1204
1205
1206 template <typename MemorySpace>
1208 SparsityPattern<MemorySpace>::begin(const size_type r) const
1209 {
1210 AssertIndexRange(r, n_rows());
1211 if (row_length(r) > 0)
1212 return const_iterator(this, r, 0);
1213 else
1214 return end(r);
1215 }
1216
1217
1218
1219 template <typename MemorySpace>
1221 SparsityPattern<MemorySpace>::end(const size_type r) const
1222 {
1223 AssertIndexRange(r, n_rows());
1224
1225 // place the iterator on the first entry
1226 // past this line, or at the end of the
1227 // matrix
1228 for (size_type i = r + 1; i < n_rows(); ++i)
1229 if (row_length(i) > 0)
1230 return const_iterator(this, i, 0);
1231
1232 // if there is no such line, then take the
1233 // end iterator of the matrix
1234 return end();
1235 }
1236
1237
1238
1239 template <typename MemorySpace>
1240 inline bool
1241 SparsityPattern<MemorySpace>::in_local_range(const size_type index) const
1242 {
1244 graph->getRowMap()->getMinGlobalIndex();
1246 graph->getRowMap()->getMaxGlobalIndex() + 1;
1247
1248 return ((index >= static_cast<size_type>(begin)) &&
1249 (index < static_cast<size_type>(end)));
1250 }
1251
1252
1253
1254 template <typename MemorySpace>
1255 inline bool
1257 {
1258 return graph->isFillComplete();
1259 }
1260
1261
1262
1263 template <typename MemorySpace>
1264 inline bool
1266 {
1267 return ((n_rows() == 0) && (n_cols() == 0));
1268 }
1269
1270
1271
1272 template <typename MemorySpace>
1273 inline void
1274 SparsityPattern<MemorySpace>::add(const size_type i, const size_type j)
1275 {
1276 add_entries(i, &j, &j + 1);
1277 }
1278
1279
1280
1281 template <typename MemorySpace>
1282 template <typename ForwardIterator>
1283 inline void
1285 ForwardIterator begin,
1286 ForwardIterator end,
1287 const bool /*indices_are_sorted*/)
1288 {
1289 if (begin == end)
1290 return;
1291
1292 // verify that the size of the data type Trilinos expects matches that the
1293 // iterator points to. we allow for some slippage between signed and
1294 // unsigned and only compare that they are both either 32 or 64 bit. to
1295 // write this test properly, not that we cannot compare the size of
1296 // '*begin' because 'begin' may be an iterator and '*begin' may be an
1297 // accessor class. consequently, we need to somehow get an actual value
1298 // from it which we can by evaluating an expression such as when
1299 // multiplying the value produced by 2
1300 Assert(sizeof(TrilinosWrappers::types::int_type) == sizeof((*begin) * 2),
1302
1303 const TrilinosWrappers::types::int_type *col_index_ptr_begin =
1304 reinterpret_cast<TrilinosWrappers::types::int_type *>(
1305 const_cast<std::decay_t<decltype(*begin)> *>(&*begin));
1306
1307 const TrilinosWrappers::types::int_type *col_index_ptr_end =
1308 reinterpret_cast<TrilinosWrappers::types::int_type *>(
1309 const_cast<std::decay_t<decltype(*end)> *>(&*end));
1310
1311 // Check at least for the first index that the conversion actually works
1312 AssertDimension(*col_index_ptr_begin, *begin);
1313 AssertDimension(*col_index_ptr_end, *end);
1314 TrilinosWrappers::types::int_type trilinos_row_index = row;
1315
1316 // TODO: The following line creates an array by copying the entries.
1317 // Perhaps there is a way to only create a 'view' of these arrays
1318 // and pass that to Tpetra?
1319 Teuchos::Array<TrilinosWrappers::types::int_type> array(
1320 col_index_ptr_begin, col_index_ptr_end);
1321
1322 if (row_is_stored_locally(row))
1323 graph->insertGlobalIndices(trilinos_row_index, array());
1324 else if (nonlocal_graph.get() != nullptr)
1325 {
1326 // this is the case when we have explicitly set the off-processor rows
1327 // and want to create a separate matrix object for them (to retain
1328 // thread-safety)
1329 Assert(nonlocal_graph->getRowMap()->getLocalElement(row) !=
1330 Teuchos::OrdinalTraits<
1332 ExcMessage("Attempted to write into off-processor matrix row "
1333 "that has not be specified as being writable upon "
1334 "initialization"));
1335 nonlocal_graph->insertGlobalIndices(trilinos_row_index, array);
1336 }
1337 else
1338 graph->insertGlobalIndices(trilinos_row_index, array);
1339 }
1340
1341
1342
1343 template <typename MemorySpace>
1344 inline Teuchos::RCP<Tpetra::CrsGraph<int,
1346 TpetraTypes::NodeType<MemorySpace>>>
1348 {
1349 return graph;
1350 }
1351
1352
1353
1354 template <typename MemorySpace>
1355 inline IndexSet
1357 {
1358 return IndexSet(graph->getDomainMap().getConst());
1359 }
1360
1361
1362
1363 template <typename MemorySpace>
1364 inline IndexSet
1366 {
1367 return IndexSet(graph->getRangeMap().getConst());
1368 }
1369
1370# endif // DOXYGEN
1371 } // namespace TpetraWrappers
1372
1373} // namespace LinearAlgebra
1374
1375
1377
1378
1379#endif // DEAL_II_TRILINOS_WITH_TPETRA
1380
1381#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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define DeclException0(Exception0)
Definition exceptions.h:466
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:580
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:534
#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:557
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:511
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
int signed_global_dof_index
Definition types.h:92
unsigned int global_dof_index
Definition types.h:81