Reference documentation for deal.II version GIT c1ddcf411d 2023-11-30 16:30:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
sparse_matrix.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2023 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 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/mpi_stub.h>
24 
25 #include <deal.II/lac/exceptions.h>
29 
30 #include <iterator>
31 #include <memory>
32 
33 
35 
36 // Forward declarations
37 #ifndef DOXYGEN
38 template <typename number>
39 class Vector;
40 template <typename number>
41 class FullMatrix;
42 template <typename Matrix>
43 class BlockMatrixBase;
44 template <typename number>
45 class SparseILU;
46 # ifdef DEAL_II_WITH_MPI
47 namespace Utilities
48 {
49  namespace MPI
50  {
51  template <typename Number>
52  void
53  sum(const SparseMatrix<Number> &, const MPI_Comm, SparseMatrix<Number> &);
54  }
55 } // namespace Utilities
56 # endif
57 
58 # ifdef DEAL_II_WITH_TRILINOS
59 namespace TrilinosWrappers
60 {
61  class SparseMatrix;
62 }
63 # endif
64 #endif
65 
76 {
81 
82  // forward declaration
83  template <typename number, bool Constness>
84  class Iterator;
85 
96  template <typename number, bool Constness>
98  {
99  public:
103  number
104  value() const;
105 
109  number &
110  value();
111 
116  const SparseMatrix<number> &
117  get_matrix() const;
118  };
119 
120 
121 
128  template <typename number>
129  class Accessor<number, true> : public SparsityPatternIterators::Accessor
130  {
131  public:
137 
141  Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
142 
147 
152 
156  number
157  value() const;
158 
163  const MatrixType &
164  get_matrix() const;
165 
166  private:
171 
176 
177  // Make iterator class a friend.
178  template <typename, bool>
179  friend class Iterator;
180  };
181 
182 
189  template <typename number>
190  class Accessor<number, false> : public SparsityPatternIterators::Accessor
191  {
192  private:
217  class Reference
218  {
219  public:
224  Reference(const Accessor *accessor, const bool dummy);
225 
229  operator number() const;
230 
234  const Reference &
235  operator=(const number n) const;
236 
240  const Reference &
241  operator+=(const number n) const;
242 
246  const Reference &
247  operator-=(const number n) const;
248 
252  const Reference &
253  operator*=(const number n) const;
254 
258  const Reference &
259  operator/=(const number n) const;
260 
261  private:
267  };
268 
269  public:
275 
279  Accessor(MatrixType *matrix, const std::size_t index);
280 
285 
289  Reference
290  value() const;
291 
296  MatrixType &
297  get_matrix() const;
298 
299  private:
304 
309 
310  // Make iterator class a friend.
311  template <typename, bool>
312  friend class Iterator;
313  };
314 
315 
316 
346  template <typename number, bool Constness>
347  class Iterator
348  {
349  public:
354 
360 
366 
371  Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
372 
377 
383 
389 
393  Iterator &
395 
399  Iterator
401 
406  operator*() const;
407 
412  operator->() const;
413 
417  bool
418  operator==(const Iterator &) const;
419 
423  bool
424  operator!=(const Iterator &) const;
425 
433  bool
434  operator<(const Iterator &) const;
435 
440  bool
441  operator>(const Iterator &) const;
442 
449  int
450  operator-(const Iterator &p) const;
451 
455  Iterator
456  operator+(const size_type n) const;
457 
458  private:
463  };
464 
465 } // namespace SparseMatrixIterators
466 
468 
469 namespace std
470 {
471  template <typename number, bool Constness>
473  ::SparseMatrixIterators::Iterator<number, Constness>>
474  {
475  using iterator_category = forward_iterator_tag;
476  using value_type =
477  typename ::SparseMatrixIterators::Iterator<number,
478  Constness>::value_type;
479  using difference_type = typename ::SparseMatrixIterators::
480  Iterator<number, Constness>::difference_type;
481  };
482 } // namespace std
483 
485 
491 // TODO: Add multithreading to the other vmult functions.
492 
519 template <typename number>
520 class SparseMatrix : public virtual Subscriptor
521 {
522 public:
527 
532  using value_type = number;
533 
544 
550 
558 
565  struct Traits
566  {
571  static const bool zero_addition_can_be_elided = true;
572  };
573 
589 
599 
608 
622  explicit SparseMatrix(const SparsityPattern &sparsity);
623 
630  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
631 
636  virtual ~SparseMatrix() override;
637 
649 
656 
665 
677  SparseMatrix &
678  operator=(const double d);
679 
693  virtual void
694  reinit(const SparsityPattern &sparsity);
695 
702  template <typename number2>
703  void
704  reinit(const SparseMatrix<number2> &sparse_matrix);
705 
711  virtual void
712  clear();
722  bool
723  empty() const;
724 
729  size_type
730  m() const;
731 
736  size_type
737  n() const;
738 
742  size_type
743  get_row_length(const size_type row) const;
744 
750  std::size_t
752 
762  std::size_t
763  n_actually_nonzero_elements(const double threshold = 0.) const;
764 
773  const SparsityPattern &
775 
780  std::size_t
782 
787 
798  void
799  set(const size_type i, const size_type j, const number value);
800 
816  template <typename number2>
817  void
818  set(const std::vector<size_type> &indices,
819  const FullMatrix<number2> &full_matrix,
820  const bool elide_zero_values = false);
821 
827  template <typename number2>
828  void
829  set(const std::vector<size_type> &row_indices,
830  const std::vector<size_type> &col_indices,
831  const FullMatrix<number2> &full_matrix,
832  const bool elide_zero_values = false);
833 
844  template <typename number2>
845  void
846  set(const size_type row,
847  const std::vector<size_type> &col_indices,
848  const std::vector<number2> &values,
849  const bool elide_zero_values = false);
850 
860  template <typename number2>
861  void
862  set(const size_type row,
863  const size_type n_cols,
864  const size_type *col_indices,
865  const number2 *values,
866  const bool elide_zero_values = false);
867 
873  void
874  add(const size_type i, const size_type j, const number value);
875 
890  template <typename number2>
891  void
892  add(const std::vector<size_type> &indices,
893  const FullMatrix<number2> &full_matrix,
894  const bool elide_zero_values = true);
895 
901  template <typename number2>
902  void
903  add(const std::vector<size_type> &row_indices,
904  const std::vector<size_type> &col_indices,
905  const FullMatrix<number2> &full_matrix,
906  const bool elide_zero_values = true);
907 
917  template <typename number2>
918  void
919  add(const size_type row,
920  const std::vector<size_type> &col_indices,
921  const std::vector<number2> &values,
922  const bool elide_zero_values = true);
923 
933  template <typename number2>
934  void
935  add(const size_type row,
936  const size_type n_cols,
937  const size_type *col_indices,
938  const number2 *values,
939  const bool elide_zero_values = true,
940  const bool col_indices_are_sorted = false);
941 
945  SparseMatrix &
946  operator*=(const number factor);
947 
951  SparseMatrix &
952  operator/=(const number factor);
953 
966  void
968 
985  template <typename somenumber>
988 
1005  template <typename ForwardIterator>
1006  void
1007  copy_from(const ForwardIterator begin, const ForwardIterator end);
1008 
1018  template <typename somenumber>
1019  void
1021 
1022 #ifdef DEAL_II_WITH_TRILINOS
1034 #endif
1035 
1047  template <typename somenumber>
1048  void
1049  add(const number factor, const SparseMatrix<somenumber> &matrix);
1050 
1070  const number &
1071  operator()(const size_type i, const size_type j) const;
1072 
1076  number &
1077  operator()(const size_type i, const size_type j);
1078 
1091  number
1092  el(const size_type i, const size_type j) const;
1093 
1103  number
1104  diag_element(const size_type i) const;
1105 
1110  number &
1112 
1134  template <class OutVector, class InVector>
1135  void
1136  vmult(OutVector &dst, const InVector &src) const;
1137 
1153  template <class OutVector, class InVector>
1154  void
1155  Tvmult(OutVector &dst, const InVector &src) const;
1156 
1173  template <class OutVector, class InVector>
1174  void
1175  vmult_add(OutVector &dst, const InVector &src) const;
1176 
1192  template <class OutVector, class InVector>
1193  void
1194  Tvmult_add(OutVector &dst, const InVector &src) const;
1195 
1213  template <typename somenumber>
1214  somenumber
1216 
1222  template <typename somenumber>
1223  somenumber
1225  const Vector<somenumber> &v) const;
1226 
1236  template <typename somenumber>
1237  somenumber
1239  const Vector<somenumber> &x,
1240  const Vector<somenumber> &b) const;
1241 
1277  template <typename numberB, typename numberC>
1278  void
1280  const SparseMatrix<numberB> &B,
1281  const Vector<number> &V = Vector<number>(),
1282  const bool rebuild_sparsity_pattern = true) const;
1283 
1308  template <typename numberB, typename numberC>
1309  void
1311  const SparseMatrix<numberB> &B,
1312  const Vector<number> &V = Vector<number>(),
1313  const bool rebuild_sparsity_pattern = true) const;
1314 
1328  real_type
1329  l1_norm() const;
1330 
1338  real_type
1339  linfty_norm() const;
1340 
1345  real_type
1358  template <typename somenumber>
1359  void
1361  const Vector<somenumber> &src,
1362  const number omega = 1.) const;
1363 
1370  template <typename somenumber>
1371  void
1373  const Vector<somenumber> &src,
1374  const number omega = 1.,
1375  const std::vector<std::size_t> &pos_right_of_diagonal =
1376  std::vector<std::size_t>()) const;
1377 
1381  template <typename somenumber>
1382  void
1384  const Vector<somenumber> &src,
1385  const number omega = 1.) const;
1386 
1390  template <typename somenumber>
1391  void
1393  const Vector<somenumber> &src,
1394  const number omega = 1.) const;
1395 
1401  template <typename somenumber>
1402  void
1403  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1404 
1409  template <typename somenumber>
1410  void
1411  SOR(Vector<somenumber> &v, const number omega = 1.) const;
1412 
1417  template <typename somenumber>
1418  void
1419  TSOR(Vector<somenumber> &v, const number omega = 1.) const;
1420 
1431  template <typename somenumber>
1432  void
1434  const std::vector<size_type> &permutation,
1435  const std::vector<size_type> &inverse_permutation,
1436  const number omega = 1.) const;
1437 
1448  template <typename somenumber>
1449  void
1451  const std::vector<size_type> &permutation,
1452  const std::vector<size_type> &inverse_permutation,
1453  const number omega = 1.) const;
1454 
1460  template <typename somenumber>
1461  void
1463  const Vector<somenumber> &b,
1464  const number omega = 1.) const;
1465 
1470  template <typename somenumber>
1471  void
1473  const Vector<somenumber> &b,
1474  const number omega = 1.) const;
1475 
1480  template <typename somenumber>
1481  void
1483  const Vector<somenumber> &b,
1484  const number omega = 1.) const;
1485 
1490  template <typename somenumber>
1491  void
1493  const Vector<somenumber> &b,
1494  const number omega = 1.) const;
1508  begin() const;
1509 
1513  iterator
1515 
1520  end() const;
1521 
1525  iterator
1526  end();
1527 
1538  begin(const size_type r) const;
1539 
1543  iterator
1544  begin(const size_type r);
1545 
1556  end(const size_type r) const;
1557 
1561  iterator
1562  end(const size_type r);
1580  template <typename StreamType>
1581  void
1582  print(StreamType &out,
1583  const bool across = false,
1584  const bool diagonal_first = true) const;
1585 
1608  void
1609  print_formatted(std::ostream &out,
1610  const unsigned int precision = 3,
1611  const bool scientific = true,
1612  const unsigned int width = 0,
1613  const char *zero_string = " ",
1614  const double denominator = 1.,
1615  const char *separator = " ") const;
1616 
1622  void
1623  print_pattern(std::ostream &out, const double threshold = 0.) const;
1624 
1633  void
1634  print_as_numpy_arrays(std::ostream &out,
1635  const unsigned int precision = 9) const;
1636 
1647  void
1648  block_write(std::ostream &out) const;
1649 
1666  void
1667  block_read(std::istream &in);
1678  int,
1679  int,
1680  << "You are trying to access the matrix entry with index <"
1681  << arg1 << ',' << arg2
1682  << ">, but this entry does not exist in the sparsity pattern "
1683  "of this matrix."
1684  "\n\n"
1685  "The most common cause for this problem is that you used "
1686  "a method to build the sparsity pattern that did not "
1687  "(completely) take into account all of the entries you "
1688  "will later try to write into. An example would be "
1689  "building a sparsity pattern that does not include "
1690  "the entries you will write into due to constraints "
1691  "on degrees of freedom such as hanging nodes or periodic "
1692  "boundary conditions. In such cases, building the "
1693  "sparsity pattern will succeed, but you will get errors "
1694  "such as the current one at one point or other when "
1695  "trying to write into the entries of the matrix.");
1700  "When copying one sparse matrix into another, "
1701  "or when adding one sparse matrix to another, "
1702  "both matrices need to refer to the same "
1703  "sparsity pattern.");
1708  int,
1709  int,
1710  << "The iterators denote a range of " << arg1
1711  << " elements, but the given number of rows was " << arg2);
1716  "You are attempting an operation on two vectors that "
1717  "are the same object, but the operation requires that the "
1718  "two objects are in fact different.");
1721 protected:
1732  void
1734 
1739  void
1741 
1742 private:
1749 
1757  std::unique_ptr<number[]> val;
1758 
1765  std::size_t max_len;
1766 
1767  // make all other sparse matrices friends
1768  template <typename somenumber>
1769  friend class SparseMatrix;
1770  template <typename somenumber>
1772  template <typename>
1773  friend class SparseILU;
1774 
1775  // To allow it calling private prepare_add() and prepare_set().
1776  template <typename>
1777  friend class BlockMatrixBase;
1778 
1779  // Also give access to internal details to the iterator/accessor classes.
1780  template <typename, bool>
1782  template <typename, bool>
1784 
1785 #ifdef DEAL_II_WITH_MPI
1786  // Give access to internal datastructures to perform MPI operations.
1787  template <typename Number>
1788  friend void
1790  const MPI_Comm,
1792 #endif
1793 };
1794 
1795 #ifndef DOXYGEN
1796 /*---------------------- Inline functions -----------------------------------*/
1797 
1798 
1799 
1800 template <typename number>
1801 template <typename number2>
1802 void
1804 {
1805  this->reinit(sparse_matrix.get_sparsity_pattern());
1806 }
1807 
1808 
1809 
1810 template <typename number>
1811 inline typename SparseMatrix<number>::size_type
1813 {
1814  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1815  return cols->rows;
1816 }
1817 
1818 
1819 template <typename number>
1820 inline typename SparseMatrix<number>::size_type
1822 {
1823  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1824  return cols->cols;
1825 }
1826 
1827 
1828 // Inline the set() and add() functions, since they will be called frequently.
1829 template <typename number>
1830 inline void
1832  const size_type j,
1833  const number value)
1834 {
1836 
1837  const size_type index = cols->operator()(i, j);
1838 
1839  // it is allowed to set elements of the matrix that are not part of the
1840  // sparsity pattern, if the value to which we set it is zero
1842  {
1843  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1844  ExcInvalidIndex(i, j));
1845  return;
1846  }
1847 
1848  val[index] = value;
1849 }
1850 
1851 
1852 
1853 template <typename number>
1854 template <typename number2>
1855 inline void
1856 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1857  const FullMatrix<number2> &values,
1858  const bool elide_zero_values)
1859 {
1860  Assert(indices.size() == values.m(),
1861  ExcDimensionMismatch(indices.size(), values.m()));
1862  Assert(values.m() == values.n(), ExcNotQuadratic());
1863 
1864  for (size_type i = 0; i < indices.size(); ++i)
1865  set(indices[i],
1866  indices.size(),
1867  indices.data(),
1868  &values(i, 0),
1869  elide_zero_values);
1870 }
1871 
1872 
1873 
1874 template <typename number>
1875 template <typename number2>
1876 inline void
1877 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1878  const std::vector<size_type> &col_indices,
1879  const FullMatrix<number2> &values,
1880  const bool elide_zero_values)
1881 {
1882  Assert(row_indices.size() == values.m(),
1883  ExcDimensionMismatch(row_indices.size(), values.m()));
1884  Assert(col_indices.size() == values.n(),
1885  ExcDimensionMismatch(col_indices.size(), values.n()));
1886 
1887  for (size_type i = 0; i < row_indices.size(); ++i)
1888  set(row_indices[i],
1889  col_indices.size(),
1890  col_indices.data(),
1891  &values(i, 0),
1892  elide_zero_values);
1893 }
1894 
1895 
1896 
1897 template <typename number>
1898 template <typename number2>
1899 inline void
1901  const std::vector<size_type> &col_indices,
1902  const std::vector<number2> &values,
1903  const bool elide_zero_values)
1904 {
1905  Assert(col_indices.size() == values.size(),
1906  ExcDimensionMismatch(col_indices.size(), values.size()));
1907 
1908  set(row,
1909  col_indices.size(),
1910  col_indices.data(),
1911  values.data(),
1912  elide_zero_values);
1913 }
1914 
1915 
1916 
1917 template <typename number>
1918 inline void
1920  const size_type j,
1921  const number value)
1922 {
1924 
1925  if (value == number())
1926  return;
1927 
1928  const size_type index = cols->operator()(i, j);
1929 
1930  // it is allowed to add elements to the matrix that are not part of the
1931  // sparsity pattern, if the value to which we set it is zero
1933  {
1934  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1935  ExcInvalidIndex(i, j));
1936  return;
1937  }
1938 
1939  val[index] += value;
1940 }
1941 
1942 
1943 
1944 template <typename number>
1945 template <typename number2>
1946 inline void
1947 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1948  const FullMatrix<number2> &values,
1949  const bool elide_zero_values)
1950 {
1951  Assert(indices.size() == values.m(),
1952  ExcDimensionMismatch(indices.size(), values.m()));
1953  Assert(values.m() == values.n(), ExcNotQuadratic());
1954 
1955  for (size_type i = 0; i < indices.size(); ++i)
1956  add(indices[i],
1957  indices.size(),
1958  indices.data(),
1959  &values(i, 0),
1960  elide_zero_values);
1961 }
1962 
1963 
1964 
1965 template <typename number>
1966 template <typename number2>
1967 inline void
1968 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1969  const std::vector<size_type> &col_indices,
1970  const FullMatrix<number2> &values,
1971  const bool elide_zero_values)
1972 {
1973  Assert(row_indices.size() == values.m(),
1974  ExcDimensionMismatch(row_indices.size(), values.m()));
1975  Assert(col_indices.size() == values.n(),
1976  ExcDimensionMismatch(col_indices.size(), values.n()));
1977 
1978  for (size_type i = 0; i < row_indices.size(); ++i)
1979  add(row_indices[i],
1980  col_indices.size(),
1981  col_indices.data(),
1982  &values(i, 0),
1983  elide_zero_values);
1984 }
1985 
1986 
1987 
1988 template <typename number>
1989 template <typename number2>
1990 inline void
1992  const std::vector<size_type> &col_indices,
1993  const std::vector<number2> &values,
1994  const bool elide_zero_values)
1995 {
1996  Assert(col_indices.size() == values.size(),
1997  ExcDimensionMismatch(col_indices.size(), values.size()));
1998 
1999  add(row,
2000  col_indices.size(),
2001  col_indices.data(),
2002  values.data(),
2003  elide_zero_values,
2004  std::is_sorted(col_indices.begin(), col_indices.end()));
2005 }
2006 
2007 
2008 
2009 template <typename number>
2010 inline SparseMatrix<number> &
2011 SparseMatrix<number>::operator*=(const number factor)
2012 {
2013  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2014  Assert(val != nullptr, ExcNotInitialized());
2015 
2016  number *val_ptr = val.get();
2017  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2018 
2019  while (val_ptr != end_ptr)
2020  *val_ptr++ *= factor;
2021 
2022  return *this;
2023 }
2024 
2025 
2026 
2027 template <typename number>
2028 inline SparseMatrix<number> &
2029 SparseMatrix<number>::operator/=(const number factor)
2030 {
2031  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2032  Assert(val != nullptr, ExcNotInitialized());
2033  Assert(factor != number(), ExcDivideByZero());
2034 
2035  const number factor_inv = number(1.) / factor;
2036 
2037  number *val_ptr = val.get();
2038  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2039 
2040  while (val_ptr != end_ptr)
2041  *val_ptr++ *= factor_inv;
2042 
2043  return *this;
2044 }
2045 
2046 
2047 
2048 template <typename number>
2049 inline const number &
2051 {
2052  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2053  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2054  ExcInvalidIndex(i, j));
2055  return val[cols->operator()(i, j)];
2056 }
2057 
2058 
2059 
2060 template <typename number>
2061 inline number &
2063 {
2064  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2065  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2066  ExcInvalidIndex(i, j));
2067  return val[cols->operator()(i, j)];
2068 }
2069 
2070 
2071 
2072 template <typename number>
2073 inline number
2074 SparseMatrix<number>::el(const size_type i, const size_type j) const
2075 {
2076  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2077  const size_type index = cols->operator()(i, j);
2078 
2080  return val[index];
2081  else
2082  return 0;
2083 }
2084 
2085 
2086 
2087 template <typename number>
2088 inline number
2090 {
2091  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2092  Assert(m() == n(), ExcNotQuadratic());
2093  AssertIndexRange(i, m());
2094 
2095  // Use that the first element in each row of a quadratic matrix is the main
2096  // diagonal
2097  return val[cols->rowstart[i]];
2098 }
2099 
2100 
2101 
2102 template <typename number>
2103 inline number &
2105 {
2106  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2107  Assert(m() == n(), ExcNotQuadratic());
2108  AssertIndexRange(i, m());
2109 
2110  // Use that the first element in each row of a quadratic matrix is the main
2111  // diagonal
2112  return val[cols->rowstart[i]];
2113 }
2114 
2115 
2116 
2117 template <typename number>
2118 template <typename ForwardIterator>
2119 void
2120 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2121  const ForwardIterator end)
2122 {
2123  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2124  ExcIteratorRange(std::distance(begin, end), m()));
2125 
2126  // for use in the inner loop, we define an alias to the type of the inner
2127  // iterators
2128  using inner_iterator =
2129  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2130  size_type row = 0;
2131  for (ForwardIterator i = begin; i != end; ++i, ++row)
2132  {
2133  const inner_iterator end_of_row = i->end();
2134  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2135  // write entries
2136  set(row, j->first, j->second);
2137  };
2138 }
2139 
2140 
2141 //---------------------------------------------------------------------------
2142 
2143 
2144 namespace SparseMatrixIterators
2145 {
2146  template <typename number>
2147  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2148  const std::size_t index_within_matrix)
2149  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2150  index_within_matrix)
2151  , matrix(matrix)
2152  {}
2153 
2154 
2155 
2156  template <typename number>
2157  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2158  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2159  , matrix(matrix)
2160  {}
2161 
2162 
2163 
2164  template <typename number>
2165  inline Accessor<number, true>::Accessor(
2167  : SparsityPatternIterators::Accessor(a)
2168  , matrix(&a.get_matrix())
2169  {}
2170 
2171 
2172 
2173  template <typename number>
2174  inline number
2175  Accessor<number, true>::value() const
2176  {
2177  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2178  return matrix->val[linear_index];
2179  }
2180 
2181 
2182 
2183  template <typename number>
2184  inline const typename Accessor<number, true>::MatrixType &
2185  Accessor<number, true>::get_matrix() const
2186  {
2187  return *matrix;
2188  }
2189 
2190 
2191 
2192  template <typename number>
2193  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2194  const bool)
2195  : accessor(accessor)
2196  {}
2197 
2198 
2199  template <typename number>
2200  inline Accessor<number, false>::Reference::operator number() const
2201  {
2202  AssertIndexRange(accessor->linear_index,
2203  accessor->matrix->n_nonzero_elements());
2204  return accessor->matrix->val[accessor->linear_index];
2205  }
2206 
2207 
2208 
2209  template <typename number>
2210  inline const typename Accessor<number, false>::Reference &
2211  Accessor<number, false>::Reference::operator=(const number n) const
2212  {
2213  AssertIndexRange(accessor->linear_index,
2214  accessor->matrix->n_nonzero_elements());
2215  accessor->matrix->val[accessor->linear_index] = n;
2216  return *this;
2217  }
2218 
2219 
2220 
2221  template <typename number>
2222  inline const typename Accessor<number, false>::Reference &
2223  Accessor<number, false>::Reference::operator+=(const number n) const
2224  {
2225  AssertIndexRange(accessor->linear_index,
2226  accessor->matrix->n_nonzero_elements());
2227  accessor->matrix->val[accessor->linear_index] += n;
2228  return *this;
2229  }
2230 
2231 
2232 
2233  template <typename number>
2234  inline const typename Accessor<number, false>::Reference &
2235  Accessor<number, false>::Reference::operator-=(const number n) const
2236  {
2237  AssertIndexRange(accessor->linear_index,
2238  accessor->matrix->n_nonzero_elements());
2239  accessor->matrix->val[accessor->linear_index] -= n;
2240  return *this;
2241  }
2242 
2243 
2244 
2245  template <typename number>
2246  inline const typename Accessor<number, false>::Reference &
2247  Accessor<number, false>::Reference::operator*=(const number n) const
2248  {
2249  AssertIndexRange(accessor->linear_index,
2250  accessor->matrix->n_nonzero_elements());
2251  accessor->matrix->val[accessor->linear_index] *= n;
2252  return *this;
2253  }
2254 
2255 
2256 
2257  template <typename number>
2258  inline const typename Accessor<number, false>::Reference &
2259  Accessor<number, false>::Reference::operator/=(const number n) const
2260  {
2261  AssertIndexRange(accessor->linear_index,
2262  accessor->matrix->n_nonzero_elements());
2263  accessor->matrix->val[accessor->linear_index] /= n;
2264  return *this;
2265  }
2266 
2267 
2268 
2269  template <typename number>
2270  inline Accessor<number, false>::Accessor(MatrixType *matrix,
2271  const std::size_t index)
2272  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2273  , matrix(matrix)
2274  {}
2275 
2276 
2277 
2278  template <typename number>
2279  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2280  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2281  , matrix(matrix)
2282  {}
2283 
2284 
2285 
2286  template <typename number>
2287  inline typename Accessor<number, false>::Reference
2288  Accessor<number, false>::value() const
2289  {
2290  return Reference(this, true);
2291  }
2292 
2293 
2294 
2295  template <typename number>
2296  inline typename Accessor<number, false>::MatrixType &
2297  Accessor<number, false>::get_matrix() const
2298  {
2299  return *matrix;
2300  }
2301 
2302 
2303 
2304  template <typename number, bool Constness>
2305  inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
2306  const std::size_t index)
2307  : accessor(matrix, index)
2308  {}
2309 
2310 
2311 
2312  template <typename number, bool Constness>
2313  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2314  : accessor(matrix)
2315  {}
2316 
2317 
2318 
2319  template <typename number, bool Constness>
2320  inline Iterator<number, Constness>::Iterator(
2322  : accessor(*i)
2323  {}
2324 
2325 
2326 
2327  template <typename number, bool Constness>
2328  inline const Iterator<number, Constness> &
2329  Iterator<number, Constness>::operator=(
2331  {
2332  accessor = *i;
2333  return *this;
2334  }
2335 
2336 
2337 
2338  template <typename number, bool Constness>
2339  inline Iterator<number, Constness> &
2341  {
2342  accessor.advance();
2343  return *this;
2344  }
2345 
2346 
2347  template <typename number, bool Constness>
2348  inline Iterator<number, Constness>
2350  {
2351  const Iterator iter = *this;
2352  accessor.advance();
2353  return iter;
2354  }
2355 
2356 
2357  template <typename number, bool Constness>
2358  inline const Accessor<number, Constness> &
2360  {
2361  return accessor;
2362  }
2363 
2364 
2365  template <typename number, bool Constness>
2366  inline const Accessor<number, Constness> *
2367  Iterator<number, Constness>::operator->() const
2368  {
2369  return &accessor;
2370  }
2371 
2372 
2373  template <typename number, bool Constness>
2374  inline bool
2375  Iterator<number, Constness>::operator==(const Iterator &other) const
2376  {
2377  return (accessor == other.accessor);
2378  }
2379 
2380 
2381  template <typename number, bool Constness>
2382  inline bool
2383  Iterator<number, Constness>::operator!=(const Iterator &other) const
2384  {
2385  return !(*this == other);
2386  }
2387 
2388 
2389  template <typename number, bool Constness>
2390  inline bool
2391  Iterator<number, Constness>::operator<(const Iterator &other) const
2392  {
2393  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2394  ExcInternalError());
2395 
2396  return (accessor < other.accessor);
2397  }
2398 
2399 
2400  template <typename number, bool Constness>
2401  inline bool
2402  Iterator<number, Constness>::operator>(const Iterator &other) const
2403  {
2404  return (other < *this);
2405  }
2406 
2407 
2408  template <typename number, bool Constness>
2409  inline int
2410  Iterator<number, Constness>::operator-(const Iterator &other) const
2411  {
2412  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2413  ExcInternalError());
2414 
2415  return (*this)->linear_index - other->linear_index;
2416  }
2417 
2418 
2419 
2420  template <typename number, bool Constness>
2421  inline Iterator<number, Constness>
2423  {
2424  Iterator x = *this;
2425  for (size_type i = 0; i < n; ++i)
2426  ++x;
2427 
2428  return x;
2429  }
2430 
2431 } // namespace SparseMatrixIterators
2432 
2433 
2434 
2435 template <typename number>
2438 {
2439  return const_iterator(this, 0);
2440 }
2441 
2442 
2443 template <typename number>
2446 {
2447  return const_iterator(this);
2448 }
2449 
2450 
2451 template <typename number>
2452 inline typename SparseMatrix<number>::iterator
2454 {
2455  return iterator(this, 0);
2456 }
2457 
2458 
2459 template <typename number>
2460 inline typename SparseMatrix<number>::iterator
2462 {
2463  return iterator(this, cols->rowstart[cols->rows]);
2464 }
2465 
2466 
2467 template <typename number>
2470 {
2471  AssertIndexRange(r, m());
2472 
2473  return const_iterator(this, cols->rowstart[r]);
2474 }
2475 
2476 
2477 
2478 template <typename number>
2480 SparseMatrix<number>::end(const size_type r) const
2481 {
2482  AssertIndexRange(r, m());
2483 
2484  return const_iterator(this, cols->rowstart[r + 1]);
2485 }
2486 
2487 
2488 
2489 template <typename number>
2490 inline typename SparseMatrix<number>::iterator
2492 {
2493  AssertIndexRange(r, m());
2494 
2495  return iterator(this, cols->rowstart[r]);
2496 }
2497 
2498 
2499 
2500 template <typename number>
2501 inline typename SparseMatrix<number>::iterator
2503 {
2504  AssertIndexRange(r, m());
2505 
2506  return iterator(this, cols->rowstart[r + 1]);
2507 }
2508 
2509 
2510 
2511 template <typename number>
2512 template <typename StreamType>
2513 inline void
2514 SparseMatrix<number>::print(StreamType &out,
2515  const bool across,
2516  const bool diagonal_first) const
2517 {
2518  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2519  Assert(val != nullptr, ExcNotInitialized());
2520 
2521  bool hanging_diagonal = false;
2522  number diagonal = number();
2523 
2524  for (size_type i = 0; i < cols->rows; ++i)
2525  {
2526  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2527  {
2528  if (!diagonal_first && i == cols->colnums[j])
2529  {
2530  diagonal = val[j];
2531  hanging_diagonal = true;
2532  }
2533  else
2534  {
2535  if (hanging_diagonal && cols->colnums[j] > i)
2536  {
2537  if (across)
2538  out << ' ' << i << ',' << i << ':' << diagonal;
2539  else
2540  out << '(' << i << ',' << i << ") " << diagonal
2541  << std::endl;
2542  hanging_diagonal = false;
2543  }
2544  if (across)
2545  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2546  else
2547  out << '(' << i << ',' << cols->colnums[j] << ") " << val[j]
2548  << std::endl;
2549  }
2550  }
2551  if (hanging_diagonal)
2552  {
2553  if (across)
2554  out << ' ' << i << ',' << i << ':' << diagonal;
2555  else
2556  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2557  hanging_diagonal = false;
2558  }
2559  }
2560  if (across)
2561  out << std::endl;
2562 }
2563 
2564 
2565 template <typename number>
2566 inline void
2568 {
2569  // nothing to do here
2570 }
2571 
2572 
2573 
2574 template <typename number>
2575 inline void
2577 {
2578  // nothing to do here
2579 }
2580 
2581 #endif // DOXYGEN
2582 
2583 
2585 
2586 #endif
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
const Reference & operator+=(const number n) const
Reference(const Accessor *accessor, const bool dummy)
const Reference & operator=(const number n) const
const Reference & operator/=(const number n) const
const Reference & operator*=(const number n) const
const Reference & operator-=(const number n) const
Accessor(MatrixType *matrix, const std::size_t index)
Accessor(const SparseMatrixIterators::Accessor< number, false > &a)
Accessor(MatrixType *matrix, const std::size_t index_within_matrix)
const SparseMatrix< number > & get_matrix() const
bool operator>(const Iterator &) const
bool operator==(const Iterator &) const
int operator-(const Iterator &p) const
bool operator<(const Iterator &) const
Iterator(MatrixType *matrix)
const Accessor< number, Constness > & value_type
Iterator operator+(const size_type n) const
Accessor< number, Constness > accessor
Iterator(const SparseMatrixIterators::Iterator< number, false > &i)
const Accessor< number, Constness > * operator->() const
Iterator(MatrixType *matrix, const std::size_t index_within_matrix)
const Iterator< number, Constness > & operator=(const SparseMatrixIterators::Iterator< number, false > &i)
typename Accessor< number, Constness >::MatrixType MatrixType
bool operator!=(const Iterator &) const
const Accessor< number, Constness > & operator*() const
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number2 > &values, const bool elide_zero_values=false)
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
SparseMatrix(const SparsityPattern &sparsity)
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
void add(const number factor, const SparseMatrix< somenumber > &matrix)
number & diag_element(const size_type i)
std::size_t n_nonzero_elements() const
size_type get_row_length(const size_type row) const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
iterator end()
void prepare_add()
void Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
void Tvmult(OutVector &dst, const InVector &src) const
const_iterator begin(const size_type r) const
const_iterator end() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
void set(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
void print_as_numpy_arrays(std::ostream &out, const unsigned int precision=9) const
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
const number & operator()(const size_type i, const size_type j) const
void set(const size_type i, const size_type j, const number value)
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 char *separator=" ") const
const_iterator begin() const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
virtual void clear()
virtual ~SparseMatrix() override
void vmult_add(OutVector &dst, const InVector &src) const
number diag_element(const size_type i) const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=true)
void SSOR(Vector< somenumber > &v, const number omega=1.) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
SparseMatrix(SparseMatrix< number > &&m) noexcept
void TSOR(Vector< somenumber > &v, const number omega=1.) const
SparseMatrix< number > & copy_from(const TrilinosWrappers::SparseMatrix &matrix)
void symmetrize()
SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id)
void print_pattern(std::ostream &out, const double threshold=0.) const
void vmult(OutVector &dst, const InVector &src) const
std::size_t memory_consumption() const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number omega=1.) const
SparseMatrix & operator=(const double d)
void block_write(std::ostream &out) const
void copy_from(const ForwardIterator begin, const ForwardIterator end)
SparseMatrix & operator/=(const number factor)
void block_read(std::istream &in)
void prepare_set()
number el(const size_type i, const size_type j) const
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
const SparsityPattern & get_sparsity_pattern() const
const_iterator end(const size_type r) const
iterator begin()
void SOR(Vector< somenumber > &v, const number omega=1.) const
std::unique_ptr< number[]> val
iterator begin(const size_type r)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
typename numbers::NumberTraits< number >::real_type real_type
types::global_dof_index size_type
void reinit(const SparseMatrix< number2 > &sparse_matrix)
number & operator()(const size_type i, const size_type j)
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
void add(const size_type i, const size_type j, const number value)
size_type n() const
size_type m() const
void copy_from(const FullMatrix< somenumber > &matrix)
std::size_t max_len
iterator end(const size_type r)
SparseMatrix< number > & operator=(SparseMatrix< number > &&m) noexcept
void compress(VectorOperation::values)
SparseMatrix(const SparseMatrix &)
real_type frobenius_norm() const
void Tvmult_add(OutVector &dst, const InVector &src) const
real_type linfty_norm() const
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number2 > &values, const bool elide_zero_values=true)
number value_type
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
real_type l1_norm() const
SparseMatrix & operator*=(const number factor)
SparseMatrix< number > & operator=(const IdentityMatrix &id)
std::size_t n_actually_nonzero_elements(const double threshold=0.) const
bool empty() const
void add(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=true)
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number omega=1.) const
virtual void reinit(const SparsityPattern &sparsity)
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=false)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
static constexpr size_type invalid_entry
Definition: vector.h:110
std::enable_if_t< std::is_floating_point_v< T > &&std::is_floating_point_v< U >, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNeedsSparsityPattern()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcDifferentSparsityPatterns()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
#define AssertIsFinite(number)
Definition: exceptions.h:1915
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:540
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:495
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
static ::ExceptionBase & ExcNotQuadratic()
Expression operator>(const Expression &lhs, const Expression &rhs)
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
static const char V
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
types::global_dof_index size_type
Definition: sparse_matrix.h:80
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
unsigned int global_dof_index
Definition: types.h:82
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
static const bool zero_addition_can_be_elided
typename ::SparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
typename ::SparseMatrixIterators::Iterator< number, Constness >::value_type value_type
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)