Reference documentation for deal.II version Git 50c3491829 2021-08-01 13:40:40 +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\}}\)
sparse_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_sparse_matrix_h
17 # define dealii_sparse_matrix_h
18 
19 
20 # include <deal.II/base/config.h>
21 
24 
25 # include <deal.II/lac/exceptions.h>
29 # ifdef DEAL_II_WITH_MPI
30 # include <mpi.h>
31 # endif
32 
33 # include <memory>
34 
35 
37 
38 // Forward declarations
39 # ifndef DOXYGEN
40 template <typename number>
41 class Vector;
42 template <typename number>
43 class FullMatrix;
44 template <typename Matrix>
45 class BlockMatrixBase;
46 template <typename number>
47 class SparseILU;
48 # ifdef DEAL_II_WITH_MPI
49 namespace Utilities
50 {
51  namespace MPI
52  {
53  template <typename Number>
54  void
56  }
57 } // namespace Utilities
58 # endif
59 
60 # ifdef DEAL_II_WITH_TRILINOS
61 namespace TrilinosWrappers
62 {
63  class SparseMatrix;
64 }
65 # endif
66 # endif
67 
78 {
83 
84  // forward declaration
85  template <typename number, bool Constness>
86  class Iterator;
87 
98  template <typename number, bool Constness>
100  {
101  public:
105  number
106  value() const;
107 
111  number &
112  value();
113 
118  const SparseMatrix<number> &
119  get_matrix() const;
120  };
121 
122 
123 
130  template <typename number>
131  class Accessor<number, true> : public SparsityPatternIterators::Accessor
132  {
133  public:
139 
143  Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
144 
149 
154 
158  number
159  value() const;
160 
165  const MatrixType &
166  get_matrix() const;
167 
168  private:
173 
178 
179  // Make iterator class a friend.
180  template <typename, bool>
181  friend class Iterator;
182  };
183 
184 
191  template <typename number>
192  class Accessor<number, false> : public SparsityPatternIterators::Accessor
193  {
194  private:
219  class Reference
220  {
221  public:
226  Reference(const Accessor *accessor, const bool dummy);
227 
231  operator number() const;
232 
236  const Reference &
237  operator=(const number n) const;
238 
242  const Reference &
243  operator+=(const number n) const;
244 
248  const Reference &
249  operator-=(const number n) const;
250 
254  const Reference &
255  operator*=(const number n) const;
256 
260  const Reference &
261  operator/=(const number n) const;
262 
263  private:
269  };
270 
271  public:
277 
281  Accessor(MatrixType *matrix, const std::size_t index);
282 
287 
291  Reference
292  value() const;
293 
298  MatrixType &
299  get_matrix() const;
300 
301  private:
306 
311 
312  // Make iterator class a friend.
313  template <typename, bool>
314  friend class Iterator;
315  };
316 
317 
318 
348  template <typename number, bool Constness>
349  class Iterator
350  {
351  public:
356 
362 
367  Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
368 
373 
379 
385 
389  Iterator &
390  operator++();
391 
395  Iterator
396  operator++(int);
397 
402  operator*() const;
403 
408  operator->() const;
409 
413  bool
414  operator==(const Iterator &) const;
415 
419  bool
420  operator!=(const Iterator &) const;
421 
429  bool
430  operator<(const Iterator &) const;
431 
436  bool
437  operator>(const Iterator &) const;
438 
445  int
446  operator-(const Iterator &p) const;
447 
451  Iterator
452  operator+(const size_type n) const;
453 
454  private:
459  };
460 
461 } // namespace SparseMatrixIterators
462 
468 // TODO: Add multithreading to the other vmult functions.
469 
496 template <typename number>
497 class SparseMatrix : public virtual Subscriptor
498 {
499 public:
504 
509  using value_type = number;
510 
521 
527 
535 
542  struct Traits
543  {
548  static const bool zero_addition_can_be_elided = true;
549  };
550 
565  SparseMatrix();
566 
575  SparseMatrix(const SparseMatrix &);
576 
584  SparseMatrix(SparseMatrix<number> &&m) noexcept;
585 
599  explicit SparseMatrix(const SparsityPattern &sparsity);
600 
607  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
608 
613  virtual ~SparseMatrix() override;
614 
625  operator=(const SparseMatrix<number> &);
626 
632  operator=(SparseMatrix<number> &&m) noexcept;
633 
641  operator=(const IdentityMatrix &id);
642 
654  SparseMatrix &
655  operator=(const double d);
656 
670  virtual void
671  reinit(const SparsityPattern &sparsity);
672 
678  virtual void
679  clear();
681 
689  bool
690  empty() const;
691 
696  size_type
697  m() const;
698 
703  size_type
704  n() const;
705 
709  size_type
710  get_row_length(const size_type row) const;
711 
717  std::size_t
718  n_nonzero_elements() const;
719 
729  std::size_t
730  n_actually_nonzero_elements(const double threshold = 0.) const;
731 
740  const SparsityPattern &
741  get_sparsity_pattern() const;
742 
747  std::size_t
748  memory_consumption() const;
749 
754 
756 
765  void
766  set(const size_type i, const size_type j, const number value);
767 
783  template <typename number2>
784  void
785  set(const std::vector<size_type> &indices,
786  const FullMatrix<number2> & full_matrix,
787  const bool elide_zero_values = false);
788 
794  template <typename number2>
795  void
796  set(const std::vector<size_type> &row_indices,
797  const std::vector<size_type> &col_indices,
798  const FullMatrix<number2> & full_matrix,
799  const bool elide_zero_values = false);
800 
811  template <typename number2>
812  void
813  set(const size_type row,
814  const std::vector<size_type> &col_indices,
815  const std::vector<number2> & values,
816  const bool elide_zero_values = false);
817 
827  template <typename number2>
828  void
829  set(const size_type row,
830  const size_type n_cols,
831  const size_type *col_indices,
832  const number2 * values,
833  const bool elide_zero_values = false);
834 
840  void
841  add(const size_type i, const size_type j, const number value);
842 
857  template <typename number2>
858  void
859  add(const std::vector<size_type> &indices,
860  const FullMatrix<number2> & full_matrix,
861  const bool elide_zero_values = true);
862 
868  template <typename number2>
869  void
870  add(const std::vector<size_type> &row_indices,
871  const std::vector<size_type> &col_indices,
872  const FullMatrix<number2> & full_matrix,
873  const bool elide_zero_values = true);
874 
884  template <typename number2>
885  void
886  add(const size_type row,
887  const std::vector<size_type> &col_indices,
888  const std::vector<number2> & values,
889  const bool elide_zero_values = true);
890 
900  template <typename number2>
901  void
902  add(const size_type row,
903  const size_type n_cols,
904  const size_type *col_indices,
905  const number2 * values,
906  const bool elide_zero_values = true,
907  const bool col_indices_are_sorted = false);
908 
912  SparseMatrix &
913  operator*=(const number factor);
914 
918  SparseMatrix &
919  operator/=(const number factor);
920 
933  void
934  symmetrize();
935 
952  template <typename somenumber>
954  copy_from(const SparseMatrix<somenumber> &source);
955 
972  template <typename ForwardIterator>
973  void
974  copy_from(const ForwardIterator begin, const ForwardIterator end);
975 
985  template <typename somenumber>
986  void
987  copy_from(const FullMatrix<somenumber> &matrix);
988 
989 # ifdef DEAL_II_WITH_TRILINOS
990 
1000  copy_from(const TrilinosWrappers::SparseMatrix &matrix);
1001 # endif
1002 
1014  template <typename somenumber>
1015  void
1016  add(const number factor, const SparseMatrix<somenumber> &matrix);
1017 
1019 
1023 
1037  const number &
1038  operator()(const size_type i, const size_type j) const;
1039 
1043  number &
1044  operator()(const size_type i, const size_type j);
1045 
1058  number
1059  el(const size_type i, const size_type j) const;
1060 
1070  number
1071  diag_element(const size_type i) const;
1072 
1077  number &
1078  diag_element(const size_type i);
1079 
1081 
1101  template <class OutVector, class InVector>
1102  void
1103  vmult(OutVector &dst, const InVector &src) const;
1104 
1120  template <class OutVector, class InVector>
1121  void
1122  Tvmult(OutVector &dst, const InVector &src) const;
1123 
1140  template <class OutVector, class InVector>
1141  void
1142  vmult_add(OutVector &dst, const InVector &src) const;
1143 
1159  template <class OutVector, class InVector>
1160  void
1161  Tvmult_add(OutVector &dst, const InVector &src) const;
1162 
1180  template <typename somenumber>
1181  somenumber
1182  matrix_norm_square(const Vector<somenumber> &v) const;
1183 
1189  template <typename somenumber>
1190  somenumber
1191  matrix_scalar_product(const Vector<somenumber> &u,
1192  const Vector<somenumber> &v) const;
1193 
1203  template <typename somenumber>
1204  somenumber
1205  residual(Vector<somenumber> & dst,
1206  const Vector<somenumber> &x,
1207  const Vector<somenumber> &b) const;
1208 
1244  template <typename numberB, typename numberC>
1245  void
1246  mmult(SparseMatrix<numberC> & C,
1247  const SparseMatrix<numberB> &B,
1248  const Vector<number> & V = Vector<number>(),
1249  const bool rebuild_sparsity_pattern = true) const;
1250 
1275  template <typename numberB, typename numberC>
1276  void
1277  Tmmult(SparseMatrix<numberC> & C,
1278  const SparseMatrix<numberB> &B,
1279  const Vector<number> & V = Vector<number>(),
1280  const bool rebuild_sparsity_pattern = true) const;
1281 
1283 
1287 
1295  real_type
1296  l1_norm() const;
1297 
1305  real_type
1306  linfty_norm() const;
1307 
1312  real_type
1313  frobenius_norm() const;
1315 
1319 
1325  template <typename somenumber>
1326  void
1327  precondition_Jacobi(Vector<somenumber> & dst,
1328  const Vector<somenumber> &src,
1329  const number omega = 1.) const;
1330 
1337  template <typename somenumber>
1338  void
1339  precondition_SSOR(Vector<somenumber> & dst,
1340  const Vector<somenumber> & src,
1341  const number omega = 1.,
1342  const std::vector<std::size_t> &pos_right_of_diagonal =
1343  std::vector<std::size_t>()) const;
1344 
1348  template <typename somenumber>
1349  void
1350  precondition_SOR(Vector<somenumber> & dst,
1351  const Vector<somenumber> &src,
1352  const number om = 1.) const;
1353 
1357  template <typename somenumber>
1358  void
1359  precondition_TSOR(Vector<somenumber> & dst,
1360  const Vector<somenumber> &src,
1361  const number om = 1.) const;
1362 
1368  template <typename somenumber>
1369  void
1370  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1371 
1376  template <typename somenumber>
1377  void
1378  SOR(Vector<somenumber> &v, const number om = 1.) const;
1379 
1384  template <typename somenumber>
1385  void
1386  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1387 
1398  template <typename somenumber>
1399  void
1400  PSOR(Vector<somenumber> & v,
1401  const std::vector<size_type> &permutation,
1402  const std::vector<size_type> &inverse_permutation,
1403  const number om = 1.) const;
1404 
1415  template <typename somenumber>
1416  void
1417  TPSOR(Vector<somenumber> & v,
1418  const std::vector<size_type> &permutation,
1419  const std::vector<size_type> &inverse_permutation,
1420  const number om = 1.) const;
1421 
1427  template <typename somenumber>
1428  void
1429  Jacobi_step(Vector<somenumber> & v,
1430  const Vector<somenumber> &b,
1431  const number om = 1.) const;
1432 
1437  template <typename somenumber>
1438  void
1439  SOR_step(Vector<somenumber> & v,
1440  const Vector<somenumber> &b,
1441  const number om = 1.) const;
1442 
1447  template <typename somenumber>
1448  void
1449  TSOR_step(Vector<somenumber> & v,
1450  const Vector<somenumber> &b,
1451  const number om = 1.) const;
1452 
1457  template <typename somenumber>
1458  void
1459  SSOR_step(Vector<somenumber> & v,
1460  const Vector<somenumber> &b,
1461  const number om = 1.) const;
1463 
1467 
1475  begin() const;
1476 
1480  iterator
1481  begin();
1482 
1487  end() const;
1488 
1492  iterator
1493  end();
1494 
1505  begin(const size_type r) const;
1506 
1510  iterator
1511  begin(const size_type r);
1512 
1523  end(const size_type r) const;
1524 
1528  iterator
1529  end(const size_type r);
1531 
1535 
1547  template <class StreamType>
1548  void
1549  print(StreamType &out,
1550  const bool across = false,
1551  const bool diagonal_first = true) const;
1552 
1573  void
1574  print_formatted(std::ostream & out,
1575  const unsigned int precision = 3,
1576  const bool scientific = true,
1577  const unsigned int width = 0,
1578  const char * zero_string = " ",
1579  const double denominator = 1.) const;
1580 
1586  void
1587  print_pattern(std::ostream &out, const double threshold = 0.) const;
1588 
1597  void
1598  print_as_numpy_arrays(std::ostream & out,
1599  const unsigned int precision = 9) const;
1600 
1611  void
1612  block_write(std::ostream &out) const;
1613 
1630  void
1631  block_read(std::istream &in);
1633 
1641  DeclException2(ExcInvalidIndex,
1642  int,
1643  int,
1644  << "You are trying to access the matrix entry with index <"
1645  << arg1 << ',' << arg2
1646  << ">, but this entry does not exist in the sparsity pattern "
1647  "of this matrix."
1648  "\n\n"
1649  "The most common cause for this problem is that you used "
1650  "a method to build the sparsity pattern that did not "
1651  "(completely) take into account all of the entries you "
1652  "will later try to write into. An example would be "
1653  "building a sparsity pattern that does not include "
1654  "the entries you will write into due to constraints "
1655  "on degrees of freedom such as hanging nodes or periodic "
1656  "boundary conditions. In such cases, building the "
1657  "sparsity pattern will succeed, but you will get errors "
1658  "such as the current one at one point or other when "
1659  "trying to write into the entries of the matrix.");
1663  DeclExceptionMsg(ExcDifferentSparsityPatterns,
1664  "When copying one sparse matrix into another, "
1665  "or when adding one sparse matrix to another, "
1666  "both matrices need to refer to the same "
1667  "sparsity pattern.");
1671  DeclException2(ExcIteratorRange,
1672  int,
1673  int,
1674  << "The iterators denote a range of " << arg1
1675  << " elements, but the given number of rows was " << arg2);
1679  DeclExceptionMsg(ExcSourceEqualsDestination,
1680  "You are attempting an operation on two matrices that "
1681  "are the same object, but the operation requires that the "
1682  "two objects are in fact different.");
1684 
1685 protected:
1696  void
1697  prepare_add();
1698 
1703  void
1704  prepare_set();
1705 
1706 private:
1713 
1721  std::unique_ptr<number[]> val;
1722 
1729  std::size_t max_len;
1730 
1731  // make all other sparse matrices friends
1732  template <typename somenumber>
1733  friend class SparseMatrix;
1734  template <typename somenumber>
1736  template <typename>
1737  friend class SparseILU;
1738 
1739  // To allow it calling private prepare_add() and prepare_set().
1740  template <typename>
1741  friend class BlockMatrixBase;
1742 
1743  // Also give access to internal details to the iterator/accessor classes.
1744  template <typename, bool>
1746  template <typename, bool>
1748 
1749 # ifdef DEAL_II_WITH_MPI
1750  // Give access to internal datastructures to perform MPI operations.
1751  template <typename Number>
1752  friend void
1754  const MPI_Comm &,
1756 # endif
1757 };
1758 
1759 # ifndef DOXYGEN
1760 /*---------------------- Inline functions -----------------------------------*/
1761 
1762 
1763 
1764 template <typename number>
1765 inline typename SparseMatrix<number>::size_type
1767 {
1768  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1769  return cols->rows;
1770 }
1771 
1772 
1773 template <typename number>
1774 inline typename SparseMatrix<number>::size_type
1776 {
1777  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1778  return cols->cols;
1779 }
1780 
1781 
1782 // Inline the set() and add() functions, since they will be called frequently.
1783 template <typename number>
1784 inline void
1786  const size_type j,
1787  const number value)
1788 {
1789  AssertIsFinite(value);
1790 
1791  const size_type index = cols->operator()(i, j);
1792 
1793  // it is allowed to set elements of the matrix that are not part of the
1794  // sparsity pattern, if the value to which we set it is zero
1795  if (index == SparsityPattern::invalid_entry)
1796  {
1797  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1798  ExcInvalidIndex(i, j));
1799  return;
1800  }
1801 
1802  val[index] = value;
1803 }
1804 
1805 
1806 
1807 template <typename number>
1808 template <typename number2>
1809 inline void
1810 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1811  const FullMatrix<number2> & values,
1812  const bool elide_zero_values)
1813 {
1814  Assert(indices.size() == values.m(),
1815  ExcDimensionMismatch(indices.size(), values.m()));
1816  Assert(values.m() == values.n(), ExcNotQuadratic());
1817 
1818  for (size_type i = 0; i < indices.size(); ++i)
1819  set(indices[i],
1820  indices.size(),
1821  indices.data(),
1822  &values(i, 0),
1823  elide_zero_values);
1824 }
1825 
1826 
1827 
1828 template <typename number>
1829 template <typename number2>
1830 inline void
1831 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1832  const std::vector<size_type> &col_indices,
1833  const FullMatrix<number2> & values,
1834  const bool elide_zero_values)
1835 {
1836  Assert(row_indices.size() == values.m(),
1837  ExcDimensionMismatch(row_indices.size(), values.m()));
1838  Assert(col_indices.size() == values.n(),
1839  ExcDimensionMismatch(col_indices.size(), values.n()));
1840 
1841  for (size_type i = 0; i < row_indices.size(); ++i)
1842  set(row_indices[i],
1843  col_indices.size(),
1844  col_indices.data(),
1845  &values(i, 0),
1846  elide_zero_values);
1847 }
1848 
1849 
1850 
1851 template <typename number>
1852 template <typename number2>
1853 inline void
1855  const std::vector<size_type> &col_indices,
1856  const std::vector<number2> & values,
1857  const bool elide_zero_values)
1858 {
1859  Assert(col_indices.size() == values.size(),
1860  ExcDimensionMismatch(col_indices.size(), values.size()));
1861 
1862  set(row,
1863  col_indices.size(),
1864  col_indices.data(),
1865  values.data(),
1866  elide_zero_values);
1867 }
1868 
1869 
1870 
1871 template <typename number>
1872 inline void
1874  const size_type j,
1875  const number value)
1876 {
1877  AssertIsFinite(value);
1878 
1879  if (value == number())
1880  return;
1881 
1882  const size_type index = cols->operator()(i, j);
1883 
1884  // it is allowed to add elements to the matrix that are not part of the
1885  // sparsity pattern, if the value to which we set it is zero
1886  if (index == SparsityPattern::invalid_entry)
1887  {
1888  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1889  ExcInvalidIndex(i, j));
1890  return;
1891  }
1892 
1893  val[index] += value;
1894 }
1895 
1896 
1897 
1898 template <typename number>
1899 template <typename number2>
1900 inline void
1901 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1902  const FullMatrix<number2> & values,
1903  const bool elide_zero_values)
1904 {
1905  Assert(indices.size() == values.m(),
1906  ExcDimensionMismatch(indices.size(), values.m()));
1907  Assert(values.m() == values.n(), ExcNotQuadratic());
1908 
1909  for (size_type i = 0; i < indices.size(); ++i)
1910  add(indices[i],
1911  indices.size(),
1912  indices.data(),
1913  &values(i, 0),
1914  elide_zero_values);
1915 }
1916 
1917 
1918 
1919 template <typename number>
1920 template <typename number2>
1921 inline void
1922 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1923  const std::vector<size_type> &col_indices,
1924  const FullMatrix<number2> & values,
1925  const bool elide_zero_values)
1926 {
1927  Assert(row_indices.size() == values.m(),
1928  ExcDimensionMismatch(row_indices.size(), values.m()));
1929  Assert(col_indices.size() == values.n(),
1930  ExcDimensionMismatch(col_indices.size(), values.n()));
1931 
1932  for (size_type i = 0; i < row_indices.size(); ++i)
1933  add(row_indices[i],
1934  col_indices.size(),
1935  col_indices.data(),
1936  &values(i, 0),
1937  elide_zero_values);
1938 }
1939 
1940 
1941 
1942 template <typename number>
1943 template <typename number2>
1944 inline void
1946  const std::vector<size_type> &col_indices,
1947  const std::vector<number2> & values,
1948  const bool elide_zero_values)
1949 {
1950  Assert(col_indices.size() == values.size(),
1951  ExcDimensionMismatch(col_indices.size(), values.size()));
1952 
1953  add(row,
1954  col_indices.size(),
1955  col_indices.data(),
1956  values.data(),
1957  elide_zero_values);
1958 }
1959 
1960 
1961 
1962 template <typename number>
1963 inline SparseMatrix<number> &
1964 SparseMatrix<number>::operator*=(const number factor)
1965 {
1966  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1967  Assert(val != nullptr, ExcNotInitialized());
1968 
1969  number * val_ptr = val.get();
1970  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1971 
1972  while (val_ptr != end_ptr)
1973  *val_ptr++ *= factor;
1974 
1975  return *this;
1976 }
1977 
1978 
1979 
1980 template <typename number>
1981 inline SparseMatrix<number> &
1982 SparseMatrix<number>::operator/=(const number factor)
1983 {
1984  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1985  Assert(val != nullptr, ExcNotInitialized());
1986  Assert(factor != number(), ExcDivideByZero());
1987 
1988  const number factor_inv = number(1.) / factor;
1989 
1990  number * val_ptr = val.get();
1991  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1992 
1993  while (val_ptr != end_ptr)
1994  *val_ptr++ *= factor_inv;
1995 
1996  return *this;
1997 }
1998 
1999 
2000 
2001 template <typename number>
2002 inline const number &
2004 {
2005  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2006  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2007  ExcInvalidIndex(i, j));
2008  return val[cols->operator()(i, j)];
2009 }
2010 
2011 
2012 
2013 template <typename number>
2014 inline number &
2016 {
2017  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2018  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2019  ExcInvalidIndex(i, j));
2020  return val[cols->operator()(i, j)];
2021 }
2022 
2023 
2024 
2025 template <typename number>
2026 inline number
2027 SparseMatrix<number>::el(const size_type i, const size_type j) const
2028 {
2029  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2030  const size_type index = cols->operator()(i, j);
2031 
2032  if (index != SparsityPattern::invalid_entry)
2033  return val[index];
2034  else
2035  return 0;
2036 }
2037 
2038 
2039 
2040 template <typename number>
2041 inline number
2043 {
2044  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2045  Assert(m() == n(), ExcNotQuadratic());
2046  AssertIndexRange(i, m());
2047 
2048  // Use that the first element in each row of a quadratic matrix is the main
2049  // diagonal
2050  return val[cols->rowstart[i]];
2051 }
2052 
2053 
2054 
2055 template <typename number>
2056 inline number &
2058 {
2059  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2060  Assert(m() == n(), ExcNotQuadratic());
2061  AssertIndexRange(i, m());
2062 
2063  // Use that the first element in each row of a quadratic matrix is the main
2064  // diagonal
2065  return val[cols->rowstart[i]];
2066 }
2067 
2068 
2069 
2070 template <typename number>
2071 template <typename ForwardIterator>
2072 void
2073 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2074  const ForwardIterator end)
2075 {
2076  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2077  ExcIteratorRange(std::distance(begin, end), m()));
2078 
2079  // for use in the inner loop, we define an alias to the type of the inner
2080  // iterators
2081  using inner_iterator =
2082  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2083  size_type row = 0;
2084  for (ForwardIterator i = begin; i != end; ++i, ++row)
2085  {
2086  const inner_iterator end_of_row = i->end();
2087  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2088  // write entries
2089  set(row, j->first, j->second);
2090  };
2091 }
2092 
2093 
2094 //---------------------------------------------------------------------------
2095 
2096 
2097 namespace SparseMatrixIterators
2098 {
2099  template <typename number>
2100  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2101  const std::size_t index_within_matrix)
2102  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2103  index_within_matrix)
2104  , matrix(matrix)
2105  {}
2106 
2107 
2108 
2109  template <typename number>
2110  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2111  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2112  , matrix(matrix)
2113  {}
2114 
2115 
2116 
2117  template <typename number>
2118  inline Accessor<number, true>::Accessor(
2120  : SparsityPatternIterators::Accessor(a)
2121  , matrix(&a.get_matrix())
2122  {}
2123 
2124 
2125 
2126  template <typename number>
2127  inline number
2128  Accessor<number, true>::value() const
2129  {
2130  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2131  return matrix->val[linear_index];
2132  }
2133 
2134 
2135 
2136  template <typename number>
2137  inline const typename Accessor<number, true>::MatrixType &
2138  Accessor<number, true>::get_matrix() const
2139  {
2140  return *matrix;
2141  }
2142 
2143 
2144 
2145  template <typename number>
2146  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2147  const bool)
2148  : accessor(accessor)
2149  {}
2150 
2151 
2152  template <typename number>
2153  inline Accessor<number, false>::Reference::operator number() const
2154  {
2155  AssertIndexRange(accessor->linear_index,
2156  accessor->matrix->n_nonzero_elements());
2157  return accessor->matrix->val[accessor->linear_index];
2158  }
2159 
2160 
2161 
2162  template <typename number>
2163  inline const typename Accessor<number, false>::Reference &
2164  Accessor<number, false>::Reference::operator=(const number n) const
2165  {
2166  AssertIndexRange(accessor->linear_index,
2167  accessor->matrix->n_nonzero_elements());
2168  accessor->matrix->val[accessor->linear_index] = n;
2169  return *this;
2170  }
2171 
2172 
2173 
2174  template <typename number>
2175  inline const typename Accessor<number, false>::Reference &
2176  Accessor<number, false>::Reference::operator+=(const number n) const
2177  {
2178  AssertIndexRange(accessor->linear_index,
2179  accessor->matrix->n_nonzero_elements());
2180  accessor->matrix->val[accessor->linear_index] += n;
2181  return *this;
2182  }
2183 
2184 
2185 
2186  template <typename number>
2187  inline const typename Accessor<number, false>::Reference &
2188  Accessor<number, false>::Reference::operator-=(const number n) const
2189  {
2190  AssertIndexRange(accessor->linear_index,
2191  accessor->matrix->n_nonzero_elements());
2192  accessor->matrix->val[accessor->linear_index] -= n;
2193  return *this;
2194  }
2195 
2196 
2197 
2198  template <typename number>
2199  inline const typename Accessor<number, false>::Reference &
2200  Accessor<number, false>::Reference::operator*=(const number n) const
2201  {
2202  AssertIndexRange(accessor->linear_index,
2203  accessor->matrix->n_nonzero_elements());
2204  accessor->matrix->val[accessor->linear_index] *= n;
2205  return *this;
2206  }
2207 
2208 
2209 
2210  template <typename number>
2211  inline const typename Accessor<number, false>::Reference &
2212  Accessor<number, false>::Reference::operator/=(const number n) const
2213  {
2214  AssertIndexRange(accessor->linear_index,
2215  accessor->matrix->n_nonzero_elements());
2216  accessor->matrix->val[accessor->linear_index] /= n;
2217  return *this;
2218  }
2219 
2220 
2221 
2222  template <typename number>
2223  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2224  const std::size_t index)
2225  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2226  , matrix(matrix)
2227  {}
2228 
2229 
2230 
2231  template <typename number>
2232  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2233  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2234  , matrix(matrix)
2235  {}
2236 
2237 
2238 
2239  template <typename number>
2240  inline typename Accessor<number, false>::Reference
2241  Accessor<number, false>::value() const
2242  {
2243  return Reference(this, true);
2244  }
2245 
2246 
2247 
2248  template <typename number>
2249  inline typename Accessor<number, false>::MatrixType &
2250  Accessor<number, false>::get_matrix() const
2251  {
2252  return *matrix;
2253  }
2254 
2255 
2256 
2257  template <typename number, bool Constness>
2258  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2259  const std::size_t index)
2260  : accessor(matrix, index)
2261  {}
2262 
2263 
2264 
2265  template <typename number, bool Constness>
2266  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2267  : accessor(matrix)
2268  {}
2269 
2270 
2271 
2272  template <typename number, bool Constness>
2273  inline Iterator<number, Constness>::Iterator(
2275  : accessor(*i)
2276  {}
2277 
2278 
2279 
2280  template <typename number, bool Constness>
2281  inline const Iterator<number, Constness> &
2282  Iterator<number, Constness>::operator=(
2284  {
2285  accessor = *i;
2286  return *this;
2287  }
2288 
2289 
2290 
2291  template <typename number, bool Constness>
2292  inline Iterator<number, Constness> &
2294  {
2295  accessor.advance();
2296  return *this;
2297  }
2298 
2299 
2300  template <typename number, bool Constness>
2301  inline Iterator<number, Constness>
2303  {
2304  const Iterator iter = *this;
2305  accessor.advance();
2306  return iter;
2307  }
2308 
2309 
2310  template <typename number, bool Constness>
2311  inline const Accessor<number, Constness> &
2313  {
2314  return accessor;
2315  }
2316 
2317 
2318  template <typename number, bool Constness>
2319  inline const Accessor<number, Constness> *
2320  Iterator<number, Constness>::operator->() const
2321  {
2322  return &accessor;
2323  }
2324 
2325 
2326  template <typename number, bool Constness>
2327  inline bool
2328  Iterator<number, Constness>::operator==(const Iterator &other) const
2329  {
2330  return (accessor == other.accessor);
2331  }
2332 
2333 
2334  template <typename number, bool Constness>
2335  inline bool
2336  Iterator<number, Constness>::operator!=(const Iterator &other) const
2337  {
2338  return !(*this == other);
2339  }
2340 
2341 
2342  template <typename number, bool Constness>
2343  inline bool
2344  Iterator<number, Constness>::operator<(const Iterator &other) const
2345  {
2346  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2347  ExcInternalError());
2348 
2349  return (accessor < other.accessor);
2350  }
2351 
2352 
2353  template <typename number, bool Constness>
2354  inline bool
2355  Iterator<number, Constness>::operator>(const Iterator &other) const
2356  {
2357  return (other < *this);
2358  }
2359 
2360 
2361  template <typename number, bool Constness>
2362  inline int
2363  Iterator<number, Constness>::operator-(const Iterator &other) const
2364  {
2365  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2366  ExcInternalError());
2367 
2368  return (*this)->linear_index - other->linear_index;
2369  }
2370 
2371 
2372 
2373  template <typename number, bool Constness>
2374  inline Iterator<number, Constness>
2376  {
2377  Iterator x = *this;
2378  for (size_type i = 0; i < n; ++i)
2379  ++x;
2380 
2381  return x;
2382  }
2383 
2384 } // namespace SparseMatrixIterators
2385 
2386 
2387 
2388 template <typename number>
2391 {
2392  return const_iterator(this, 0);
2393 }
2394 
2395 
2396 template <typename number>
2399 {
2400  return const_iterator(this);
2401 }
2402 
2403 
2404 template <typename number>
2405 inline typename SparseMatrix<number>::iterator
2407 {
2408  return iterator(this, 0);
2409 }
2410 
2411 
2412 template <typename number>
2413 inline typename SparseMatrix<number>::iterator
2415 {
2416  return iterator(this, cols->rowstart[cols->rows]);
2417 }
2418 
2419 
2420 template <typename number>
2423 {
2424  AssertIndexRange(r, m());
2425 
2426  return const_iterator(this, cols->rowstart[r]);
2427 }
2428 
2429 
2430 
2431 template <typename number>
2433 SparseMatrix<number>::end(const size_type r) const
2434 {
2435  AssertIndexRange(r, m());
2436 
2437  return const_iterator(this, cols->rowstart[r + 1]);
2438 }
2439 
2440 
2441 
2442 template <typename number>
2443 inline typename SparseMatrix<number>::iterator
2445 {
2446  AssertIndexRange(r, m());
2447 
2448  return iterator(this, cols->rowstart[r]);
2449 }
2450 
2451 
2452 
2453 template <typename number>
2454 inline typename SparseMatrix<number>::iterator
2456 {
2457  AssertIndexRange(r, m());
2458 
2459  return iterator(this, cols->rowstart[r + 1]);
2460 }
2461 
2462 
2463 
2464 template <typename number>
2465 template <class StreamType>
2466 inline void
2467 SparseMatrix<number>::print(StreamType &out,
2468  const bool across,
2469  const bool diagonal_first) const
2470 {
2471  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2472  Assert(val != nullptr, ExcNotInitialized());
2473 
2474  bool hanging_diagonal = false;
2475  number diagonal = number();
2476 
2477  for (size_type i = 0; i < cols->rows; ++i)
2478  {
2479  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2480  {
2481  if (!diagonal_first && i == cols->colnums[j])
2482  {
2483  diagonal = val[j];
2484  hanging_diagonal = true;
2485  }
2486  else
2487  {
2488  if (hanging_diagonal && cols->colnums[j] > i)
2489  {
2490  if (across)
2491  out << ' ' << i << ',' << i << ':' << diagonal;
2492  else
2493  out << '(' << i << ',' << i << ") " << diagonal
2494  << std::endl;
2495  hanging_diagonal = false;
2496  }
2497  if (across)
2498  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2499  else
2500  out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2501  << std::endl;
2502  }
2503  }
2504  if (hanging_diagonal)
2505  {
2506  if (across)
2507  out << ' ' << i << ',' << i << ':' << diagonal;
2508  else
2509  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2510  hanging_diagonal = false;
2511  }
2512  }
2513  if (across)
2514  out << std::endl;
2515 }
2516 
2517 
2518 template <typename number>
2519 inline void
2521 {
2522  // nothing to do here
2523 }
2524 
2525 
2526 
2527 template <typename number>
2528 inline void
2530 {
2531  // nothing to do here
2532 }
2533 
2534 # endif // DOXYGEN
2535 
2536 
2537 /*---------------------------- sparse_matrix.h ---------------------------*/
2538 
2540 
2541 #endif
2542 /*---------------------------- sparse_matrix.h ---------------------------*/
typename Accessor< number, Constness >::MatrixType MatrixType
size_type m() const
SparseMatrix & operator/=(const number factor)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
typename numbers::NumberTraits< number >::real_type real_type
void prepare_add()
const_iterator end() const
Contents is actually a matrix.
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
std::unique_ptr< number[]> val
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void set(const size_type i, const size_type j, const number value)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1724
number value_type
static const char V
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcNotInitialized()
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcDivideByZero()
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
std::unique_ptr< std::size_t[]> rowstart
Matrix is diagonal.
size_type n() const
std::string compress(const std::string &input)
Definition: utilities.cc:392
size_type n() const
T sum(const T &t, const MPI_Comm &mpi_communicator)
number diag_element(const size_type i) const
#define Assert(cond, exc)
Definition: exceptions.h:1467
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
SparseMatrixIterators::Iterator< number, false > iterator
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2944
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
Accessor< number, Constness > accessor
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
VectorType::value_type * end(VectorType &V)
void prepare_set()
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
Expression operator>(const Expression &lhs, const Expression &rhs)
void add(const size_type i, const size_type j, const number value)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparseMatrix< double > SparseMatrix
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
types::global_dof_index size_type
std::unique_ptr< size_type[]> colnums
static ::ExceptionBase & ExcNotQuadratic()
unsigned int global_dof_index
Definition: types.h:76
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
number el(const size_type i, const size_type j) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2918
size_type m() const
static const size_type invalid_entry
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)
const SparsityPattern & get_sparsity_pattern() const
const_iterator begin() const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
const number & operator()(const size_type i, const size_type j) const
static ::ExceptionBase & ExcNeedsSparsityPattern()
SparseMatrixIterators::Iterator< number, true > const_iterator
#define AssertIsFinite(number)
Definition: exceptions.h:1750
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::size_t max_len
const Accessor< number, Constness > & value_type
SparseMatrix & operator*=(const number factor)
static ::ExceptionBase & ExcInternalError()