Reference documentation for deal.II version Git e8f99f7e31 2020-07-08 12:50:43 +0200
\(\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 <memory>
30 
31 
33 
34 // Forward declarations
35 # ifndef DOXYGEN
36 template <typename number>
37 class Vector;
38 template <typename number>
39 class FullMatrix;
40 # endif
41 
51 {
52  // forward declaration
53  template <typename number, bool Constness>
54  class Iterator;
55 
66  template <typename number, bool Constness>
68  {
69  public:
73  number
74  value() const;
75 
79  number &
80  value();
81 
87  get_matrix() const;
88  };
89 
90 
91 
98  template <typename number>
100  {
101  public:
107 
111  Accessor(MatrixType *matrix, const unsigned int row);
112 
117 
122 
126  number
127  value() const;
128 
133  const MatrixType &
134  get_matrix() const;
135 
136  private:
141 
146 
147  // Make iterator class a friend.
148  template <typename, bool>
149  friend class Iterator;
150  };
151 
152 
159  template <typename number>
160  class Accessor<number, false> : public ChunkSparsityPatternIterators::Accessor
161  {
162  private:
185  class Reference
186  {
187  public:
192  Reference(const Accessor *accessor, const bool dummy);
193 
197  operator number() const;
198 
202  const Reference &
203  operator=(const number n) const;
204 
208  const Reference &
209  operator+=(const number n) const;
210 
214  const Reference &
215  operator-=(const number n) const;
216 
220  const Reference &
221  operator*=(const number n) const;
222 
226  const Reference &
227  operator/=(const number n) const;
228 
229  private:
235  };
236 
237  public:
243 
247  Accessor(MatrixType *matrix, const unsigned int row);
248 
253 
257  Reference
258  value() const;
259 
264  MatrixType &
265  get_matrix() const;
266 
267  private:
272 
277 
278  // Make iterator class a friend.
279  template <typename, bool>
280  friend class Iterator;
281  };
282 
283 
284 
295  template <typename number, bool Constness>
296  class Iterator
297  {
298  public:
303 
309 
314  Iterator(MatrixType *matrix, const unsigned int row);
315 
320 
326 
330  Iterator &
331  operator++();
332 
336  Iterator
337  operator++(int);
338 
342  const Accessor<number, Constness> &operator*() const;
343 
347  const Accessor<number, Constness> *operator->() const;
348 
352  bool
353  operator==(const Iterator &) const;
354 
358  bool
359  operator!=(const Iterator &) const;
360 
368  bool
369  operator<(const Iterator &) const;
370 
375  bool
376  operator>(const Iterator &) const;
377 
384  int
385  operator-(const Iterator &p) const;
386 
390  Iterator
391  operator+(const unsigned int n) const;
392 
393  private:
398  };
399 
400 } // namespace ChunkSparseMatrixIterators
401 
402 
403 
420 template <typename number>
421 class ChunkSparseMatrix : public virtual Subscriptor
422 {
423 public:
428 
433  using value_type = number;
434 
445 
451 
459 
466  struct Traits
467  {
472  static const bool zero_addition_can_be_elided = true;
473  };
474 
490 
500 
514  explicit ChunkSparseMatrix(const ChunkSparsityPattern &sparsity);
515 
522  ChunkSparseMatrix(const ChunkSparsityPattern &sparsity,
523  const IdentityMatrix & id);
524 
529  virtual ~ChunkSparseMatrix() override;
530 
541  operator=(const ChunkSparseMatrix<number> &);
542 
550  operator=(const IdentityMatrix &id);
551 
562  operator=(const double d);
563 
577  virtual void
578  reinit(const ChunkSparsityPattern &sparsity);
579 
585  virtual void
586  clear();
588 
596  bool
597  empty() const;
598 
603  size_type
604  m() const;
605 
610  size_type
611  n() const;
612 
618  size_type
619  n_nonzero_elements() const;
620 
628  size_type
629  n_actually_nonzero_elements() const;
630 
639  const ChunkSparsityPattern &
640  get_sparsity_pattern() const;
641 
646  std::size_t
647  memory_consumption() const;
648 
650 
659  void
660  set(const size_type i, const size_type j, const number value);
661 
667  void
668  add(const size_type i, const size_type j, const number value);
669 
679  template <typename number2>
680  void
681  add(const size_type row,
682  const size_type n_cols,
683  const size_type *col_indices,
684  const number2 * values,
685  const bool elide_zero_values = true,
686  const bool col_indices_are_sorted = false);
687 
692  operator*=(const number factor);
693 
698  operator/=(const number factor);
699 
712  void
713  symmetrize();
714 
731  template <typename somenumber>
733  copy_from(const ChunkSparseMatrix<somenumber> &source);
734 
752  template <typename ForwardIterator>
753  void
754  copy_from(const ForwardIterator begin, const ForwardIterator end);
755 
761  template <typename somenumber>
762  void
763  copy_from(const FullMatrix<somenumber> &matrix);
764 
776  template <typename somenumber>
777  void
778  add(const number factor, const ChunkSparseMatrix<somenumber> &matrix);
779 
781 
785 
799  number
800  operator()(const size_type i, const size_type j) const;
801 
814  number
815  el(const size_type i, const size_type j) const;
816 
826  number
827  diag_element(const size_type i) const;
828 
838  void
839  extract_row_copy(const size_type row,
840  const size_type array_length,
841  size_type & row_length,
842  size_type * column_indices,
843  number * values) const;
844 
846 
864  template <class OutVector, class InVector>
865  void
866  vmult(OutVector &dst, const InVector &src) const;
867 
883  template <class OutVector, class InVector>
884  void
885  Tvmult(OutVector &dst, const InVector &src) const;
886 
901  template <class OutVector, class InVector>
902  void
903  vmult_add(OutVector &dst, const InVector &src) const;
904 
920  template <class OutVector, class InVector>
921  void
922  Tvmult_add(OutVector &dst, const InVector &src) const;
923 
939  template <typename somenumber>
940  somenumber
941  matrix_norm_square(const Vector<somenumber> &v) const;
942 
946  template <typename somenumber>
947  somenumber
948  matrix_scalar_product(const Vector<somenumber> &u,
949  const Vector<somenumber> &v) const;
957  template <typename somenumber>
958  somenumber
959  residual(Vector<somenumber> & dst,
960  const Vector<somenumber> &x,
961  const Vector<somenumber> &b) const;
962 
964 
968 
976  real_type
977  l1_norm() const;
978 
986  real_type
987  linfty_norm() const;
988 
993  real_type
994  frobenius_norm() const;
996 
1000 
1006  template <typename somenumber>
1007  void
1008  precondition_Jacobi(Vector<somenumber> & dst,
1009  const Vector<somenumber> &src,
1010  const number omega = 1.) const;
1011 
1015  template <typename somenumber>
1016  void
1017  precondition_SSOR(Vector<somenumber> & dst,
1018  const Vector<somenumber> &src,
1019  const number om = 1.) const;
1020 
1024  template <typename somenumber>
1025  void
1026  precondition_SOR(Vector<somenumber> & dst,
1027  const Vector<somenumber> &src,
1028  const number om = 1.) const;
1029 
1033  template <typename somenumber>
1034  void
1035  precondition_TSOR(Vector<somenumber> & dst,
1036  const Vector<somenumber> &src,
1037  const number om = 1.) const;
1038 
1044  template <typename somenumber>
1045  void
1046  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1047 
1052  template <typename somenumber>
1053  void
1054  SOR(Vector<somenumber> &v, const number om = 1.) const;
1055 
1060  template <typename somenumber>
1061  void
1062  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1063 
1074  template <typename somenumber>
1075  void
1076  PSOR(Vector<somenumber> & v,
1077  const std::vector<size_type> &permutation,
1078  const std::vector<size_type> &inverse_permutation,
1079  const number om = 1.) const;
1080 
1091  template <typename somenumber>
1092  void
1093  TPSOR(Vector<somenumber> & v,
1094  const std::vector<size_type> &permutation,
1095  const std::vector<size_type> &inverse_permutation,
1096  const number om = 1.) const;
1097 
1102  template <typename somenumber>
1103  void
1104  SOR_step(Vector<somenumber> & v,
1105  const Vector<somenumber> &b,
1106  const number om = 1.) const;
1107 
1112  template <typename somenumber>
1113  void
1114  TSOR_step(Vector<somenumber> & v,
1115  const Vector<somenumber> &b,
1116  const number om = 1.) const;
1117 
1122  template <typename somenumber>
1123  void
1124  SSOR_step(Vector<somenumber> & v,
1125  const Vector<somenumber> &b,
1126  const number om = 1.) const;
1128 
1132 
1143  begin() const;
1144 
1154  end() const;
1155 
1165  iterator
1166  begin();
1167 
1176  iterator
1177  end();
1178 
1194  begin(const unsigned int r) const;
1195 
1211  end(const unsigned int r) const;
1212 
1227  iterator
1228  begin(const unsigned int r);
1229 
1244  iterator
1245  end(const unsigned int r);
1247 
1251 
1256  void
1257  print(std::ostream &out) const;
1258 
1279  void
1280  print_formatted(std::ostream & out,
1281  const unsigned int precision = 3,
1282  const bool scientific = true,
1283  const unsigned int width = 0,
1284  const char * zero_string = " ",
1285  const double denominator = 1.) const;
1286 
1292  void
1293  print_pattern(std::ostream &out, const double threshold = 0.) const;
1294 
1305  void
1306  block_write(std::ostream &out) const;
1307 
1324  void
1325  block_read(std::istream &in);
1327 
1335  DeclException2(ExcInvalidIndex,
1336  int,
1337  int,
1338  << "You are trying to access the matrix entry with index <"
1339  << arg1 << ',' << arg2
1340  << ">, but this entry does not exist in the sparsity pattern "
1341  "of this matrix."
1342  "\n\n"
1343  "The most common cause for this problem is that you used "
1344  "a method to build the sparsity pattern that did not "
1345  "(completely) take into account all of the entries you "
1346  "will later try to write into. An example would be "
1347  "building a sparsity pattern that does not include "
1348  "the entries you will write into due to constraints "
1349  "on degrees of freedom such as hanging nodes or periodic "
1350  "boundary conditions. In such cases, building the "
1351  "sparsity pattern will succeed, but you will get errors "
1352  "such as the current one at one point or other when "
1353  "trying to write into the entries of the matrix.");
1357  DeclException0(ExcDifferentChunkSparsityPatterns);
1361  DeclException2(ExcIteratorRange,
1362  int,
1363  int,
1364  << "The iterators denote a range of " << arg1
1365  << " elements, but the given number of rows was " << arg2);
1369  DeclExceptionMsg(ExcSourceEqualsDestination,
1370  "You are attempting an operation on two matrices that "
1371  "are the same object, but the operation requires that the "
1372  "two objects are in fact different.");
1374 private:
1381 
1389  std::unique_ptr<number[]> val;
1390 
1398 
1402  size_type
1403  compute_location(const size_type i, const size_type j) const;
1404 
1405  // make all other sparse matrices friends
1406  template <typename somenumber>
1407  friend class ChunkSparseMatrix;
1408 
1409  // Also give access to internal details to the iterator/accessor classes.
1410  template <typename, bool>
1412  template <typename, bool>
1414 };
1415 
1418 # ifndef DOXYGEN
1419 /*---------------------- Inline functions -----------------------------------*/
1420 
1421 
1422 
1423 template <typename number>
1426 {
1427  Assert(cols != nullptr, ExcNotInitialized());
1428  return cols->rows;
1429 }
1430 
1431 
1432 template <typename number>
1435 {
1436  Assert(cols != nullptr, ExcNotInitialized());
1437  return cols->cols;
1438 }
1439 
1440 
1441 
1442 template <typename number>
1443 inline const ChunkSparsityPattern &
1445 {
1446  Assert(cols != nullptr, ExcNotInitialized());
1447  return *cols;
1448 }
1449 
1450 
1451 
1452 template <typename number>
1455  const size_type j) const
1456 {
1457  const size_type chunk_size = cols->get_chunk_size();
1458  const size_type chunk_index =
1459  cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1460 
1461  if (chunk_index == ChunkSparsityPattern::invalid_entry)
1463  else
1464  {
1465  return (chunk_index * chunk_size * chunk_size +
1466  (i % chunk_size) * chunk_size + (j % chunk_size));
1467  }
1468 }
1469 
1470 
1471 template <typename number>
1472 inline void
1474  const size_type j,
1475  const number value)
1476 {
1477  AssertIsFinite(value);
1478 
1479  Assert(cols != nullptr, ExcNotInitialized());
1480  // it is allowed to set elements of the matrix that are not part of the
1481  // sparsity pattern, if the value to which we set it is zero
1482  const size_type index = compute_location(i, j);
1483  Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1484  ExcInvalidIndex(i, j));
1485 
1486  if (index != SparsityPattern::invalid_entry)
1487  val[index] = value;
1488 }
1489 
1490 
1491 
1492 template <typename number>
1493 inline void
1495  const size_type j,
1496  const number value)
1497 {
1498  AssertIsFinite(value);
1499 
1500  Assert(cols != nullptr, ExcNotInitialized());
1501 
1502  if (std::abs(value) != 0.)
1503  {
1504  const size_type index = compute_location(i, j);
1506  ExcInvalidIndex(i, j));
1507 
1508  val[index] += value;
1509  }
1510 }
1511 
1512 
1513 
1514 template <typename number>
1515 template <typename number2>
1516 inline void
1518  const size_type n_cols,
1519  const size_type *col_indices,
1520  const number2 * values,
1521  const bool /*elide_zero_values*/,
1522  const bool /*col_indices_are_sorted*/)
1523 {
1524  // TODO: could be done more efficiently...
1525  for (size_type col = 0; col < n_cols; ++col)
1526  add(row, col_indices[col], static_cast<number>(values[col]));
1527 }
1528 
1529 
1530 
1531 template <typename number>
1533 ChunkSparseMatrix<number>::operator*=(const number factor)
1534 {
1535  Assert(cols != nullptr, ExcNotInitialized());
1536  Assert(val != nullptr, ExcNotInitialized());
1537 
1538  const size_type chunk_size = cols->get_chunk_size();
1539 
1540  // multiply all elements of the matrix with the given factor. this includes
1541  // the padding elements in chunks that overlap the boundaries of the actual
1542  // matrix -- but since multiplication with a number does not violate the
1543  // invariant of keeping these elements at zero nothing can happen
1544  number * val_ptr = val.get();
1545  const number *const end_ptr =
1546  val.get() +
1547  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1548  while (val_ptr != end_ptr)
1549  *val_ptr++ *= factor;
1550 
1551  return *this;
1552 }
1553 
1554 
1555 
1556 template <typename number>
1558 ChunkSparseMatrix<number>::operator/=(const number factor)
1559 {
1560  Assert(cols != nullptr, ExcNotInitialized());
1561  Assert(val != nullptr, ExcNotInitialized());
1562  Assert(std::abs(factor) != 0, ExcDivideByZero());
1563 
1564  const number factor_inv = 1. / factor;
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 
1577  while (val_ptr != end_ptr)
1578  *val_ptr++ *= factor_inv;
1579 
1580  return *this;
1581 }
1582 
1583 
1584 
1585 template <typename number>
1586 inline number
1588  const size_type j) const
1589 {
1590  Assert(cols != nullptr, ExcNotInitialized());
1591  AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1592  ExcInvalidIndex(i, j));
1593  return val[compute_location(i, j)];
1594 }
1595 
1596 
1597 
1598 template <typename number>
1599 inline number
1600 ChunkSparseMatrix<number>::el(const size_type i, const size_type j) const
1601 {
1602  Assert(cols != nullptr, ExcNotInitialized());
1603  const size_type index = compute_location(i, j);
1604 
1606  return val[index];
1607  else
1608  return 0;
1609 }
1610 
1611 
1612 
1613 template <typename number>
1614 inline number
1616 {
1617  Assert(cols != nullptr, ExcNotInitialized());
1618  Assert(m() == n(), ExcNotQuadratic());
1619  AssertIndexRange(i, m());
1620 
1621  // Use that the first element in each row of a quadratic matrix is the main
1622  // diagonal of the chunk sparsity pattern
1623  const size_type chunk_size = cols->get_chunk_size();
1624  return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1625  chunk_size +
1626  (i % chunk_size) * chunk_size + (i % chunk_size)];
1627 }
1628 
1629 
1630 
1631 template <typename number>
1632 template <typename ForwardIterator>
1633 inline void
1634 ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1635  const ForwardIterator end)
1636 {
1637  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1638  ExcIteratorRange(std::distance(begin, end), m()));
1639 
1640  // for use in the inner loop, we define an alias to the type of the inner
1641  // iterators
1642  using inner_iterator =
1643  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1644  size_type row = 0;
1645  for (ForwardIterator i = begin; i != end; ++i, ++row)
1646  {
1647  const inner_iterator end_of_row = i->end();
1648  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1649  // write entries
1650  set(row, j->first, j->second);
1651  }
1652 }
1653 
1654 
1655 
1656 //---------------------------------------------------------------------------
1657 
1658 
1660 {
1661  template <typename number>
1662  inline Accessor<number, true>::Accessor(const MatrixType * matrix,
1663  const unsigned int row)
1664  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1665  row)
1666  , matrix(matrix)
1667  {}
1668 
1669 
1670 
1671  template <typename number>
1672  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1673  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1674  , matrix(matrix)
1675  {}
1676 
1677 
1678 
1679  template <typename number>
1680  inline Accessor<number, true>::Accessor(
1682  : ChunkSparsityPatternIterators::Accessor(a)
1683  , matrix(&a.get_matrix())
1684  {}
1685 
1686 
1687 
1688  template <typename number>
1689  inline number
1690  Accessor<number, true>::value() const
1691  {
1692  const unsigned int chunk_size =
1693  matrix->get_sparsity_pattern().get_chunk_size();
1694  return matrix->val[reduced_index() * chunk_size * chunk_size +
1695  chunk_row * chunk_size + chunk_col];
1696  }
1697 
1698 
1699 
1700  template <typename number>
1701  inline const typename Accessor<number, true>::MatrixType &
1702  Accessor<number, true>::get_matrix() const
1703  {
1704  return *matrix;
1705  }
1706 
1707 
1708 
1709  template <typename number>
1710  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1711  const bool)
1712  : accessor(accessor)
1713  {}
1714 
1715 
1716  template <typename number>
1717  inline Accessor<number, false>::Reference::operator number() const
1718  {
1719  const unsigned int chunk_size =
1720  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1721  return accessor->matrix
1722  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1723  accessor->chunk_row * chunk_size + accessor->chunk_col];
1724  }
1725 
1726 
1727 
1728  template <typename number>
1729  inline const typename Accessor<number, false>::Reference &
1730  Accessor<number, false>::Reference::operator=(const number n) const
1731  {
1732  const unsigned int chunk_size =
1733  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1734  accessor->matrix
1735  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1736  accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1737  return *this;
1738  }
1739 
1740 
1741 
1742  template <typename number>
1743  inline const typename Accessor<number, false>::Reference &
1744  Accessor<number, false>::Reference::operator+=(const number n) const
1745  {
1746  const unsigned int chunk_size =
1747  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1748  accessor->matrix
1749  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1750  accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1751  return *this;
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 Accessor<number, false>::Accessor(MatrixType * matrix,
1800  const unsigned int row)
1801  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1802  row)
1803  , matrix(matrix)
1804  {}
1805 
1806 
1807 
1808  template <typename number>
1809  inline Accessor<number, false>::Accessor(MatrixType *matrix)
1810  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1811  , matrix(matrix)
1812  {}
1813 
1814 
1815 
1816  template <typename number>
1817  inline typename Accessor<number, false>::Reference
1818  Accessor<number, false>::value() const
1819  {
1820  return Reference(this, true);
1821  }
1822 
1823 
1824 
1825  template <typename number>
1826  inline typename Accessor<number, false>::MatrixType &
1827  Accessor<number, false>::get_matrix() const
1828  {
1829  return *matrix;
1830  }
1831 
1832 
1833 
1834  template <typename number, bool Constness>
1835  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
1836  const unsigned int row)
1837  : accessor(matrix, row)
1838  {}
1839 
1840 
1841 
1842  template <typename number, bool Constness>
1843  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1844  : accessor(matrix)
1845  {}
1846 
1847 
1848 
1849  template <typename number, bool Constness>
1850  inline Iterator<number, Constness>::Iterator(
1852  : accessor(*i)
1853  {}
1854 
1855 
1856 
1857  template <typename number, bool Constness>
1858  inline Iterator<number, Constness> &
1860  {
1861  accessor.advance();
1862  return *this;
1863  }
1864 
1865 
1866  template <typename number, bool Constness>
1867  inline Iterator<number, Constness>
1869  {
1870  const Iterator iter = *this;
1871  accessor.advance();
1872  return iter;
1873  }
1874 
1875 
1876  template <typename number, bool Constness>
1877  inline const Accessor<number, Constness> &Iterator<number, Constness>::
1878  operator*() const
1879  {
1880  return accessor;
1881  }
1882 
1883 
1884  template <typename number, bool Constness>
1885  inline const Accessor<number, Constness> *Iterator<number, Constness>::
1886  operator->() const
1887  {
1888  return &accessor;
1889  }
1890 
1891 
1892  template <typename number, bool Constness>
1893  inline bool
1894  Iterator<number, Constness>::operator==(const Iterator &other) const
1895  {
1896  return (accessor == other.accessor);
1897  }
1898 
1899 
1900  template <typename number, bool Constness>
1901  inline bool
1902  Iterator<number, Constness>::operator!=(const Iterator &other) const
1903  {
1904  return !(*this == other);
1905  }
1906 
1907 
1908  template <typename number, bool Constness>
1909  inline bool
1910  Iterator<number, Constness>::operator<(const Iterator &other) const
1911  {
1912  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1913  ExcInternalError());
1914 
1915  return (accessor < other.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 (other < *this);
1924  }
1925 
1926 
1927  template <typename number, bool Constness>
1928  inline int
1929  Iterator<number, Constness>::operator-(const Iterator &other) const
1930  {
1931  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1932  ExcInternalError());
1933 
1934  // TODO: can be optimized
1935  int difference = 0;
1936  if (*this < other)
1937  {
1938  Iterator copy = *this;
1939  while (copy != other)
1940  {
1941  ++copy;
1942  --difference;
1943  }
1944  }
1945  else
1946  {
1947  Iterator copy = other;
1948  while (copy != *this)
1949  {
1950  ++copy;
1951  ++difference;
1952  }
1953  }
1954  return difference;
1955  }
1956 
1957 
1958 
1959  template <typename number, bool Constness>
1960  inline Iterator<number, Constness>
1961  Iterator<number, Constness>::operator+(const unsigned int n) const
1962  {
1963  Iterator x = *this;
1964  for (unsigned int i = 0; i < n; ++i)
1965  ++x;
1966 
1967  return x;
1968  }
1969 
1970 } // namespace ChunkSparseMatrixIterators
1971 
1972 
1973 
1974 template <typename number>
1977 {
1978  return const_iterator(this, 0);
1979 }
1980 
1981 
1982 template <typename number>
1985 {
1986  return const_iterator(this);
1987 }
1988 
1989 
1990 template <typename number>
1993 {
1994  return iterator(this, 0);
1995 }
1996 
1997 
1998 template <typename number>
2001 {
2002  return iterator(this);
2003 }
2004 
2005 
2006 template <typename number>
2008 ChunkSparseMatrix<number>::begin(const unsigned int r) const
2009 {
2010  AssertIndexRange(r, m());
2011  return const_iterator(this, r);
2012 }
2013 
2014 
2015 
2016 template <typename number>
2018 ChunkSparseMatrix<number>::end(const unsigned int r) const
2019 {
2020  AssertIndexRange(r, m());
2021  return const_iterator(this, r + 1);
2022 }
2023 
2024 
2025 
2026 template <typename number>
2028 ChunkSparseMatrix<number>::begin(const unsigned int r)
2029 {
2030  AssertIndexRange(r, m());
2031  return iterator(this, r);
2032 }
2033 
2034 
2035 
2036 template <typename number>
2038 ChunkSparseMatrix<number>::end(const unsigned int r)
2039 {
2040  AssertIndexRange(r, m());
2041  return iterator(this, r + 1);
2042 }
2043 
2044 
2045 
2046 # endif // DOXYGEN
2047 
2049 
2050 #endif
2051 /*--------------------------- chunk_sparse_matrix.h -------------------------*/
typename numbers::NumberTraits< number >::real_type real_type
size_type n() const
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
Contents is actually a matrix.
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
size_type m() const
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
std::unique_ptr< number[]> val
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcDivideByZero()
const ChunkSparseMatrix< number > & get_matrix() const
ChunkSparseMatrixIterators::Iterator< number, false > iterator
Accessor< number, Constness > accessor
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
const_iterator begin() const
number operator()(const size_type i, const size_type j) const
#define Assert(cond, exc)
Definition: exceptions.h:1403
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
ChunkSparseMatrix & operator*=(const number factor)
size_type compute_location(const size_type i, const size_type j) const
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2779
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
VectorType::value_type * end(VectorType &V)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
ChunkSparseMatrixIterators::Iterator< number, true > const_iterator
Expression operator>(const Expression &lhs, const Expression &rhs)
number el(const size_type i, const size_type j) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
static ::ExceptionBase & ExcNotQuadratic()
number diag_element(const size_type i) const
unsigned int global_dof_index
Definition: types.h:76
const_iterator end() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
VectorType::value_type * begin(VectorType &V)
constexpr int chunk_size
Definition: cuda_size.h:35
const ChunkSparsityPattern & get_sparsity_pattern() const
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2753
ChunkSparseMatrix & operator/=(const number factor)
const Accessor< number, Constness > & value_type
void set(const size_type i, const size_type j, const number value)
static const size_type invalid_entry
static const size_type invalid_entry
bool operator<(const Accessor &) const
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)
void add(const size_type i, const size_type j, const number value)
typename Accessor< number, Constness >::MatrixType MatrixType
void copy(const T *begin, const T *end, U *dest)
#define AssertIsFinite(number)
Definition: exceptions.h:1659
bool operator==(const Accessor &) const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()