Reference documentation for deal.II version 9.5.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
sparsity_pattern.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_sparsity_pattern_h
17#define dealii_sparsity_pattern_h
18
19
20#include <deal.II/base/config.h>
21
26
28
29// boost::serialization::make_array used to be in array.hpp, but was
30// moved to a different file in BOOST 1.64
31#include <boost/version.hpp>
32#if BOOST_VERSION >= 106400
33# include <boost/serialization/array_wrapper.hpp>
34#else
35# include <boost/serialization/array.hpp>
36#endif
37#include <boost/serialization/split_member.hpp>
38
39#include <algorithm>
40#include <iostream>
41#include <memory>
42#include <vector>
43
45
46// Forward declarations
47#ifndef DOXYGEN
48class SparsityPattern;
51template <typename number>
52class FullMatrix;
53template <typename number>
54class SparseMatrix;
55template <typename number>
57template <typename number>
58class SparseILU;
59
61{
62 class Accessor;
63}
64#endif
65
71namespace internals
72{
73 namespace SparsityPatternTools
74 {
79
87
93 template <typename value>
95 get_column_index_from_iterator(const std::pair<size_type, value> &i);
96
102 template <typename value>
104 get_column_index_from_iterator(const std::pair<const size_type, value> &i);
105
106 } // namespace SparsityPatternTools
107} // namespace internals
108
109
114{
115 // forward declaration
116 class Iterator;
117
122
134 {
135 public:
140
144 Accessor(const SparsityPattern *matrix, const std::size_t linear_index);
145
149 Accessor(const SparsityPattern *matrix);
150
157
163 row() const;
164
171 index() const;
172
184
190 column() const;
191
201 bool
203
207 bool
208 operator==(const Accessor &) const;
209
217 bool
218 operator<(const Accessor &) const;
219
220 protected:
222 "The instance of this class was initialized"
223 " without SparsityPattern object, which"
224 " means that it is a dummy accessor that can"
225 " not do any operations.");
226
231
239 std::size_t linear_index;
240
244 void
246
247 // Grant access to iterator class.
249
250 // Grant access to accessor class of ChunkSparsityPattern.
252 };
253
254
255
280 class Iterator : public LinearIndexIterator<Iterator, Accessor>
281 {
282 public:
287
292
298 Iterator(const SparsityPattern *sp, const std::size_t linear_index);
299
305 };
306} // namespace SparsityPatternIterators
307
344{
345public:
350
366
372
381
397 begin() const;
398
403 end() const;
404
417 begin(const size_type r) const;
418
428 end(const size_type r) const;
429
450
467
476 const size_type n,
477 const unsigned int max_per_row);
478
479
489 const size_type n,
490 const std::vector<unsigned int> &row_lengths);
491
500 SparsityPattern(const size_type m, const unsigned int max_per_row);
501
510 const std::vector<unsigned int> &row_lengths);
511
534 SparsityPattern(const SparsityPattern &original,
535 const unsigned int max_per_row,
536 const size_type extra_off_diagonals);
537
544 operator=(const SparsityPattern &);
545
556 void
557 reinit(const size_type m,
558 const size_type n,
559 const ArrayView<const unsigned int> &row_lengths);
560
564 void
565 reinit(const size_type m,
566 const size_type n,
567 const std::vector<unsigned int> &row_lengths);
568
577 void
578 reinit(const size_type m, const size_type n, const unsigned int max_per_row);
579
592 void
593 compress();
594
595
672 template <typename ForwardIterator>
673 void
675 const size_type n_cols,
676 const ForwardIterator begin,
677 const ForwardIterator end);
678
683 void
685
690 void
691 copy_from(const SparsityPattern &sp);
692
707 template <typename number>
708 void
709 copy_from(const FullMatrix<number> &matrix);
710
717 void
718 add(const size_type i, const size_type j);
719
726 template <typename ForwardIterator>
727 void
728 add_entries(const size_type row,
729 ForwardIterator begin,
730 ForwardIterator end,
731 const bool indices_are_sorted = false);
732
733 virtual void
734 add_row_entries(const size_type & row,
735 const ArrayView<const size_type> &columns,
736 const bool indices_are_sorted = false) override;
737
739
747 void
748 symmetrize();
749
770 std::size_t
772
777 bool
778 empty() const;
779
783 bool
784 exists(const size_type i, const size_type j) const;
785
797 std::pair<size_type, size_type>
798 matrix_position(const std::size_t global_index) const;
799
810 bandwidth() const;
811
815 bool
817
824 max_entries_per_row() const;
825
829 bool
831
844 bool
846
850 unsigned int
851 row_length(const size_type row) const;
852
857 std::size_t
858 memory_consumption() const;
859
893 operator()(const size_type i, const size_type j) const;
894
906 column_number(const size_type row, const unsigned int index) const;
907
916 row_position(const size_type i, const size_type j) const;
917
934 void
935 print(std::ostream &out) const;
936
950 void
951 print_gnuplot(std::ostream &out) const;
952
960 void
961 print_svg(std::ostream &out) const;
962
973 void
974 block_write(std::ostream &out) const;
975
989 void
990 block_read(std::istream &in);
991
996 template <class Archive>
997 void
998 save(Archive &ar, const unsigned int version) const;
999
1004 template <class Archive>
1005 void
1006 load(Archive &ar, const unsigned int version);
1007
1008#ifdef DOXYGEN
1013 template <class Archive>
1014 void
1015 serialize(Archive &archive, const unsigned int version);
1016#else
1017 // This macro defines the serialize() method that is compatible with
1018 // the templated save() and load() method that have been implemented.
1019 BOOST_SERIALIZATION_SPLIT_MEMBER()
1020#endif
1021
1035 int,
1036 int,
1037 << "The iterators denote a range of " << arg1
1038 << " elements, but the given number of rows was " << arg2);
1043 int,
1044 << "The number of partitions you gave is " << arg1
1045 << ", but must be greater than zero.");
1046
1051 int,
1052 int,
1053 << "Upon entering a new entry to row " << arg1
1054 << ": there was no free entry any more. " << std::endl
1055 << "(Maximum number of entries for this row: " << arg2
1056 << "; maybe the matrix is already compressed?)");
1057
1064 "The operation you attempted changes the structure of the SparsityPattern "
1065 "and is not possible after compress() has been called.");
1066
1073 "The operation you attempted is only allowed after the SparsityPattern "
1074 "has been set up and compress() was called.");
1075
1079private:
1084
1093
1099 std::size_t max_vec_len;
1100
1108 unsigned int max_row_length;
1109
1123 std::unique_ptr<std::size_t[]> rowstart;
1124
1147 std::unique_ptr<size_type[]> colnums;
1148
1153
1154 // Make all sparse matrices friends of this class.
1155 template <typename number>
1156 friend class SparseMatrix;
1157 template <typename number>
1159 template <typename number>
1160 friend class SparseILU;
1161 template <typename number>
1162 friend class ChunkSparseMatrix;
1163
1166
1167 // Also give access to internal details to the iterator/accessor classes.
1171};
1172
1173
1177/*---------------------- Inline functions -----------------------------------*/
1178
1179#ifndef DOXYGEN
1180
1181
1183{
1184 inline Accessor::Accessor(const SparsityPattern *sparsity_pattern,
1185 const std::size_t i)
1186 : container(sparsity_pattern)
1187 , linear_index(i)
1188 {}
1189
1190
1191
1192 inline Accessor::Accessor(const SparsityPattern *sparsity_pattern)
1193 : container(sparsity_pattern)
1194 , linear_index(container->rowstart[container->rows])
1195 {}
1196
1197
1198
1199 inline Accessor::Accessor()
1200 : container(nullptr)
1201 , linear_index(numbers::invalid_size_type)
1202 {}
1203
1204
1205
1206 inline bool
1207 Accessor::is_valid_entry() const
1208 {
1209 Assert(container != nullptr, DummyAccessor());
1210 return (linear_index < container->rowstart[container->rows] &&
1211 container->colnums[linear_index] != SparsityPattern::invalid_entry);
1212 }
1213
1214
1215
1216 inline size_type
1217 Accessor::row() const
1218 {
1219 Assert(is_valid_entry() == true, ExcInvalidIterator());
1220
1221 const std::size_t *insert_point =
1222 std::upper_bound(container->rowstart.get(),
1223 container->rowstart.get() + container->rows + 1,
1224 linear_index);
1225 return insert_point - container->rowstart.get() - 1;
1226 }
1227
1228
1229
1230 inline size_type
1231 Accessor::column() const
1232 {
1233 Assert(is_valid_entry() == true, ExcInvalidIterator());
1234
1235 return (container->colnums[linear_index]);
1236 }
1237
1238
1239
1240 inline size_type
1241 Accessor::index() const
1242 {
1243 Assert(is_valid_entry() == true, ExcInvalidIterator());
1244
1245 return linear_index - container->rowstart[row()];
1246 }
1247
1248
1249
1250 inline size_type
1251 Accessor::global_index() const
1252 {
1253 Assert(is_valid_entry() == true, ExcInvalidIterator());
1254
1255 return linear_index;
1256 }
1257
1258
1259
1260 inline bool
1261 Accessor::operator==(const Accessor &other) const
1262 {
1263 return (container == other.container && linear_index == other.linear_index);
1264 }
1265
1266
1267
1268 inline bool
1269 Accessor::operator<(const Accessor &other) const
1270 {
1271 Assert(container != nullptr, DummyAccessor());
1272 Assert(other.container != nullptr, DummyAccessor());
1273 Assert(container == other.container, ExcInternalError());
1274
1275 return linear_index < other.linear_index;
1276 }
1277
1278
1279
1280 inline void
1281 Accessor::advance()
1282 {
1283 Assert(container != nullptr, DummyAccessor());
1284 Assert(linear_index < container->rowstart[container->rows],
1286 ++linear_index;
1287 }
1288
1289
1290 inline Iterator::Iterator(const SparsityPattern *sp,
1291 const std::size_t linear_index)
1292 : LinearIndexIterator<Iterator, Accessor>(Accessor(sp, linear_index))
1293 {}
1294
1295
1296 inline Iterator::Iterator(const Accessor &accessor)
1297 : LinearIndexIterator<Iterator, Accessor>(accessor)
1298 {}
1299
1300
1301} // namespace SparsityPatternIterators
1302
1303
1304
1305inline std::size_t
1307{
1309
1310 if ((rowstart != nullptr) && (colnums != nullptr))
1311 return rowstart[rows] - rowstart[0];
1312 else
1313 // the object is empty or has zero size
1314 return 0;
1315}
1316
1317
1318
1319inline bool
1321{
1322 return compressed;
1323}
1324
1325
1326
1327inline bool
1329{
1330 return (store_diagonal_first_in_row == false);
1331}
1332
1333
1334
1335inline unsigned int
1336SparsityPattern::row_length(const size_type row) const
1337{
1338 AssertIndexRange(row, rows);
1339 return rowstart[row + 1] - rowstart[row];
1340}
1341
1342
1343
1345SparsityPattern::column_number(const size_type row,
1346 const unsigned int index) const
1347{
1348 AssertIndexRange(row, rows);
1349 AssertIndexRange(index, row_length(row));
1350
1351 return colnums[rowstart[row] + index];
1352}
1353
1354
1355
1358{
1359 if (n_rows() > 0)
1360 return {this, rowstart[0]};
1361 else
1362 return end();
1363}
1364
1365
1366
1369{
1370 if (n_rows() > 0)
1371 return {this, rowstart[rows]};
1372 else
1373 return {nullptr, 0};
1374}
1375
1376
1377
1379SparsityPattern::begin(const size_type r) const
1380{
1382
1383 return {this, rowstart[r]};
1384}
1385
1386
1387
1389SparsityPattern::end(const size_type r) const
1390{
1392
1393 return {this, rowstart[r + 1]};
1394}
1395
1396
1397
1398inline bool
1400{
1402 return false;
1403
1404 // it isn't quite necessary to compare *all* member variables. by only
1405 // comparing the essential ones, we can say that two sparsity patterns are
1406 // equal even if one is compressed and the other is not (in which case some
1407 // of the member variables are not yet set correctly)
1408 if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed)
1409 return false;
1410
1411 if (rows > 0)
1412 for (size_type i = 0; i < rows + 1; ++i)
1413 if (rowstart[i] != sp2.rowstart[i])
1414 return false;
1415
1416 if (rows > 0)
1417 for (size_type i = 0; i < rowstart[rows]; ++i)
1418 if (colnums[i] != sp2.colnums[i])
1419 return false;
1420
1421 return true;
1422}
1423
1424
1425
1426namespace internal
1427{
1428 namespace SparsityPatternTools
1429 {
1434
1435 inline size_type
1436 get_column_index_from_iterator(const size_type i)
1437 {
1438 return i;
1439 }
1440
1441
1442
1443 template <typename value>
1444 inline size_type
1445 get_column_index_from_iterator(const std::pair<size_type, value> &i)
1446 {
1447 return i.first;
1448 }
1449
1450
1451
1452 template <typename value>
1453 inline size_type
1454 get_column_index_from_iterator(const std::pair<const size_type, value> &i)
1455 {
1456 return i.first;
1457 }
1458 } // namespace SparsityPatternTools
1459} // namespace internal
1460
1461
1462
1463template <typename ForwardIterator>
1464void
1465SparsityPattern::copy_from(const size_type n_rows,
1466 const size_type n_cols,
1467 const ForwardIterator begin,
1468 const ForwardIterator end)
1469{
1470 Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1471 ExcIteratorRange(std::distance(begin, end), n_rows));
1472
1473 // first determine row lengths for each row. if the matrix is quadratic,
1474 // then we might have to add an additional entry for the diagonal, if that
1475 // is not yet present. as we have to call compress anyway later on, don't
1476 // bother to check whether that diagonal entry is in a certain row or not
1477 const bool is_square = (n_rows == n_cols);
1478 std::vector<unsigned int> row_lengths;
1479 row_lengths.reserve(n_rows);
1480 for (ForwardIterator i = begin; i != end; ++i)
1481 row_lengths.push_back(std::distance(i->begin(), i->end()) +
1482 (is_square ? 1 : 0));
1483 reinit(n_rows, n_cols, row_lengths);
1484
1485 // now enter all the elements into the matrix. note that if the matrix is
1486 // quadratic, then we already have the diagonal element preallocated
1487 //
1488 // for use in the inner loop, we define an alias to the type of the inner
1489 // iterators
1490 size_type row = 0;
1491 using inner_iterator =
1492 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1493 for (ForwardIterator i = begin; i != end; ++i, ++row)
1494 {
1495 size_type * cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
1496 const inner_iterator end_of_row = i->end();
1497 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1498 {
1499 const size_type col =
1500 internal::SparsityPatternTools::get_column_index_from_iterator(*j);
1502
1503 if ((col != row) || !is_square)
1504 *cols++ = col;
1505 }
1506 }
1507
1508 // finally compress everything. this also sorts the entries within each row
1509 compress();
1510}
1511
1512
1513
1514template <class Archive>
1515inline void
1516SparsityPattern::save(Archive &ar, const unsigned int) const
1517{
1518 // forward to serialization function in the base class.
1519 ar &boost::serialization::base_object<const Subscriptor>(*this);
1520
1522
1523 if (max_dim != 0)
1524 ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1525 else
1526 Assert(rowstart.get() == nullptr, ExcInternalError());
1527
1528 if (max_vec_len != 0)
1529 ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1530 else
1531 Assert(colnums.get() == nullptr, ExcInternalError());
1533}
1534
1535
1536
1537template <class Archive>
1538inline void
1539SparsityPattern::load(Archive &ar, const unsigned int)
1540{
1541 // forward to serialization function in the base class.
1542 ar &boost::serialization::base_object<Subscriptor>(*this);
1543
1545
1546 if (max_dim != 0)
1547 {
1548 rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
1549 ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1550 }
1551 else
1552 rowstart.reset();
1553
1554 if (max_vec_len != 0)
1555 {
1556 colnums = std::make_unique<size_type[]>(max_vec_len);
1557 ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1558 }
1559 else
1560 colnums.reset();
1562}
1563
1564
1565#endif // DOXYGEN
1566
1568
1569#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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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:533
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ExcMatrixIsCompressed()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
types::global_dof_index size_type
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:222
unsigned int global_dof_index
Definition types.h:82