Reference documentation for deal.II version GIT relicensing-437-g81ec864850 2024-04-19 07: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
sparsity_pattern.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2023 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_sparsity_pattern_h
16#define dealii_sparsity_pattern_h
17
18
19#include <deal.II/base/config.h>
20
25
27
28// boost::serialization::make_array used to be in array.hpp, but was
29// moved to a different file in BOOST 1.64
30#include <boost/version.hpp>
31#if BOOST_VERSION >= 106400
32# include <boost/serialization/array_wrapper.hpp>
33#else
34# include <boost/serialization/array.hpp>
35#endif
36#include <boost/serialization/split_member.hpp>
37
38#include <algorithm>
39#include <iostream>
40#include <memory>
41#include <vector>
42
44
45// Forward declarations
46#ifndef DOXYGEN
47class SparsityPattern;
50template <typename number>
51class FullMatrix;
52template <typename number>
53class SparseMatrix;
54template <typename number>
56template <typename number>
57class SparseILU;
58
60{
61 class Accessor;
62}
63#endif
64
70namespace internals
71{
72 namespace SparsityPatternTools
73 {
78
86
92 template <typename value>
94 get_column_index_from_iterator(const std::pair<size_type, value> &i);
95
101 template <typename value>
103 get_column_index_from_iterator(const std::pair<const size_type, value> &i);
104
105 } // namespace SparsityPatternTools
106} // namespace internals
107
108
113{
114 // forward declaration
115 class Iterator;
116
121
133 {
134 public:
139
143 Accessor(const SparsityPattern *matrix, const std::size_t linear_index);
144
148 Accessor(const SparsityPattern *matrix);
149
156
162 row() const;
163
170 index() const;
171
183
189 column() const;
190
200 bool
202
206 bool
207 operator==(const Accessor &) const;
208
216 bool
217 operator<(const Accessor &) const;
218
219 protected:
221 "The instance of this class was initialized"
222 " without SparsityPattern object, which"
223 " means that it is a dummy accessor that can"
224 " not do any operations.");
225
230
238 std::size_t linear_index;
239
243 void
245
246 // Grant access to iterator class.
248
249 // Grant access to accessor class of ChunkSparsityPattern.
251 };
252
253
254
279 class Iterator : public LinearIndexIterator<Iterator, Accessor>
280 {
281 public:
286
291
297 Iterator(const SparsityPattern *sp, const std::size_t linear_index);
298
304 };
305} // namespace SparsityPatternIterators
306
343{
344public:
349
365
371
380
396 begin() const;
397
402 end() const;
403
416 begin(const size_type r) const;
417
427 end(const size_type r) const;
428
449
466
475 const size_type n,
476 const unsigned int max_per_row);
477
478
488 const size_type n,
489 const std::vector<unsigned int> &row_lengths);
490
499 SparsityPattern(const size_type m, const unsigned int max_per_row);
500
509 const std::vector<unsigned int> &row_lengths);
510
533 SparsityPattern(const SparsityPattern &original,
534 const unsigned int max_per_row,
535 const size_type extra_off_diagonals);
536
543 operator=(const SparsityPattern &);
544
555 void
556 reinit(const size_type m,
557 const size_type n,
558 const ArrayView<const unsigned int> &row_lengths);
559
563 void
564 reinit(const size_type m,
565 const size_type n,
566 const std::vector<unsigned int> &row_lengths);
567
576 void
577 reinit(const size_type m, const size_type n, const unsigned int max_per_row);
578
591 void
592 compress();
593
594
671 template <typename ForwardIterator>
672 void
674 const size_type n_cols,
675 const ForwardIterator begin,
676 const ForwardIterator end);
677
682 void
684
689 void
690 copy_from(const SparsityPattern &sp);
691
706 template <typename number>
707 void
708 copy_from(const FullMatrix<number> &matrix);
709
716 void
717 add(const size_type i, const size_type j);
718
725 template <typename ForwardIterator>
726 void
727 add_entries(const size_type row,
728 ForwardIterator begin,
729 ForwardIterator end,
730 const bool indices_are_sorted = false);
731
732 virtual void
733 add_row_entries(const size_type &row,
734 const ArrayView<const size_type> &columns,
735 const bool indices_are_sorted = false) override;
736
738
746 void
747 symmetrize();
748
769 std::size_t
771
776 bool
777 empty() const;
778
782 bool
783 exists(const size_type i, const size_type j) const;
784
796 std::pair<size_type, size_type>
797 matrix_position(const std::size_t global_index) const;
798
809 bandwidth() const;
810
814 bool
816
823 max_entries_per_row() const;
824
828 bool
830
843 bool
845
849 unsigned int
850 row_length(const size_type row) const;
851
856 std::size_t
857 memory_consumption() const;
858
892 operator()(const size_type i, const size_type j) const;
893
905 column_number(const size_type row, const unsigned int index) const;
906
915 row_position(const size_type i, const size_type j) const;
916
933 void
934 print(std::ostream &out) const;
935
949 void
950 print_gnuplot(std::ostream &out) const;
951
959 void
960 print_svg(std::ostream &out) const;
961
972 void
973 block_write(std::ostream &out) const;
974
988 void
989 block_read(std::istream &in);
990
995 template <class Archive>
996 void
997 save(Archive &ar, const unsigned int version) const;
998
1003 template <class Archive>
1004 void
1005 load(Archive &ar, const unsigned int version);
1006
1007#ifdef DOXYGEN
1012 template <class Archive>
1013 void
1014 serialize(Archive &archive, const unsigned int version);
1015#else
1016 // This macro defines the serialize() method that is compatible with
1017 // the templated save() and load() method that have been implemented.
1018 BOOST_SERIALIZATION_SPLIT_MEMBER()
1019#endif
1020
1034 int,
1035 int,
1036 << "The iterators denote a range of " << arg1
1037 << " elements, but the given number of rows was " << arg2);
1042 int,
1043 << "The number of partitions you gave is " << arg1
1044 << ", but must be greater than zero.");
1045
1050 int,
1051 int,
1052 << "Upon entering a new entry to row " << arg1
1053 << ": there was no free entry any more. " << std::endl
1054 << "(Maximum number of entries for this row: " << arg2
1055 << "; maybe the matrix is already compressed?)");
1056
1063 "The operation you attempted changes the structure of the SparsityPattern "
1064 "and is not possible after compress() has been called.");
1065
1072 "The operation you attempted is only allowed after the SparsityPattern "
1073 "has been set up and compress() was called.");
1074
1078private:
1083
1092
1098 std::size_t max_vec_len;
1099
1107 unsigned int max_row_length;
1108
1122 std::unique_ptr<std::size_t[]> rowstart;
1123
1146 std::unique_ptr<size_type[]> colnums;
1147
1152
1153 // Make all sparse matrices friends of this class.
1154 template <typename number>
1155 friend class SparseMatrix;
1156 template <typename number>
1158 template <typename number>
1159 friend class SparseILU;
1160 template <typename number>
1161 friend class ChunkSparseMatrix;
1162
1165
1166 // Also give access to internal details to the iterator/accessor classes.
1170};
1171
1172
1176/*---------------------- Inline functions -----------------------------------*/
1177
1178#ifndef DOXYGEN
1179
1180
1182{
1183 inline Accessor::Accessor(const SparsityPattern *sparsity_pattern,
1184 const std::size_t i)
1185 : container(sparsity_pattern)
1186 , linear_index(i)
1187 {}
1188
1189
1190
1191 inline Accessor::Accessor(const SparsityPattern *sparsity_pattern)
1192 : container(sparsity_pattern)
1193 , linear_index(container->rowstart[container->rows])
1194 {}
1195
1196
1197
1198 inline Accessor::Accessor()
1199 : container(nullptr)
1200 , linear_index(numbers::invalid_size_type)
1201 {}
1202
1203
1204
1205 inline bool
1206 Accessor::is_valid_entry() const
1207 {
1208 Assert(container != nullptr, DummyAccessor());
1209 return (linear_index < container->rowstart[container->rows] &&
1210 container->colnums[linear_index] != SparsityPattern::invalid_entry);
1211 }
1212
1213
1214
1215 inline size_type
1216 Accessor::row() const
1217 {
1218 Assert(is_valid_entry() == true, ExcInvalidIterator());
1219
1220 const std::size_t *insert_point =
1221 std::upper_bound(container->rowstart.get(),
1222 container->rowstart.get() + container->rows + 1,
1223 linear_index);
1224 return insert_point - container->rowstart.get() - 1;
1225 }
1226
1227
1228
1229 inline size_type
1230 Accessor::column() const
1231 {
1232 Assert(is_valid_entry() == true, ExcInvalidIterator());
1233
1234 return (container->colnums[linear_index]);
1235 }
1236
1237
1238
1239 inline size_type
1240 Accessor::index() const
1241 {
1242 Assert(is_valid_entry() == true, ExcInvalidIterator());
1243
1244 return linear_index - container->rowstart[row()];
1245 }
1246
1247
1248
1249 inline size_type
1250 Accessor::global_index() const
1251 {
1252 Assert(is_valid_entry() == true, ExcInvalidIterator());
1253
1254 return linear_index;
1255 }
1256
1257
1258
1259 inline bool
1260 Accessor::operator==(const Accessor &other) const
1261 {
1262 return (container == other.container && linear_index == other.linear_index);
1263 }
1264
1265
1266
1267 inline bool
1268 Accessor::operator<(const Accessor &other) const
1269 {
1270 Assert(container != nullptr, DummyAccessor());
1271 Assert(other.container != nullptr, DummyAccessor());
1272 Assert(container == other.container, ExcInternalError());
1273
1274 return linear_index < other.linear_index;
1275 }
1276
1277
1278
1279 inline void
1280 Accessor::advance()
1281 {
1282 Assert(container != nullptr, DummyAccessor());
1283 Assert(linear_index < container->rowstart[container->rows],
1285 ++linear_index;
1286 }
1287
1288
1289 inline Iterator::Iterator(const SparsityPattern *sp,
1290 const std::size_t linear_index)
1291 : LinearIndexIterator<Iterator, Accessor>(Accessor(sp, linear_index))
1292 {}
1293
1294
1295 inline Iterator::Iterator(const Accessor &accessor)
1296 : LinearIndexIterator<Iterator, Accessor>(accessor)
1297 {}
1298
1299
1300} // namespace SparsityPatternIterators
1301
1302
1303
1304inline std::size_t
1306{
1308
1309 if ((rowstart != nullptr) && (colnums != nullptr))
1310 return rowstart[rows] - rowstart[0];
1311 else
1312 // the object is empty or has zero size
1313 return 0;
1314}
1315
1316
1317
1318inline bool
1320{
1321 return compressed;
1322}
1323
1324
1325
1326inline bool
1328{
1329 return (store_diagonal_first_in_row == false);
1330}
1331
1332
1333
1334inline unsigned int
1335SparsityPattern::row_length(const size_type row) const
1336{
1337 AssertIndexRange(row, rows);
1338 return rowstart[row + 1] - rowstart[row];
1339}
1340
1341
1342
1344SparsityPattern::column_number(const size_type row,
1345 const unsigned int index) const
1346{
1347 AssertIndexRange(row, rows);
1348 AssertIndexRange(index, row_length(row));
1349
1350 return colnums[rowstart[row] + index];
1351}
1352
1353
1354
1357{
1358 if (n_rows() > 0)
1359 return {this, rowstart[0]};
1360 else
1361 return end();
1362}
1363
1364
1365
1368{
1369 if (n_rows() > 0)
1370 return {this, rowstart[rows]};
1371 else
1372 return {nullptr, 0};
1373}
1374
1375
1376
1378SparsityPattern::begin(const size_type r) const
1379{
1381
1382 return {this, rowstart[r]};
1383}
1384
1385
1386
1388SparsityPattern::end(const size_type r) const
1389{
1391
1392 return {this, rowstart[r + 1]};
1393}
1394
1395
1396
1397inline bool
1399{
1401 return false;
1402
1403 // it isn't quite necessary to compare *all* member variables. by only
1404 // comparing the essential ones, we can say that two sparsity patterns are
1405 // equal even if one is compressed and the other is not (in which case some
1406 // of the member variables are not yet set correctly)
1407 if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed)
1408 return false;
1409
1410 if (rows > 0)
1411 for (size_type i = 0; i < rows + 1; ++i)
1412 if (rowstart[i] != sp2.rowstart[i])
1413 return false;
1414
1415 if (rows > 0)
1416 for (size_type i = 0; i < rowstart[rows]; ++i)
1417 if (colnums[i] != sp2.colnums[i])
1418 return false;
1419
1420 return true;
1421}
1422
1423
1424
1425namespace internal
1426{
1427 namespace SparsityPatternTools
1428 {
1433
1434 inline size_type
1435 get_column_index_from_iterator(const size_type i)
1436 {
1437 return i;
1438 }
1439
1440
1441
1442 template <typename value>
1443 inline size_type
1444 get_column_index_from_iterator(const std::pair<size_type, value> &i)
1445 {
1446 return i.first;
1447 }
1448
1449
1450
1451 template <typename value>
1452 inline size_type
1453 get_column_index_from_iterator(const std::pair<const size_type, value> &i)
1454 {
1455 return i.first;
1456 }
1457 } // namespace SparsityPatternTools
1458} // namespace internal
1459
1460
1461
1462template <typename ForwardIterator>
1463void
1464SparsityPattern::copy_from(const size_type n_rows,
1465 const size_type n_cols,
1466 const ForwardIterator begin,
1467 const ForwardIterator end)
1468{
1469 Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1470 ExcIteratorRange(std::distance(begin, end), n_rows));
1471
1472 // first determine row lengths for each row. if the matrix is quadratic,
1473 // then we might have to add an additional entry for the diagonal, if that
1474 // is not yet present. as we have to call compress anyway later on, don't
1475 // bother to check whether that diagonal entry is in a certain row or not
1476 const bool is_square = (n_rows == n_cols);
1477 std::vector<unsigned int> row_lengths;
1478 row_lengths.reserve(n_rows);
1479 for (ForwardIterator i = begin; i != end; ++i)
1480 row_lengths.push_back(std::distance(i->begin(), i->end()) +
1481 (is_square ? 1 : 0));
1482 reinit(n_rows, n_cols, row_lengths);
1483
1484 // now enter all the elements into the matrix. note that if the matrix is
1485 // quadratic, then we already have the diagonal element preallocated
1486 //
1487 // for use in the inner loop, we define an alias to the type of the inner
1488 // iterators
1489 size_type row = 0;
1490 using inner_iterator =
1491 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1492 for (ForwardIterator i = begin; i != end; ++i, ++row)
1493 {
1494 size_type *cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
1495 const inner_iterator end_of_row = i->end();
1496 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1497 {
1498 const size_type col =
1499 internal::SparsityPatternTools::get_column_index_from_iterator(*j);
1501
1502 if ((col != row) || !is_square)
1503 *cols++ = col;
1504 }
1505 }
1506
1507 // finally compress everything. this also sorts the entries within each row
1508 compress();
1509}
1510
1511
1512
1513template <class Archive>
1514inline void
1515SparsityPattern::save(Archive &ar, const unsigned int) const
1516{
1517 // forward to serialization function in the base class.
1518 ar &boost::serialization::base_object<const Subscriptor>(*this);
1519
1521
1522 if (max_dim != 0)
1523 ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1524 else
1525 Assert(rowstart.get() == nullptr, ExcInternalError());
1526
1527 if (max_vec_len != 0)
1528 ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1529 else
1530 Assert(colnums.get() == nullptr, ExcInternalError());
1532}
1533
1534
1535
1536template <class Archive>
1537inline void
1538SparsityPattern::load(Archive &ar, const unsigned int)
1539{
1540 // forward to serialization function in the base class.
1541 ar &boost::serialization::base_object<Subscriptor>(*this);
1542
1544
1545 if (max_dim != 0)
1546 {
1547 rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
1548 ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1549 }
1550 else
1551 rowstart.reset();
1552
1553 if (max_vec_len != 0)
1554 {
1555 colnums = std::make_unique<size_type[]>(max_vec_len);
1556 ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1557 }
1558 else
1559 colnums.reset();
1561}
1562
1563
1564#endif // DOXYGEN
1565
1567
1568#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
friend class ChunkSparsityPatternIterators::Accessor
bool operator==(const Accessor &) const
Accessor(const SparsityPattern *matrix)
Accessor(const SparsityPattern *matrix, const std::size_t linear_index)
bool operator<(const Accessor &) const
Iterator(const SparsityPattern *sp, const std::size_t linear_index)
Iterator(const Accessor &accessor)
void load(Archive &ar, const unsigned int version)
std::pair< size_type, size_type > matrix_position(const std::size_t global_index) const
iterator begin(const size_type r) const
void block_write(std::ostream &out) const
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
void print_svg(std::ostream &out) const
bool stores_only_added_elements() const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
void serialize(Archive &archive, const unsigned int version)
size_type bandwidth() const
void print_gnuplot(std::ostream &out) const
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
bool is_compressed() const
SparsityPattern & operator=(const SparsityPattern &)
std::size_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
iterator begin() const
std::unique_ptr< size_type[]> colnums
std::size_t max_vec_len
size_type column_number(const size_type row, const unsigned int index) const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
static constexpr size_type invalid_entry
size_type row_position(const size_type i, const size_type j) const
std::unique_ptr< std::size_t[]> rowstart
void print(std::ostream &out) const
iterator end(const size_type r) const
void add(const size_type i, const size_type j)
size_type max_entries_per_row() const
unsigned int max_row_length
iterator end() const
size_type operator()(const size_type i, const size_type j) const
void save(Archive &ar, const unsigned int version) const
types::global_dof_index size_type
unsigned int row_length(const size_type row) const
std::size_t memory_consumption() const
bool operator==(const SparsityPattern &) const
void block_read(std::istream &in)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcInvalidIterator()
static ::ExceptionBase & DummyAccessor()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotEnoughSpace(int arg1, int arg2)
static ::ExceptionBase & ExcIteratorPastEnd()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcMatrixIsCompressed()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
types::global_dof_index size_type
types::global_dof_index size_type
size_type get_column_index_from_iterator(const size_type i)
const types::global_dof_index invalid_size_type
Definition types.h:233
unsigned int global_dof_index
Definition types.h:81