Reference documentation for deal.II version 9.4.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\}}\)
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
27// boost::serialization::make_array used to be in array.hpp, but was
28// moved to a different file in BOOST 1.64
29#include <boost/version.hpp>
30#if BOOST_VERSION >= 106400
31# include <boost/serialization/array_wrapper.hpp>
32#else
33# include <boost/serialization/array.hpp>
34#endif
35#include <boost/serialization/split_member.hpp>
36
37#include <algorithm>
38#include <iostream>
39#include <memory>
40#include <vector>
41
43
44// Forward declarations
45#ifndef DOXYGEN
46class 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
69namespace internals
70{
71 namespace SparsityPatternTools
72 {
77
85
91 template <typename value>
93 get_column_index_from_iterator(const std::pair<size_type, value> &i);
94
100 template <typename value>
102 get_column_index_from_iterator(const std::pair<const size_type, value> &i);
103
104 } // namespace SparsityPatternTools
105} // namespace internals
106
107
112{
113 // forward declaration
114 class Iterator;
115
120
132 {
133 public:
138
143
148
155
161 row() const;
162
169 index() const;
170
182
188 column() const;
189
199 bool
201
205 bool
206 operator==(const Accessor &) const;
207
215 bool
216 operator<(const Accessor &) const;
217
218 protected:
220 "The instance of this class was initialized"
221 " without SparsityPattern object, which"
222 " means that it is a dummy accessor that can"
223 " not do any operations.");
224
229
237 std::size_t linear_index;
238
242 void
244
245 // Grant access to iterator class.
247
248 // Grant access to accessor class of ChunkSparsityPattern.
250 };
251
252
253
278 class Iterator : public LinearIndexIterator<Iterator, Accessor>
279 {
280 public:
285
290
296 Iterator(const SparsityPatternBase *sp, const std::size_t linear_index);
297
303 };
304} // namespace SparsityPatternIterators
305
306
307
326{
327public:
332
338
347
363
379
383 ~SparsityPatternBase() override = default;
384
393 void
394 reinit(const size_type m, const size_type n, const unsigned int max_per_row);
395
406 void
407 reinit(const size_type m,
408 const size_type n,
409 const std::vector<unsigned int> &row_lengths);
410
416 virtual void
418 const size_type n,
419 const ArrayView<const unsigned int> &row_lengths) = 0;
420
428 void
429 symmetrize();
430
437 void
438 add(const size_type i, const size_type j);
439
459 begin() const;
460
465 end() const;
466
479 begin(const size_type r) const;
480
490 end(const size_type r) const;
491
492
506 bool
508
513 bool
514 empty() const;
515
519 bool
520 exists(const size_type i, const size_type j) const;
521
528 max_entries_per_row() const;
529
540 bandwidth() const;
541
550 std::size_t
552
556 bool
558
564 n_rows() const;
565
571 n_cols() const;
572
576 unsigned int
577 row_length(const size_type row) const;
578
583 std::size_t
584 memory_consumption() const;
585
607 column_number(const size_type row, const unsigned int index) const;
608
617 row_position(const size_type i, const size_type j) const;
618
630 std::pair<size_type, size_type>
631 matrix_position(const std::size_t global_index) const;
632
648 void
649 print(std::ostream &out) const;
650
664 void
665 print_gnuplot(std::ostream &out) const;
666
674 void
675 print_svg(std::ostream &out) const;
676
682 template <class Archive>
683 void
684 save(Archive &ar, const unsigned int version) const;
685
691 template <class Archive>
692 void
693 load(Archive &ar, const unsigned int version);
694
695#ifdef DOXYGEN
701 template <class Archive>
702 void
703 serialize(Archive &archive, const unsigned int version);
704#else
705 // This macro defines the serialize() method that is compatible with
706 // the templated save() and load() method that have been implemented.
707 BOOST_SERIALIZATION_SPLIT_MEMBER()
708#endif
709
726 "The operation you attempted is only allowed after the SparsityPattern "
727 "has been set up and compress() was called.");
728
733 int,
734 int,
735 << "Upon entering a new entry to row " << arg1
736 << ": there was no free entry any more. " << std::endl
737 << "(Maximum number of entries for this row: " << arg2
738 << "; maybe the matrix is already compressed?)");
739
746 "The operation you attempted changes the structure of the SparsityPattern "
747 "and is not possible after compress() has been called.");
748
754protected:
763
768
773
779 std::size_t max_vec_len;
780
788 unsigned int max_row_length;
789
803 std::unique_ptr<std::size_t[]> rowstart;
804
827 std::unique_ptr<size_type[]> colnums;
828
833
834 // Make all sparse matrices friends of this class.
835 template <typename number>
836 friend class SparseMatrix;
837 template <typename number>
839 template <typename number>
840 friend class SparseILU;
841 template <typename number>
842 friend class ChunkSparseMatrix;
843
846
847 // Also give access to internal details to the iterator/accessor classes.
851};
852
889{
890public:
895
901
910
917
934
951
960 const size_type n,
961 const unsigned int max_per_row);
962
963
973 const size_type n,
974 const std::vector<unsigned int> &row_lengths);
975
984 SparsityPattern(const size_type m, const unsigned int max_per_row);
985
994 const std::vector<unsigned int> &row_lengths);
995
1018 SparsityPattern(const SparsityPattern &original,
1019 const unsigned int max_per_row,
1020 const size_type extra_off_diagonals);
1021
1025 ~SparsityPattern() override = default;
1026
1033 operator=(const SparsityPattern &);
1034
1040 virtual void
1041 reinit(const size_type m,
1042 const size_type n,
1043 const ArrayView<const unsigned int> &row_lengths) override;
1044
1057 void
1058 compress();
1059
1060
1137 template <typename ForwardIterator>
1138 void
1140 const size_type n_cols,
1141 const ForwardIterator begin,
1142 const ForwardIterator end);
1143
1148 void
1150
1155 void
1156 copy_from(const SparsityPattern &sp);
1157
1172 template <typename number>
1173 void
1175
1182 template <typename ForwardIterator>
1183 void
1184 add_entries(const size_type row,
1185 ForwardIterator begin,
1186 ForwardIterator end,
1187 const bool indices_are_sorted = false);
1188
1204 bool
1206
1219 bool
1221
1226 std::size_t
1227 memory_consumption() const;
1228
1261 size_type
1262 operator()(const size_type i, const size_type j) const;
1263
1285 void
1286 block_write(std::ostream &out) const;
1287
1301 void
1302 block_read(std::istream &in);
1303
1308 template <class Archive>
1309 void
1310 save(Archive &ar, const unsigned int version) const;
1311
1316 template <class Archive>
1317 void
1318 load(Archive &ar, const unsigned int version);
1319
1320#ifdef DOXYGEN
1325 template <class Archive>
1326 void
1327 serialize(Archive &archive, const unsigned int version);
1328#else
1329 // This macro defines the serialize() method that is compatible with
1330 // the templated save() and load() method that have been implemented.
1331 BOOST_SERIALIZATION_SPLIT_MEMBER()
1332#endif
1333
1347 int,
1348 int,
1349 << "The iterators denote a range of " << arg1
1350 << " elements, but the given number of rows was " << arg2);
1355 int,
1356 << "The number of partitions you gave is " << arg1
1357 << ", but must be greater than zero.");
1361private:
1366
1367 // Make all sparse matrices friends of this class.
1368 template <typename number>
1369 friend class SparseMatrix;
1370 template <typename number>
1372 template <typename number>
1373 friend class SparseILU;
1374 template <typename number>
1375 friend class ChunkSparseMatrix;
1376
1379
1380 // Also give access to internal details to the iterator/accessor classes.
1384};
1385
1386
1390/*---------------------- Inline functions -----------------------------------*/
1391
1392#ifndef DOXYGEN
1393
1394
1396{
1397 inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern,
1398 const std::size_t i)
1399 : container(sparsity_pattern)
1400 , linear_index(i)
1401 {}
1402
1403
1404
1405 inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern)
1406 : container(sparsity_pattern)
1407 , linear_index(container->rowstart[container->rows])
1408 {}
1409
1410
1411
1412 inline Accessor::Accessor()
1413 : container(nullptr)
1414 , linear_index(numbers::invalid_size_type)
1415 {}
1416
1417
1418
1419 inline bool
1420 Accessor::is_valid_entry() const
1421 {
1422 Assert(container != nullptr, DummyAccessor());
1423 return (linear_index < container->rowstart[container->rows] &&
1424 container->colnums[linear_index] != SparsityPattern::invalid_entry);
1425 }
1426
1427
1428
1429 inline size_type
1430 Accessor::row() const
1431 {
1432 Assert(is_valid_entry() == true, ExcInvalidIterator());
1433
1434 const std::size_t *insert_point =
1435 std::upper_bound(container->rowstart.get(),
1436 container->rowstart.get() + container->rows + 1,
1437 linear_index);
1438 return insert_point - container->rowstart.get() - 1;
1439 }
1440
1441
1442
1443 inline size_type
1444 Accessor::column() const
1445 {
1446 Assert(is_valid_entry() == true, ExcInvalidIterator());
1447
1448 return (container->colnums[linear_index]);
1449 }
1450
1451
1452
1453 inline size_type
1454 Accessor::index() const
1455 {
1456 Assert(is_valid_entry() == true, ExcInvalidIterator());
1457
1458 return linear_index - container->rowstart[row()];
1459 }
1460
1461
1462
1463 inline size_type
1465 {
1466 Assert(is_valid_entry() == true, ExcInvalidIterator());
1467
1468 return linear_index;
1469 }
1470
1471
1472
1473 inline bool
1474 Accessor::operator==(const Accessor &other) const
1475 {
1476 return (container == other.container && linear_index == other.linear_index);
1477 }
1478
1479
1480
1481 inline bool
1482 Accessor::operator<(const Accessor &other) const
1483 {
1484 Assert(container != nullptr, DummyAccessor());
1485 Assert(other.container != nullptr, DummyAccessor());
1486 Assert(container == other.container, ExcInternalError());
1487
1488 return linear_index < other.linear_index;
1489 }
1490
1491
1492
1493 inline void
1495 {
1496 Assert(container != nullptr, DummyAccessor());
1497 Assert(linear_index < container->rowstart[container->rows],
1499 ++linear_index;
1500 }
1501
1502
1503 inline Iterator::Iterator(const SparsityPatternBase *sp,
1504 const std::size_t linear_index)
1505 : LinearIndexIterator<Iterator, Accessor>(Accessor(sp, linear_index))
1506 {}
1507
1508
1509 inline Iterator::Iterator(const Accessor &accessor)
1510 : LinearIndexIterator<Iterator, Accessor>(accessor)
1511 {}
1512
1513
1514} // namespace SparsityPatternIterators
1515
1516
1517
1520{
1521 if (n_rows() > 0)
1522 return {this, rowstart[0]};
1523 else
1524 return end();
1525}
1526
1527
1528
1531{
1532 if (n_rows() > 0)
1533 return {this, rowstart[rows]};
1534 else
1535 return {nullptr, 0};
1536}
1537
1538
1539
1542{
1544
1545 return {this, rowstart[r]};
1546}
1547
1548
1549
1552{
1554
1555 return {this, rowstart[r + 1]};
1556}
1557
1558
1559
1562{
1563 return rows;
1564}
1565
1566
1567
1570{
1571 return cols;
1572}
1573
1574
1575
1576inline bool
1578{
1579 return compressed;
1580}
1581
1582
1583
1584inline bool
1586{
1587 return (store_diagonal_first_in_row == false);
1588}
1589
1590
1591
1592inline unsigned int
1594{
1595 AssertIndexRange(row, rows);
1596 return rowstart[row + 1] - rowstart[row];
1597}
1598
1599
1600
1603 const unsigned int index) const
1604{
1605 AssertIndexRange(row, rows);
1606 AssertIndexRange(index, row_length(row));
1607
1608 return colnums[rowstart[row] + index];
1609}
1610
1611
1612
1613inline std::size_t
1615{
1617
1618 if ((rowstart != nullptr) && (colnums != nullptr))
1619 return rowstart[rows] - rowstart[0];
1620 else
1621 // the object is empty or has zero size
1622 return 0;
1623}
1624
1625
1626
1627template <class Archive>
1628inline void
1629SparsityPatternBase::save(Archive &ar, const unsigned int) const
1630{
1631 // forward to serialization function in the base class.
1632 ar &boost::serialization::base_object<const Subscriptor>(*this);
1633
1635
1636 if (max_dim != 0)
1637 ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1638 else
1639 Assert(rowstart.get() == nullptr, ExcInternalError());
1640
1641 if (max_vec_len != 0)
1642 ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1643 else
1644 Assert(colnums.get() == nullptr, ExcInternalError());
1645}
1646
1647
1648
1649template <class Archive>
1650inline void
1651SparsityPatternBase::load(Archive &ar, const unsigned int)
1652{
1653 // forward to serialization function in the base class.
1654 ar &boost::serialization::base_object<Subscriptor>(*this);
1655
1657
1658 if (max_dim != 0)
1659 {
1660 rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
1661 ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1662 }
1663 else
1664 rowstart.reset();
1665
1666 if (max_vec_len != 0)
1667 {
1668 colnums = std::make_unique<size_type[]>(max_vec_len);
1669 ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1670 }
1671 else
1672 colnums.reset();
1673}
1674
1675
1676
1677template <class Archive>
1678inline void
1679SparsityPattern::save(Archive &ar, const unsigned int) const
1680{
1681 // forward to serialization function in the base class.
1682 ar &boost::serialization::base_object<const SparsityPatternBase>(*this);
1684}
1685
1686
1687
1688template <class Archive>
1689inline void
1690SparsityPattern::load(Archive &ar, const unsigned int)
1691{
1692 // forward to serialization function in the base class.
1693 ar &boost::serialization::base_object<SparsityPatternBase>(*this);
1695}
1696
1697
1698
1699inline bool
1701{
1702 // it isn't quite necessary to compare *all* member variables. by only
1703 // comparing the essential ones, we can say that two sparsity patterns are
1704 // equal even if one is compressed and the other is not (in which case some
1705 // of the member variables are not yet set correctly)
1706 if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed)
1707 return false;
1708
1709 if (rows > 0)
1710 for (size_type i = 0; i < rows + 1; ++i)
1711 if (rowstart[i] != sp2.rowstart[i])
1712 return false;
1713
1714 if (rows > 0)
1715 for (size_type i = 0; i < rowstart[rows]; ++i)
1716 if (colnums[i] != sp2.colnums[i])
1717 return false;
1718
1719 return true;
1720}
1721
1722
1723
1724inline bool
1726{
1727 return (static_cast<const SparsityPatternBase &>(*this) == sp2) &&
1729}
1730
1731
1732
1733namespace internal
1734{
1735 namespace SparsityPatternTools
1736 {
1741
1742 inline size_type
1744 {
1745 return i;
1746 }
1747
1748
1749
1750 template <typename value>
1751 inline size_type
1752 get_column_index_from_iterator(const std::pair<size_type, value> &i)
1753 {
1754 return i.first;
1755 }
1756
1757
1758
1759 template <typename value>
1760 inline size_type
1761 get_column_index_from_iterator(const std::pair<const size_type, value> &i)
1762 {
1763 return i.first;
1764 }
1765 } // namespace SparsityPatternTools
1766} // namespace internal
1767
1768
1769
1770template <typename ForwardIterator>
1771void
1773 const size_type n_cols,
1774 const ForwardIterator begin,
1775 const ForwardIterator end)
1776{
1777 Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1778 ExcIteratorRange(std::distance(begin, end), n_rows));
1779
1780 // first determine row lengths for each row. if the matrix is quadratic,
1781 // then we might have to add an additional entry for the diagonal, if that
1782 // is not yet present. as we have to call compress anyway later on, don't
1783 // bother to check whether that diagonal entry is in a certain row or not
1784 const bool is_square = (n_rows == n_cols);
1785 std::vector<unsigned int> row_lengths;
1786 row_lengths.reserve(n_rows);
1787 for (ForwardIterator i = begin; i != end; ++i)
1788 row_lengths.push_back(std::distance(i->begin(), i->end()) +
1789 (is_square ? 1 : 0));
1790 reinit(n_rows, n_cols, row_lengths);
1791
1792 // now enter all the elements into the matrix. note that if the matrix is
1793 // quadratic, then we already have the diagonal element preallocated
1794 //
1795 // for use in the inner loop, we define an alias to the type of the inner
1796 // iterators
1797 size_type row = 0;
1798 using inner_iterator =
1799 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1800 for (ForwardIterator i = begin; i != end; ++i, ++row)
1801 {
1802 size_type * cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
1803 const inner_iterator end_of_row = i->end();
1804 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1805 {
1806 const size_type col =
1809
1810 if ((col != row) || !is_square)
1811 *cols++ = col;
1812 }
1813 }
1814
1815 // finally compress everything. this also sorts the entries within each row
1816 compress();
1817}
1818
1819
1820#endif // DOXYGEN
1821
1823
1824#endif
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void serialize(Archive &archive, const unsigned int version)
void save(Archive &ar, const unsigned int version) const
std::size_t n_nonzero_elements() const
bool operator==(const SparsityPatternBase &) const
unsigned int max_row_length
void print_svg(std::ostream &out) const
void add(const size_type i, const size_type j)
SparsityPatternIterators::Iterator iterator
types::global_dof_index size_type
void reinit(const size_type m, const size_type n, const unsigned int max_per_row)
unsigned int row_length(const size_type row) const
size_type n_rows() const
std::unique_ptr< size_type[]> colnums
std::unique_ptr< std::size_t[]> rowstart
void print(std::ostream &out) const
bool exists(const size_type i, const size_type j) const
iterator end() const
size_type column_number(const size_type row, const unsigned int index) const
iterator begin(const size_type r) const
virtual void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)=0
void print_gnuplot(std::ostream &out) const
static const size_type invalid_entry
size_type bandwidth() const
bool is_compressed() const
size_type max_entries_per_row() const
void load(Archive &ar, const unsigned int version)
std::pair< size_type, size_type > matrix_position(const std::size_t global_index) const
~SparsityPatternBase() override=default
std::size_t memory_consumption() const
SparsityPatternIterators::Iterator const_iterator
iterator begin() const
size_type row_position(const size_type i, const size_type j) const
iterator end(const size_type r) const
size_type n_cols() const
bool operator==(const Accessor &) const
const SparsityPatternBase * container
bool operator<(const Accessor &) const
Accessor(const SparsityPatternBase *matrix)
Accessor(const SparsityPatternBase *matrix, const std::size_t linear_index)
Iterator(const SparsityPatternBase *sp, const std::size_t linear_index)
Iterator(const Accessor &accessor)
void load(Archive &ar, const unsigned int version)
void block_write(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)
SparsityPattern & operator=(const SparsityPattern &)
SparsityPatternBase::size_type size_type
~SparsityPattern() override=default
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
virtual void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths) override
size_type operator()(const size_type i, const size_type j) const
void save(Archive &ar, const unsigned int version) 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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInvalidIterator()
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & DummyAccessor()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcNotEnoughSpace(int arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcMatrixIsCompressed()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
types::global_dof_index size_type
size_type get_column_index_from_iterator(const size_type i)
size_type get_column_index_from_iterator(const std::pair< const size_type, value > &i)
const types::global_dof_index invalid_size_type
Definition: types.h:210
unsigned int global_dof_index
Definition: types.h:76
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
void advance(std::tuple< I1, I2 > &t, const unsigned int n)