Reference documentation for deal.II version GIT c00406db16 2023-09-28 18:00: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 
1606  void
1607  print_formatted(std::ostream &out,
1608  const unsigned int precision = 3,
1609  const bool scientific = true,
1610  const unsigned int width = 0,
1611  const char *zero_string = " ",
1612  const double denominator = 1.) const;
1613 
1619  void
1620  print_pattern(std::ostream &out, const double threshold = 0.) const;
1621 
1630  void
1631  print_as_numpy_arrays(std::ostream &out,
1632  const unsigned int precision = 9) const;
1633 
1644  void
1645  block_write(std::ostream &out) const;
1646 
1663  void
1664  block_read(std::istream &in);
1675  int,
1676  int,
1677  << "You are trying to access the matrix entry with index <"
1678  << arg1 << ',' << arg2
1679  << ">, but this entry does not exist in the sparsity pattern "
1680  "of this matrix."
1681  "\n\n"
1682  "The most common cause for this problem is that you used "
1683  "a method to build the sparsity pattern that did not "
1684  "(completely) take into account all of the entries you "
1685  "will later try to write into. An example would be "
1686  "building a sparsity pattern that does not include "
1687  "the entries you will write into due to constraints "
1688  "on degrees of freedom such as hanging nodes or periodic "
1689  "boundary conditions. In such cases, building the "
1690  "sparsity pattern will succeed, but you will get errors "
1691  "such as the current one at one point or other when "
1692  "trying to write into the entries of the matrix.");
1697  "When copying one sparse matrix into another, "
1698  "or when adding one sparse matrix to another, "
1699  "both matrices need to refer to the same "
1700  "sparsity pattern.");
1705  int,
1706  int,
1707  << "The iterators denote a range of " << arg1
1708  << " elements, but the given number of rows was " << arg2);
1713  "You are attempting an operation on two vectors that "
1714  "are the same object, but the operation requires that the "
1715  "two objects are in fact different.");
1718 protected:
1729  void
1731 
1736  void
1738 
1739 private:
1746 
1754  std::unique_ptr<number[]> val;
1755 
1762  std::size_t max_len;
1763 
1764  // make all other sparse matrices friends
1765  template <typename somenumber>
1766  friend class SparseMatrix;
1767  template <typename somenumber>
1769  template <typename>
1770  friend class SparseILU;
1771 
1772  // To allow it calling private prepare_add() and prepare_set().
1773  template <typename>
1774  friend class BlockMatrixBase;
1775 
1776  // Also give access to internal details to the iterator/accessor classes.
1777  template <typename, bool>
1779  template <typename, bool>
1781 
1782 #ifdef DEAL_II_WITH_MPI
1783  // Give access to internal datastructures to perform MPI operations.
1784  template <typename Number>
1785  friend void
1787  const MPI_Comm,
1789 #endif
1790 };
1791 
1792 #ifndef DOXYGEN
1793 /*---------------------- Inline functions -----------------------------------*/
1794 
1795 
1796 
1797 template <typename number>
1798 template <typename number2>
1799 void
1801 {
1802  this->reinit(sparse_matrix.get_sparsity_pattern());
1803 }
1804 
1805 
1806 
1807 template <typename number>
1808 inline typename SparseMatrix<number>::size_type
1810 {
1811  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1812  return cols->rows;
1813 }
1814 
1815 
1816 template <typename number>
1817 inline typename SparseMatrix<number>::size_type
1819 {
1820  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1821  return cols->cols;
1822 }
1823 
1824 
1825 // Inline the set() and add() functions, since they will be called frequently.
1826 template <typename number>
1827 inline void
1829  const size_type j,
1830  const number value)
1831 {
1833 
1834  const size_type index = cols->operator()(i, j);
1835 
1836  // it is allowed to set elements of the matrix that are not part of the
1837  // sparsity pattern, if the value to which we set it is zero
1839  {
1840  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1841  ExcInvalidIndex(i, j));
1842  return;
1843  }
1844 
1845  val[index] = value;
1846 }
1847 
1848 
1849 
1850 template <typename number>
1851 template <typename number2>
1852 inline void
1853 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1854  const FullMatrix<number2> &values,
1855  const bool elide_zero_values)
1856 {
1857  Assert(indices.size() == values.m(),
1858  ExcDimensionMismatch(indices.size(), values.m()));
1859  Assert(values.m() == values.n(), ExcNotQuadratic());
1860 
1861  for (size_type i = 0; i < indices.size(); ++i)
1862  set(indices[i],
1863  indices.size(),
1864  indices.data(),
1865  &values(i, 0),
1866  elide_zero_values);
1867 }
1868 
1869 
1870 
1871 template <typename number>
1872 template <typename number2>
1873 inline void
1874 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1875  const std::vector<size_type> &col_indices,
1876  const FullMatrix<number2> &values,
1877  const bool elide_zero_values)
1878 {
1879  Assert(row_indices.size() == values.m(),
1880  ExcDimensionMismatch(row_indices.size(), values.m()));
1881  Assert(col_indices.size() == values.n(),
1882  ExcDimensionMismatch(col_indices.size(), values.n()));
1883 
1884  for (size_type i = 0; i < row_indices.size(); ++i)
1885  set(row_indices[i],
1886  col_indices.size(),
1887  col_indices.data(),
1888  &values(i, 0),
1889  elide_zero_values);
1890 }
1891 
1892 
1893 
1894 template <typename number>
1895 template <typename number2>
1896 inline void
1898  const std::vector<size_type> &col_indices,
1899  const std::vector<number2> &values,
1900  const bool elide_zero_values)
1901 {
1902  Assert(col_indices.size() == values.size(),
1903  ExcDimensionMismatch(col_indices.size(), values.size()));
1904 
1905  set(row,
1906  col_indices.size(),
1907  col_indices.data(),
1908  values.data(),
1909  elide_zero_values);
1910 }
1911 
1912 
1913 
1914 template <typename number>
1915 inline void
1917  const size_type j,
1918  const number value)
1919 {
1921 
1922  if (value == number())
1923  return;
1924 
1925  const size_type index = cols->operator()(i, j);
1926 
1927  // it is allowed to add elements to the matrix that are not part of the
1928  // sparsity pattern, if the value to which we set it is zero
1930  {
1931  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1932  ExcInvalidIndex(i, j));
1933  return;
1934  }
1935 
1936  val[index] += value;
1937 }
1938 
1939 
1940 
1941 template <typename number>
1942 template <typename number2>
1943 inline void
1944 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1945  const FullMatrix<number2> &values,
1946  const bool elide_zero_values)
1947 {
1948  Assert(indices.size() == values.m(),
1949  ExcDimensionMismatch(indices.size(), values.m()));
1950  Assert(values.m() == values.n(), ExcNotQuadratic());
1951 
1952  for (size_type i = 0; i < indices.size(); ++i)
1953  add(indices[i],
1954  indices.size(),
1955  indices.data(),
1956  &values(i, 0),
1957  elide_zero_values);
1958 }
1959 
1960 
1961 
1962 template <typename number>
1963 template <typename number2>
1964 inline void
1965 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1966  const std::vector<size_type> &col_indices,
1967  const FullMatrix<number2> &values,
1968  const bool elide_zero_values)
1969 {
1970  Assert(row_indices.size() == values.m(),
1971  ExcDimensionMismatch(row_indices.size(), values.m()));
1972  Assert(col_indices.size() == values.n(),
1973  ExcDimensionMismatch(col_indices.size(), values.n()));
1974 
1975  for (size_type i = 0; i < row_indices.size(); ++i)
1976  add(row_indices[i],
1977  col_indices.size(),
1978  col_indices.data(),
1979  &values(i, 0),
1980  elide_zero_values);
1981 }
1982 
1983 
1984 
1985 template <typename number>
1986 template <typename number2>
1987 inline void
1989  const std::vector<size_type> &col_indices,
1990  const std::vector<number2> &values,
1991  const bool elide_zero_values)
1992 {
1993  Assert(col_indices.size() == values.size(),
1994  ExcDimensionMismatch(col_indices.size(), values.size()));
1995 
1996  add(row,
1997  col_indices.size(),
1998  col_indices.data(),
1999  values.data(),
2000  elide_zero_values,
2001  std::is_sorted(col_indices.begin(), col_indices.end()));
2002 }
2003 
2004 
2005 
2006 template <typename number>
2007 inline SparseMatrix<number> &
2008 SparseMatrix<number>::operator*=(const number factor)
2009 {
2010  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2011  Assert(val != nullptr, ExcNotInitialized());
2012 
2013  number *val_ptr = val.get();
2014  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2015 
2016  while (val_ptr != end_ptr)
2017  *val_ptr++ *= factor;
2018 
2019  return *this;
2020 }
2021 
2022 
2023 
2024 template <typename number>
2025 inline SparseMatrix<number> &
2026 SparseMatrix<number>::operator/=(const number factor)
2027 {
2028  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2029  Assert(val != nullptr, ExcNotInitialized());
2030  Assert(factor != number(), ExcDivideByZero());
2031 
2032  const number factor_inv = number(1.) / factor;
2033 
2034  number *val_ptr = val.get();
2035  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2036 
2037  while (val_ptr != end_ptr)
2038  *val_ptr++ *= factor_inv;
2039 
2040  return *this;
2041 }
2042 
2043 
2044 
2045 template <typename number>
2046 inline const number &
2048 {
2049  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2050  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2051  ExcInvalidIndex(i, j));
2052  return val[cols->operator()(i, j)];
2053 }
2054 
2055 
2056 
2057 template <typename number>
2058 inline number &
2060 {
2061  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2062  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2063  ExcInvalidIndex(i, j));
2064  return val[cols->operator()(i, j)];
2065 }
2066 
2067 
2068 
2069 template <typename number>
2070 inline number
2071 SparseMatrix<number>::el(const size_type i, const size_type j) const
2072 {
2073  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2074  const size_type index = cols->operator()(i, j);
2075 
2077  return val[index];
2078  else
2079  return 0;
2080 }
2081 
2082 
2083 
2084 template <typename number>
2085 inline number
2087 {
2088  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2089  Assert(m() == n(), ExcNotQuadratic());
2090  AssertIndexRange(i, m());
2091 
2092  // Use that the first element in each row of a quadratic matrix is the main
2093  // diagonal
2094  return val[cols->rowstart[i]];
2095 }
2096 
2097 
2098 
2099 template <typename number>
2100 inline number &
2102 {
2103  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2104  Assert(m() == n(), ExcNotQuadratic());
2105  AssertIndexRange(i, m());
2106 
2107  // Use that the first element in each row of a quadratic matrix is the main
2108  // diagonal
2109  return val[cols->rowstart[i]];
2110 }
2111 
2112 
2113 
2114 template <typename number>
2115 template <typename ForwardIterator>
2116 void
2117 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2118  const ForwardIterator end)
2119 {
2120  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2121  ExcIteratorRange(std::distance(begin, end), m()));
2122 
2123  // for use in the inner loop, we define an alias to the type of the inner
2124  // iterators
2125  using inner_iterator =
2126  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2127  size_type row = 0;
2128  for (ForwardIterator i = begin; i != end; ++i, ++row)
2129  {
2130  const inner_iterator end_of_row = i->end();
2131  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2132  // write entries
2133  set(row, j->first, j->second);
2134  };
2135 }
2136 
2137 
2138 //---------------------------------------------------------------------------
2139 
2140 
2141 namespace SparseMatrixIterators
2142 {
2143  template <typename number>
2144  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2145  const std::size_t index_within_matrix)
2146  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2147  index_within_matrix)
2148  , matrix(matrix)
2149  {}
2150 
2151 
2152 
2153  template <typename number>
2154  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2155  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2156  , matrix(matrix)
2157  {}
2158 
2159 
2160 
2161  template <typename number>
2162  inline Accessor<number, true>::Accessor(
2164  : SparsityPatternIterators::Accessor(a)
2165  , matrix(&a.get_matrix())
2166  {}
2167 
2168 
2169 
2170  template <typename number>
2171  inline number
2172  Accessor<number, true>::value() const
2173  {
2174  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2175  return matrix->val[linear_index];
2176  }
2177 
2178 
2179 
2180  template <typename number>
2181  inline const typename Accessor<number, true>::MatrixType &
2182  Accessor<number, true>::get_matrix() const
2183  {
2184  return *matrix;
2185  }
2186 
2187 
2188 
2189  template <typename number>
2190  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2191  const bool)
2192  : accessor(accessor)
2193  {}
2194 
2195 
2196  template <typename number>
2197  inline Accessor<number, false>::Reference::operator number() const
2198  {
2199  AssertIndexRange(accessor->linear_index,
2200  accessor->matrix->n_nonzero_elements());
2201  return accessor->matrix->val[accessor->linear_index];
2202  }
2203 
2204 
2205 
2206  template <typename number>
2207  inline const typename Accessor<number, false>::Reference &
2208  Accessor<number, false>::Reference::operator=(const number n) const
2209  {
2210  AssertIndexRange(accessor->linear_index,
2211  accessor->matrix->n_nonzero_elements());
2212  accessor->matrix->val[accessor->linear_index] = n;
2213  return *this;
2214  }
2215 
2216 
2217 
2218  template <typename number>
2219  inline const typename Accessor<number, false>::Reference &
2220  Accessor<number, false>::Reference::operator+=(const number n) const
2221  {
2222  AssertIndexRange(accessor->linear_index,
2223  accessor->matrix->n_nonzero_elements());
2224  accessor->matrix->val[accessor->linear_index] += n;
2225  return *this;
2226  }
2227 
2228 
2229 
2230  template <typename number>
2231  inline const typename Accessor<number, false>::Reference &
2232  Accessor<number, false>::Reference::operator-=(const number n) const
2233  {
2234  AssertIndexRange(accessor->linear_index,
2235  accessor->matrix->n_nonzero_elements());
2236  accessor->matrix->val[accessor->linear_index] -= n;
2237  return *this;
2238  }
2239 
2240 
2241 
2242  template <typename number>
2243  inline const typename Accessor<number, false>::Reference &
2244  Accessor<number, false>::Reference::operator*=(const number n) const
2245  {
2246  AssertIndexRange(accessor->linear_index,
2247  accessor->matrix->n_nonzero_elements());
2248  accessor->matrix->val[accessor->linear_index] *= n;
2249  return *this;
2250  }
2251 
2252 
2253 
2254  template <typename number>
2255  inline const typename Accessor<number, false>::Reference &
2256  Accessor<number, false>::Reference::operator/=(const number n) const
2257  {
2258  AssertIndexRange(accessor->linear_index,
2259  accessor->matrix->n_nonzero_elements());
2260  accessor->matrix->val[accessor->linear_index] /= n;
2261  return *this;
2262  }
2263 
2264 
2265 
2266  template <typename number>
2267  inline Accessor<number, false>::Accessor(MatrixType *matrix,
2268  const std::size_t index)
2269  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2270  , matrix(matrix)
2271  {}
2272 
2273 
2274 
2275  template <typename number>
2276  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2277  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2278  , matrix(matrix)
2279  {}
2280 
2281 
2282 
2283  template <typename number>
2284  inline typename Accessor<number, false>::Reference
2285  Accessor<number, false>::value() const
2286  {
2287  return Reference(this, true);
2288  }
2289 
2290 
2291 
2292  template <typename number>
2293  inline typename Accessor<number, false>::MatrixType &
2294  Accessor<number, false>::get_matrix() const
2295  {
2296  return *matrix;
2297  }
2298 
2299 
2300 
2301  template <typename number, bool Constness>
2302  inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
2303  const std::size_t index)
2304  : accessor(matrix, index)
2305  {}
2306 
2307 
2308 
2309  template <typename number, bool Constness>
2310  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2311  : accessor(matrix)
2312  {}
2313 
2314 
2315 
2316  template <typename number, bool Constness>
2317  inline Iterator<number, Constness>::Iterator(
2319  : accessor(*i)
2320  {}
2321 
2322 
2323 
2324  template <typename number, bool Constness>
2325  inline const Iterator<number, Constness> &
2326  Iterator<number, Constness>::operator=(
2328  {
2329  accessor = *i;
2330  return *this;
2331  }
2332 
2333 
2334 
2335  template <typename number, bool Constness>
2336  inline Iterator<number, Constness> &
2338  {
2339  accessor.advance();
2340  return *this;
2341  }
2342 
2343 
2344  template <typename number, bool Constness>
2345  inline Iterator<number, Constness>
2347  {
2348  const Iterator iter = *this;
2349  accessor.advance();
2350  return iter;
2351  }
2352 
2353 
2354  template <typename number, bool Constness>
2355  inline const Accessor<number, Constness> &
2357  {
2358  return accessor;
2359  }
2360 
2361 
2362  template <typename number, bool Constness>
2363  inline const Accessor<number, Constness> *
2364  Iterator<number, Constness>::operator->() const
2365  {
2366  return &accessor;
2367  }
2368 
2369 
2370  template <typename number, bool Constness>
2371  inline bool
2372  Iterator<number, Constness>::operator==(const Iterator &other) const
2373  {
2374  return (accessor == other.accessor);
2375  }
2376 
2377 
2378  template <typename number, bool Constness>
2379  inline bool
2380  Iterator<number, Constness>::operator!=(const Iterator &other) const
2381  {
2382  return !(*this == other);
2383  }
2384 
2385 
2386  template <typename number, bool Constness>
2387  inline bool
2388  Iterator<number, Constness>::operator<(const Iterator &other) const
2389  {
2390  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2391  ExcInternalError());
2392 
2393  return (accessor < other.accessor);
2394  }
2395 
2396 
2397  template <typename number, bool Constness>
2398  inline bool
2399  Iterator<number, Constness>::operator>(const Iterator &other) const
2400  {
2401  return (other < *this);
2402  }
2403 
2404 
2405  template <typename number, bool Constness>
2406  inline int
2407  Iterator<number, Constness>::operator-(const Iterator &other) const
2408  {
2409  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2410  ExcInternalError());
2411 
2412  return (*this)->linear_index - other->linear_index;
2413  }
2414 
2415 
2416 
2417  template <typename number, bool Constness>
2418  inline Iterator<number, Constness>
2420  {
2421  Iterator x = *this;
2422  for (size_type i = 0; i < n; ++i)
2423  ++x;
2424 
2425  return x;
2426  }
2427 
2428 } // namespace SparseMatrixIterators
2429 
2430 
2431 
2432 template <typename number>
2435 {
2436  return const_iterator(this, 0);
2437 }
2438 
2439 
2440 template <typename number>
2443 {
2444  return const_iterator(this);
2445 }
2446 
2447 
2448 template <typename number>
2449 inline typename SparseMatrix<number>::iterator
2451 {
2452  return iterator(this, 0);
2453 }
2454 
2455 
2456 template <typename number>
2457 inline typename SparseMatrix<number>::iterator
2459 {
2460  return iterator(this, cols->rowstart[cols->rows]);
2461 }
2462 
2463 
2464 template <typename number>
2467 {
2468  AssertIndexRange(r, m());
2469 
2470  return const_iterator(this, cols->rowstart[r]);
2471 }
2472 
2473 
2474 
2475 template <typename number>
2477 SparseMatrix<number>::end(const size_type r) const
2478 {
2479  AssertIndexRange(r, m());
2480 
2481  return const_iterator(this, cols->rowstart[r + 1]);
2482 }
2483 
2484 
2485 
2486 template <typename number>
2487 inline typename SparseMatrix<number>::iterator
2489 {
2490  AssertIndexRange(r, m());
2491 
2492  return iterator(this, cols->rowstart[r]);
2493 }
2494 
2495 
2496 
2497 template <typename number>
2498 inline typename SparseMatrix<number>::iterator
2500 {
2501  AssertIndexRange(r, m());
2502 
2503  return iterator(this, cols->rowstart[r + 1]);
2504 }
2505 
2506 
2507 
2508 template <typename number>
2509 template <typename StreamType>
2510 inline void
2511 SparseMatrix<number>::print(StreamType &out,
2512  const bool across,
2513  const bool diagonal_first) const
2514 {
2515  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2516  Assert(val != nullptr, ExcNotInitialized());
2517 
2518  bool hanging_diagonal = false;
2519  number diagonal = number();
2520 
2521  for (size_type i = 0; i < cols->rows; ++i)
2522  {
2523  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2524  {
2525  if (!diagonal_first && i == cols->colnums[j])
2526  {
2527  diagonal = val[j];
2528  hanging_diagonal = true;
2529  }
2530  else
2531  {
2532  if (hanging_diagonal && cols->colnums[j] > i)
2533  {
2534  if (across)
2535  out << ' ' << i << ',' << i << ':' << diagonal;
2536  else
2537  out << '(' << i << ',' << i << ") " << diagonal
2538  << std::endl;
2539  hanging_diagonal = false;
2540  }
2541  if (across)
2542  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2543  else
2544  out << '(' << i << ',' << cols->colnums[j] << ") " << val[j]
2545  << std::endl;
2546  }
2547  }
2548  if (hanging_diagonal)
2549  {
2550  if (across)
2551  out << ' ' << i << ',' << i << ':' << diagonal;
2552  else
2553  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2554  hanging_diagonal = false;
2555  }
2556  }
2557  if (across)
2558  out << std::endl;
2559 }
2560 
2561 
2562 template <typename number>
2563 inline void
2565 {
2566  // nothing to do here
2567 }
2568 
2569 
2570 
2571 template <typename number>
2572 inline void
2574 {
2575  // nothing to do here
2576 }
2577 
2578 #endif // DOXYGEN
2579 
2580 
2582 
2583 #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)
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
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
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:1616
static ::ExceptionBase & ExcDifferentSparsityPatterns()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
#define AssertIsFinite(number)
Definition: exceptions.h:1884
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:535
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:490
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)