Reference documentation for deal.II version GIT 69313620a0 2022-11-28 13:00: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\}}\)
chunk_sparse_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2021 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_chunk_sparse_matrix_h
17 #define dealii_chunk_sparse_matrix_h
18 
19 
20 #include <deal.II/base/config.h>
21 
24 
26 #include <deal.II/lac/exceptions.h>
28 
29 #include <iterator>
30 #include <memory>
31 
32 
34 
35 // Forward declarations
36 #ifndef DOXYGEN
37 template <typename number>
38 class Vector;
39 template <typename number>
40 class FullMatrix;
41 #endif
42 
53 {
54  // forward declaration
55  template <typename number, bool Constness>
56  class Iterator;
57 
68  template <typename number, bool Constness>
70  {
71  public:
75  number
76  value() const;
77 
81  number &
82  value();
83 
89  get_matrix() const;
90  };
91 
92 
93 
100  template <typename number>
102  {
103  public:
109 
113  Accessor(MatrixType *matrix, const unsigned int row);
114 
119 
124 
128  number
129  value() const;
130 
135  const MatrixType &
136  get_matrix() const;
137 
138  private:
143 
148 
149  // Make iterator class a friend.
150  template <typename, bool>
151  friend class Iterator;
152  };
153 
154 
161  template <typename number>
162  class Accessor<number, false> : public ChunkSparsityPatternIterators::Accessor
163  {
164  private:
187  class Reference
188  {
189  public:
194  Reference(const Accessor *accessor, const bool dummy);
195 
199  operator number() const;
200 
204  const Reference &
205  operator=(const number n) const;
206 
210  const Reference &
211  operator+=(const number n) const;
212 
216  const Reference &
217  operator-=(const number n) const;
218 
222  const Reference &
223  operator*=(const number n) const;
224 
228  const Reference &
229  operator/=(const number n) const;
230 
231  private:
237  };
238 
239  public:
245 
249  Accessor(MatrixType *matrix, const unsigned int row);
250 
255 
259  Reference
260  value() const;
261 
266  MatrixType &
267  get_matrix() const;
268 
269  private:
274 
279 
280  // Make iterator class a friend.
281  template <typename, bool>
282  friend class Iterator;
283  };
284 
285 
286 
297  template <typename number, bool Constness>
298  class Iterator
299  {
300  public:
305 
311 
317 
322  Iterator(MatrixType *matrix, const unsigned int row);
323 
328 
334 
338  Iterator &
340 
344  Iterator
346 
351  operator*() const;
352 
357  operator->() const;
358 
362  bool
363  operator==(const Iterator &) const;
364 
368  bool
369  operator!=(const Iterator &) const;
370 
378  bool
379  operator<(const Iterator &) const;
380 
385  bool
386  operator>(const Iterator &) const;
387 
394  int
395  operator-(const Iterator &p) const;
396 
400  Iterator
401  operator+(const unsigned int n) const;
402 
403  private:
408  };
409 
410 } // namespace ChunkSparseMatrixIterators
411 
413 
414 namespace std
415 {
416  template <typename number, bool Constness>
417  struct iterator_traits<
418  ::ChunkSparseMatrixIterators::Iterator<number, Constness>>
419  {
420  using iterator_category = forward_iterator_tag;
421  using value_type = typename ::ChunkSparseMatrixIterators::
422  Iterator<number, Constness>::value_type;
423  using difference_type = typename ::ChunkSparseMatrixIterators::
424  Iterator<number, Constness>::difference_type;
425  };
426 } // namespace std
427 
429 
430 
431 
448 template <typename number>
449 class ChunkSparseMatrix : public virtual Subscriptor
450 {
451 public:
456 
461  using value_type = number;
462 
473 
479 
487 
494  struct Traits
495  {
500  static const bool zero_addition_can_be_elided = true;
501  };
502 
518 
528 
542  explicit ChunkSparseMatrix(const ChunkSparsityPattern &sparsity);
543 
551  const IdentityMatrix & id);
552 
557  virtual ~ChunkSparseMatrix() override;
558 
570 
579 
590  operator=(const double d);
591 
605  virtual void
606  reinit(const ChunkSparsityPattern &sparsity);
607 
613  virtual void
614  clear();
624  bool
625  empty() const;
626 
631  size_type
632  m() const;
633 
638  size_type
639  n() const;
640 
646  size_type
648 
656  size_type
658 
667  const ChunkSparsityPattern &
669 
674  std::size_t
676 
687  void
688  set(const size_type i, const size_type j, const number value);
689 
695  void
696  add(const size_type i, const size_type j, const number value);
697 
707  template <typename number2>
708  void
709  add(const size_type row,
710  const size_type n_cols,
711  const size_type *col_indices,
712  const number2 * values,
713  const bool elide_zero_values = true,
714  const bool col_indices_are_sorted = false);
715 
720  operator*=(const number factor);
721 
726  operator/=(const number factor);
727 
740  void
742 
759  template <typename somenumber>
762 
780  template <typename ForwardIterator>
781  void
782  copy_from(const ForwardIterator begin, const ForwardIterator end);
783 
789  template <typename somenumber>
790  void
792 
804  template <typename somenumber>
805  void
806  add(const number factor, const ChunkSparseMatrix<somenumber> &matrix);
807 
827  number
828  operator()(const size_type i, const size_type j) const;
829 
842  number
843  el(const size_type i, const size_type j) const;
844 
854  number
855  diag_element(const size_type i) const;
856 
866  void
868  const size_type array_length,
869  size_type & row_length,
870  size_type * column_indices,
871  number * values) const;
872 
892  template <class OutVector, class InVector>
893  void
894  vmult(OutVector &dst, const InVector &src) const;
895 
911  template <class OutVector, class InVector>
912  void
913  Tvmult(OutVector &dst, const InVector &src) const;
914 
929  template <class OutVector, class InVector>
930  void
931  vmult_add(OutVector &dst, const InVector &src) const;
932 
948  template <class OutVector, class InVector>
949  void
950  Tvmult_add(OutVector &dst, const InVector &src) const;
951 
967  template <typename somenumber>
968  somenumber
970 
974  template <typename somenumber>
975  somenumber
977  const Vector<somenumber> &v) const;
985  template <typename somenumber>
986  somenumber
988  const Vector<somenumber> &x,
989  const Vector<somenumber> &b) const;
990 
1004  real_type
1005  l1_norm() const;
1006 
1014  real_type
1015  linfty_norm() const;
1016 
1021  real_type
1034  template <typename somenumber>
1035  void
1037  const Vector<somenumber> &src,
1038  const number omega = 1.) const;
1039 
1043  template <typename somenumber>
1044  void
1046  const Vector<somenumber> &src,
1047  const number om = 1.) const;
1048 
1052  template <typename somenumber>
1053  void
1055  const Vector<somenumber> &src,
1056  const number om = 1.) const;
1057 
1061  template <typename somenumber>
1062  void
1064  const Vector<somenumber> &src,
1065  const number om = 1.) const;
1066 
1072  template <typename somenumber>
1073  void
1074  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1075 
1080  template <typename somenumber>
1081  void
1082  SOR(Vector<somenumber> &v, const number om = 1.) const;
1083 
1088  template <typename somenumber>
1089  void
1090  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1091 
1102  template <typename somenumber>
1103  void
1105  const std::vector<size_type> &permutation,
1106  const std::vector<size_type> &inverse_permutation,
1107  const number om = 1.) const;
1108 
1119  template <typename somenumber>
1120  void
1122  const std::vector<size_type> &permutation,
1123  const std::vector<size_type> &inverse_permutation,
1124  const number om = 1.) const;
1125 
1130  template <typename somenumber>
1131  void
1133  const Vector<somenumber> &b,
1134  const number om = 1.) const;
1135 
1140  template <typename somenumber>
1141  void
1143  const Vector<somenumber> &b,
1144  const number om = 1.) const;
1145 
1150  template <typename somenumber>
1151  void
1153  const Vector<somenumber> &b,
1154  const number om = 1.) const;
1171  begin() const;
1172 
1182  end() const;
1183 
1193  iterator
1195 
1204  iterator
1205  end();
1206 
1222  begin(const unsigned int r) const;
1223 
1239  end(const unsigned int r) const;
1240 
1255  iterator
1256  begin(const unsigned int r);
1257 
1272  iterator
1273  end(const unsigned int r);
1284  void
1285  print(std::ostream &out) const;
1286 
1307  void
1308  print_formatted(std::ostream & out,
1309  const unsigned int precision = 3,
1310  const bool scientific = true,
1311  const unsigned int width = 0,
1312  const char * zero_string = " ",
1313  const double denominator = 1.) const;
1314 
1320  void
1321  print_pattern(std::ostream &out, const double threshold = 0.) const;
1322 
1333  void
1334  block_write(std::ostream &out) const;
1335 
1352  void
1353  block_read(std::istream &in);
1364  int,
1365  int,
1366  << "You are trying to access the matrix entry with index <"
1367  << arg1 << ',' << arg2
1368  << ">, but this entry does not exist in the sparsity pattern "
1369  "of this matrix."
1370  "\n\n"
1371  "The most common cause for this problem is that you used "
1372  "a method to build the sparsity pattern that did not "
1373  "(completely) take into account all of the entries you "
1374  "will later try to write into. An example would be "
1375  "building a sparsity pattern that does not include "
1376  "the entries you will write into due to constraints "
1377  "on degrees of freedom such as hanging nodes or periodic "
1378  "boundary conditions. In such cases, building the "
1379  "sparsity pattern will succeed, but you will get errors "
1380  "such as the current one at one point or other when "
1381  "trying to write into the entries of the matrix.");
1390  int,
1391  int,
1392  << "The iterators denote a range of " << arg1
1393  << " elements, but the given number of rows was " << arg2);
1398  "You are attempting an operation on two vectors that "
1399  "are the same object, but the operation requires that the "
1400  "two objects are in fact different.");
1402 private:
1409 
1417  std::unique_ptr<number[]> val;
1418 
1426 
1430  size_type
1431  compute_location(const size_type i, const size_type j) const;
1432 
1433  // make all other sparse matrices friends
1434  template <typename somenumber>
1435  friend class ChunkSparseMatrix;
1436 
1437  // Also give access to internal details to the iterator/accessor classes.
1438  template <typename, bool>
1440  template <typename, bool>
1442 };
1443 
1446 #ifndef DOXYGEN
1447 /*---------------------- Inline functions -----------------------------------*/
1448 
1449 
1450 
1451 template <typename number>
1454 {
1455  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1456  return cols->rows;
1457 }
1458 
1459 
1460 template <typename number>
1463 {
1464  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1465  return cols->cols;
1466 }
1467 
1468 
1469 
1470 template <typename number>
1471 inline const ChunkSparsityPattern &
1473 {
1474  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1475  return *cols;
1476 }
1477 
1478 
1479 
1480 template <typename number>
1483  const size_type j) const
1484 {
1485  const size_type chunk_size = cols->get_chunk_size();
1486  const size_type chunk_index =
1487  cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1488 
1489  if (chunk_index == ChunkSparsityPattern::invalid_entry)
1491  else
1492  {
1493  return (chunk_index * chunk_size * chunk_size +
1494  (i % chunk_size) * chunk_size + (j % chunk_size));
1495  }
1496 }
1497 
1498 
1499 template <typename number>
1500 inline void
1502  const size_type j,
1503  const number value)
1504 {
1506 
1507  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1508  // it is allowed to set elements of the matrix that are not part of the
1509  // sparsity pattern, if the value to which we set it is zero
1510  const size_type index = compute_location(i, j);
1511  Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1512  ExcInvalidIndex(i, j));
1513 
1514  if (index != SparsityPattern::invalid_entry)
1515  val[index] = value;
1516 }
1517 
1518 
1519 
1520 template <typename number>
1521 inline void
1523  const size_type j,
1524  const number value)
1525 {
1527 
1528  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1529 
1530  if (std::abs(value) != 0.)
1531  {
1532  const size_type index = compute_location(i, j);
1534  ExcInvalidIndex(i, j));
1535 
1536  val[index] += value;
1537  }
1538 }
1539 
1540 
1541 
1542 template <typename number>
1543 template <typename number2>
1544 inline void
1546  const size_type n_cols,
1547  const size_type *col_indices,
1548  const number2 * values,
1549  const bool /*elide_zero_values*/,
1550  const bool /*col_indices_are_sorted*/)
1551 {
1552  // TODO: could be done more efficiently...
1553  for (size_type col = 0; col < n_cols; ++col)
1554  add(row, col_indices[col], static_cast<number>(values[col]));
1555 }
1556 
1557 
1558 
1559 template <typename number>
1561 ChunkSparseMatrix<number>::operator*=(const number factor)
1562 {
1563  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1564  Assert(val != nullptr, ExcNotInitialized());
1565 
1566  const size_type chunk_size = cols->get_chunk_size();
1567 
1568  // multiply all elements of the matrix with the given factor. this includes
1569  // the padding elements in chunks that overlap the boundaries of the actual
1570  // matrix -- but since multiplication with a number does not violate the
1571  // invariant of keeping these elements at zero nothing can happen
1572  number * val_ptr = val.get();
1573  const number *const end_ptr =
1574  val.get() +
1575  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1576  while (val_ptr != end_ptr)
1577  *val_ptr++ *= factor;
1578 
1579  return *this;
1580 }
1581 
1582 
1583 
1584 template <typename number>
1586 ChunkSparseMatrix<number>::operator/=(const number factor)
1587 {
1588  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1589  Assert(val != nullptr, ExcNotInitialized());
1590  Assert(std::abs(factor) != 0, ExcDivideByZero());
1591 
1592  const number factor_inv = 1. / factor;
1593 
1594  const size_type chunk_size = cols->get_chunk_size();
1595 
1596  // multiply all elements of the matrix with the given factor. this includes
1597  // the padding elements in chunks that overlap the boundaries of the actual
1598  // matrix -- but since multiplication with a number does not violate the
1599  // invariant of keeping these elements at zero nothing can happen
1600  number * val_ptr = val.get();
1601  const number *const end_ptr =
1602  val.get() +
1603  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1604 
1605  while (val_ptr != end_ptr)
1606  *val_ptr++ *= factor_inv;
1607 
1608  return *this;
1609 }
1610 
1611 
1612 
1613 template <typename number>
1614 inline number
1616  const size_type j) const
1617 {
1618  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1619  AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1620  ExcInvalidIndex(i, j));
1621  return val[compute_location(i, j)];
1622 }
1623 
1624 
1625 
1626 template <typename number>
1627 inline number
1628 ChunkSparseMatrix<number>::el(const size_type i, const size_type j) const
1629 {
1630  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1631  const size_type index = compute_location(i, j);
1632 
1634  return val[index];
1635  else
1636  return 0;
1637 }
1638 
1639 
1640 
1641 template <typename number>
1642 inline number
1644 {
1645  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1646  Assert(m() == n(), ExcNotQuadratic());
1647  AssertIndexRange(i, m());
1648 
1649  // Use that the first element in each row of a quadratic matrix is the main
1650  // diagonal of the chunk sparsity pattern
1651  const size_type chunk_size = cols->get_chunk_size();
1652  return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1653  chunk_size +
1654  (i % chunk_size) * chunk_size + (i % chunk_size)];
1655 }
1656 
1657 
1658 
1659 template <typename number>
1660 template <typename ForwardIterator>
1661 inline void
1662 ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1663  const ForwardIterator end)
1664 {
1665  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1666  ExcIteratorRange(std::distance(begin, end), m()));
1667 
1668  // for use in the inner loop, we define an alias to the type of the inner
1669  // iterators
1670  using inner_iterator =
1671  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1672  size_type row = 0;
1673  for (ForwardIterator i = begin; i != end; ++i, ++row)
1674  {
1675  const inner_iterator end_of_row = i->end();
1676  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1677  // write entries
1678  set(row, j->first, j->second);
1679  }
1680 }
1681 
1682 
1683 
1684 //---------------------------------------------------------------------------
1685 
1686 
1688 {
1689  template <typename number>
1690  inline Accessor<number, true>::Accessor(const MatrixType * matrix,
1691  const unsigned int row)
1692  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1693  row)
1694  , matrix(matrix)
1695  {}
1696 
1697 
1698 
1699  template <typename number>
1700  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1701  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1702  , matrix(matrix)
1703  {}
1704 
1705 
1706 
1707  template <typename number>
1708  inline Accessor<number, true>::Accessor(
1710  : ChunkSparsityPatternIterators::Accessor(a)
1711  , matrix(&a.get_matrix())
1712  {}
1713 
1714 
1715 
1716  template <typename number>
1717  inline number
1718  Accessor<number, true>::value() const
1719  {
1720  const unsigned int chunk_size =
1721  matrix->get_sparsity_pattern().get_chunk_size();
1722  return matrix->val[reduced_index() * chunk_size * chunk_size +
1723  chunk_row * chunk_size + chunk_col];
1724  }
1725 
1726 
1727 
1728  template <typename number>
1729  inline const typename Accessor<number, true>::MatrixType &
1730  Accessor<number, true>::get_matrix() const
1731  {
1732  return *matrix;
1733  }
1734 
1735 
1736 
1737  template <typename number>
1738  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1739  const bool)
1740  : accessor(accessor)
1741  {}
1742 
1743 
1744  template <typename number>
1745  inline Accessor<number, false>::Reference::operator number() const
1746  {
1747  const unsigned int chunk_size =
1748  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1749  return accessor->matrix
1750  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1751  accessor->chunk_row * chunk_size + accessor->chunk_col];
1752  }
1753 
1754 
1755 
1756  template <typename number>
1757  inline const typename Accessor<number, false>::Reference &
1758  Accessor<number, false>::Reference::operator=(const number n) const
1759  {
1760  const unsigned int chunk_size =
1761  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1762  accessor->matrix
1763  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1764  accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1765  return *this;
1766  }
1767 
1768 
1769 
1770  template <typename number>
1771  inline const typename Accessor<number, false>::Reference &
1772  Accessor<number, false>::Reference::operator+=(const number n) const
1773  {
1774  const unsigned int chunk_size =
1775  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1776  accessor->matrix
1777  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1778  accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1779  return *this;
1780  }
1781 
1782 
1783 
1784  template <typename number>
1785  inline const typename Accessor<number, false>::Reference &
1786  Accessor<number, false>::Reference::operator-=(const number n) const
1787  {
1788  const unsigned int chunk_size =
1789  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1790  accessor->matrix
1791  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1792  accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1793  return *this;
1794  }
1795 
1796 
1797 
1798  template <typename number>
1799  inline const typename Accessor<number, false>::Reference &
1800  Accessor<number, false>::Reference::operator*=(const number n) const
1801  {
1802  const unsigned int chunk_size =
1803  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1804  accessor->matrix
1805  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1806  accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1807  return *this;
1808  }
1809 
1810 
1811 
1812  template <typename number>
1813  inline const typename Accessor<number, false>::Reference &
1814  Accessor<number, false>::Reference::operator/=(const number n) const
1815  {
1816  const unsigned int chunk_size =
1817  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1818  accessor->matrix
1819  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1820  accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1821  return *this;
1822  }
1823 
1824 
1825 
1826  template <typename number>
1827  inline Accessor<number, false>::Accessor(MatrixType * matrix,
1828  const unsigned int row)
1829  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1830  row)
1831  , matrix(matrix)
1832  {}
1833 
1834 
1835 
1836  template <typename number>
1837  inline Accessor<number, false>::Accessor(MatrixType *matrix)
1838  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1839  , matrix(matrix)
1840  {}
1841 
1842 
1843 
1844  template <typename number>
1845  inline typename Accessor<number, false>::Reference
1846  Accessor<number, false>::value() const
1847  {
1848  return Reference(this, true);
1849  }
1850 
1851 
1852 
1853  template <typename number>
1854  inline typename Accessor<number, false>::MatrixType &
1855  Accessor<number, false>::get_matrix() const
1856  {
1857  return *matrix;
1858  }
1859 
1860 
1861 
1862  template <typename number, bool Constness>
1863  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
1864  const unsigned int row)
1865  : accessor(matrix, row)
1866  {}
1867 
1868 
1869 
1870  template <typename number, bool Constness>
1871  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1872  : accessor(matrix)
1873  {}
1874 
1875 
1876 
1877  template <typename number, bool Constness>
1878  inline Iterator<number, Constness>::Iterator(
1880  : accessor(*i)
1881  {}
1882 
1883 
1884 
1885  template <typename number, bool Constness>
1886  inline Iterator<number, Constness> &
1888  {
1889  accessor.advance();
1890  return *this;
1891  }
1892 
1893 
1894  template <typename number, bool Constness>
1895  inline Iterator<number, Constness>
1897  {
1898  const Iterator iter = *this;
1899  accessor.advance();
1900  return iter;
1901  }
1902 
1903 
1904  template <typename number, bool Constness>
1905  inline const Accessor<number, Constness> &
1907  {
1908  return accessor;
1909  }
1910 
1911 
1912  template <typename number, bool Constness>
1913  inline const Accessor<number, Constness> *
1914  Iterator<number, Constness>::operator->() const
1915  {
1916  return &accessor;
1917  }
1918 
1919 
1920  template <typename number, bool Constness>
1921  inline bool
1922  Iterator<number, Constness>::operator==(const Iterator &other) const
1923  {
1924  return (accessor == other.accessor);
1925  }
1926 
1927 
1928  template <typename number, bool Constness>
1929  inline bool
1930  Iterator<number, Constness>::operator!=(const Iterator &other) const
1931  {
1932  return !(*this == other);
1933  }
1934 
1935 
1936  template <typename number, bool Constness>
1937  inline bool
1938  Iterator<number, Constness>::operator<(const Iterator &other) const
1939  {
1940  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1941  ExcInternalError());
1942 
1943  return (accessor < other.accessor);
1944  }
1945 
1946 
1947  template <typename number, bool Constness>
1948  inline bool
1949  Iterator<number, Constness>::operator>(const Iterator &other) const
1950  {
1951  return (other < *this);
1952  }
1953 
1954 
1955  template <typename number, bool Constness>
1956  inline int
1957  Iterator<number, Constness>::operator-(const Iterator &other) const
1958  {
1959  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1960  ExcInternalError());
1961 
1962  // TODO: can be optimized
1963  int difference = 0;
1964  if (*this < other)
1965  {
1966  Iterator copy = *this;
1967  while (copy != other)
1968  {
1969  ++copy;
1970  --difference;
1971  }
1972  }
1973  else
1974  {
1975  Iterator copy = other;
1976  while (copy != *this)
1977  {
1978  ++copy;
1979  ++difference;
1980  }
1981  }
1982  return difference;
1983  }
1984 
1985 
1986 
1987  template <typename number, bool Constness>
1988  inline Iterator<number, Constness>
1989  Iterator<number, Constness>::operator+(const unsigned int n) const
1990  {
1991  Iterator x = *this;
1992  for (unsigned int i = 0; i < n; ++i)
1993  ++x;
1994 
1995  return x;
1996  }
1997 
1998 } // namespace ChunkSparseMatrixIterators
1999 
2000 
2001 
2002 template <typename number>
2005 {
2006  return const_iterator(this, 0);
2007 }
2008 
2009 
2010 template <typename number>
2013 {
2014  return const_iterator(this);
2015 }
2016 
2017 
2018 template <typename number>
2021 {
2022  return iterator(this, 0);
2023 }
2024 
2025 
2026 template <typename number>
2029 {
2030  return iterator(this);
2031 }
2032 
2033 
2034 template <typename number>
2036 ChunkSparseMatrix<number>::begin(const unsigned int r) const
2037 {
2038  AssertIndexRange(r, m());
2039  return const_iterator(this, r);
2040 }
2041 
2042 
2043 
2044 template <typename number>
2046 ChunkSparseMatrix<number>::end(const unsigned int r) const
2047 {
2048  AssertIndexRange(r, m());
2049  return const_iterator(this, r + 1);
2050 }
2051 
2052 
2053 
2054 template <typename number>
2056 ChunkSparseMatrix<number>::begin(const unsigned int r)
2057 {
2058  AssertIndexRange(r, m());
2059  return iterator(this, r);
2060 }
2061 
2062 
2063 
2064 template <typename number>
2066 ChunkSparseMatrix<number>::end(const unsigned int r)
2067 {
2068  AssertIndexRange(r, m());
2069  return iterator(this, r + 1);
2070 }
2071 
2072 
2073 
2074 #endif // DOXYGEN
2075 
2077 
2078 #endif
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
Reference(const Accessor *accessor, const bool dummy)
Accessor(MatrixType *matrix, const unsigned int row)
Accessor(MatrixType *matrix, const unsigned int row)
Accessor(const ChunkSparseMatrixIterators::Accessor< number, false > &a)
const ChunkSparseMatrix< number > & get_matrix() const
bool operator>(const Iterator &) const
Iterator(const ChunkSparseMatrixIterators::Iterator< number, false > &i)
const Accessor< number, Constness > & value_type
const Accessor< number, Constness > * operator->() const
bool operator<(const Iterator &) const
Iterator(MatrixType *matrix, const unsigned int row)
bool operator==(const Iterator &) const
Accessor< number, Constness > accessor
int operator-(const Iterator &p) const
const Accessor< number, Constness > & operator*() const
typename Accessor< number, Constness >::MatrixType MatrixType
bool operator!=(const Iterator &) const
Iterator operator+(const unsigned int n) const
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
ChunkSparseMatrix & operator/=(const number factor)
ChunkSparseMatrix & operator*=(const number factor)
ChunkSparseMatrix< number > & operator=(const IdentityMatrix &id)
real_type linfty_norm() const
size_type compute_location(const size_type i, const size_type j) const
ChunkSparseMatrix(const ChunkSparsityPattern &sparsity)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
void Tvmult(OutVector &dst, const InVector &src) const
ChunkSparseMatrix< number > & operator=(const ChunkSparseMatrix< number > &)
const ChunkSparsityPattern & get_sparsity_pattern() const
const_iterator begin() const
void add(const size_type i, const size_type j, const number value)
void set(const size_type i, const size_type j, const number value)
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const_iterator begin(const unsigned int r) const
void print(std::ostream &out) const
bool empty() const
iterator begin()
virtual void reinit(const ChunkSparsityPattern &sparsity)
void extract_row_copy(const size_type row, const size_type array_length, size_type &row_length, size_type *column_indices, number *values) const
iterator end(const unsigned int r)
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
void print_pattern(std::ostream &out, const double threshold=0.) const
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void copy_from(const ForwardIterator begin, const ForwardIterator end)
real_type l1_norm() const
const_iterator end() const
void vmult_add(OutVector &dst, const InVector &src) const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
virtual ~ChunkSparseMatrix() override
std::unique_ptr< number[]> val
void block_write(std::ostream &out) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type m() const
size_type n_actually_nonzero_elements() const
ChunkSparseMatrix(const ChunkSparsityPattern &sparsity, const IdentityMatrix &id)
void vmult(OutVector &dst, const InVector &src) const
void TSOR(Vector< somenumber > &v, const number om=1.) const
void add(const number factor, const ChunkSparseMatrix< somenumber > &matrix)
size_type n_nonzero_elements() const
void copy_from(const FullMatrix< somenumber > &matrix)
virtual void clear()
void block_read(std::istream &in)
number operator()(const size_type i, const size_type j) const
void Tvmult_add(OutVector &dst, const InVector &src) const
number el(const size_type i, const size_type j) const
types::global_dof_index size_type
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
typename numbers::NumberTraits< number >::real_type real_type
ChunkSparseMatrix(const ChunkSparseMatrix &)
number diag_element(const size_type i) const
void SOR(Vector< somenumber > &v, const number om=1.) const
iterator begin(const unsigned int r)
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
std::size_t memory_consumption() const
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
real_type frobenius_norm() const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
size_type n() const
ChunkSparseMatrix & operator=(const double d)
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
const_iterator end(const unsigned int r) const
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
static const size_type invalid_entry
static constexpr size_type invalid_entry
Definition: vector.h:109
std::enable_if_t< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define AssertIsFinite(number)
Definition: exceptions.h:1786
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcDifferentChunkSparsityPatterns()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
constexpr int chunk_size
Definition: cuda_size.h:35
Expression operator>(const Expression &lhs, const Expression &rhs)
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
void copy(const T *begin, const T *end, U *dest)
unsigned int global_dof_index
Definition: types.h:81
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
static const bool zero_addition_can_be_elided
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::value_type value_type
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)