Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2019 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 
22 #include <deal.II/base/array_view.h>
23 #include <deal.II/base/exceptions.h>
24 #include <deal.II/base/linear_index_iterator.h>
25 #include <deal.II/base/std_cxx14/memory.h>
26 #include <deal.II/base/subscriptor.h>
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 <vector>
41 
42 DEAL_II_NAMESPACE_OPEN
43 
44 class SparsityPattern;
48 template <typename number>
49 class FullMatrix;
50 template <typename number>
51 class SparseMatrix;
52 template <typename number>
54 template <typename number>
55 class SparseILU;
56 
58 {
59  class Accessor;
60 }
61 
62 
67 namespace internals
68 {
69  namespace SparsityPatternTools
70  {
74  using size_type = types::global_dof_index;
75 
81  size_type
82  get_column_index_from_iterator(const size_type i);
83 
89  template <typename value>
90  size_type
91  get_column_index_from_iterator(const std::pair<size_type, value> &i);
92 
98  template <typename value>
99  size_type
100  get_column_index_from_iterator(const std::pair<const size_type, value> &i);
101 
102  } // namespace SparsityPatternTools
103 } // namespace internals
104 
105 
110 {
111  // forward declaration
112  class Iterator;
113 
118 
132  class Accessor
133  {
134  public:
139 
143  Accessor(const SparsityPatternBase *matrix, const std::size_t linear_index);
144 
148  Accessor(const SparsityPatternBase *matrix);
149 
155  Accessor();
156 
161  size_type
162  row() const;
163 
169  size_type
170  index() const;
171 
181  size_type
182  global_index() const;
183 
188  size_type
189  column() const;
190 
200  bool
201  is_valid_entry() const;
202 
206  bool
207  operator==(const Accessor &) const;
208 
216  bool
217  operator<(const Accessor &) const;
218 
219  protected:
220  DeclExceptionMsg(DummyAccessor,
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
244  advance();
245 
250 
255  };
256 
257 
258 
283  class Iterator : public LinearIndexIterator<Iterator, Accessor>
284  {
285  public:
290 
295 
301  Iterator(const SparsityPatternBase *sp, const std::size_t linear_index);
302 
307  Iterator(const Accessor &accessor);
308  };
309 } // namespace SparsityPatternIterators
310 
311 
312 
332 {
333 public:
338 
344 
353 
368  static const size_type invalid_entry = numbers::invalid_size_type;
369 
376  // @{
383 
387  ~SparsityPatternBase() override = default;
388 
397  void
398  reinit(const size_type m, const size_type n, const unsigned int max_per_row);
399 
410  void
411  reinit(const size_type m,
412  const size_type n,
413  const std::vector<unsigned int> &row_lengths);
414 
420  virtual void
421  reinit(const size_type m,
422  const size_type n,
423  const ArrayView<const unsigned int> &row_lengths) = 0;
424 
432  void
433  symmetrize();
434 
441  void
442  add(const size_type i, const size_type j);
443 
444  // @}
445 
449  // @{
450 
459  iterator
460  begin() const;
461 
465  iterator
466  end() const;
467 
479  iterator
480  begin(const size_type r) const;
481 
490  iterator
491  end(const size_type r) const;
492 
493 
494  // @}
495 
499  // @{
500 
504  bool
505  operator==(const SparsityPatternBase &) const;
506 
511  bool
512  empty() const;
513 
517  bool
518  exists(const size_type i, const size_type j) const;
519 
525  size_type
526  max_entries_per_row() const;
527 
537  size_type
538  bandwidth() const;
539 
548  std::size_t
549  n_nonzero_elements() const;
550 
554  bool
555  is_compressed() const;
556 
561  size_type
562  n_rows() const;
563 
568  size_type
569  n_cols() const;
570 
574  unsigned int
575  row_length(const size_type row) const;
576 
581  std::size_t
582  memory_consumption() const;
583 
584  // @}
585 
589  // @{
590 
601  size_type
602  column_number(const size_type row, const unsigned int index) const;
603 
611  size_type
612  row_position(const size_type i, const size_type j) const;
613 
625  std::pair<size_type, size_type>
626  matrix_position(const std::size_t global_index) const;
627 
628  // @}
629 
633  // @{
634 
640  void
641  print(std::ostream &out) const;
642 
656  void
657  print_gnuplot(std::ostream &out) const;
658 
666  void
667  print_svg(std::ostream &out) const;
668 
673  template <class Archive>
674  void
675  save(Archive &ar, const unsigned int version) const;
676 
681  template <class Archive>
682  void
683  load(Archive &ar, const unsigned int version);
684 
685  BOOST_SERIALIZATION_SPLIT_MEMBER()
686 
687  // @}
688 
689 
699  ExcNotCompressed,
700  "The operation you attempted is only allowed after the SparsityPattern "
701  "has been set up and compress() was called.");
702 
706  DeclException2(ExcNotEnoughSpace,
707  int,
708  int,
709  << "Upon entering a new entry to row " << arg1
710  << ": there was no free entry any more. " << std::endl
711  << "(Maximum number of entries for this row: " << arg2
712  << "; maybe the matrix is already compressed?)");
713 
719  ExcMatrixIsCompressed,
720  "The operation you attempted changes the structure of the SparsityPattern "
721  "and is not possible after compress() has been called.");
722 
723  // @}
724 
725 
726 protected:
734  size_type max_dim;
735 
739  size_type rows;
740 
744  size_type cols;
745 
751  std::size_t max_vec_len;
752 
760  unsigned int max_row_length;
761 
775  std::unique_ptr<std::size_t[]> rowstart;
776 
799  std::unique_ptr<size_type[]> colnums;
800 
804  bool compressed;
805 
809  template <typename number>
810  friend class SparseMatrix;
811  template <typename number>
812  friend class SparseLUDecomposition;
813  template <typename number>
814  friend class SparseILU;
815  template <typename number>
816  friend class ChunkSparseMatrix;
817 
818  friend class ChunkSparsityPattern;
819  friend class DynamicSparsityPattern;
820 
824  friend class SparsityPatternIterators::Iterator;
825  friend class SparsityPatternIterators::Accessor;
826  friend class ChunkSparsityPatternIterators::Accessor;
827 };
828 
860 {
861 public:
866 
872 
881 
888 
895  // @{
901  SparsityPattern();
902 
919 
927  SparsityPattern(const size_type m,
928  const size_type n,
929  const unsigned int max_per_row);
930 
931 
940  SparsityPattern(const size_type m,
941  const size_type n,
942  const std::vector<unsigned int> &row_lengths);
943 
952  SparsityPattern(const size_type m, const unsigned int max_per_row);
953 
961  SparsityPattern(const size_type m,
962  const std::vector<unsigned int> &row_lengths);
963 
986  SparsityPattern(const SparsityPattern &original,
987  const unsigned int max_per_row,
988  const size_type extra_off_diagonals);
989 
993  ~SparsityPattern() override = default;
994 
1000  SparsityPattern &
1001  operator=(const SparsityPattern &);
1002 
1008  virtual void
1009  reinit(const size_type m,
1010  const size_type n,
1011  const ArrayView<const unsigned int> &row_lengths) override;
1012 
1025  void
1026  compress();
1027 
1028 
1105  template <typename ForwardIterator>
1106  void
1107  copy_from(const size_type n_rows,
1108  const size_type n_cols,
1109  const ForwardIterator begin,
1110  const ForwardIterator end);
1111 
1116  void
1117  copy_from(const DynamicSparsityPattern &dsp);
1118 
1123  void
1124  copy_from(const SparsityPattern &sp);
1125 
1133  template <typename number>
1134  void
1135  copy_from(const FullMatrix<number> &matrix);
1136 
1143  template <typename ForwardIterator>
1144  void
1145  add_entries(const size_type row,
1146  ForwardIterator begin,
1147  ForwardIterator end,
1148  const bool indices_are_sorted = false);
1149 
1150  // @}
1151 
1152 
1156  // @{
1160  bool
1161  operator==(const SparsityPattern &) const;
1162 
1175  bool
1176  stores_only_added_elements() const;
1177 
1182  std::size_t
1183  memory_consumption() const;
1184 
1185  // @}
1189  // @{
1212  size_type
1213  operator()(const size_type i, const size_type j) const;
1214 
1215  // @}
1219  // @{
1220 
1231  void
1232  block_write(std::ostream &out) const;
1233 
1247  void
1248  block_read(std::istream &in);
1249 
1254  template <class Archive>
1255  void
1256  save(Archive &ar, const unsigned int version) const;
1257 
1262  template <class Archive>
1263  void
1264  load(Archive &ar, const unsigned int version);
1265 
1266  BOOST_SERIALIZATION_SPLIT_MEMBER()
1267 
1268  // @}
1269 
1270 
1277  DeclException2(ExcIteratorRange,
1278  int,
1279  int,
1280  << "The iterators denote a range of " << arg1
1281  << " elements, but the given number of rows was " << arg2);
1286  int,
1287  << "The number of partitions you gave is " << arg1
1288  << ", but must be greater than zero.");
1290 private:
1294  bool store_diagonal_first_in_row;
1295 
1299  template <typename number>
1300  friend class SparseMatrix;
1301  template <typename number>
1302  friend class SparseLUDecomposition;
1303  template <typename number>
1304  friend class SparseILU;
1305  template <typename number>
1306  friend class ChunkSparseMatrix;
1307 
1308  friend class ChunkSparsityPattern;
1309  friend class DynamicSparsityPattern;
1310 
1314  friend class SparsityPatternIterators::Iterator;
1315  friend class SparsityPatternIterators::Accessor;
1316  friend class ChunkSparsityPatternIterators::Accessor;
1317 };
1318 
1319 
1321 /*---------------------- Inline functions -----------------------------------*/
1322 
1323 #ifndef DOXYGEN
1324 
1325 
1326 namespace SparsityPatternIterators
1327 {
1328  inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern,
1329  const std::size_t i)
1330  : container(sparsity_pattern)
1331  , linear_index(i)
1332  {}
1333 
1334 
1335 
1336  inline Accessor::Accessor(const SparsityPatternBase *sparsity_pattern)
1337  : container(sparsity_pattern)
1338  , linear_index(container->rowstart[container->rows])
1339  {}
1340 
1341 
1342 
1343  inline Accessor::Accessor()
1344  : container(nullptr)
1345  , linear_index(numbers::invalid_size_type)
1346  {}
1347 
1348 
1349 
1350  inline bool
1351  Accessor::is_valid_entry() const
1352  {
1353  Assert(container != nullptr, DummyAccessor());
1354  return (linear_index < container->rowstart[container->rows] &&
1355  container->colnums[linear_index] != SparsityPattern::invalid_entry);
1356  }
1357 
1358 
1359 
1360  inline size_type
1361  Accessor::row() const
1362  {
1363  Assert(is_valid_entry() == true, ExcInvalidIterator());
1364 
1365  const std::size_t *insert_point =
1366  std::upper_bound(container->rowstart.get(),
1367  container->rowstart.get() + container->rows + 1,
1368  linear_index);
1369  return insert_point - container->rowstart.get() - 1;
1370  }
1371 
1372 
1373 
1374  inline size_type
1375  Accessor::column() const
1376  {
1377  Assert(is_valid_entry() == true, ExcInvalidIterator());
1378 
1379  return (container->colnums[linear_index]);
1380  }
1381 
1382 
1383 
1384  inline size_type
1385  Accessor::index() const
1386  {
1387  Assert(is_valid_entry() == true, ExcInvalidIterator());
1388 
1389  return linear_index - container->rowstart[row()];
1390  }
1391 
1392 
1393 
1394  inline size_type
1395  Accessor::global_index() const
1396  {
1397  Assert(is_valid_entry() == true, ExcInvalidIterator());
1398 
1399  return linear_index;
1400  }
1401 
1402 
1403 
1404  inline bool
1405  Accessor::operator==(const Accessor &other) const
1406  {
1407  Assert(container != nullptr, DummyAccessor());
1408  Assert(other.container != nullptr, DummyAccessor());
1409  return (container == other.container && linear_index == other.linear_index);
1410  }
1411 
1412 
1413 
1414  inline bool
1415  Accessor::operator<(const Accessor &other) const
1416  {
1417  Assert(container != nullptr, DummyAccessor());
1418  Assert(other.container != nullptr, DummyAccessor());
1419  Assert(container == other.container, ExcInternalError());
1420 
1421  return linear_index < other.linear_index;
1422  }
1423 
1424 
1425 
1426  inline void
1428  {
1429  Assert(container != nullptr, DummyAccessor());
1430  Assert(linear_index < container->rowstart[container->rows],
1431  ExcIteratorPastEnd());
1432  ++linear_index;
1433  }
1434 
1435 
1436  inline Iterator::Iterator(const SparsityPatternBase *sp,
1437  const std::size_t linear_index)
1438  : LinearIndexIterator<Iterator, Accessor>(Accessor(sp, linear_index))
1439  {}
1440 
1441 
1442  inline Iterator::Iterator(const Accessor &accessor)
1444  {}
1445 
1446 
1447 } // namespace SparsityPatternIterators
1448 
1449 
1450 
1453 {
1454  if (n_rows() > 0)
1455  return {this, rowstart[0]};
1456  else
1457  return end();
1458 }
1459 
1460 
1461 
1464 {
1465  if (n_rows() > 0)
1466  return {this, rowstart[rows]};
1467  else
1468  return {nullptr, 0};
1469 }
1470 
1471 
1472 
1475 {
1476  Assert(r < n_rows(), ExcIndexRangeType<size_type>(r, 0, n_rows()));
1477 
1478  return {this, rowstart[r]};
1479 }
1480 
1481 
1482 
1484 SparsityPatternBase::end(const size_type r) const
1485 {
1486  Assert(r < n_rows(), ExcIndexRangeType<size_type>(r, 0, n_rows()));
1487 
1488  return {this, rowstart[r + 1]};
1489 }
1490 
1491 
1492 
1495 {
1496  return rows;
1497 }
1498 
1499 
1500 
1503 {
1504  return cols;
1505 }
1506 
1507 
1508 
1509 inline bool
1511 {
1512  return compressed;
1513 }
1514 
1515 
1516 
1517 inline bool
1519 {
1520  return (store_diagonal_first_in_row == false);
1521 }
1522 
1523 
1524 
1525 inline unsigned int
1527 {
1528  Assert(row < rows, ExcIndexRangeType<size_type>(row, 0, rows));
1529  return rowstart[row + 1] - rowstart[row];
1530 }
1531 
1532 
1533 
1536  const unsigned int index) const
1537 {
1538  Assert(row < rows, ExcIndexRangeType<size_type>(row, 0, rows));
1539  Assert(index < row_length(row), ExcIndexRange(index, 0, row_length(row)));
1540 
1541  return colnums[rowstart[row] + index];
1542 }
1543 
1544 
1545 
1546 inline std::size_t
1548 {
1549  Assert(compressed, ExcNotCompressed());
1550 
1551  if ((rowstart != nullptr) && (colnums != nullptr))
1552  return rowstart[rows] - rowstart[0];
1553  else
1554  // the object is empty or has zero size
1555  return 0;
1556 }
1557 
1558 
1559 
1560 template <class Archive>
1561 inline void
1562 SparsityPatternBase::save(Archive &ar, const unsigned int) const
1563 {
1564  // forward to serialization function in the base class.
1565  ar &boost::serialization::base_object<const Subscriptor>(*this);
1566 
1567  ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed;
1568 
1569  ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1570  ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1571 }
1572 
1573 
1574 
1575 template <class Archive>
1576 inline void
1577 SparsityPatternBase::load(Archive &ar, const unsigned int)
1578 {
1579  // forward to serialization function in the base class.
1580  ar &boost::serialization::base_object<Subscriptor>(*this);
1581 
1582  ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed;
1583 
1584  rowstart = std_cxx14::make_unique<std::size_t[]>(max_dim + 1);
1585  colnums = std_cxx14::make_unique<size_type[]>(max_vec_len);
1586 
1587  ar &boost::serialization::make_array(rowstart.get(), max_dim + 1);
1588  ar &boost::serialization::make_array(colnums.get(), max_vec_len);
1589 }
1590 
1591 
1592 
1593 template <class Archive>
1594 inline void
1595 SparsityPattern::save(Archive &ar, const unsigned int) const
1596 {
1597  // forward to serialization function in the base class.
1598  ar &boost::serialization::base_object<const SparsityPatternBase>(*this);
1599  ar &store_diagonal_first_in_row;
1600 }
1601 
1602 
1603 
1604 template <class Archive>
1605 inline void
1606 SparsityPattern::load(Archive &ar, const unsigned int)
1607 {
1608  // forward to serialization function in the base class.
1609  ar &boost::serialization::base_object<SparsityPatternBase>(*this);
1610  ar &store_diagonal_first_in_row;
1611 }
1612 
1613 
1614 
1615 inline bool
1616 SparsityPatternBase::operator==(const SparsityPatternBase &sp2) const
1617 {
1618  // it isn't quite necessary to compare *all* member variables. by only
1619  // comparing the essential ones, we can say that two sparsity patterns are
1620  // equal even if one is compressed and the other is not (in which case some
1621  // of the member variables are not yet set correctly)
1622  if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed)
1623  return false;
1624 
1625  for (size_type i = 0; i < rows + 1; ++i)
1626  if (rowstart[i] != sp2.rowstart[i])
1627  return false;
1628 
1629  for (size_type i = 0; i < rowstart[rows]; ++i)
1630  if (colnums[i] != sp2.colnums[i])
1631  return false;
1632 
1633  return true;
1634 }
1635 
1636 
1637 
1638 inline bool
1640 {
1641  return (static_cast<const SparsityPatternBase &>(*this) == sp2) &&
1642  (store_diagonal_first_in_row == sp2.store_diagonal_first_in_row);
1643 }
1644 
1645 
1646 
1647 namespace internal
1648 {
1649  namespace SparsityPatternTools
1650  {
1655 
1656  inline size_type
1657  get_column_index_from_iterator(const size_type i)
1658  {
1659  return i;
1660  }
1661 
1662 
1663 
1664  template <typename value>
1665  inline size_type
1666  get_column_index_from_iterator(const std::pair<size_type, value> &i)
1667  {
1668  return i.first;
1669  }
1670 
1671 
1672 
1673  template <typename value>
1674  inline size_type
1675  get_column_index_from_iterator(const std::pair<const size_type, value> &i)
1676  {
1677  return i.first;
1678  }
1679  } // namespace SparsityPatternTools
1680 } // namespace internal
1681 
1682 
1683 
1684 template <typename ForwardIterator>
1685 void
1687  const size_type n_cols,
1688  const ForwardIterator begin,
1689  const ForwardIterator end)
1690 {
1691  Assert(static_cast<size_type>(std::distance(begin, end)) == n_rows,
1692  ExcIteratorRange(std::distance(begin, end), n_rows));
1693 
1694  // first determine row lengths for each row. if the matrix is quadratic,
1695  // then we might have to add an additional entry for the diagonal, if that
1696  // is not yet present. as we have to call compress anyway later on, don't
1697  // bother to check whether that diagonal entry is in a certain row or not
1698  const bool is_square = (n_rows == n_cols);
1699  std::vector<unsigned int> row_lengths;
1700  row_lengths.reserve(n_rows);
1701  for (ForwardIterator i = begin; i != end; ++i)
1702  row_lengths.push_back(std::distance(i->begin(), i->end()) +
1703  (is_square ? 1 : 0));
1704  reinit(n_rows, n_cols, row_lengths);
1705 
1706  // now enter all the elements into the matrix. note that if the matrix is
1707  // quadratic, then we already have the diagonal element preallocated
1708  //
1709  // for use in the inner loop, we define an alias to the type of the inner
1710  // iterators
1711  size_type row = 0;
1712  using inner_iterator =
1713  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1714  for (ForwardIterator i = begin; i != end; ++i, ++row)
1715  {
1716  size_type * cols = &colnums[rowstart[row]] + (is_square ? 1 : 0);
1717  const inner_iterator end_of_row = i->end();
1718  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1719  {
1720  const size_type col =
1721  internal::SparsityPatternTools::get_column_index_from_iterator(*j);
1722  Assert(col < n_cols, ExcIndexRange(col, 0, n_cols));
1723 
1724  if ((col != row) || !is_square)
1725  *cols++ = col;
1726  }
1727  }
1728 
1729  // finally compress everything. this also sorts the entries within each row
1730  compress();
1731 }
1732 
1733 
1734 #endif // DOXYGEN
1735 
1736 DEAL_II_NAMESPACE_CLOSE
1737 
1738 #endif
SparsityPatternIterators::size_type size_type
const types::global_dof_index invalid_size_type
Definition: types.h:182
iterator end() const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
iterator begin() const
void load(Archive &ar, const unsigned int version)
types::global_dof_index size_type
Iterator(const ChunkSparsityPattern *sp, const size_type row)
STL namespace.
bool operator==(const SparsityPatternBase &) const
types::global_dof_index size_type
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void reinit(const size_type m, const size_type n, const unsigned int max_per_row)
std::unique_ptr< std::size_t[]> rowstart
static ::ExceptionBase & ExcInvalidIterator()
static ::ExceptionBase & ExcInvalidNumberOfPartitions(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
SparsityPatternBase::size_type size_type
#define Assert(cond, exc)
Definition: exceptions.h:1407
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:496
SparsityPatternIterators::Iterator const_iterator
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
void load(Archive &ar, const unsigned int version)
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
std::size_t n_nonzero_elements() const
static ::ExceptionBase & ExcIteratorPastEnd()
std::unique_ptr< size_type[]> colnums
unsigned int global_dof_index
Definition: types.h:89
__global__ void set(Number *val, const Number s, const size_type N)
bool is_compressed() const
bool operator==(const SparsityPattern &) const
bool stores_only_added_elements() const
void save(Archive &ar, const unsigned int version) const
static const size_type invalid_entry
types::global_dof_index size_type
bool operator<(const Accessor &) const
size_type column_number(const size_type row, const unsigned int index) const
void save(Archive &ar, const unsigned int version) const
unsigned int row_length(const size_type row) const
size_type n_rows() const
bool operator==(const Accessor &) const
const SparsityPatternBase * container
static ::ExceptionBase & ExcInternalError()
SparsityPatternIterators::Iterator iterator
size_type n_cols() const