Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
chunk_sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2018 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 
22 # include <deal.II/base/smartpointer.h>
23 # include <deal.II/base/subscriptor.h>
24 
25 # include <deal.II/lac/chunk_sparsity_pattern.h>
26 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/identity_matrix.h>
28 
29 # include <memory>
30 
31 
32 DEAL_II_NAMESPACE_OPEN
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 
116  Accessor(MatrixType *matrix);
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 
252  Accessor(MatrixType *matrix);
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 
319  Iterator(MatrixType *matrix);
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 
422 template <typename number>
423 class ChunkSparseMatrix : public virtual Subscriptor
424 {
425 public:
430 
435  using value_type = number;
436 
447 
453 
461 
468  struct Traits
469  {
474  static const bool zero_addition_can_be_elided = true;
475  };
476 
492 
502 
516  explicit ChunkSparseMatrix(const ChunkSparsityPattern &sparsity);
517 
524  ChunkSparseMatrix(const ChunkSparsityPattern &sparsity,
525  const IdentityMatrix & id);
526 
531  virtual ~ChunkSparseMatrix() override;
532 
543  operator=(const ChunkSparseMatrix<number> &);
544 
552  operator=(const IdentityMatrix &id);
553 
564  operator=(const double d);
565 
579  virtual void
580  reinit(const ChunkSparsityPattern &sparsity);
581 
587  virtual void
588  clear();
590 
598  bool
599  empty() const;
600 
605  size_type
606  m() const;
607 
612  size_type
613  n() const;
614 
620  size_type
621  n_nonzero_elements() const;
622 
630  size_type
631  n_actually_nonzero_elements() const;
632 
641  const ChunkSparsityPattern &
642  get_sparsity_pattern() const;
643 
648  std::size_t
649  memory_consumption() const;
650 
652 
661  void
662  set(const size_type i, const size_type j, const number value);
663 
669  void
670  add(const size_type i, const size_type j, const number value);
671 
681  template <typename number2>
682  void
683  add(const size_type row,
684  const size_type n_cols,
685  const size_type *col_indices,
686  const number2 * values,
687  const bool elide_zero_values = true,
688  const bool col_indices_are_sorted = false);
689 
694  operator*=(const number factor);
695 
700  operator/=(const number factor);
701 
714  void
715  symmetrize();
716 
733  template <typename somenumber>
735  copy_from(const ChunkSparseMatrix<somenumber> &source);
736 
754  template <typename ForwardIterator>
755  void
756  copy_from(const ForwardIterator begin, const ForwardIterator end);
757 
763  template <typename somenumber>
764  void
765  copy_from(const FullMatrix<somenumber> &matrix);
766 
778  template <typename somenumber>
779  void
780  add(const number factor, const ChunkSparseMatrix<somenumber> &matrix);
781 
783 
787 
801  number
802  operator()(const size_type i, const size_type j) const;
803 
816  number
817  el(const size_type i, const size_type j) const;
818 
828  number
829  diag_element(const size_type i) const;
830 
840  void
841  extract_row_copy(const size_type row,
842  const size_type array_length,
843  size_type & row_length,
844  size_type * column_indices,
845  number * values) const;
846 
848 
866  template <class OutVector, class InVector>
867  void
868  vmult(OutVector &dst, const InVector &src) const;
869 
885  template <class OutVector, class InVector>
886  void
887  Tvmult(OutVector &dst, const InVector &src) const;
888 
903  template <class OutVector, class InVector>
904  void
905  vmult_add(OutVector &dst, const InVector &src) const;
906 
922  template <class OutVector, class InVector>
923  void
924  Tvmult_add(OutVector &dst, const InVector &src) const;
925 
941  template <typename somenumber>
942  somenumber
943  matrix_norm_square(const Vector<somenumber> &v) const;
944 
948  template <typename somenumber>
949  somenumber
950  matrix_scalar_product(const Vector<somenumber> &u,
951  const Vector<somenumber> &v) const;
959  template <typename somenumber>
960  somenumber
961  residual(Vector<somenumber> & dst,
962  const Vector<somenumber> &x,
963  const Vector<somenumber> &b) const;
964 
966 
970 
978  real_type
979  l1_norm() const;
980 
988  real_type
989  linfty_norm() const;
990 
995  real_type
996  frobenius_norm() const;
998 
1002 
1008  template <typename somenumber>
1009  void
1010  precondition_Jacobi(Vector<somenumber> & dst,
1011  const Vector<somenumber> &src,
1012  const number omega = 1.) const;
1013 
1017  template <typename somenumber>
1018  void
1019  precondition_SSOR(Vector<somenumber> & dst,
1020  const Vector<somenumber> &src,
1021  const number om = 1.) const;
1022 
1026  template <typename somenumber>
1027  void
1028  precondition_SOR(Vector<somenumber> & dst,
1029  const Vector<somenumber> &src,
1030  const number om = 1.) const;
1031 
1035  template <typename somenumber>
1036  void
1037  precondition_TSOR(Vector<somenumber> & dst,
1038  const Vector<somenumber> &src,
1039  const number om = 1.) const;
1040 
1046  template <typename somenumber>
1047  void
1048  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1049 
1054  template <typename somenumber>
1055  void
1056  SOR(Vector<somenumber> &v, const number om = 1.) const;
1057 
1062  template <typename somenumber>
1063  void
1064  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1065 
1076  template <typename somenumber>
1077  void
1078  PSOR(Vector<somenumber> & v,
1079  const std::vector<size_type> &permutation,
1080  const std::vector<size_type> &inverse_permutation,
1081  const number om = 1.) const;
1082 
1093  template <typename somenumber>
1094  void
1095  TPSOR(Vector<somenumber> & v,
1096  const std::vector<size_type> &permutation,
1097  const std::vector<size_type> &inverse_permutation,
1098  const number om = 1.) const;
1099 
1104  template <typename somenumber>
1105  void
1106  SOR_step(Vector<somenumber> & v,
1107  const Vector<somenumber> &b,
1108  const number om = 1.) const;
1109 
1114  template <typename somenumber>
1115  void
1116  TSOR_step(Vector<somenumber> & v,
1117  const Vector<somenumber> &b,
1118  const number om = 1.) const;
1119 
1124  template <typename somenumber>
1125  void
1126  SSOR_step(Vector<somenumber> & v,
1127  const Vector<somenumber> &b,
1128  const number om = 1.) const;
1130 
1134 
1145  begin() const;
1146 
1156  end() const;
1157 
1167  iterator
1168  begin();
1169 
1178  iterator
1179  end();
1180 
1196  begin(const unsigned int r) const;
1197 
1213  end(const unsigned int r) const;
1214 
1229  iterator
1230  begin(const unsigned int r);
1231 
1246  iterator
1247  end(const unsigned int r);
1249 
1253 
1258  void
1259  print(std::ostream &out) const;
1260 
1281  void
1282  print_formatted(std::ostream & out,
1283  const unsigned int precision = 3,
1284  const bool scientific = true,
1285  const unsigned int width = 0,
1286  const char * zero_string = " ",
1287  const double denominator = 1.) const;
1288 
1294  void
1295  print_pattern(std::ostream &out, const double threshold = 0.) const;
1296 
1307  void
1308  block_write(std::ostream &out) const;
1309 
1326  void
1327  block_read(std::istream &in);
1329 
1337  DeclException2(ExcInvalidIndex,
1338  int,
1339  int,
1340  << "You are trying to access the matrix entry with index <"
1341  << arg1 << ',' << arg2
1342  << ">, but this entry does not exist in the sparsity pattern "
1343  "of this matrix."
1344  "\n\n"
1345  "The most common cause for this problem is that you used "
1346  "a method to build the sparsity pattern that did not "
1347  "(completely) take into account all of the entries you "
1348  "will later try to write into. An example would be "
1349  "building a sparsity pattern that does not include "
1350  "the entries you will write into due to constraints "
1351  "on degrees of freedom such as hanging nodes or periodic "
1352  "boundary conditions. In such cases, building the "
1353  "sparsity pattern will succeed, but you will get errors "
1354  "such as the current one at one point or other when "
1355  "trying to write into the entries of the matrix.");
1359  DeclException0(ExcDifferentChunkSparsityPatterns);
1363  DeclException2(ExcIteratorRange,
1364  int,
1365  int,
1366  << "The iterators denote a range of " << arg1
1367  << " elements, but the given number of rows was " << arg2);
1371  DeclExceptionMsg(ExcSourceEqualsDestination,
1372  "You are attempting an operation on two matrices that "
1373  "are the same object, but the operation requires that the "
1374  "two objects are in fact different.");
1376 private:
1383 
1391  std::unique_ptr<number[]> val;
1392 
1400 
1404  size_type
1405  compute_location(const size_type i, const size_type j) const;
1406 
1407  // make all other sparse matrices friends
1408  template <typename somenumber>
1409  friend class ChunkSparseMatrix;
1410 
1411  // Also give access to internal details to the iterator/accessor classes.
1412  template <typename, bool>
1414  template <typename, bool>
1416 };
1417 
1420 # ifndef DOXYGEN
1421 /*---------------------- Inline functions -----------------------------------*/
1422 
1423 
1424 
1425 template <typename number>
1428 {
1429  Assert(cols != nullptr, ExcNotInitialized());
1430  return cols->rows;
1431 }
1432 
1433 
1434 template <typename number>
1437 {
1438  Assert(cols != nullptr, ExcNotInitialized());
1439  return cols->cols;
1440 }
1441 
1442 
1443 
1444 template <typename number>
1445 inline const ChunkSparsityPattern &
1447 {
1448  Assert(cols != nullptr, ExcNotInitialized());
1449  return *cols;
1450 }
1451 
1452 
1453 
1454 template <typename number>
1457  const size_type j) const
1458 {
1459  const size_type chunk_size = cols->get_chunk_size();
1460  const size_type chunk_index =
1461  cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1462 
1463  if (chunk_index == ChunkSparsityPattern::invalid_entry)
1465  else
1466  {
1467  return (chunk_index * chunk_size * chunk_size +
1468  (i % chunk_size) * chunk_size + (j % chunk_size));
1469  }
1470 }
1471 
1472 
1473 template <typename number>
1474 inline void
1476  const size_type j,
1477  const number value)
1478 {
1479  AssertIsFinite(value);
1480 
1481  Assert(cols != nullptr, ExcNotInitialized());
1482  // it is allowed to set elements of the matrix that are not part of the
1483  // sparsity pattern, if the value to which we set it is zero
1484  const size_type index = compute_location(i, j);
1485  Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1486  ExcInvalidIndex(i, j));
1487 
1488  if (index != SparsityPattern::invalid_entry)
1489  val[index] = value;
1490 }
1491 
1492 
1493 
1494 template <typename number>
1495 inline void
1497  const size_type j,
1498  const number value)
1499 {
1500  AssertIsFinite(value);
1501 
1502  Assert(cols != nullptr, ExcNotInitialized());
1503 
1504  if (std::abs(value) != 0.)
1505  {
1506  const size_type index = compute_location(i, j);
1508  ExcInvalidIndex(i, j));
1509 
1510  val[index] += value;
1511  }
1512 }
1513 
1514 
1515 
1516 template <typename number>
1517 template <typename number2>
1518 inline void
1520  const size_type n_cols,
1521  const size_type *col_indices,
1522  const number2 * values,
1523  const bool /*elide_zero_values*/,
1524  const bool /*col_indices_are_sorted*/)
1525 {
1526  // TODO: could be done more efficiently...
1527  for (size_type col = 0; col < n_cols; ++col)
1528  add(row, col_indices[col], static_cast<number>(values[col]));
1529 }
1530 
1531 
1532 
1533 template <typename number>
1535 ChunkSparseMatrix<number>::operator*=(const number factor)
1536 {
1537  Assert(cols != nullptr, ExcNotInitialized());
1538  Assert(val != nullptr, ExcNotInitialized());
1539 
1540  const size_type chunk_size = cols->get_chunk_size();
1541 
1542  // multiply all elements of the matrix with the given factor. this includes
1543  // the padding elements in chunks that overlap the boundaries of the actual
1544  // matrix -- but since multiplication with a number does not violate the
1545  // invariant of keeping these elements at zero nothing can happen
1546  number * val_ptr = val.get();
1547  const number *const end_ptr =
1548  val.get() +
1549  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1550  while (val_ptr != end_ptr)
1551  *val_ptr++ *= factor;
1552 
1553  return *this;
1554 }
1555 
1556 
1557 
1558 template <typename number>
1560 ChunkSparseMatrix<number>::operator/=(const number factor)
1561 {
1562  Assert(cols != nullptr, ExcNotInitialized());
1563  Assert(val != nullptr, ExcNotInitialized());
1564  Assert(std::abs(factor) != 0, ExcDivideByZero());
1565 
1566  const number factor_inv = 1. / factor;
1567 
1568  const size_type chunk_size = cols->get_chunk_size();
1569 
1570  // multiply all elements of the matrix with the given factor. this includes
1571  // the padding elements in chunks that overlap the boundaries of the actual
1572  // matrix -- but since multiplication with a number does not violate the
1573  // invariant of keeping these elements at zero nothing can happen
1574  number * val_ptr = val.get();
1575  const number *const end_ptr =
1576  val.get() +
1577  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1578 
1579  while (val_ptr != end_ptr)
1580  *val_ptr++ *= factor_inv;
1581 
1582  return *this;
1583 }
1584 
1585 
1586 
1587 template <typename number>
1588 inline number
1590  const size_type j) const
1591 {
1592  Assert(cols != nullptr, ExcNotInitialized());
1593  AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1594  ExcInvalidIndex(i, j));
1595  return val[compute_location(i, j)];
1596 }
1597 
1598 
1599 
1600 template <typename number>
1601 inline number
1602 ChunkSparseMatrix<number>::el(const size_type i, const size_type j) const
1603 {
1604  Assert(cols != nullptr, ExcNotInitialized());
1605  const size_type index = compute_location(i, j);
1606 
1608  return val[index];
1609  else
1610  return 0;
1611 }
1612 
1613 
1614 
1615 template <typename number>
1616 inline number
1618 {
1619  Assert(cols != nullptr, ExcNotInitialized());
1620  Assert(m() == n(), ExcNotQuadratic());
1621  AssertIndexRange(i, m());
1622 
1623  // Use that the first element in each row of a quadratic matrix is the main
1624  // diagonal of the chunk sparsity pattern
1625  const size_type chunk_size = cols->get_chunk_size();
1626  return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1627  chunk_size +
1628  (i % chunk_size) * chunk_size + (i % chunk_size)];
1629 }
1630 
1631 
1632 
1633 template <typename number>
1634 template <typename ForwardIterator>
1635 inline void
1636 ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1637  const ForwardIterator end)
1638 {
1639  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1640  ExcIteratorRange(std::distance(begin, end), m()));
1641 
1642  // for use in the inner loop, we define an alias to the type of the inner
1643  // iterators
1644  using inner_iterator =
1645  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1646  size_type row = 0;
1647  for (ForwardIterator i = begin; i != end; ++i, ++row)
1648  {
1649  const inner_iterator end_of_row = i->end();
1650  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1651  // write entries
1652  set(row, j->first, j->second);
1653  }
1654 }
1655 
1656 
1657 
1658 //---------------------------------------------------------------------------
1659 
1660 
1662 {
1663  template <typename number>
1664  inline Accessor<number, true>::Accessor(const MatrixType * matrix,
1665  const unsigned int row)
1666  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1667  row)
1668  , matrix(matrix)
1669  {}
1670 
1671 
1672 
1673  template <typename number>
1674  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1675  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1676  , matrix(matrix)
1677  {}
1678 
1679 
1680 
1681  template <typename number>
1682  inline Accessor<number, true>::Accessor(
1684  : ChunkSparsityPatternIterators::Accessor(a)
1685  , matrix(&a.get_matrix())
1686  {}
1687 
1688 
1689 
1690  template <typename number>
1691  inline number
1692  Accessor<number, true>::value() const
1693  {
1694  const unsigned int chunk_size =
1695  matrix->get_sparsity_pattern().get_chunk_size();
1696  return matrix->val[reduced_index() * chunk_size * chunk_size +
1697  chunk_row * chunk_size + chunk_col];
1698  }
1699 
1700 
1701 
1702  template <typename number>
1703  inline const typename Accessor<number, true>::MatrixType &
1704  Accessor<number, true>::get_matrix() const
1705  {
1706  return *matrix;
1707  }
1708 
1709 
1710 
1711  template <typename number>
1712  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1713  const bool)
1714  : accessor(accessor)
1715  {}
1716 
1717 
1718  template <typename number>
1719  inline Accessor<number, false>::Reference::operator number() const
1720  {
1721  const unsigned int chunk_size =
1722  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1723  return accessor->matrix
1724  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1725  accessor->chunk_row * chunk_size + accessor->chunk_col];
1726  }
1727 
1728 
1729 
1730  template <typename number>
1731  inline const typename Accessor<number, false>::Reference &
1732  Accessor<number, false>::Reference::operator=(const number n) const
1733  {
1734  const unsigned int chunk_size =
1735  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1736  accessor->matrix
1737  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1738  accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1739  return *this;
1740  }
1741 
1742 
1743 
1744  template <typename number>
1745  inline const typename Accessor<number, false>::Reference &
1746  Accessor<number, false>::Reference::operator+=(const number n) const
1747  {
1748  const unsigned int chunk_size =
1749  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1750  accessor->matrix
1751  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1752  accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1753  return *this;
1754  }
1755 
1756 
1757 
1758  template <typename number>
1759  inline const typename Accessor<number, false>::Reference &
1760  Accessor<number, false>::Reference::operator-=(const number n) const
1761  {
1762  const unsigned int chunk_size =
1763  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1764  accessor->matrix
1765  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1766  accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1767  return *this;
1768  }
1769 
1770 
1771 
1772  template <typename number>
1773  inline const typename Accessor<number, false>::Reference &
1774  Accessor<number, false>::Reference::operator*=(const number n) const
1775  {
1776  const unsigned int chunk_size =
1777  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1778  accessor->matrix
1779  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1780  accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1781  return *this;
1782  }
1783 
1784 
1785 
1786  template <typename number>
1787  inline const typename Accessor<number, false>::Reference &
1788  Accessor<number, false>::Reference::operator/=(const number n) const
1789  {
1790  const unsigned int chunk_size =
1791  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1792  accessor->matrix
1793  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1794  accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1795  return *this;
1796  }
1797 
1798 
1799 
1800  template <typename number>
1801  inline Accessor<number, false>::Accessor(MatrixType * matrix,
1802  const unsigned int row)
1803  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1804  row)
1805  , matrix(matrix)
1806  {}
1807 
1808 
1809 
1810  template <typename number>
1811  inline Accessor<number, false>::Accessor(MatrixType *matrix)
1812  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1813  , matrix(matrix)
1814  {}
1815 
1816 
1817 
1818  template <typename number>
1819  inline typename Accessor<number, false>::Reference
1820  Accessor<number, false>::value() const
1821  {
1822  return Reference(this, true);
1823  }
1824 
1825 
1826 
1827  template <typename number>
1828  inline typename Accessor<number, false>::MatrixType &
1829  Accessor<number, false>::get_matrix() const
1830  {
1831  return *matrix;
1832  }
1833 
1834 
1835 
1836  template <typename number, bool Constness>
1837  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
1838  const unsigned int row)
1839  : accessor(matrix, row)
1840  {}
1841 
1842 
1843 
1844  template <typename number, bool Constness>
1845  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1846  : accessor(matrix)
1847  {}
1848 
1849 
1850 
1851  template <typename number, bool Constness>
1852  inline Iterator<number, Constness>::Iterator(
1854  : accessor(*i)
1855  {}
1856 
1857 
1858 
1859  template <typename number, bool Constness>
1860  inline Iterator<number, Constness> &
1861  Iterator<number, Constness>::operator++()
1862  {
1863  accessor.advance();
1864  return *this;
1865  }
1866 
1867 
1868  template <typename number, bool Constness>
1869  inline Iterator<number, Constness>
1870  Iterator<number, Constness>::operator++(int)
1871  {
1872  const Iterator iter = *this;
1873  accessor.advance();
1874  return iter;
1875  }
1876 
1877 
1878  template <typename number, bool Constness>
1879  inline const Accessor<number, Constness> &Iterator<number, Constness>::
1880  operator*() const
1881  {
1882  return accessor;
1883  }
1884 
1885 
1886  template <typename number, bool Constness>
1887  inline const Accessor<number, Constness> *Iterator<number, Constness>::
1888  operator->() const
1889  {
1890  return &accessor;
1891  }
1892 
1893 
1894  template <typename number, bool Constness>
1895  inline bool
1896  Iterator<number, Constness>::operator==(const Iterator &other) const
1897  {
1898  return (accessor == other.accessor);
1899  }
1900 
1901 
1902  template <typename number, bool Constness>
1903  inline bool
1904  Iterator<number, Constness>::operator!=(const Iterator &other) const
1905  {
1906  return !(*this == other);
1907  }
1908 
1909 
1910  template <typename number, bool Constness>
1911  inline bool
1912  Iterator<number, Constness>::operator<(const Iterator &other) const
1913  {
1914  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1915  ExcInternalError());
1916 
1917  return (accessor < other.accessor);
1918  }
1919 
1920 
1921  template <typename number, bool Constness>
1922  inline bool
1923  Iterator<number, Constness>::operator>(const Iterator &other) const
1924  {
1925  return (other < *this);
1926  }
1927 
1928 
1929  template <typename number, bool Constness>
1930  inline int
1931  Iterator<number, Constness>::operator-(const Iterator &other) const
1932  {
1933  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1934  ExcInternalError());
1935 
1936  // TODO: can be optimized
1937  int difference = 0;
1938  if (*this < other)
1939  {
1940  Iterator copy = *this;
1941  while (copy != other)
1942  {
1943  ++copy;
1944  --difference;
1945  }
1946  }
1947  else
1948  {
1949  Iterator copy = other;
1950  while (copy != *this)
1951  {
1952  ++copy;
1953  ++difference;
1954  }
1955  }
1956  return difference;
1957  }
1958 
1959 
1960 
1961  template <typename number, bool Constness>
1962  inline Iterator<number, Constness>
1963  Iterator<number, Constness>::operator+(const unsigned int n) const
1964  {
1965  Iterator x = *this;
1966  for (unsigned int i = 0; i < n; ++i)
1967  ++x;
1968 
1969  return x;
1970  }
1971 
1972 } // namespace ChunkSparseMatrixIterators
1973 
1974 
1975 
1976 template <typename number>
1979 {
1980  return const_iterator(this, 0);
1981 }
1982 
1983 
1984 template <typename number>
1987 {
1988  return const_iterator(this);
1989 }
1990 
1991 
1992 template <typename number>
1995 {
1996  return iterator(this, 0);
1997 }
1998 
1999 
2000 template <typename number>
2003 {
2004  return iterator(this);
2005 }
2006 
2007 
2008 template <typename number>
2010 ChunkSparseMatrix<number>::begin(const unsigned int r) const
2011 {
2012  Assert(r < m(), ExcIndexRange(r, 0, m()));
2013  return const_iterator(this, r);
2014 }
2015 
2016 
2017 
2018 template <typename number>
2020 ChunkSparseMatrix<number>::end(const unsigned int r) const
2021 {
2022  Assert(r < m(), ExcIndexRange(r, 0, m()));
2023  return const_iterator(this, r + 1);
2024 }
2025 
2026 
2027 
2028 template <typename number>
2030 ChunkSparseMatrix<number>::begin(const unsigned int r)
2031 {
2032  Assert(r < m(), ExcIndexRange(r, 0, m()));
2033  return iterator(this, r);
2034 }
2035 
2036 
2037 
2038 template <typename number>
2040 ChunkSparseMatrix<number>::end(const unsigned int r)
2041 {
2042  Assert(r < m(), ExcIndexRange(r, 0, m()));
2043  return iterator(this, r + 1);
2044 }
2045 
2046 
2047 
2048 # endif // DOXYGEN
2049 
2050 DEAL_II_NAMESPACE_CLOSE
2051 
2052 #endif
2053 /*--------------------------- chunk_sparse_matrix.h -------------------------*/
typename numbers::NumberTraits< number >::real_type real_type
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
size_type n() const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
size_type m() const
static ::ExceptionBase & ExcNotInitialized()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
std::unique_ptr< number[]> val
static ::ExceptionBase & ExcDivideByZero()
LinearAlgebra::distributed::Vector< Number > Vector
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:1407
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:496
#define DeclException0(Exception0)
Definition: exceptions.h:473
ChunkSparseMatrixIterators::Iterator< number, true > const_iterator
number el(const size_type i, const size_type j) const
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
static ::ExceptionBase & ExcNotQuadratic()
number diag_element(const size_type i) const
unsigned int global_dof_index
Definition: types.h:89
const_iterator end() const
const ChunkSparsityPattern & get_sparsity_pattern() const
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
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
types::global_dof_index size_type
void add(const size_type i, const size_type j, const number value)
typename Accessor< number, Constness >::MatrixType MatrixType
#define AssertIsFinite(number)
Definition: exceptions.h:1669
bool operator==(const Accessor &) const
static ::ExceptionBase & ExcInternalError()