Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18:55: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 
20 # include <deal.II/base/config.h>
21 
24 
25 # include <deal.II/lac/exceptions.h>
29 # ifdef DEAL_II_WITH_MPI
30 # include <mpi.h>
31 # endif
32 
33 # include <iterator>
34 # include <memory>
35 
36 
38 
39 // Forward declarations
40 # ifndef DOXYGEN
41 template <typename number>
42 class Vector;
43 template <typename number>
44 class FullMatrix;
45 template <typename Matrix>
46 class BlockMatrixBase;
47 template <typename number>
48 class SparseILU;
49 # ifdef DEAL_II_WITH_MPI
50 namespace Utilities
51 {
52  namespace MPI
53  {
54  template <typename Number>
55  void
56  sum(const SparseMatrix<Number> &, const MPI_Comm &, SparseMatrix<Number> &);
57  }
58 } // namespace Utilities
59 # endif
60 
61 # ifdef DEAL_II_WITH_TRILINOS
62 namespace TrilinosWrappers
63 {
64  class SparseMatrix;
65 }
66 # endif
67 # endif
68 
79 {
84 
85  // forward declaration
86  template <typename number, bool Constness>
87  class Iterator;
88 
99  template <typename number, bool Constness>
101  {
102  public:
106  number
107  value() const;
108 
112  number &
113  value();
114 
119  const SparseMatrix<number> &
120  get_matrix() const;
121  };
122 
123 
124 
131  template <typename number>
132  class Accessor<number, true> : public SparsityPatternIterators::Accessor
133  {
134  public:
140 
144  Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
145 
150 
155 
159  number
160  value() const;
161 
166  const MatrixType &
167  get_matrix() const;
168 
169  private:
174 
179 
180  // Make iterator class a friend.
181  template <typename, bool>
182  friend class Iterator;
183  };
184 
185 
192  template <typename number>
193  class Accessor<number, false> : public SparsityPatternIterators::Accessor
194  {
195  private:
220  class Reference
221  {
222  public:
227  Reference(const Accessor *accessor, const bool dummy);
228 
232  operator number() const;
233 
237  const Reference &
238  operator=(const number n) const;
239 
243  const Reference &
244  operator+=(const number n) const;
245 
249  const Reference &
250  operator-=(const number n) const;
251 
255  const Reference &
256  operator*=(const number n) const;
257 
261  const Reference &
262  operator/=(const number n) const;
263 
264  private:
270  };
271 
272  public:
278 
282  Accessor(MatrixType *matrix, const std::size_t index);
283 
288 
292  Reference
293  value() const;
294 
299  MatrixType &
300  get_matrix() const;
301 
302  private:
307 
312 
313  // Make iterator class a friend.
314  template <typename, bool>
315  friend class Iterator;
316  };
317 
318 
319 
349  template <typename number, bool Constness>
350  class Iterator
351  {
352  public:
357 
363 
369 
374  Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
375 
380 
386 
392 
396  Iterator &
398 
402  Iterator
404 
409  operator*() const;
410 
415  operator->() const;
416 
420  bool
421  operator==(const Iterator &) const;
422 
426  bool
427  operator!=(const Iterator &) const;
428 
436  bool
437  operator<(const Iterator &) const;
438 
443  bool
444  operator>(const Iterator &) const;
445 
452  int
453  operator-(const Iterator &p) const;
454 
458  Iterator
459  operator+(const size_type n) const;
460 
461  private:
466  };
467 
468 } // namespace SparseMatrixIterators
469 
471 
472 namespace std
473 {
474  template <typename number, bool Constness>
475  struct iterator_traits<
476  ::SparseMatrixIterators::Iterator<number, Constness>>
477  {
478  using iterator_category = forward_iterator_tag;
479  using value_type =
480  typename ::SparseMatrixIterators::Iterator<number,
481  Constness>::value_type;
482  using difference_type = typename ::SparseMatrixIterators::
483  Iterator<number, Constness>::difference_type;
484  };
485 } // namespace std
486 
488 
494 // TODO: Add multithreading to the other vmult functions.
495 
522 template <typename number>
523 class SparseMatrix : public virtual Subscriptor
524 {
525 public:
530 
535  using value_type = number;
536 
547 
553 
561 
568  struct Traits
569  {
574  static const bool zero_addition_can_be_elided = true;
575  };
576 
592 
602 
611 
625  explicit SparseMatrix(const SparsityPattern &sparsity);
626 
633  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
634 
639  virtual ~SparseMatrix() override;
640 
652 
659 
668 
680  SparseMatrix &
681  operator=(const double d);
682 
696  virtual void
697  reinit(const SparsityPattern &sparsity);
698 
705  template <typename number2>
706  void
707  reinit(const SparseMatrix<number2> &sparse_matrix);
708 
714  virtual void
715  clear();
717 
725  bool
726  empty() const;
727 
732  size_type
733  m() const;
734 
739  size_type
740  n() const;
741 
745  size_type
746  get_row_length(const size_type row) const;
747 
753  std::size_t
755 
765  std::size_t
766  n_actually_nonzero_elements(const double threshold = 0.) const;
767 
776  const SparsityPattern &
778 
783  std::size_t
785 
790 
792 
801  void
802  set(const size_type i, const size_type j, const number value);
803 
819  template <typename number2>
820  void
821  set(const std::vector<size_type> &indices,
822  const FullMatrix<number2> & full_matrix,
823  const bool elide_zero_values = false);
824 
830  template <typename number2>
831  void
832  set(const std::vector<size_type> &row_indices,
833  const std::vector<size_type> &col_indices,
834  const FullMatrix<number2> & full_matrix,
835  const bool elide_zero_values = false);
836 
847  template <typename number2>
848  void
849  set(const size_type row,
850  const std::vector<size_type> &col_indices,
851  const std::vector<number2> & values,
852  const bool elide_zero_values = false);
853 
863  template <typename number2>
864  void
865  set(const size_type row,
866  const size_type n_cols,
867  const size_type *col_indices,
868  const number2 * values,
869  const bool elide_zero_values = false);
870 
876  void
877  add(const size_type i, const size_type j, const number value);
878 
893  template <typename number2>
894  void
895  add(const std::vector<size_type> &indices,
896  const FullMatrix<number2> & full_matrix,
897  const bool elide_zero_values = true);
898 
904  template <typename number2>
905  void
906  add(const std::vector<size_type> &row_indices,
907  const std::vector<size_type> &col_indices,
908  const FullMatrix<number2> & full_matrix,
909  const bool elide_zero_values = true);
910 
920  template <typename number2>
921  void
922  add(const size_type row,
923  const std::vector<size_type> &col_indices,
924  const std::vector<number2> & values,
925  const bool elide_zero_values = true);
926 
936  template <typename number2>
937  void
938  add(const size_type row,
939  const size_type n_cols,
940  const size_type *col_indices,
941  const number2 * values,
942  const bool elide_zero_values = true,
943  const bool col_indices_are_sorted = false);
944 
948  SparseMatrix &
949  operator*=(const number factor);
950 
954  SparseMatrix &
955  operator/=(const number factor);
956 
969  void
971 
988  template <typename somenumber>
991 
1008  template <typename ForwardIterator>
1009  void
1010  copy_from(const ForwardIterator begin, const ForwardIterator end);
1011 
1021  template <typename somenumber>
1022  void
1024 
1025 # ifdef DEAL_II_WITH_TRILINOS
1037 # endif
1038 
1050  template <typename somenumber>
1051  void
1052  add(const number factor, const SparseMatrix<somenumber> &matrix);
1053 
1055 
1059 
1073  const number &
1074  operator()(const size_type i, const size_type j) const;
1075 
1079  number &
1080  operator()(const size_type i, const size_type j);
1081 
1094  number
1095  el(const size_type i, const size_type j) const;
1096 
1106  number
1107  diag_element(const size_type i) const;
1108 
1113  number &
1115 
1117 
1137  template <class OutVector, class InVector>
1138  void
1139  vmult(OutVector &dst, const InVector &src) const;
1140 
1156  template <class OutVector, class InVector>
1157  void
1158  Tvmult(OutVector &dst, const InVector &src) const;
1159 
1176  template <class OutVector, class InVector>
1177  void
1178  vmult_add(OutVector &dst, const InVector &src) const;
1179 
1195  template <class OutVector, class InVector>
1196  void
1197  Tvmult_add(OutVector &dst, const InVector &src) const;
1198 
1216  template <typename somenumber>
1217  somenumber
1219 
1225  template <typename somenumber>
1226  somenumber
1228  const Vector<somenumber> &v) const;
1229 
1239  template <typename somenumber>
1240  somenumber
1242  const Vector<somenumber> &x,
1243  const Vector<somenumber> &b) const;
1244 
1280  template <typename numberB, typename numberC>
1281  void
1283  const SparseMatrix<numberB> &B,
1284  const Vector<number> & V = Vector<number>(),
1285  const bool rebuild_sparsity_pattern = true) const;
1286 
1311  template <typename numberB, typename numberC>
1312  void
1314  const SparseMatrix<numberB> &B,
1315  const Vector<number> & V = Vector<number>(),
1316  const bool rebuild_sparsity_pattern = true) const;
1317 
1319 
1323 
1331  real_type
1332  l1_norm() const;
1333 
1341  real_type
1342  linfty_norm() const;
1343 
1348  real_type
1351 
1355 
1361  template <typename somenumber>
1362  void
1364  const Vector<somenumber> &src,
1365  const number omega = 1.) const;
1366 
1373  template <typename somenumber>
1374  void
1376  const Vector<somenumber> & src,
1377  const number omega = 1.,
1378  const std::vector<std::size_t> &pos_right_of_diagonal =
1379  std::vector<std::size_t>()) const;
1380 
1384  template <typename somenumber>
1385  void
1387  const Vector<somenumber> &src,
1388  const number omega = 1.) const;
1389 
1393  template <typename somenumber>
1394  void
1396  const Vector<somenumber> &src,
1397  const number omega = 1.) const;
1398 
1404  template <typename somenumber>
1405  void
1406  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1407 
1412  template <typename somenumber>
1413  void
1414  SOR(Vector<somenumber> &v, const number omega = 1.) const;
1415 
1420  template <typename somenumber>
1421  void
1422  TSOR(Vector<somenumber> &v, const number omega = 1.) const;
1423 
1434  template <typename somenumber>
1435  void
1437  const std::vector<size_type> &permutation,
1438  const std::vector<size_type> &inverse_permutation,
1439  const number omega = 1.) const;
1440 
1451  template <typename somenumber>
1452  void
1454  const std::vector<size_type> &permutation,
1455  const std::vector<size_type> &inverse_permutation,
1456  const number omega = 1.) const;
1457 
1463  template <typename somenumber>
1464  void
1466  const Vector<somenumber> &b,
1467  const number omega = 1.) const;
1468 
1473  template <typename somenumber>
1474  void
1476  const Vector<somenumber> &b,
1477  const number omega = 1.) const;
1478 
1483  template <typename somenumber>
1484  void
1486  const Vector<somenumber> &b,
1487  const number omega = 1.) const;
1488 
1493  template <typename somenumber>
1494  void
1496  const Vector<somenumber> &b,
1497  const number omega = 1.) const;
1499 
1503 
1511  begin() const;
1512 
1516  iterator
1518 
1523  end() const;
1524 
1528  iterator
1529  end();
1530 
1541  begin(const size_type r) const;
1542 
1546  iterator
1547  begin(const size_type r);
1548 
1559  end(const size_type r) const;
1560 
1564  iterator
1565  end(const size_type r);
1567 
1571 
1583  template <class StreamType>
1584  void
1585  print(StreamType &out,
1586  const bool across = false,
1587  const bool diagonal_first = true) const;
1588 
1609  void
1610  print_formatted(std::ostream & out,
1611  const unsigned int precision = 3,
1612  const bool scientific = true,
1613  const unsigned int width = 0,
1614  const char * zero_string = " ",
1615  const double denominator = 1.) 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);
1669 
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 matrices that "
1717  "are the same object, but the operation requires that the "
1718  "two objects are in fact different.");
1720 
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 }
2005 
2006 
2007 
2008 template <typename number>
2009 inline SparseMatrix<number> &
2010 SparseMatrix<number>::operator*=(const number factor)
2011 {
2012  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2013  Assert(val != nullptr, ExcNotInitialized());
2014 
2015  number * val_ptr = val.get();
2016  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2017 
2018  while (val_ptr != end_ptr)
2019  *val_ptr++ *= factor;
2020 
2021  return *this;
2022 }
2023 
2024 
2025 
2026 template <typename number>
2027 inline SparseMatrix<number> &
2028 SparseMatrix<number>::operator/=(const number factor)
2029 {
2030  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2031  Assert(val != nullptr, ExcNotInitialized());
2032  Assert(factor != number(), ExcDivideByZero());
2033 
2034  const number factor_inv = number(1.) / factor;
2035 
2036  number * val_ptr = val.get();
2037  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2038 
2039  while (val_ptr != end_ptr)
2040  *val_ptr++ *= factor_inv;
2041 
2042  return *this;
2043 }
2044 
2045 
2046 
2047 template <typename number>
2048 inline const number &
2050 {
2051  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2052  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2053  ExcInvalidIndex(i, j));
2054  return val[cols->operator()(i, j)];
2055 }
2056 
2057 
2058 
2059 template <typename number>
2060 inline number &
2062 {
2063  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2064  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2065  ExcInvalidIndex(i, j));
2066  return val[cols->operator()(i, j)];
2067 }
2068 
2069 
2070 
2071 template <typename number>
2072 inline number
2073 SparseMatrix<number>::el(const size_type i, const size_type j) const
2074 {
2075  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2076  const size_type index = cols->operator()(i, j);
2077 
2079  return val[index];
2080  else
2081  return 0;
2082 }
2083 
2084 
2085 
2086 template <typename number>
2087 inline number
2089 {
2090  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2091  Assert(m() == n(), ExcNotQuadratic());
2092  AssertIndexRange(i, m());
2093 
2094  // Use that the first element in each row of a quadratic matrix is the main
2095  // diagonal
2096  return val[cols->rowstart[i]];
2097 }
2098 
2099 
2100 
2101 template <typename number>
2102 inline number &
2104 {
2105  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2106  Assert(m() == n(), ExcNotQuadratic());
2107  AssertIndexRange(i, m());
2108 
2109  // Use that the first element in each row of a quadratic matrix is the main
2110  // diagonal
2111  return val[cols->rowstart[i]];
2112 }
2113 
2114 
2115 
2116 template <typename number>
2117 template <typename ForwardIterator>
2118 void
2119 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2120  const ForwardIterator end)
2121 {
2122  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2123  ExcIteratorRange(std::distance(begin, end), m()));
2124 
2125  // for use in the inner loop, we define an alias to the type of the inner
2126  // iterators
2127  using inner_iterator =
2128  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2129  size_type row = 0;
2130  for (ForwardIterator i = begin; i != end; ++i, ++row)
2131  {
2132  const inner_iterator end_of_row = i->end();
2133  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2134  // write entries
2135  set(row, j->first, j->second);
2136  };
2137 }
2138 
2139 
2140 //---------------------------------------------------------------------------
2141 
2142 
2143 namespace SparseMatrixIterators
2144 {
2145  template <typename number>
2146  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2147  const std::size_t index_within_matrix)
2148  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2149  index_within_matrix)
2150  , matrix(matrix)
2151  {}
2152 
2153 
2154 
2155  template <typename number>
2156  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2157  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2158  , matrix(matrix)
2159  {}
2160 
2161 
2162 
2163  template <typename number>
2164  inline Accessor<number, true>::Accessor(
2166  : SparsityPatternIterators::Accessor(a)
2167  , matrix(&a.get_matrix())
2168  {}
2169 
2170 
2171 
2172  template <typename number>
2173  inline number
2174  Accessor<number, true>::value() const
2175  {
2176  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2177  return matrix->val[linear_index];
2178  }
2179 
2180 
2181 
2182  template <typename number>
2183  inline const typename Accessor<number, true>::MatrixType &
2184  Accessor<number, true>::get_matrix() const
2185  {
2186  return *matrix;
2187  }
2188 
2189 
2190 
2191  template <typename number>
2192  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2193  const bool)
2194  : accessor(accessor)
2195  {}
2196 
2197 
2198  template <typename number>
2199  inline Accessor<number, false>::Reference::operator number() const
2200  {
2201  AssertIndexRange(accessor->linear_index,
2202  accessor->matrix->n_nonzero_elements());
2203  return accessor->matrix->val[accessor->linear_index];
2204  }
2205 
2206 
2207 
2208  template <typename number>
2209  inline const typename Accessor<number, false>::Reference &
2210  Accessor<number, false>::Reference::operator=(const number n) const
2211  {
2212  AssertIndexRange(accessor->linear_index,
2213  accessor->matrix->n_nonzero_elements());
2214  accessor->matrix->val[accessor->linear_index] = n;
2215  return *this;
2216  }
2217 
2218 
2219 
2220  template <typename number>
2221  inline const typename Accessor<number, false>::Reference &
2222  Accessor<number, false>::Reference::operator+=(const number n) const
2223  {
2224  AssertIndexRange(accessor->linear_index,
2225  accessor->matrix->n_nonzero_elements());
2226  accessor->matrix->val[accessor->linear_index] += n;
2227  return *this;
2228  }
2229 
2230 
2231 
2232  template <typename number>
2233  inline const typename Accessor<number, false>::Reference &
2234  Accessor<number, false>::Reference::operator-=(const number n) const
2235  {
2236  AssertIndexRange(accessor->linear_index,
2237  accessor->matrix->n_nonzero_elements());
2238  accessor->matrix->val[accessor->linear_index] -= n;
2239  return *this;
2240  }
2241 
2242 
2243 
2244  template <typename number>
2245  inline const typename Accessor<number, false>::Reference &
2246  Accessor<number, false>::Reference::operator*=(const number n) const
2247  {
2248  AssertIndexRange(accessor->linear_index,
2249  accessor->matrix->n_nonzero_elements());
2250  accessor->matrix->val[accessor->linear_index] *= n;
2251  return *this;
2252  }
2253 
2254 
2255 
2256  template <typename number>
2257  inline const typename Accessor<number, false>::Reference &
2258  Accessor<number, false>::Reference::operator/=(const number n) const
2259  {
2260  AssertIndexRange(accessor->linear_index,
2261  accessor->matrix->n_nonzero_elements());
2262  accessor->matrix->val[accessor->linear_index] /= n;
2263  return *this;
2264  }
2265 
2266 
2267 
2268  template <typename number>
2269  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2270  const std::size_t index)
2271  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2272  , matrix(matrix)
2273  {}
2274 
2275 
2276 
2277  template <typename number>
2278  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2279  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2280  , matrix(matrix)
2281  {}
2282 
2283 
2284 
2285  template <typename number>
2286  inline typename Accessor<number, false>::Reference
2287  Accessor<number, false>::value() const
2288  {
2289  return Reference(this, true);
2290  }
2291 
2292 
2293 
2294  template <typename number>
2295  inline typename Accessor<number, false>::MatrixType &
2296  Accessor<number, false>::get_matrix() const
2297  {
2298  return *matrix;
2299  }
2300 
2301 
2302 
2303  template <typename number, bool Constness>
2304  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2305  const std::size_t index)
2306  : accessor(matrix, index)
2307  {}
2308 
2309 
2310 
2311  template <typename number, bool Constness>
2312  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2313  : accessor(matrix)
2314  {}
2315 
2316 
2317 
2318  template <typename number, bool Constness>
2319  inline Iterator<number, Constness>::Iterator(
2321  : accessor(*i)
2322  {}
2323 
2324 
2325 
2326  template <typename number, bool Constness>
2327  inline const Iterator<number, Constness> &
2328  Iterator<number, Constness>::operator=(
2330  {
2331  accessor = *i;
2332  return *this;
2333  }
2334 
2335 
2336 
2337  template <typename number, bool Constness>
2338  inline Iterator<number, Constness> &
2340  {
2341  accessor.advance();
2342  return *this;
2343  }
2344 
2345 
2346  template <typename number, bool Constness>
2347  inline Iterator<number, Constness>
2349  {
2350  const Iterator iter = *this;
2351  accessor.advance();
2352  return iter;
2353  }
2354 
2355 
2356  template <typename number, bool Constness>
2357  inline const Accessor<number, Constness> &
2359  {
2360  return accessor;
2361  }
2362 
2363 
2364  template <typename number, bool Constness>
2365  inline const Accessor<number, Constness> *
2366  Iterator<number, Constness>::operator->() const
2367  {
2368  return &accessor;
2369  }
2370 
2371 
2372  template <typename number, bool Constness>
2373  inline bool
2374  Iterator<number, Constness>::operator==(const Iterator &other) const
2375  {
2376  return (accessor == other.accessor);
2377  }
2378 
2379 
2380  template <typename number, bool Constness>
2381  inline bool
2382  Iterator<number, Constness>::operator!=(const Iterator &other) const
2383  {
2384  return !(*this == other);
2385  }
2386 
2387 
2388  template <typename number, bool Constness>
2389  inline bool
2390  Iterator<number, Constness>::operator<(const Iterator &other) const
2391  {
2392  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2393  ExcInternalError());
2394 
2395  return (accessor < other.accessor);
2396  }
2397 
2398 
2399  template <typename number, bool Constness>
2400  inline bool
2401  Iterator<number, Constness>::operator>(const Iterator &other) const
2402  {
2403  return (other < *this);
2404  }
2405 
2406 
2407  template <typename number, bool Constness>
2408  inline int
2409  Iterator<number, Constness>::operator-(const Iterator &other) const
2410  {
2411  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2412  ExcInternalError());
2413 
2414  return (*this)->linear_index - other->linear_index;
2415  }
2416 
2417 
2418 
2419  template <typename number, bool Constness>
2420  inline Iterator<number, Constness>
2422  {
2423  Iterator x = *this;
2424  for (size_type i = 0; i < n; ++i)
2425  ++x;
2426 
2427  return x;
2428  }
2429 
2430 } // namespace SparseMatrixIterators
2431 
2432 
2433 
2434 template <typename number>
2437 {
2438  return const_iterator(this, 0);
2439 }
2440 
2441 
2442 template <typename number>
2445 {
2446  return const_iterator(this);
2447 }
2448 
2449 
2450 template <typename number>
2451 inline typename SparseMatrix<number>::iterator
2453 {
2454  return iterator(this, 0);
2455 }
2456 
2457 
2458 template <typename number>
2459 inline typename SparseMatrix<number>::iterator
2461 {
2462  return iterator(this, cols->rowstart[cols->rows]);
2463 }
2464 
2465 
2466 template <typename number>
2469 {
2470  AssertIndexRange(r, m());
2471 
2472  return const_iterator(this, cols->rowstart[r]);
2473 }
2474 
2475 
2476 
2477 template <typename number>
2479 SparseMatrix<number>::end(const size_type r) const
2480 {
2481  AssertIndexRange(r, m());
2482 
2483  return const_iterator(this, cols->rowstart[r + 1]);
2484 }
2485 
2486 
2487 
2488 template <typename number>
2489 inline typename SparseMatrix<number>::iterator
2491 {
2492  AssertIndexRange(r, m());
2493 
2494  return iterator(this, cols->rowstart[r]);
2495 }
2496 
2497 
2498 
2499 template <typename number>
2500 inline typename SparseMatrix<number>::iterator
2502 {
2503  AssertIndexRange(r, m());
2504 
2505  return iterator(this, cols->rowstart[r + 1]);
2506 }
2507 
2508 
2509 
2510 template <typename number>
2511 template <class StreamType>
2512 inline void
2513 SparseMatrix<number>::print(StreamType &out,
2514  const bool across,
2515  const bool diagonal_first) const
2516 {
2517  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2518  Assert(val != nullptr, ExcNotInitialized());
2519 
2520  bool hanging_diagonal = false;
2521  number diagonal = number();
2522 
2523  for (size_type i = 0; i < cols->rows; ++i)
2524  {
2525  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2526  {
2527  if (!diagonal_first && i == cols->colnums[j])
2528  {
2529  diagonal = val[j];
2530  hanging_diagonal = true;
2531  }
2532  else
2533  {
2534  if (hanging_diagonal && cols->colnums[j] > i)
2535  {
2536  if (across)
2537  out << ' ' << i << ',' << i << ':' << diagonal;
2538  else
2539  out << '(' << i << ',' << i << ") " << diagonal
2540  << std::endl;
2541  hanging_diagonal = false;
2542  }
2543  if (across)
2544  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2545  else
2546  out << '(' << i << ',' << cols->colnums[j] << ") " << val[j]
2547  << std::endl;
2548  }
2549  }
2550  if (hanging_diagonal)
2551  {
2552  if (across)
2553  out << ' ' << i << ',' << i << ':' << diagonal;
2554  else
2555  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2556  hanging_diagonal = false;
2557  }
2558  }
2559  if (across)
2560  out << std::endl;
2561 }
2562 
2563 
2564 template <typename number>
2565 inline void
2567 {
2568  // nothing to do here
2569 }
2570 
2571 
2572 
2573 template <typename number>
2574 inline void
2576 {
2577  // nothing to do here
2578 }
2579 
2580 # endif // DOXYGEN
2581 
2582 
2583 /*---------------------------- sparse_matrix.h ---------------------------*/
2584 
2586 
2587 #endif
2588 /*---------------------------- sparse_matrix.h ---------------------------*/
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 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 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
SparseMatrix< number > & copy_from(const TrilinosWrappers::SparseMatrix &matrix)
void symmetrize()
SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id)
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 copy_from(const ForwardIterator begin, const ForwardIterator end)
SparseMatrix & operator/=(const number factor)
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
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
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)
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 const size_type invalid_entry
Definition: vector.h:109
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
__global__ void set(Number *val, const Number s, const size_type N)
void prepare_add()
static ::ExceptionBase & ExcSourceEqualsDestination()
static ::ExceptionBase & ExcDivideByZero()
void print_as_numpy_arrays(std::ostream &out, const unsigned int precision=9) const
static ::ExceptionBase & ExcInternalError()
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
static const bool zero_addition_can_be_elided
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNeedsSparsityPattern()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
#define Assert(cond, exc)
Definition: exceptions.h:1473
void print_pattern(std::ostream &out, const double threshold=0.) const
static ::ExceptionBase & ExcDifferentSparsityPatterns()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
#define AssertIsFinite(number)
Definition: exceptions.h:1758
void block_write(std::ostream &out) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
void block_read(std::istream &in)
void prepare_set()
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
Expression operator>(const Expression &lhs, const Expression &rhs)
std::unique_ptr< number[]> val
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
std::size_t max_len
static ::ExceptionBase & ExcNotQuadratic()
@ 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:83
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:76
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
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)