Reference documentation for deal.II version GIT eef19498cf 2022-10-05 14:45: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 - 2022 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>
472  struct iterator_traits<
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 <class 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 matrices 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 }
2002 
2003 
2004 
2005 template <typename number>
2006 inline SparseMatrix<number> &
2007 SparseMatrix<number>::operator*=(const number factor)
2008 {
2009  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2010  Assert(val != nullptr, ExcNotInitialized());
2011 
2012  number * val_ptr = val.get();
2013  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2014 
2015  while (val_ptr != end_ptr)
2016  *val_ptr++ *= factor;
2017 
2018  return *this;
2019 }
2020 
2021 
2022 
2023 template <typename number>
2024 inline SparseMatrix<number> &
2025 SparseMatrix<number>::operator/=(const number factor)
2026 {
2027  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2028  Assert(val != nullptr, ExcNotInitialized());
2029  Assert(factor != number(), ExcDivideByZero());
2030 
2031  const number factor_inv = number(1.) / factor;
2032 
2033  number * val_ptr = val.get();
2034  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2035 
2036  while (val_ptr != end_ptr)
2037  *val_ptr++ *= factor_inv;
2038 
2039  return *this;
2040 }
2041 
2042 
2043 
2044 template <typename number>
2045 inline const number &
2047 {
2048  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2049  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2050  ExcInvalidIndex(i, j));
2051  return val[cols->operator()(i, j)];
2052 }
2053 
2054 
2055 
2056 template <typename number>
2057 inline number &
2059 {
2060  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2061  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2062  ExcInvalidIndex(i, j));
2063  return val[cols->operator()(i, j)];
2064 }
2065 
2066 
2067 
2068 template <typename number>
2069 inline number
2070 SparseMatrix<number>::el(const size_type i, const size_type j) const
2071 {
2072  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2073  const size_type index = cols->operator()(i, j);
2074 
2076  return val[index];
2077  else
2078  return 0;
2079 }
2080 
2081 
2082 
2083 template <typename number>
2084 inline number
2086 {
2087  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2088  Assert(m() == n(), ExcNotQuadratic());
2089  AssertIndexRange(i, m());
2090 
2091  // Use that the first element in each row of a quadratic matrix is the main
2092  // diagonal
2093  return val[cols->rowstart[i]];
2094 }
2095 
2096 
2097 
2098 template <typename number>
2099 inline number &
2101 {
2102  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2103  Assert(m() == n(), ExcNotQuadratic());
2104  AssertIndexRange(i, m());
2105 
2106  // Use that the first element in each row of a quadratic matrix is the main
2107  // diagonal
2108  return val[cols->rowstart[i]];
2109 }
2110 
2111 
2112 
2113 template <typename number>
2114 template <typename ForwardIterator>
2115 void
2116 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2117  const ForwardIterator end)
2118 {
2119  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2120  ExcIteratorRange(std::distance(begin, end), m()));
2121 
2122  // for use in the inner loop, we define an alias to the type of the inner
2123  // iterators
2124  using inner_iterator =
2125  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2126  size_type row = 0;
2127  for (ForwardIterator i = begin; i != end; ++i, ++row)
2128  {
2129  const inner_iterator end_of_row = i->end();
2130  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2131  // write entries
2132  set(row, j->first, j->second);
2133  };
2134 }
2135 
2136 
2137 //---------------------------------------------------------------------------
2138 
2139 
2140 namespace SparseMatrixIterators
2141 {
2142  template <typename number>
2143  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2144  const std::size_t index_within_matrix)
2145  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2146  index_within_matrix)
2147  , matrix(matrix)
2148  {}
2149 
2150 
2151 
2152  template <typename number>
2153  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2154  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2155  , matrix(matrix)
2156  {}
2157 
2158 
2159 
2160  template <typename number>
2161  inline Accessor<number, true>::Accessor(
2163  : SparsityPatternIterators::Accessor(a)
2164  , matrix(&a.get_matrix())
2165  {}
2166 
2167 
2168 
2169  template <typename number>
2170  inline number
2171  Accessor<number, true>::value() const
2172  {
2173  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2174  return matrix->val[linear_index];
2175  }
2176 
2177 
2178 
2179  template <typename number>
2180  inline const typename Accessor<number, true>::MatrixType &
2181  Accessor<number, true>::get_matrix() const
2182  {
2183  return *matrix;
2184  }
2185 
2186 
2187 
2188  template <typename number>
2189  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2190  const bool)
2191  : accessor(accessor)
2192  {}
2193 
2194 
2195  template <typename number>
2196  inline Accessor<number, false>::Reference::operator number() const
2197  {
2198  AssertIndexRange(accessor->linear_index,
2199  accessor->matrix->n_nonzero_elements());
2200  return accessor->matrix->val[accessor->linear_index];
2201  }
2202 
2203 
2204 
2205  template <typename number>
2206  inline const typename Accessor<number, false>::Reference &
2207  Accessor<number, false>::Reference::operator=(const number n) const
2208  {
2209  AssertIndexRange(accessor->linear_index,
2210  accessor->matrix->n_nonzero_elements());
2211  accessor->matrix->val[accessor->linear_index] = n;
2212  return *this;
2213  }
2214 
2215 
2216 
2217  template <typename number>
2218  inline const typename Accessor<number, false>::Reference &
2219  Accessor<number, false>::Reference::operator+=(const number n) const
2220  {
2221  AssertIndexRange(accessor->linear_index,
2222  accessor->matrix->n_nonzero_elements());
2223  accessor->matrix->val[accessor->linear_index] += n;
2224  return *this;
2225  }
2226 
2227 
2228 
2229  template <typename number>
2230  inline const typename Accessor<number, false>::Reference &
2231  Accessor<number, false>::Reference::operator-=(const number n) const
2232  {
2233  AssertIndexRange(accessor->linear_index,
2234  accessor->matrix->n_nonzero_elements());
2235  accessor->matrix->val[accessor->linear_index] -= n;
2236  return *this;
2237  }
2238 
2239 
2240 
2241  template <typename number>
2242  inline const typename Accessor<number, false>::Reference &
2243  Accessor<number, false>::Reference::operator*=(const number n) const
2244  {
2245  AssertIndexRange(accessor->linear_index,
2246  accessor->matrix->n_nonzero_elements());
2247  accessor->matrix->val[accessor->linear_index] *= n;
2248  return *this;
2249  }
2250 
2251 
2252 
2253  template <typename number>
2254  inline const typename Accessor<number, false>::Reference &
2255  Accessor<number, false>::Reference::operator/=(const number n) const
2256  {
2257  AssertIndexRange(accessor->linear_index,
2258  accessor->matrix->n_nonzero_elements());
2259  accessor->matrix->val[accessor->linear_index] /= n;
2260  return *this;
2261  }
2262 
2263 
2264 
2265  template <typename number>
2266  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2267  const std::size_t index)
2268  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2269  , matrix(matrix)
2270  {}
2271 
2272 
2273 
2274  template <typename number>
2275  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2276  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2277  , matrix(matrix)
2278  {}
2279 
2280 
2281 
2282  template <typename number>
2283  inline typename Accessor<number, false>::Reference
2284  Accessor<number, false>::value() const
2285  {
2286  return Reference(this, true);
2287  }
2288 
2289 
2290 
2291  template <typename number>
2292  inline typename Accessor<number, false>::MatrixType &
2293  Accessor<number, false>::get_matrix() const
2294  {
2295  return *matrix;
2296  }
2297 
2298 
2299 
2300  template <typename number, bool Constness>
2301  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2302  const std::size_t index)
2303  : accessor(matrix, index)
2304  {}
2305 
2306 
2307 
2308  template <typename number, bool Constness>
2309  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2310  : accessor(matrix)
2311  {}
2312 
2313 
2314 
2315  template <typename number, bool Constness>
2316  inline Iterator<number, Constness>::Iterator(
2318  : accessor(*i)
2319  {}
2320 
2321 
2322 
2323  template <typename number, bool Constness>
2324  inline const Iterator<number, Constness> &
2325  Iterator<number, Constness>::operator=(
2327  {
2328  accessor = *i;
2329  return *this;
2330  }
2331 
2332 
2333 
2334  template <typename number, bool Constness>
2335  inline Iterator<number, Constness> &
2337  {
2338  accessor.advance();
2339  return *this;
2340  }
2341 
2342 
2343  template <typename number, bool Constness>
2344  inline Iterator<number, Constness>
2346  {
2347  const Iterator iter = *this;
2348  accessor.advance();
2349  return iter;
2350  }
2351 
2352 
2353  template <typename number, bool Constness>
2354  inline const Accessor<number, Constness> &
2356  {
2357  return accessor;
2358  }
2359 
2360 
2361  template <typename number, bool Constness>
2362  inline const Accessor<number, Constness> *
2363  Iterator<number, Constness>::operator->() const
2364  {
2365  return &accessor;
2366  }
2367 
2368 
2369  template <typename number, bool Constness>
2370  inline bool
2371  Iterator<number, Constness>::operator==(const Iterator &other) const
2372  {
2373  return (accessor == other.accessor);
2374  }
2375 
2376 
2377  template <typename number, bool Constness>
2378  inline bool
2379  Iterator<number, Constness>::operator!=(const Iterator &other) const
2380  {
2381  return !(*this == other);
2382  }
2383 
2384 
2385  template <typename number, bool Constness>
2386  inline bool
2387  Iterator<number, Constness>::operator<(const Iterator &other) const
2388  {
2389  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2390  ExcInternalError());
2391 
2392  return (accessor < other.accessor);
2393  }
2394 
2395 
2396  template <typename number, bool Constness>
2397  inline bool
2398  Iterator<number, Constness>::operator>(const Iterator &other) const
2399  {
2400  return (other < *this);
2401  }
2402 
2403 
2404  template <typename number, bool Constness>
2405  inline int
2406  Iterator<number, Constness>::operator-(const Iterator &other) const
2407  {
2408  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2409  ExcInternalError());
2410 
2411  return (*this)->linear_index - other->linear_index;
2412  }
2413 
2414 
2415 
2416  template <typename number, bool Constness>
2417  inline Iterator<number, Constness>
2419  {
2420  Iterator x = *this;
2421  for (size_type i = 0; i < n; ++i)
2422  ++x;
2423 
2424  return x;
2425  }
2426 
2427 } // namespace SparseMatrixIterators
2428 
2429 
2430 
2431 template <typename number>
2434 {
2435  return const_iterator(this, 0);
2436 }
2437 
2438 
2439 template <typename number>
2442 {
2443  return const_iterator(this);
2444 }
2445 
2446 
2447 template <typename number>
2448 inline typename SparseMatrix<number>::iterator
2450 {
2451  return iterator(this, 0);
2452 }
2453 
2454 
2455 template <typename number>
2456 inline typename SparseMatrix<number>::iterator
2458 {
2459  return iterator(this, cols->rowstart[cols->rows]);
2460 }
2461 
2462 
2463 template <typename number>
2466 {
2467  AssertIndexRange(r, m());
2468 
2469  return const_iterator(this, cols->rowstart[r]);
2470 }
2471 
2472 
2473 
2474 template <typename number>
2476 SparseMatrix<number>::end(const size_type r) const
2477 {
2478  AssertIndexRange(r, m());
2479 
2480  return const_iterator(this, cols->rowstart[r + 1]);
2481 }
2482 
2483 
2484 
2485 template <typename number>
2486 inline typename SparseMatrix<number>::iterator
2488 {
2489  AssertIndexRange(r, m());
2490 
2491  return iterator(this, cols->rowstart[r]);
2492 }
2493 
2494 
2495 
2496 template <typename number>
2497 inline typename SparseMatrix<number>::iterator
2499 {
2500  AssertIndexRange(r, m());
2501 
2502  return iterator(this, cols->rowstart[r + 1]);
2503 }
2504 
2505 
2506 
2507 template <typename number>
2508 template <class StreamType>
2509 inline void
2510 SparseMatrix<number>::print(StreamType &out,
2511  const bool across,
2512  const bool diagonal_first) const
2513 {
2514  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2515  Assert(val != nullptr, ExcNotInitialized());
2516 
2517  bool hanging_diagonal = false;
2518  number diagonal = number();
2519 
2520  for (size_type i = 0; i < cols->rows; ++i)
2521  {
2522  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2523  {
2524  if (!diagonal_first && i == cols->colnums[j])
2525  {
2526  diagonal = val[j];
2527  hanging_diagonal = true;
2528  }
2529  else
2530  {
2531  if (hanging_diagonal && cols->colnums[j] > i)
2532  {
2533  if (across)
2534  out << ' ' << i << ',' << i << ':' << diagonal;
2535  else
2536  out << '(' << i << ',' << i << ") " << diagonal
2537  << std::endl;
2538  hanging_diagonal = false;
2539  }
2540  if (across)
2541  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2542  else
2543  out << '(' << i << ',' << cols->colnums[j] << ") " << val[j]
2544  << std::endl;
2545  }
2546  }
2547  if (hanging_diagonal)
2548  {
2549  if (across)
2550  out << ' ' << i << ',' << i << ':' << diagonal;
2551  else
2552  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2553  hanging_diagonal = false;
2554  }
2555  }
2556  if (across)
2557  out << std::endl;
2558 }
2559 
2560 
2561 template <typename number>
2562 inline void
2564 {
2565  // nothing to do here
2566 }
2567 
2568 
2569 
2570 template <typename number>
2571 inline void
2573 {
2574  // nothing to do here
2575 }
2576 
2577 #endif // DOXYGEN
2578 
2579 
2581 
2582 #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
void compress(::VectorOperation::values)
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
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:109
std::enable_if_t< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, 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:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
__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:1501
static ::ExceptionBase & ExcDifferentSparsityPatterns()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
#define AssertIsFinite(number)
Definition: exceptions.h:1786
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
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:81
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)