Reference documentation for deal.II version GIT 821ac203c4 2022-05-20 14:50: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 - 2020 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 
52 {
53  // forward declaration
54  template <typename number, bool Constness>
55  class Iterator;
56 
67  template <typename number, bool Constness>
69  {
70  public:
74  number
75  value() const;
76 
80  number &
81  value();
82 
88  get_matrix() const;
89  };
90 
91 
92 
99  template <typename number>
101  {
102  public:
108 
112  Accessor(MatrixType *matrix, const unsigned int row);
113 
118 
123 
127  number
128  value() const;
129 
134  const MatrixType &
135  get_matrix() const;
136 
137  private:
142 
147 
148  // Make iterator class a friend.
149  template <typename, bool>
150  friend class Iterator;
151  };
152 
153 
160  template <typename number>
161  class Accessor<number, false> : public ChunkSparsityPatternIterators::Accessor
162  {
163  private:
186  class Reference
187  {
188  public:
193  Reference(const Accessor *accessor, const bool dummy);
194 
198  operator number() const;
199 
203  const Reference &
204  operator=(const number n) const;
205 
209  const Reference &
210  operator+=(const number n) const;
211 
215  const Reference &
216  operator-=(const number n) const;
217 
221  const Reference &
222  operator*=(const number n) const;
223 
227  const Reference &
228  operator/=(const number n) const;
229 
230  private:
236  };
237 
238  public:
244 
248  Accessor(MatrixType *matrix, const unsigned int row);
249 
254 
258  Reference
259  value() const;
260 
265  MatrixType &
266  get_matrix() const;
267 
268  private:
273 
278 
279  // Make iterator class a friend.
280  template <typename, bool>
281  friend class Iterator;
282  };
283 
284 
285 
296  template <typename number, bool Constness>
297  class Iterator
298  {
299  public:
304 
310 
316 
321  Iterator(MatrixType *matrix, const unsigned int row);
322 
327 
333 
337  Iterator &
339 
343  Iterator
345 
350  operator*() const;
351 
356  operator->() const;
357 
361  bool
362  operator==(const Iterator &) const;
363 
367  bool
368  operator!=(const Iterator &) const;
369 
377  bool
378  operator<(const Iterator &) const;
379 
384  bool
385  operator>(const Iterator &) const;
386 
393  int
394  operator-(const Iterator &p) const;
395 
399  Iterator
400  operator+(const unsigned int n) const;
401 
402  private:
407  };
408 
409 } // namespace ChunkSparseMatrixIterators
410 
412 
413 namespace std
414 {
415  template <typename number, bool Constness>
416  struct iterator_traits<
417  ::ChunkSparseMatrixIterators::Iterator<number, Constness>>
418  {
419  using iterator_category = forward_iterator_tag;
420  using value_type = typename ::ChunkSparseMatrixIterators::
421  Iterator<number, Constness>::value_type;
422  using difference_type = typename ::ChunkSparseMatrixIterators::
423  Iterator<number, Constness>::difference_type;
424  };
425 } // namespace std
426 
428 
429 
430 
447 template <typename number>
448 class ChunkSparseMatrix : public virtual Subscriptor
449 {
450 public:
455 
460  using value_type = number;
461 
472 
478 
486 
493  struct Traits
494  {
499  static const bool zero_addition_can_be_elided = true;
500  };
501 
517 
527 
541  explicit ChunkSparseMatrix(const ChunkSparsityPattern &sparsity);
542 
550  const IdentityMatrix & id);
551 
556  virtual ~ChunkSparseMatrix() override;
557 
569 
578 
589  operator=(const double d);
590 
604  virtual void
605  reinit(const ChunkSparsityPattern &sparsity);
606 
612  virtual void
613  clear();
615 
623  bool
624  empty() const;
625 
630  size_type
631  m() const;
632 
637  size_type
638  n() const;
639 
645  size_type
647 
655  size_type
657 
666  const ChunkSparsityPattern &
668 
673  std::size_t
675 
677 
686  void
687  set(const size_type i, const size_type j, const number value);
688 
694  void
695  add(const size_type i, const size_type j, const number value);
696 
706  template <typename number2>
707  void
708  add(const size_type row,
709  const size_type n_cols,
710  const size_type *col_indices,
711  const number2 * values,
712  const bool elide_zero_values = true,
713  const bool col_indices_are_sorted = false);
714 
719  operator*=(const number factor);
720 
725  operator/=(const number factor);
726 
739  void
741 
758  template <typename somenumber>
761 
779  template <typename ForwardIterator>
780  void
781  copy_from(const ForwardIterator begin, const ForwardIterator end);
782 
788  template <typename somenumber>
789  void
791 
803  template <typename somenumber>
804  void
805  add(const number factor, const ChunkSparseMatrix<somenumber> &matrix);
806 
808 
812 
826  number
827  operator()(const size_type i, const size_type j) const;
828 
841  number
842  el(const size_type i, const size_type j) const;
843 
853  number
854  diag_element(const size_type i) const;
855 
865  void
867  const size_type array_length,
868  size_type & row_length,
869  size_type * column_indices,
870  number * values) const;
871 
873 
891  template <class OutVector, class InVector>
892  void
893  vmult(OutVector &dst, const InVector &src) const;
894 
910  template <class OutVector, class InVector>
911  void
912  Tvmult(OutVector &dst, const InVector &src) const;
913 
928  template <class OutVector, class InVector>
929  void
930  vmult_add(OutVector &dst, const InVector &src) const;
931 
947  template <class OutVector, class InVector>
948  void
949  Tvmult_add(OutVector &dst, const InVector &src) const;
950 
966  template <typename somenumber>
967  somenumber
969 
973  template <typename somenumber>
974  somenumber
976  const Vector<somenumber> &v) const;
984  template <typename somenumber>
985  somenumber
987  const Vector<somenumber> &x,
988  const Vector<somenumber> &b) const;
989 
991 
995 
1003  real_type
1004  l1_norm() const;
1005 
1013  real_type
1014  linfty_norm() const;
1015 
1020  real_type
1023 
1027 
1033  template <typename somenumber>
1034  void
1036  const Vector<somenumber> &src,
1037  const number omega = 1.) const;
1038 
1042  template <typename somenumber>
1043  void
1045  const Vector<somenumber> &src,
1046  const number om = 1.) const;
1047 
1051  template <typename somenumber>
1052  void
1054  const Vector<somenumber> &src,
1055  const number om = 1.) const;
1056 
1060  template <typename somenumber>
1061  void
1063  const Vector<somenumber> &src,
1064  const number om = 1.) const;
1065 
1071  template <typename somenumber>
1072  void
1073  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1074 
1079  template <typename somenumber>
1080  void
1081  SOR(Vector<somenumber> &v, const number om = 1.) const;
1082 
1087  template <typename somenumber>
1088  void
1089  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1090 
1101  template <typename somenumber>
1102  void
1104  const std::vector<size_type> &permutation,
1105  const std::vector<size_type> &inverse_permutation,
1106  const number om = 1.) const;
1107 
1118  template <typename somenumber>
1119  void
1121  const std::vector<size_type> &permutation,
1122  const std::vector<size_type> &inverse_permutation,
1123  const number om = 1.) const;
1124 
1129  template <typename somenumber>
1130  void
1132  const Vector<somenumber> &b,
1133  const number om = 1.) const;
1134 
1139  template <typename somenumber>
1140  void
1142  const Vector<somenumber> &b,
1143  const number om = 1.) const;
1144 
1149  template <typename somenumber>
1150  void
1152  const Vector<somenumber> &b,
1153  const number om = 1.) const;
1155 
1159 
1170  begin() const;
1171 
1181  end() const;
1182 
1192  iterator
1194 
1203  iterator
1204  end();
1205 
1221  begin(const unsigned int r) const;
1222 
1238  end(const unsigned int r) const;
1239 
1254  iterator
1255  begin(const unsigned int r);
1256 
1271  iterator
1272  end(const unsigned int r);
1274 
1278 
1283  void
1284  print(std::ostream &out) const;
1285 
1306  void
1307  print_formatted(std::ostream & out,
1308  const unsigned int precision = 3,
1309  const bool scientific = true,
1310  const unsigned int width = 0,
1311  const char * zero_string = " ",
1312  const double denominator = 1.) const;
1313 
1319  void
1320  print_pattern(std::ostream &out, const double threshold = 0.) const;
1321 
1332  void
1333  block_write(std::ostream &out) const;
1334 
1351  void
1352  block_read(std::istream &in);
1354 
1363  int,
1364  int,
1365  << "You are trying to access the matrix entry with index <"
1366  << arg1 << ',' << arg2
1367  << ">, but this entry does not exist in the sparsity pattern "
1368  "of this matrix."
1369  "\n\n"
1370  "The most common cause for this problem is that you used "
1371  "a method to build the sparsity pattern that did not "
1372  "(completely) take into account all of the entries you "
1373  "will later try to write into. An example would be "
1374  "building a sparsity pattern that does not include "
1375  "the entries you will write into due to constraints "
1376  "on degrees of freedom such as hanging nodes or periodic "
1377  "boundary conditions. In such cases, building the "
1378  "sparsity pattern will succeed, but you will get errors "
1379  "such as the current one at one point or other when "
1380  "trying to write into the entries of the matrix.");
1389  int,
1390  int,
1391  << "The iterators denote a range of " << arg1
1392  << " elements, but the given number of rows was " << arg2);
1397  "You are attempting an operation on two matrices that "
1398  "are the same object, but the operation requires that the "
1399  "two objects are in fact different.");
1401 private:
1408 
1416  std::unique_ptr<number[]> val;
1417 
1425 
1429  size_type
1430  compute_location(const size_type i, const size_type j) const;
1431 
1432  // make all other sparse matrices friends
1433  template <typename somenumber>
1434  friend class ChunkSparseMatrix;
1435 
1436  // Also give access to internal details to the iterator/accessor classes.
1437  template <typename, bool>
1439  template <typename, bool>
1441 };
1442 
1445 # ifndef DOXYGEN
1446 /*---------------------- Inline functions -----------------------------------*/
1447 
1448 
1449 
1450 template <typename number>
1453 {
1454  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1455  return cols->rows;
1456 }
1457 
1458 
1459 template <typename number>
1462 {
1463  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1464  return cols->cols;
1465 }
1466 
1467 
1468 
1469 template <typename number>
1470 inline const ChunkSparsityPattern &
1472 {
1473  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1474  return *cols;
1475 }
1476 
1477 
1478 
1479 template <typename number>
1482  const size_type j) const
1483 {
1484  const size_type chunk_size = cols->get_chunk_size();
1485  const size_type chunk_index =
1486  cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1487 
1488  if (chunk_index == ChunkSparsityPattern::invalid_entry)
1490  else
1491  {
1492  return (chunk_index * chunk_size * chunk_size +
1493  (i % chunk_size) * chunk_size + (j % chunk_size));
1494  }
1495 }
1496 
1497 
1498 template <typename number>
1499 inline void
1501  const size_type j,
1502  const number value)
1503 {
1505 
1506  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1507  // it is allowed to set elements of the matrix that are not part of the
1508  // sparsity pattern, if the value to which we set it is zero
1509  const size_type index = compute_location(i, j);
1510  Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1511  ExcInvalidIndex(i, j));
1512 
1513  if (index != SparsityPattern::invalid_entry)
1514  val[index] = value;
1515 }
1516 
1517 
1518 
1519 template <typename number>
1520 inline void
1522  const size_type j,
1523  const number value)
1524 {
1526 
1527  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1528 
1529  if (std::abs(value) != 0.)
1530  {
1531  const size_type index = compute_location(i, j);
1533  ExcInvalidIndex(i, j));
1534 
1535  val[index] += value;
1536  }
1537 }
1538 
1539 
1540 
1541 template <typename number>
1542 template <typename number2>
1543 inline void
1545  const size_type n_cols,
1546  const size_type *col_indices,
1547  const number2 * values,
1548  const bool /*elide_zero_values*/,
1549  const bool /*col_indices_are_sorted*/)
1550 {
1551  // TODO: could be done more efficiently...
1552  for (size_type col = 0; col < n_cols; ++col)
1553  add(row, col_indices[col], static_cast<number>(values[col]));
1554 }
1555 
1556 
1557 
1558 template <typename number>
1560 ChunkSparseMatrix<number>::operator*=(const number factor)
1561 {
1562  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1563  Assert(val != nullptr, ExcNotInitialized());
1564 
1565  const size_type chunk_size = cols->get_chunk_size();
1566 
1567  // multiply all elements of the matrix with the given factor. this includes
1568  // the padding elements in chunks that overlap the boundaries of the actual
1569  // matrix -- but since multiplication with a number does not violate the
1570  // invariant of keeping these elements at zero nothing can happen
1571  number * val_ptr = val.get();
1572  const number *const end_ptr =
1573  val.get() +
1574  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1575  while (val_ptr != end_ptr)
1576  *val_ptr++ *= factor;
1577 
1578  return *this;
1579 }
1580 
1581 
1582 
1583 template <typename number>
1585 ChunkSparseMatrix<number>::operator/=(const number factor)
1586 {
1587  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1588  Assert(val != nullptr, ExcNotInitialized());
1589  Assert(std::abs(factor) != 0, ExcDivideByZero());
1590 
1591  const number factor_inv = 1. / factor;
1592 
1593  const size_type chunk_size = cols->get_chunk_size();
1594 
1595  // multiply all elements of the matrix with the given factor. this includes
1596  // the padding elements in chunks that overlap the boundaries of the actual
1597  // matrix -- but since multiplication with a number does not violate the
1598  // invariant of keeping these elements at zero nothing can happen
1599  number * val_ptr = val.get();
1600  const number *const end_ptr =
1601  val.get() +
1602  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1603 
1604  while (val_ptr != end_ptr)
1605  *val_ptr++ *= factor_inv;
1606 
1607  return *this;
1608 }
1609 
1610 
1611 
1612 template <typename number>
1613 inline number
1615  const size_type j) const
1616 {
1617  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1618  AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1619  ExcInvalidIndex(i, j));
1620  return val[compute_location(i, j)];
1621 }
1622 
1623 
1624 
1625 template <typename number>
1626 inline number
1627 ChunkSparseMatrix<number>::el(const size_type i, const size_type j) const
1628 {
1629  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1630  const size_type index = compute_location(i, j);
1631 
1633  return val[index];
1634  else
1635  return 0;
1636 }
1637 
1638 
1639 
1640 template <typename number>
1641 inline number
1643 {
1644  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1645  Assert(m() == n(), ExcNotQuadratic());
1646  AssertIndexRange(i, m());
1647 
1648  // Use that the first element in each row of a quadratic matrix is the main
1649  // diagonal of the chunk sparsity pattern
1650  const size_type chunk_size = cols->get_chunk_size();
1651  return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1652  chunk_size +
1653  (i % chunk_size) * chunk_size + (i % chunk_size)];
1654 }
1655 
1656 
1657 
1658 template <typename number>
1659 template <typename ForwardIterator>
1660 inline void
1661 ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1662  const ForwardIterator end)
1663 {
1664  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1665  ExcIteratorRange(std::distance(begin, end), m()));
1666 
1667  // for use in the inner loop, we define an alias to the type of the inner
1668  // iterators
1669  using inner_iterator =
1670  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1671  size_type row = 0;
1672  for (ForwardIterator i = begin; i != end; ++i, ++row)
1673  {
1674  const inner_iterator end_of_row = i->end();
1675  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1676  // write entries
1677  set(row, j->first, j->second);
1678  }
1679 }
1680 
1681 
1682 
1683 //---------------------------------------------------------------------------
1684 
1685 
1687 {
1688  template <typename number>
1689  inline Accessor<number, true>::Accessor(const MatrixType * matrix,
1690  const unsigned int row)
1691  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1692  row)
1693  , matrix(matrix)
1694  {}
1695 
1696 
1697 
1698  template <typename number>
1699  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1700  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1701  , matrix(matrix)
1702  {}
1703 
1704 
1705 
1706  template <typename number>
1707  inline Accessor<number, true>::Accessor(
1709  : ChunkSparsityPatternIterators::Accessor(a)
1710  , matrix(&a.get_matrix())
1711  {}
1712 
1713 
1714 
1715  template <typename number>
1716  inline number
1717  Accessor<number, true>::value() const
1718  {
1719  const unsigned int chunk_size =
1720  matrix->get_sparsity_pattern().get_chunk_size();
1721  return matrix->val[reduced_index() * chunk_size * chunk_size +
1722  chunk_row * chunk_size + chunk_col];
1723  }
1724 
1725 
1726 
1727  template <typename number>
1728  inline const typename Accessor<number, true>::MatrixType &
1729  Accessor<number, true>::get_matrix() const
1730  {
1731  return *matrix;
1732  }
1733 
1734 
1735 
1736  template <typename number>
1737  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1738  const bool)
1739  : accessor(accessor)
1740  {}
1741 
1742 
1743  template <typename number>
1744  inline Accessor<number, false>::Reference::operator number() const
1745  {
1746  const unsigned int chunk_size =
1747  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1748  return accessor->matrix
1749  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1750  accessor->chunk_row * chunk_size + accessor->chunk_col];
1751  }
1752 
1753 
1754 
1755  template <typename number>
1756  inline const typename Accessor<number, false>::Reference &
1757  Accessor<number, false>::Reference::operator=(const number n) const
1758  {
1759  const unsigned int chunk_size =
1760  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1761  accessor->matrix
1762  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1763  accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1764  return *this;
1765  }
1766 
1767 
1768 
1769  template <typename number>
1770  inline const typename Accessor<number, false>::Reference &
1771  Accessor<number, false>::Reference::operator+=(const number n) const
1772  {
1773  const unsigned int chunk_size =
1774  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1775  accessor->matrix
1776  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1777  accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1778  return *this;
1779  }
1780 
1781 
1782 
1783  template <typename number>
1784  inline const typename Accessor<number, false>::Reference &
1785  Accessor<number, false>::Reference::operator-=(const number n) const
1786  {
1787  const unsigned int chunk_size =
1788  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1789  accessor->matrix
1790  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1791  accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1792  return *this;
1793  }
1794 
1795 
1796 
1797  template <typename number>
1798  inline const typename Accessor<number, false>::Reference &
1799  Accessor<number, false>::Reference::operator*=(const number n) const
1800  {
1801  const unsigned int chunk_size =
1802  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1803  accessor->matrix
1804  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1805  accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1806  return *this;
1807  }
1808 
1809 
1810 
1811  template <typename number>
1812  inline const typename Accessor<number, false>::Reference &
1813  Accessor<number, false>::Reference::operator/=(const number n) const
1814  {
1815  const unsigned int chunk_size =
1816  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1817  accessor->matrix
1818  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1819  accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1820  return *this;
1821  }
1822 
1823 
1824 
1825  template <typename number>
1826  inline Accessor<number, false>::Accessor(MatrixType * matrix,
1827  const unsigned int row)
1828  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1829  row)
1830  , matrix(matrix)
1831  {}
1832 
1833 
1834 
1835  template <typename number>
1836  inline Accessor<number, false>::Accessor(MatrixType *matrix)
1837  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1838  , matrix(matrix)
1839  {}
1840 
1841 
1842 
1843  template <typename number>
1844  inline typename Accessor<number, false>::Reference
1845  Accessor<number, false>::value() const
1846  {
1847  return Reference(this, true);
1848  }
1849 
1850 
1851 
1852  template <typename number>
1853  inline typename Accessor<number, false>::MatrixType &
1854  Accessor<number, false>::get_matrix() const
1855  {
1856  return *matrix;
1857  }
1858 
1859 
1860 
1861  template <typename number, bool Constness>
1862  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
1863  const unsigned int row)
1864  : accessor(matrix, row)
1865  {}
1866 
1867 
1868 
1869  template <typename number, bool Constness>
1870  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1871  : accessor(matrix)
1872  {}
1873 
1874 
1875 
1876  template <typename number, bool Constness>
1877  inline Iterator<number, Constness>::Iterator(
1879  : accessor(*i)
1880  {}
1881 
1882 
1883 
1884  template <typename number, bool Constness>
1885  inline Iterator<number, Constness> &
1887  {
1888  accessor.advance();
1889  return *this;
1890  }
1891 
1892 
1893  template <typename number, bool Constness>
1894  inline Iterator<number, Constness>
1896  {
1897  const Iterator iter = *this;
1898  accessor.advance();
1899  return iter;
1900  }
1901 
1902 
1903  template <typename number, bool Constness>
1904  inline const Accessor<number, Constness> &
1906  {
1907  return accessor;
1908  }
1909 
1910 
1911  template <typename number, bool Constness>
1912  inline const Accessor<number, Constness> *
1913  Iterator<number, Constness>::operator->() const
1914  {
1915  return &accessor;
1916  }
1917 
1918 
1919  template <typename number, bool Constness>
1920  inline bool
1921  Iterator<number, Constness>::operator==(const Iterator &other) const
1922  {
1923  return (accessor == other.accessor);
1924  }
1925 
1926 
1927  template <typename number, bool Constness>
1928  inline bool
1929  Iterator<number, Constness>::operator!=(const Iterator &other) const
1930  {
1931  return !(*this == other);
1932  }
1933 
1934 
1935  template <typename number, bool Constness>
1936  inline bool
1937  Iterator<number, Constness>::operator<(const Iterator &other) const
1938  {
1939  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1940  ExcInternalError());
1941 
1942  return (accessor < other.accessor);
1943  }
1944 
1945 
1946  template <typename number, bool Constness>
1947  inline bool
1948  Iterator<number, Constness>::operator>(const Iterator &other) const
1949  {
1950  return (other < *this);
1951  }
1952 
1953 
1954  template <typename number, bool Constness>
1955  inline int
1956  Iterator<number, Constness>::operator-(const Iterator &other) const
1957  {
1958  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1959  ExcInternalError());
1960 
1961  // TODO: can be optimized
1962  int difference = 0;
1963  if (*this < other)
1964  {
1965  Iterator copy = *this;
1966  while (copy != other)
1967  {
1968  ++copy;
1969  --difference;
1970  }
1971  }
1972  else
1973  {
1974  Iterator copy = other;
1975  while (copy != *this)
1976  {
1977  ++copy;
1978  ++difference;
1979  }
1980  }
1981  return difference;
1982  }
1983 
1984 
1985 
1986  template <typename number, bool Constness>
1987  inline Iterator<number, Constness>
1988  Iterator<number, Constness>::operator+(const unsigned int n) const
1989  {
1990  Iterator x = *this;
1991  for (unsigned int i = 0; i < n; ++i)
1992  ++x;
1993 
1994  return x;
1995  }
1996 
1997 } // namespace ChunkSparseMatrixIterators
1998 
1999 
2000 
2001 template <typename number>
2004 {
2005  return const_iterator(this, 0);
2006 }
2007 
2008 
2009 template <typename number>
2012 {
2013  return const_iterator(this);
2014 }
2015 
2016 
2017 template <typename number>
2020 {
2021  return iterator(this, 0);
2022 }
2023 
2024 
2025 template <typename number>
2028 {
2029  return iterator(this);
2030 }
2031 
2032 
2033 template <typename number>
2035 ChunkSparseMatrix<number>::begin(const unsigned int r) const
2036 {
2037  AssertIndexRange(r, m());
2038  return const_iterator(this, r);
2039 }
2040 
2041 
2042 
2043 template <typename number>
2045 ChunkSparseMatrix<number>::end(const unsigned int r) const
2046 {
2047  AssertIndexRange(r, m());
2048  return const_iterator(this, r + 1);
2049 }
2050 
2051 
2052 
2053 template <typename number>
2055 ChunkSparseMatrix<number>::begin(const unsigned int r)
2056 {
2057  AssertIndexRange(r, m());
2058  return iterator(this, r);
2059 }
2060 
2061 
2062 
2063 template <typename number>
2065 ChunkSparseMatrix<number>::end(const unsigned int r)
2066 {
2067  AssertIndexRange(r, m());
2068  return iterator(this, r + 1);
2069 }
2070 
2071 
2072 
2073 # endif // DOXYGEN
2074 
2076 
2077 #endif
2078 /*--------------------------- chunk_sparse_matrix.h -------------------------*/
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static const size_type invalid_entry
Definition: vector.h:109
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
__global__ void set(Number *val, const Number s, const size_type N)
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
#define DeclException0(Exception0)
Definition: exceptions.h:464
size_type compute_location(const size_type i, const size_type j) const
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcInternalError()
void print(std::ostream &out) const
static ::ExceptionBase & ExcNotInitialized()
void print_pattern(std::ostream &out, const double threshold=0.) const
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::unique_ptr< number[]> val
void block_write(std::ostream &out) const
#define AssertIsFinite(number)
Definition: exceptions.h:1758
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
void block_read(std::istream &in)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcSourceEqualsDestination()
Expression operator>(const Expression &lhs, const Expression &rhs)
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcDifferentChunkSparsityPatterns()
static const bool zero_addition_can_be_elided
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
bool operator>(const Iterator &) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
Iterator(const ChunkSparseMatrixIterators::Iterator< number, false > &i)
ChunkSparseMatrix & operator/=(const number factor)
ChunkSparseMatrix & operator*=(const number factor)
ChunkSparseMatrix< number > & operator=(const IdentityMatrix &id)
Accessor(MatrixType *matrix, const unsigned int row)
const Accessor< number, Constness > & value_type
real_type linfty_norm() const
const Accessor< number, Constness > * operator->() const
ChunkSparseMatrix(const ChunkSparsityPattern &sparsity)
bool operator<(const Iterator &) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
Iterator(MatrixType *matrix, const unsigned int row)
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
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 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
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
bool operator==(const Iterator &) 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
Accessor< number, Constness > accessor
somenumber matrix_norm_square(const Vector< somenumber > &v) const
virtual ~ChunkSparseMatrix() override
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
int operator-(const Iterator &p) 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()
number operator()(const size_type i, const size_type j) const
void Tvmult_add(OutVector &dst, const InVector &src) const
Reference(const Accessor *accessor, const bool dummy)
number el(const size_type i, const size_type j) const
Accessor(const ChunkSparseMatrixIterators::Accessor< number, false > &a)
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::value_type value_type
types::global_dof_index size_type
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
typename numbers::NumberTraits< number >::real_type real_type
const ChunkSparseMatrix< number > & get_matrix() const
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
const Accessor< number, Constness > & operator*() const
typename Accessor< number, Constness >::MatrixType MatrixType
Accessor(MatrixType *matrix, const unsigned int row)
std::size_t memory_consumption() const
bool operator!=(const Iterator &) 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)
Iterator operator+(const unsigned int n) const
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
constexpr int chunk_size
Definition: cuda_size.h:35
@ 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:76
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)