Reference documentation for deal.II version 9.1.1
\(\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 template <typename number>
35 class Vector;
36 template <typename number>
37 class FullMatrix;
38 
48 {
49  // forward declaration
50  template <typename number, bool Constness>
51  class Iterator;
52 
63  template <typename number, bool Constness>
65  {
66  public:
70  number
71  value() const;
72 
76  number &
77  value();
78 
84  get_matrix() const;
85  };
86 
87 
88 
95  template <typename number>
97  {
98  public:
104 
108  Accessor(MatrixType *matrix, const unsigned int row);
109 
113  Accessor(MatrixType *matrix);
114 
119 
123  number
124  value() const;
125 
130  const MatrixType &
131  get_matrix() const;
132 
133  private:
138 
143 
147  template <typename, bool>
148  friend class Iterator;
149  };
150 
151 
158  template <typename number>
159  class Accessor<number, false> : public ChunkSparsityPatternIterators::Accessor
160  {
161  private:
184  class Reference
185  {
186  public:
191  Reference(const Accessor *accessor, const bool dummy);
192 
196  operator number() const;
197 
201  const Reference &
202  operator=(const number n) const;
203 
207  const Reference &
208  operator+=(const number n) const;
209 
213  const Reference &
214  operator-=(const number n) const;
215 
219  const Reference &
220  operator*=(const number n) const;
221 
225  const Reference &
226  operator/=(const number n) const;
227 
228  private:
234  };
235 
236  public:
242 
246  Accessor(MatrixType *matrix, const unsigned int row);
247 
251  Accessor(MatrixType *matrix);
252 
256  Reference
257  value() const;
258 
263  MatrixType &
264  get_matrix() const;
265 
266  private:
271 
276 
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 
315  Iterator(MatrixType *matrix, const unsigned int row);
316 
320  Iterator(MatrixType *matrix);
321 
327 
331  Iterator &
332  operator++();
333 
337  Iterator
338  operator++(int);
339 
343  const Accessor<number, Constness> &operator*() const;
344 
349 
353  bool
354  operator==(const Iterator &) const;
355 
359  bool
360  operator!=(const Iterator &) const;
361 
369  bool
370  operator<(const Iterator &) const;
371 
376  bool
377  operator>(const Iterator &) const;
378 
385  int
386  operator-(const Iterator &p) const;
387 
391  Iterator
392  operator+(const unsigned int n) const;
393 
394  private:
399  };
400 
401 } // namespace ChunkSparseMatrixIterators
402 
403 
404 
423 template <typename number>
424 class ChunkSparseMatrix : public virtual Subscriptor
425 {
426 public:
431 
436  using value_type = number;
437 
448 
454 
462 
469  struct Traits
470  {
475  static const bool zero_addition_can_be_elided = true;
476  };
477 
493 
503 
517  explicit ChunkSparseMatrix(const ChunkSparsityPattern &sparsity);
518 
525  ChunkSparseMatrix(const ChunkSparsityPattern &sparsity,
526  const IdentityMatrix & id);
527 
532  virtual ~ChunkSparseMatrix() override;
533 
545 
553  operator=(const IdentityMatrix &id);
554 
565  operator=(const double d);
566 
580  virtual void
581  reinit(const ChunkSparsityPattern &sparsity);
582 
588  virtual void
589  clear();
591 
599  bool
600  empty() const;
601 
606  size_type
607  m() const;
608 
613  size_type
614  n() const;
615 
621  size_type
622  n_nonzero_elements() const;
623 
631  size_type
633 
642  const ChunkSparsityPattern &
643  get_sparsity_pattern() const;
644 
649  std::size_t
650  memory_consumption() const;
651 
653 
662  void
663  set(const size_type i, const size_type j, const number value);
664 
670  void
671  add(const size_type i, const size_type j, const number value);
672 
682  template <typename number2>
683  void
684  add(const size_type row,
685  const size_type n_cols,
686  const size_type *col_indices,
687  const number2 * values,
688  const bool elide_zero_values = true,
689  const bool col_indices_are_sorted = false);
690 
695  operator*=(const number factor);
696 
701  operator/=(const number factor);
702 
715  void
716  symmetrize();
717 
734  template <typename somenumber>
737 
755  template <typename ForwardIterator>
756  void
757  copy_from(const ForwardIterator begin, const ForwardIterator end);
758 
764  template <typename somenumber>
765  void
766  copy_from(const FullMatrix<somenumber> &matrix);
767 
779  template <typename somenumber>
780  void
781  add(const number factor, const ChunkSparseMatrix<somenumber> &matrix);
782 
784 
788 
802  number
803  operator()(const size_type i, const size_type j) const;
804 
817  number
818  el(const size_type i, const size_type j) const;
819 
829  number
830  diag_element(const size_type i) const;
831 
841  void
842  extract_row_copy(const size_type row,
843  const size_type array_length,
844  size_type & row_length,
845  size_type * column_indices,
846  number * values) const;
847 
849 
867  template <class OutVector, class InVector>
868  void
869  vmult(OutVector &dst, const InVector &src) const;
870 
886  template <class OutVector, class InVector>
887  void
888  Tvmult(OutVector &dst, const InVector &src) const;
889 
904  template <class OutVector, class InVector>
905  void
906  vmult_add(OutVector &dst, const InVector &src) const;
907 
923  template <class OutVector, class InVector>
924  void
925  Tvmult_add(OutVector &dst, const InVector &src) const;
926 
942  template <typename somenumber>
943  somenumber
944  matrix_norm_square(const Vector<somenumber> &v) const;
945 
949  template <typename somenumber>
950  somenumber
951  matrix_scalar_product(const Vector<somenumber> &u,
952  const Vector<somenumber> &v) const;
960  template <typename somenumber>
961  somenumber
962  residual(Vector<somenumber> & dst,
963  const Vector<somenumber> &x,
964  const Vector<somenumber> &b) const;
965 
967 
971 
979  real_type
980  l1_norm() const;
981 
989  real_type
990  linfty_norm() const;
991 
996  real_type
997  frobenius_norm() const;
999 
1003 
1009  template <typename somenumber>
1010  void
1011  precondition_Jacobi(Vector<somenumber> & dst,
1012  const Vector<somenumber> &src,
1013  const number omega = 1.) const;
1014 
1018  template <typename somenumber>
1019  void
1020  precondition_SSOR(Vector<somenumber> & dst,
1021  const Vector<somenumber> &src,
1022  const number om = 1.) const;
1023 
1027  template <typename somenumber>
1028  void
1029  precondition_SOR(Vector<somenumber> & dst,
1030  const Vector<somenumber> &src,
1031  const number om = 1.) const;
1032 
1036  template <typename somenumber>
1037  void
1038  precondition_TSOR(Vector<somenumber> & dst,
1039  const Vector<somenumber> &src,
1040  const number om = 1.) const;
1041 
1047  template <typename somenumber>
1048  void
1049  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1050 
1055  template <typename somenumber>
1056  void
1057  SOR(Vector<somenumber> &v, const number om = 1.) const;
1058 
1063  template <typename somenumber>
1064  void
1065  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1066 
1077  template <typename somenumber>
1078  void
1079  PSOR(Vector<somenumber> & v,
1080  const std::vector<size_type> &permutation,
1081  const std::vector<size_type> &inverse_permutation,
1082  const number om = 1.) const;
1083 
1094  template <typename somenumber>
1095  void
1096  TPSOR(Vector<somenumber> & v,
1097  const std::vector<size_type> &permutation,
1098  const std::vector<size_type> &inverse_permutation,
1099  const number om = 1.) const;
1100 
1105  template <typename somenumber>
1106  void
1107  SOR_step(Vector<somenumber> & v,
1108  const Vector<somenumber> &b,
1109  const number om = 1.) const;
1110 
1115  template <typename somenumber>
1116  void
1117  TSOR_step(Vector<somenumber> & v,
1118  const Vector<somenumber> &b,
1119  const number om = 1.) const;
1120 
1125  template <typename somenumber>
1126  void
1127  SSOR_step(Vector<somenumber> & v,
1128  const Vector<somenumber> &b,
1129  const number om = 1.) const;
1131 
1135 
1146  begin() const;
1147 
1157  end() const;
1158 
1168  iterator
1169  begin();
1170 
1179  iterator
1180  end();
1181 
1197  begin(const unsigned int r) const;
1198 
1214  end(const unsigned int r) const;
1215 
1230  iterator
1231  begin(const unsigned int r);
1232 
1247  iterator
1248  end(const unsigned int r);
1250 
1254 
1259  void
1260  print(std::ostream &out) const;
1261 
1282  void
1283  print_formatted(std::ostream & out,
1284  const unsigned int precision = 3,
1285  const bool scientific = true,
1286  const unsigned int width = 0,
1287  const char * zero_string = " ",
1288  const double denominator = 1.) const;
1289 
1295  void
1296  print_pattern(std::ostream &out, const double threshold = 0.) const;
1297 
1308  void
1309  block_write(std::ostream &out) const;
1310 
1327  void
1328  block_read(std::istream &in);
1330 
1339  int,
1340  int,
1341  << "You are trying to access the matrix entry with index <"
1342  << arg1 << ',' << arg2
1343  << ">, but this entry does not exist in the sparsity pattern"
1344  "of this matrix."
1345  "\n\n"
1346  "The most common cause for this problem is that you used "
1347  "a method to build the sparsity pattern that did not "
1348  "(completely) take into account all of the entries you "
1349  "will later try to write into. An example would be "
1350  "building a sparsity pattern that does not include "
1351  "the entries you will write into due to constraints "
1352  "on degrees of freedom such as hanging nodes or periodic "
1353  "boundary conditions. In such cases, building the "
1354  "sparsity pattern will succeed, but you will get errors "
1355  "such as the current one at one point or other when "
1356  "trying to write into the entries of the matrix.");
1365  int,
1366  int,
1367  << "The iterators denote a range of " << arg1
1368  << " elements, but the given number of rows was " << arg2);
1373  "You are attempting an operation on two matrices that "
1374  "are the same object, but the operation requires that the "
1375  "two objects are in fact different.");
1377 private:
1384 
1392  std::unique_ptr<number[]> val;
1393 
1401 
1405  size_type
1406  compute_location(const size_type i, const size_type j) const;
1407 
1408  // make all other sparse matrices friends
1409  template <typename somenumber>
1410  friend class ChunkSparseMatrix;
1411 
1415  template <typename, bool>
1417  template <typename, bool>
1419 };
1420 
1423 # ifndef DOXYGEN
1424 /*---------------------- Inline functions -----------------------------------*/
1425 
1426 
1427 
1428 template <typename number>
1431 {
1432  Assert(cols != nullptr, ExcNotInitialized());
1433  return cols->rows;
1434 }
1435 
1436 
1437 template <typename number>
1440 {
1441  Assert(cols != nullptr, ExcNotInitialized());
1442  return cols->cols;
1443 }
1444 
1445 
1446 
1447 template <typename number>
1448 inline const ChunkSparsityPattern &
1450 {
1451  Assert(cols != nullptr, ExcNotInitialized());
1452  return *cols;
1453 }
1454 
1455 
1456 
1457 template <typename number>
1460  const size_type j) const
1461 {
1462  const size_type chunk_size = cols->get_chunk_size();
1463  const size_type chunk_index =
1464  cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1465 
1466  if (chunk_index == ChunkSparsityPattern::invalid_entry)
1468  else
1469  {
1470  return (chunk_index * chunk_size * chunk_size +
1471  (i % chunk_size) * chunk_size + (j % chunk_size));
1472  }
1473 }
1474 
1475 
1476 template <typename number>
1477 inline void
1479  const size_type j,
1480  const number value)
1481 {
1483 
1484  Assert(cols != nullptr, ExcNotInitialized());
1485  // it is allowed to set elements of the matrix that are not part of the
1486  // sparsity pattern, if the value to which we set it is zero
1487  const size_type index = compute_location(i, j);
1488  Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1489  ExcInvalidIndex(i, j));
1490 
1491  if (index != SparsityPattern::invalid_entry)
1492  val[index] = value;
1493 }
1494 
1495 
1496 
1497 template <typename number>
1498 inline void
1500  const size_type j,
1501  const number value)
1502 {
1504 
1505  Assert(cols != nullptr, ExcNotInitialized());
1506 
1507  if (std::abs(value) != 0.)
1508  {
1509  const size_type index = compute_location(i, j);
1511  ExcInvalidIndex(i, j));
1512 
1513  val[index] += value;
1514  }
1515 }
1516 
1517 
1518 
1519 template <typename number>
1520 template <typename number2>
1521 inline void
1523  const size_type n_cols,
1524  const size_type *col_indices,
1525  const number2 * values,
1526  const bool /*elide_zero_values*/,
1527  const bool /*col_indices_are_sorted*/)
1528 {
1529  // TODO: could be done more efficiently...
1530  for (size_type col = 0; col < n_cols; ++col)
1531  add(row, col_indices[col], static_cast<number>(values[col]));
1532 }
1533 
1534 
1535 
1536 template <typename number>
1538 ChunkSparseMatrix<number>::operator*=(const number factor)
1539 {
1540  Assert(cols != nullptr, ExcNotInitialized());
1541  Assert(val != nullptr, ExcNotInitialized());
1542 
1543  const size_type chunk_size = cols->get_chunk_size();
1544 
1545  // multiply all elements of the matrix with the given factor. this includes
1546  // the padding elements in chunks that overlap the boundaries of the actual
1547  // matrix -- but since multiplication with a number does not violate the
1548  // invariant of keeping these elements at zero nothing can happen
1549  number * val_ptr = val.get();
1550  const number *const end_ptr =
1551  val.get() +
1552  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1553  while (val_ptr != end_ptr)
1554  *val_ptr++ *= factor;
1555 
1556  return *this;
1557 }
1558 
1559 
1560 
1561 template <typename number>
1563 ChunkSparseMatrix<number>::operator/=(const number factor)
1564 {
1565  Assert(cols != nullptr, ExcNotInitialized());
1566  Assert(val != nullptr, ExcNotInitialized());
1567  Assert(std::abs(factor) != 0, ExcDivideByZero());
1568 
1569  const number factor_inv = 1. / factor;
1570 
1571  const size_type chunk_size = cols->get_chunk_size();
1572 
1573  // multiply all elements of the matrix with the given factor. this includes
1574  // the padding elements in chunks that overlap the boundaries of the actual
1575  // matrix -- but since multiplication with a number does not violate the
1576  // invariant of keeping these elements at zero nothing can happen
1577  number * val_ptr = val.get();
1578  const number *const end_ptr =
1579  val.get() +
1580  cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1581 
1582  while (val_ptr != end_ptr)
1583  *val_ptr++ *= factor_inv;
1584 
1585  return *this;
1586 }
1587 
1588 
1589 
1590 template <typename number>
1591 inline number
1593  const size_type j) const
1594 {
1595  Assert(cols != nullptr, ExcNotInitialized());
1596  AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1597  ExcInvalidIndex(i, j));
1598  return val[compute_location(i, j)];
1599 }
1600 
1601 
1602 
1603 template <typename number>
1604 inline number
1605 ChunkSparseMatrix<number>::el(const size_type i, const size_type j) const
1606 {
1607  Assert(cols != nullptr, ExcNotInitialized());
1608  const size_type index = compute_location(i, j);
1609 
1611  return val[index];
1612  else
1613  return 0;
1614 }
1615 
1616 
1617 
1618 template <typename number>
1619 inline number
1621 {
1622  Assert(cols != nullptr, ExcNotInitialized());
1623  Assert(m() == n(), ExcNotQuadratic());
1624  AssertIndexRange(i, m());
1625 
1626  // Use that the first element in each row of a quadratic matrix is the main
1627  // diagonal of the chunk sparsity pattern
1628  const size_type chunk_size = cols->get_chunk_size();
1629  return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1630  chunk_size +
1631  (i % chunk_size) * chunk_size + (i % chunk_size)];
1632 }
1633 
1634 
1635 
1636 template <typename number>
1637 template <typename ForwardIterator>
1638 inline void
1639 ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1640  const ForwardIterator end)
1641 {
1642  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1643  ExcIteratorRange(std::distance(begin, end), m()));
1644 
1645  // for use in the inner loop, we define an alias to the type of the inner
1646  // iterators
1647  using inner_iterator =
1648  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1649  size_type row = 0;
1650  for (ForwardIterator i = begin; i != end; ++i, ++row)
1651  {
1652  const inner_iterator end_of_row = i->end();
1653  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1654  // write entries
1655  set(row, j->first, j->second);
1656  }
1657 }
1658 
1659 
1660 
1661 //---------------------------------------------------------------------------
1662 
1663 
1665 {
1666  template <typename number>
1667  inline Accessor<number, true>::Accessor(const MatrixType * matrix,
1668  const unsigned int row)
1669  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1670  row)
1671  , matrix(matrix)
1672  {}
1673 
1674 
1675 
1676  template <typename number>
1677  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1678  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1679  , matrix(matrix)
1680  {}
1681 
1682 
1683 
1684  template <typename number>
1685  inline Accessor<number, true>::Accessor(
1687  : ChunkSparsityPatternIterators::Accessor(a)
1688  , matrix(&a.get_matrix())
1689  {}
1690 
1691 
1692 
1693  template <typename number>
1694  inline number
1695  Accessor<number, true>::value() const
1696  {
1697  const unsigned int chunk_size =
1698  matrix->get_sparsity_pattern().get_chunk_size();
1699  return matrix->val[reduced_index() * chunk_size * chunk_size +
1700  chunk_row * chunk_size + chunk_col];
1701  }
1702 
1703 
1704 
1705  template <typename number>
1706  inline const typename Accessor<number, true>::MatrixType &
1707  Accessor<number, true>::get_matrix() const
1708  {
1709  return *matrix;
1710  }
1711 
1712 
1713 
1714  template <typename number>
1715  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1716  const bool)
1717  : accessor(accessor)
1718  {}
1719 
1720 
1721  template <typename number>
1722  inline Accessor<number, false>::Reference::operator number() const
1723  {
1724  const unsigned int chunk_size =
1725  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1726  return accessor->matrix
1727  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1728  accessor->chunk_row * chunk_size + accessor->chunk_col];
1729  }
1730 
1731 
1732 
1733  template <typename number>
1734  inline const typename Accessor<number, false>::Reference &
1735  Accessor<number, false>::Reference::operator=(const number n) const
1736  {
1737  const unsigned int chunk_size =
1738  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1739  accessor->matrix
1740  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1741  accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1742  return *this;
1743  }
1744 
1745 
1746 
1747  template <typename number>
1748  inline const typename Accessor<number, false>::Reference &
1749  Accessor<number, false>::Reference::operator+=(const number n) const
1750  {
1751  const unsigned int chunk_size =
1752  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1753  accessor->matrix
1754  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1755  accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1756  return *this;
1757  }
1758 
1759 
1760 
1761  template <typename number>
1762  inline const typename Accessor<number, false>::Reference &
1763  Accessor<number, false>::Reference::operator-=(const number n) const
1764  {
1765  const unsigned int chunk_size =
1766  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1767  accessor->matrix
1768  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1769  accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1770  return *this;
1771  }
1772 
1773 
1774 
1775  template <typename number>
1776  inline const typename Accessor<number, false>::Reference &
1777  Accessor<number, false>::Reference::operator*=(const number n) const
1778  {
1779  const unsigned int chunk_size =
1780  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1781  accessor->matrix
1782  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1783  accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1784  return *this;
1785  }
1786 
1787 
1788 
1789  template <typename number>
1790  inline const typename Accessor<number, false>::Reference &
1791  Accessor<number, false>::Reference::operator/=(const number n) const
1792  {
1793  const unsigned int chunk_size =
1794  accessor->matrix->get_sparsity_pattern().get_chunk_size();
1795  accessor->matrix
1796  ->val[accessor->reduced_index() * chunk_size * chunk_size +
1797  accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1798  return *this;
1799  }
1800 
1801 
1802 
1803  template <typename number>
1804  inline Accessor<number, false>::Accessor(MatrixType * matrix,
1805  const unsigned int row)
1806  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1807  row)
1808  , matrix(matrix)
1809  {}
1810 
1811 
1812 
1813  template <typename number>
1814  inline Accessor<number, false>::Accessor(MatrixType *matrix)
1815  : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1816  , matrix(matrix)
1817  {}
1818 
1819 
1820 
1821  template <typename number>
1822  inline typename Accessor<number, false>::Reference
1823  Accessor<number, false>::value() const
1824  {
1825  return Reference(this, true);
1826  }
1827 
1828 
1829 
1830  template <typename number>
1831  inline typename Accessor<number, false>::MatrixType &
1832  Accessor<number, false>::get_matrix() const
1833  {
1834  return *matrix;
1835  }
1836 
1837 
1838 
1839  template <typename number, bool Constness>
1840  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
1841  const unsigned int row)
1842  : accessor(matrix, row)
1843  {}
1844 
1845 
1846 
1847  template <typename number, bool Constness>
1848  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1849  : accessor(matrix)
1850  {}
1851 
1852 
1853 
1854  template <typename number, bool Constness>
1855  inline Iterator<number, Constness>::Iterator(
1857  : accessor(*i)
1858  {}
1859 
1860 
1861 
1862  template <typename number, bool Constness>
1863  inline Iterator<number, Constness> &
1864  Iterator<number, Constness>::operator++()
1865  {
1866  accessor.advance();
1867  return *this;
1868  }
1869 
1870 
1871  template <typename number, bool Constness>
1872  inline Iterator<number, Constness>
1873  Iterator<number, Constness>::operator++(int)
1874  {
1875  const Iterator iter = *this;
1876  accessor.advance();
1877  return iter;
1878  }
1879 
1880 
1881  template <typename number, bool Constness>
1882  inline const Accessor<number, Constness> &Iterator<number, Constness>::
1883  operator*() const
1884  {
1885  return accessor;
1886  }
1887 
1888 
1889  template <typename number, bool Constness>
1890  inline const Accessor<number, Constness> *Iterator<number, Constness>::
1891  operator->() const
1892  {
1893  return &accessor;
1894  }
1895 
1896 
1897  template <typename number, bool Constness>
1898  inline bool
1899  Iterator<number, Constness>::operator==(const Iterator &other) const
1900  {
1901  return (accessor == other.accessor);
1902  }
1903 
1904 
1905  template <typename number, bool Constness>
1906  inline bool
1907  Iterator<number, Constness>::operator!=(const Iterator &other) const
1908  {
1909  return !(*this == other);
1910  }
1911 
1912 
1913  template <typename number, bool Constness>
1914  inline bool
1915  Iterator<number, Constness>::operator<(const Iterator &other) const
1916  {
1917  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1918  ExcInternalError());
1919 
1920  return (accessor < other.accessor);
1921  }
1922 
1923 
1924  template <typename number, bool Constness>
1925  inline bool
1926  Iterator<number, Constness>::operator>(const Iterator &other) const
1927  {
1928  return (other < *this);
1929  }
1930 
1931 
1932  template <typename number, bool Constness>
1933  inline int
1934  Iterator<number, Constness>::operator-(const Iterator &other) const
1935  {
1936  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1937  ExcInternalError());
1938 
1939  // TODO: can be optimized
1940  int difference = 0;
1941  if (*this < other)
1942  {
1943  Iterator copy = *this;
1944  while (copy != other)
1945  {
1946  ++copy;
1947  --difference;
1948  }
1949  }
1950  else
1951  {
1952  Iterator copy = other;
1953  while (copy != *this)
1954  {
1955  ++copy;
1956  ++difference;
1957  }
1958  }
1959  return difference;
1960  }
1961 
1962 
1963 
1964  template <typename number, bool Constness>
1965  inline Iterator<number, Constness>
1966  Iterator<number, Constness>::operator+(const unsigned int n) const
1967  {
1968  Iterator x = *this;
1969  for (unsigned int i = 0; i < n; ++i)
1970  ++x;
1971 
1972  return x;
1973  }
1974 
1975 } // namespace ChunkSparseMatrixIterators
1976 
1977 
1978 
1979 template <typename number>
1982 {
1983  return const_iterator(this, 0);
1984 }
1985 
1986 
1987 template <typename number>
1990 {
1991  return const_iterator(this);
1992 }
1993 
1994 
1995 template <typename number>
1998 {
1999  return iterator(this, 0);
2000 }
2001 
2002 
2003 template <typename number>
2006 {
2007  return iterator(this);
2008 }
2009 
2010 
2011 template <typename number>
2013 ChunkSparseMatrix<number>::begin(const unsigned int r) const
2014 {
2015  Assert(r < m(), ExcIndexRange(r, 0, m()));
2016  return const_iterator(this, r);
2017 }
2018 
2019 
2020 
2021 template <typename number>
2023 ChunkSparseMatrix<number>::end(const unsigned int r) const
2024 {
2025  Assert(r < m(), ExcIndexRange(r, 0, m()));
2026  return const_iterator(this, r + 1);
2027 }
2028 
2029 
2030 
2031 template <typename number>
2033 ChunkSparseMatrix<number>::begin(const unsigned int r)
2034 {
2035  Assert(r < m(), ExcIndexRange(r, 0, m()));
2036  return iterator(this, r);
2037 }
2038 
2039 
2040 
2041 template <typename number>
2043 ChunkSparseMatrix<number>::end(const unsigned int r)
2044 {
2045  Assert(r < m(), ExcIndexRange(r, 0, m()));
2046  return iterator(this, r + 1);
2047 }
2048 
2049 
2050 
2051 # endif // DOXYGEN
2052 
2053 DEAL_II_NAMESPACE_CLOSE
2054 
2055 #endif
2056 /*--------------------------- chunk_sparse_matrix.h -------------------------*/
typename numbers::NumberTraits< number >::real_type real_type
void block_write(std::ostream &out) const
size_type n() const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
void vmult(OutVector &dst, const InVector &src) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
bool empty() const
Contents is actually a matrix.
const Accessor< number, Constness > & operator*() const
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
bool operator==(const Iterator &) const
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void print(std::ostream &out) const
Iterator operator+(const unsigned int n) const
void print_pattern(std::ostream &out, const double threshold=0.) const
void extract_row_copy(const size_type row, const size_type array_length, size_type &row_length, size_type *column_indices, number *values) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
size_type m() const
virtual ~ChunkSparseMatrix() override
size_type n_nonzero_elements() const
static ::ExceptionBase & ExcNotInitialized()
Iterator(MatrixType *matrix, const unsigned int row)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
std::unique_ptr< number[]> val
static ::ExceptionBase & ExcDivideByZero()
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
LinearAlgebra::distributed::Vector< Number > Vector
const ChunkSparseMatrix< number > & get_matrix() const
const Accessor< number, Constness > * operator->() const
void SOR(Vector< somenumber > &v, const number om=1.) const
static ::ExceptionBase & ExcDifferentChunkSparsityPatterns()
void TSOR(Vector< somenumber > &v, const number om=1.) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
Accessor< number, Constness > accessor
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
std::size_t memory_consumption() const
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
const_iterator begin() const
number operator()(const size_type i, const size_type j) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
#define Assert(cond, exc)
Definition: exceptions.h:1407
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
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
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
virtual void reinit(const ChunkSparsityPattern &sparsity)
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
bool operator<(const Iterator &) const
number el(const size_type i, const size_type j) const
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
real_type frobenius_norm() const
void Tvmult(OutVector &dst, const InVector &src) const
static ::ExceptionBase & ExcNotQuadratic()
number diag_element(const size_type i) const
unsigned int global_dof_index
Definition: types.h:89
MatrixTableIterators::Iterator< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Iterator
Definition: table.h:1912
const_iterator end() const
void Tvmult_add(OutVector &dst, const InVector &src) const
size_type n_actually_nonzero_elements() const
void vmult_add(OutVector &dst, const InVector &src) const
static const bool zero_addition_can_be_elided
constexpr int chunk_size
Definition: cuda_size.h:35
const ChunkSparsityPattern & get_sparsity_pattern() const
void block_read(std::istream &in)
virtual void clear()
int operator-(const Iterator &p) const
ChunkSparseMatrix & operator/=(const number factor)
real_type l1_norm() const
const Accessor< number, Constness > & value_type
bool operator>(const Iterator &) const
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
types::global_dof_index size_type
void add(const size_type i, const size_type j, const number value)
bool operator!=(const Iterator &) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
typename Accessor< number, Constness >::MatrixType MatrixType
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
static ::ExceptionBase & ExcSourceEqualsDestination()
ChunkSparseMatrix< number > & operator=(const ChunkSparseMatrix< number > &)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
real_type linfty_norm() const
#define AssertIsFinite(number)
Definition: exceptions.h:1669
static ::ExceptionBase & ExcInternalError()