Reference documentation for deal.II version Git 52d5a2a9d6 2021-04-13 22:10:47 -0400
\(\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 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_sparse_matrix_h
17 # define dealii_sparse_matrix_h
18 
19 
20 # include <deal.II/base/config.h>
21 
24 
25 # include <deal.II/lac/exceptions.h>
29 # ifdef DEAL_II_WITH_MPI
30 # include <mpi.h>
31 # endif
32 
33 # include <memory>
34 
35 
37 
38 // Forward declarations
39 # ifndef DOXYGEN
40 template <typename number>
41 class Vector;
42 template <typename number>
43 class FullMatrix;
44 template <typename Matrix>
45 class BlockMatrixBase;
46 template <typename number>
47 class SparseILU;
48 # ifdef DEAL_II_WITH_MPI
49 namespace Utilities
50 {
51  namespace MPI
52  {
53  template <typename Number>
54  void
56  }
57 } // namespace Utilities
58 # endif
59 
60 # ifdef DEAL_II_WITH_TRILINOS
61 namespace TrilinosWrappers
62 {
63  class SparseMatrix;
64 }
65 # endif
66 # endif
67 
78 {
83 
84  // forward declaration
85  template <typename number, bool Constness>
86  class Iterator;
87 
98  template <typename number, bool Constness>
100  {
101  public:
105  number
106  value() const;
107 
111  number &
112  value();
113 
118  const SparseMatrix<number> &
119  get_matrix() const;
120  };
121 
122 
123 
130  template <typename number>
131  class Accessor<number, true> : public SparsityPatternIterators::Accessor
132  {
133  public:
139 
143  Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
144 
149 
154 
158  number
159  value() const;
160 
165  const MatrixType &
166  get_matrix() const;
167 
168  private:
173 
178 
179  // Make iterator class a friend.
180  template <typename, bool>
181  friend class Iterator;
182  };
183 
184 
191  template <typename number>
192  class Accessor<number, false> : public SparsityPatternIterators::Accessor
193  {
194  private:
219  class Reference
220  {
221  public:
226  Reference(const Accessor *accessor, const bool dummy);
227 
231  operator number() const;
232 
236  const Reference &
237  operator=(const number n) const;
238 
242  const Reference &
243  operator+=(const number n) const;
244 
248  const Reference &
249  operator-=(const number n) const;
250 
254  const Reference &
255  operator*=(const number n) const;
256 
260  const Reference &
261  operator/=(const number n) const;
262 
263  private:
269  };
270 
271  public:
277 
281  Accessor(MatrixType *matrix, const std::size_t index);
282 
287 
291  Reference
292  value() const;
293 
298  MatrixType &
299  get_matrix() const;
300 
301  private:
306 
311 
312  // Make iterator class a friend.
313  template <typename, bool>
314  friend class Iterator;
315  };
316 
317 
318 
348  template <typename number, bool Constness>
349  class Iterator
350  {
351  public:
356 
362 
367  Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
368 
373 
379 
385 
389  Iterator &
390  operator++();
391 
395  Iterator
396  operator++(int);
397 
401  const Accessor<number, Constness> &operator*() const;
402 
406  const Accessor<number, Constness> *operator->() const;
407 
411  bool
412  operator==(const Iterator &) const;
413 
417  bool
418  operator!=(const Iterator &) const;
419 
427  bool
428  operator<(const Iterator &) const;
429 
434  bool
435  operator>(const Iterator &) const;
436 
443  int
444  operator-(const Iterator &p) const;
445 
449  Iterator
450  operator+(const size_type n) const;
451 
452  private:
457  };
458 
459 } // namespace SparseMatrixIterators
460 
466 // TODO: Add multithreading to the other vmult functions.
467 
494 template <typename number>
495 class SparseMatrix : public virtual Subscriptor
496 {
497 public:
502 
507  using value_type = number;
508 
519 
525 
533 
540  struct Traits
541  {
546  static const bool zero_addition_can_be_elided = true;
547  };
548 
563  SparseMatrix();
564 
573  SparseMatrix(const SparseMatrix &);
574 
582  SparseMatrix(SparseMatrix<number> &&m) noexcept;
583 
597  explicit SparseMatrix(const SparsityPattern &sparsity);
598 
605  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
606 
611  virtual ~SparseMatrix() override;
612 
623  operator=(const SparseMatrix<number> &);
624 
630  operator=(SparseMatrix<number> &&m) noexcept;
631 
639  operator=(const IdentityMatrix &id);
640 
652  SparseMatrix &
653  operator=(const double d);
654 
668  virtual void
669  reinit(const SparsityPattern &sparsity);
670 
676  virtual void
677  clear();
679 
687  bool
688  empty() const;
689 
694  size_type
695  m() const;
696 
701  size_type
702  n() const;
703 
707  size_type
708  get_row_length(const size_type row) const;
709 
715  std::size_t
716  n_nonzero_elements() const;
717 
727  std::size_t
728  n_actually_nonzero_elements(const double threshold = 0.) const;
729 
738  const SparsityPattern &
739  get_sparsity_pattern() const;
740 
745  std::size_t
746  memory_consumption() const;
747 
752 
754 
763  void
764  set(const size_type i, const size_type j, const number value);
765 
781  template <typename number2>
782  void
783  set(const std::vector<size_type> &indices,
784  const FullMatrix<number2> & full_matrix,
785  const bool elide_zero_values = false);
786 
792  template <typename number2>
793  void
794  set(const std::vector<size_type> &row_indices,
795  const std::vector<size_type> &col_indices,
796  const FullMatrix<number2> & full_matrix,
797  const bool elide_zero_values = false);
798 
809  template <typename number2>
810  void
811  set(const size_type row,
812  const std::vector<size_type> &col_indices,
813  const std::vector<number2> & values,
814  const bool elide_zero_values = false);
815 
825  template <typename number2>
826  void
827  set(const size_type row,
828  const size_type n_cols,
829  const size_type *col_indices,
830  const number2 * values,
831  const bool elide_zero_values = false);
832 
838  void
839  add(const size_type i, const size_type j, const number value);
840 
855  template <typename number2>
856  void
857  add(const std::vector<size_type> &indices,
858  const FullMatrix<number2> & full_matrix,
859  const bool elide_zero_values = true);
860 
866  template <typename number2>
867  void
868  add(const std::vector<size_type> &row_indices,
869  const std::vector<size_type> &col_indices,
870  const FullMatrix<number2> & full_matrix,
871  const bool elide_zero_values = true);
872 
882  template <typename number2>
883  void
884  add(const size_type row,
885  const std::vector<size_type> &col_indices,
886  const std::vector<number2> & values,
887  const bool elide_zero_values = true);
888 
898  template <typename number2>
899  void
900  add(const size_type row,
901  const size_type n_cols,
902  const size_type *col_indices,
903  const number2 * values,
904  const bool elide_zero_values = true,
905  const bool col_indices_are_sorted = false);
906 
910  SparseMatrix &
911  operator*=(const number factor);
912 
916  SparseMatrix &
917  operator/=(const number factor);
918 
931  void
932  symmetrize();
933 
950  template <typename somenumber>
952  copy_from(const SparseMatrix<somenumber> &source);
953 
970  template <typename ForwardIterator>
971  void
972  copy_from(const ForwardIterator begin, const ForwardIterator end);
973 
983  template <typename somenumber>
984  void
985  copy_from(const FullMatrix<somenumber> &matrix);
986 
987 # ifdef DEAL_II_WITH_TRILINOS
988 
998  copy_from(const TrilinosWrappers::SparseMatrix &matrix);
999 # endif
1000 
1012  template <typename somenumber>
1013  void
1014  add(const number factor, const SparseMatrix<somenumber> &matrix);
1015 
1017 
1021 
1035  const number &
1036  operator()(const size_type i, const size_type j) const;
1037 
1041  number &
1042  operator()(const size_type i, const size_type j);
1043 
1056  number
1057  el(const size_type i, const size_type j) const;
1058 
1068  number
1069  diag_element(const size_type i) const;
1070 
1075  number &
1076  diag_element(const size_type i);
1077 
1079 
1099  template <class OutVector, class InVector>
1100  void
1101  vmult(OutVector &dst, const InVector &src) const;
1102 
1118  template <class OutVector, class InVector>
1119  void
1120  Tvmult(OutVector &dst, const InVector &src) const;
1121 
1138  template <class OutVector, class InVector>
1139  void
1140  vmult_add(OutVector &dst, const InVector &src) const;
1141 
1157  template <class OutVector, class InVector>
1158  void
1159  Tvmult_add(OutVector &dst, const InVector &src) const;
1160 
1178  template <typename somenumber>
1179  somenumber
1180  matrix_norm_square(const Vector<somenumber> &v) const;
1181 
1187  template <typename somenumber>
1188  somenumber
1189  matrix_scalar_product(const Vector<somenumber> &u,
1190  const Vector<somenumber> &v) const;
1191 
1201  template <typename somenumber>
1202  somenumber
1203  residual(Vector<somenumber> & dst,
1204  const Vector<somenumber> &x,
1205  const Vector<somenumber> &b) const;
1206 
1242  template <typename numberB, typename numberC>
1243  void
1244  mmult(SparseMatrix<numberC> & C,
1245  const SparseMatrix<numberB> &B,
1246  const Vector<number> & V = Vector<number>(),
1247  const bool rebuild_sparsity_pattern = true) const;
1248 
1273  template <typename numberB, typename numberC>
1274  void
1275  Tmmult(SparseMatrix<numberC> & C,
1276  const SparseMatrix<numberB> &B,
1277  const Vector<number> & V = Vector<number>(),
1278  const bool rebuild_sparsity_pattern = true) const;
1279 
1281 
1285 
1293  real_type
1294  l1_norm() const;
1295 
1303  real_type
1304  linfty_norm() const;
1305 
1310  real_type
1311  frobenius_norm() const;
1313 
1317 
1323  template <typename somenumber>
1324  void
1325  precondition_Jacobi(Vector<somenumber> & dst,
1326  const Vector<somenumber> &src,
1327  const number omega = 1.) const;
1328 
1335  template <typename somenumber>
1336  void
1337  precondition_SSOR(Vector<somenumber> & dst,
1338  const Vector<somenumber> & src,
1339  const number omega = 1.,
1340  const std::vector<std::size_t> &pos_right_of_diagonal =
1341  std::vector<std::size_t>()) const;
1342 
1346  template <typename somenumber>
1347  void
1348  precondition_SOR(Vector<somenumber> & dst,
1349  const Vector<somenumber> &src,
1350  const number om = 1.) const;
1351 
1355  template <typename somenumber>
1356  void
1357  precondition_TSOR(Vector<somenumber> & dst,
1358  const Vector<somenumber> &src,
1359  const number om = 1.) const;
1360 
1366  template <typename somenumber>
1367  void
1368  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1369 
1374  template <typename somenumber>
1375  void
1376  SOR(Vector<somenumber> &v, const number om = 1.) const;
1377 
1382  template <typename somenumber>
1383  void
1384  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1385 
1396  template <typename somenumber>
1397  void
1398  PSOR(Vector<somenumber> & v,
1399  const std::vector<size_type> &permutation,
1400  const std::vector<size_type> &inverse_permutation,
1401  const number om = 1.) const;
1402 
1413  template <typename somenumber>
1414  void
1415  TPSOR(Vector<somenumber> & v,
1416  const std::vector<size_type> &permutation,
1417  const std::vector<size_type> &inverse_permutation,
1418  const number om = 1.) const;
1419 
1425  template <typename somenumber>
1426  void
1427  Jacobi_step(Vector<somenumber> & v,
1428  const Vector<somenumber> &b,
1429  const number om = 1.) const;
1430 
1435  template <typename somenumber>
1436  void
1437  SOR_step(Vector<somenumber> & v,
1438  const Vector<somenumber> &b,
1439  const number om = 1.) const;
1440 
1445  template <typename somenumber>
1446  void
1447  TSOR_step(Vector<somenumber> & v,
1448  const Vector<somenumber> &b,
1449  const number om = 1.) const;
1450 
1455  template <typename somenumber>
1456  void
1457  SSOR_step(Vector<somenumber> & v,
1458  const Vector<somenumber> &b,
1459  const number om = 1.) const;
1461 
1465 
1473  begin() const;
1474 
1478  iterator
1479  begin();
1480 
1485  end() const;
1486 
1490  iterator
1491  end();
1492 
1503  begin(const size_type r) const;
1504 
1508  iterator
1509  begin(const size_type r);
1510 
1521  end(const size_type r) const;
1522 
1526  iterator
1527  end(const size_type r);
1529 
1533 
1545  template <class StreamType>
1546  void
1547  print(StreamType &out,
1548  const bool across = false,
1549  const bool diagonal_first = true) const;
1550 
1571  void
1572  print_formatted(std::ostream & out,
1573  const unsigned int precision = 3,
1574  const bool scientific = true,
1575  const unsigned int width = 0,
1576  const char * zero_string = " ",
1577  const double denominator = 1.) const;
1578 
1584  void
1585  print_pattern(std::ostream &out, const double threshold = 0.) const;
1586 
1595  void
1596  print_as_numpy_arrays(std::ostream & out,
1597  const unsigned int precision = 9) const;
1598 
1609  void
1610  block_write(std::ostream &out) const;
1611 
1628  void
1629  block_read(std::istream &in);
1631 
1639  DeclException2(ExcInvalidIndex,
1640  int,
1641  int,
1642  << "You are trying to access the matrix entry with index <"
1643  << arg1 << ',' << arg2
1644  << ">, but this entry does not exist in the sparsity pattern "
1645  "of this matrix."
1646  "\n\n"
1647  "The most common cause for this problem is that you used "
1648  "a method to build the sparsity pattern that did not "
1649  "(completely) take into account all of the entries you "
1650  "will later try to write into. An example would be "
1651  "building a sparsity pattern that does not include "
1652  "the entries you will write into due to constraints "
1653  "on degrees of freedom such as hanging nodes or periodic "
1654  "boundary conditions. In such cases, building the "
1655  "sparsity pattern will succeed, but you will get errors "
1656  "such as the current one at one point or other when "
1657  "trying to write into the entries of the matrix.");
1661  DeclExceptionMsg(ExcDifferentSparsityPatterns,
1662  "When copying one sparse matrix into another, "
1663  "or when adding one sparse matrix to another, "
1664  "both matrices need to refer to the same "
1665  "sparsity pattern.");
1669  DeclException2(ExcIteratorRange,
1670  int,
1671  int,
1672  << "The iterators denote a range of " << arg1
1673  << " elements, but the given number of rows was " << arg2);
1677  DeclExceptionMsg(ExcSourceEqualsDestination,
1678  "You are attempting an operation on two matrices that "
1679  "are the same object, but the operation requires that the "
1680  "two objects are in fact different.");
1682 
1683 protected:
1694  void
1695  prepare_add();
1696 
1701  void
1702  prepare_set();
1703 
1704 private:
1711 
1719  std::unique_ptr<number[]> val;
1720 
1727  std::size_t max_len;
1728 
1729  // make all other sparse matrices friends
1730  template <typename somenumber>
1731  friend class SparseMatrix;
1732  template <typename somenumber>
1734  template <typename>
1735  friend class SparseILU;
1736 
1737  // To allow it calling private prepare_add() and prepare_set().
1738  template <typename>
1739  friend class BlockMatrixBase;
1740 
1741  // Also give access to internal details to the iterator/accessor classes.
1742  template <typename, bool>
1744  template <typename, bool>
1746 
1747 # ifdef DEAL_II_WITH_MPI
1748  // Give access to internal datastructures to perform MPI operations.
1749  template <typename Number>
1750  friend void
1752  const MPI_Comm &,
1754 # endif
1755 };
1756 
1757 # ifndef DOXYGEN
1758 /*---------------------- Inline functions -----------------------------------*/
1759 
1760 
1761 
1762 template <typename number>
1763 inline typename SparseMatrix<number>::size_type
1765 {
1766  Assert(cols != nullptr, ExcNotInitialized());
1767  return cols->rows;
1768 }
1769 
1770 
1771 template <typename number>
1772 inline typename SparseMatrix<number>::size_type
1774 {
1775  Assert(cols != nullptr, ExcNotInitialized());
1776  return cols->cols;
1777 }
1778 
1779 
1780 // Inline the set() and add() functions, since they will be called frequently.
1781 template <typename number>
1782 inline void
1784  const size_type j,
1785  const number value)
1786 {
1787  AssertIsFinite(value);
1788 
1789  const size_type index = cols->operator()(i, j);
1790 
1791  // it is allowed to set elements of the matrix that are not part of the
1792  // sparsity pattern, if the value to which we set it is zero
1793  if (index == SparsityPattern::invalid_entry)
1794  {
1795  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1796  ExcInvalidIndex(i, j));
1797  return;
1798  }
1799 
1800  val[index] = value;
1801 }
1802 
1803 
1804 
1805 template <typename number>
1806 template <typename number2>
1807 inline void
1808 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1809  const FullMatrix<number2> & values,
1810  const bool elide_zero_values)
1811 {
1812  Assert(indices.size() == values.m(),
1813  ExcDimensionMismatch(indices.size(), values.m()));
1814  Assert(values.m() == values.n(), ExcNotQuadratic());
1815 
1816  for (size_type i = 0; i < indices.size(); ++i)
1817  set(indices[i],
1818  indices.size(),
1819  indices.data(),
1820  &values(i, 0),
1821  elide_zero_values);
1822 }
1823 
1824 
1825 
1826 template <typename number>
1827 template <typename number2>
1828 inline void
1829 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1830  const std::vector<size_type> &col_indices,
1831  const FullMatrix<number2> & values,
1832  const bool elide_zero_values)
1833 {
1834  Assert(row_indices.size() == values.m(),
1835  ExcDimensionMismatch(row_indices.size(), values.m()));
1836  Assert(col_indices.size() == values.n(),
1837  ExcDimensionMismatch(col_indices.size(), values.n()));
1838 
1839  for (size_type i = 0; i < row_indices.size(); ++i)
1840  set(row_indices[i],
1841  col_indices.size(),
1842  col_indices.data(),
1843  &values(i, 0),
1844  elide_zero_values);
1845 }
1846 
1847 
1848 
1849 template <typename number>
1850 template <typename number2>
1851 inline void
1853  const std::vector<size_type> &col_indices,
1854  const std::vector<number2> & values,
1855  const bool elide_zero_values)
1856 {
1857  Assert(col_indices.size() == values.size(),
1858  ExcDimensionMismatch(col_indices.size(), values.size()));
1859 
1860  set(row,
1861  col_indices.size(),
1862  col_indices.data(),
1863  values.data(),
1864  elide_zero_values);
1865 }
1866 
1867 
1868 
1869 template <typename number>
1870 inline void
1872  const size_type j,
1873  const number value)
1874 {
1875  AssertIsFinite(value);
1876 
1877  if (value == number())
1878  return;
1879 
1880  const size_type index = cols->operator()(i, j);
1881 
1882  // it is allowed to add elements to the matrix that are not part of the
1883  // sparsity pattern, if the value to which we set it is zero
1884  if (index == SparsityPattern::invalid_entry)
1885  {
1886  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1887  ExcInvalidIndex(i, j));
1888  return;
1889  }
1890 
1891  val[index] += value;
1892 }
1893 
1894 
1895 
1896 template <typename number>
1897 template <typename number2>
1898 inline void
1899 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1900  const FullMatrix<number2> & values,
1901  const bool elide_zero_values)
1902 {
1903  Assert(indices.size() == values.m(),
1904  ExcDimensionMismatch(indices.size(), values.m()));
1905  Assert(values.m() == values.n(), ExcNotQuadratic());
1906 
1907  for (size_type i = 0; i < indices.size(); ++i)
1908  add(indices[i],
1909  indices.size(),
1910  indices.data(),
1911  &values(i, 0),
1912  elide_zero_values);
1913 }
1914 
1915 
1916 
1917 template <typename number>
1918 template <typename number2>
1919 inline void
1920 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1921  const std::vector<size_type> &col_indices,
1922  const FullMatrix<number2> & values,
1923  const bool elide_zero_values)
1924 {
1925  Assert(row_indices.size() == values.m(),
1926  ExcDimensionMismatch(row_indices.size(), values.m()));
1927  Assert(col_indices.size() == values.n(),
1928  ExcDimensionMismatch(col_indices.size(), values.n()));
1929 
1930  for (size_type i = 0; i < row_indices.size(); ++i)
1931  add(row_indices[i],
1932  col_indices.size(),
1933  col_indices.data(),
1934  &values(i, 0),
1935  elide_zero_values);
1936 }
1937 
1938 
1939 
1940 template <typename number>
1941 template <typename number2>
1942 inline void
1944  const std::vector<size_type> &col_indices,
1945  const std::vector<number2> & values,
1946  const bool elide_zero_values)
1947 {
1948  Assert(col_indices.size() == values.size(),
1949  ExcDimensionMismatch(col_indices.size(), values.size()));
1950 
1951  add(row,
1952  col_indices.size(),
1953  col_indices.data(),
1954  values.data(),
1955  elide_zero_values);
1956 }
1957 
1958 
1959 
1960 template <typename number>
1961 inline SparseMatrix<number> &
1962 SparseMatrix<number>::operator*=(const number factor)
1963 {
1964  Assert(cols != nullptr, ExcNotInitialized());
1965  Assert(val != nullptr, ExcNotInitialized());
1966 
1967  number * val_ptr = val.get();
1968  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1969 
1970  while (val_ptr != end_ptr)
1971  *val_ptr++ *= factor;
1972 
1973  return *this;
1974 }
1975 
1976 
1977 
1978 template <typename number>
1979 inline SparseMatrix<number> &
1980 SparseMatrix<number>::operator/=(const number factor)
1981 {
1982  Assert(cols != nullptr, ExcNotInitialized());
1983  Assert(val != nullptr, ExcNotInitialized());
1984  Assert(factor != number(), ExcDivideByZero());
1985 
1986  const number factor_inv = number(1.) / factor;
1987 
1988  number * val_ptr = val.get();
1989  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1990 
1991  while (val_ptr != end_ptr)
1992  *val_ptr++ *= factor_inv;
1993 
1994  return *this;
1995 }
1996 
1997 
1998 
1999 template <typename number>
2000 inline const number &
2002 {
2003  Assert(cols != nullptr, ExcNotInitialized());
2004  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2005  ExcInvalidIndex(i, j));
2006  return val[cols->operator()(i, j)];
2007 }
2008 
2009 
2010 
2011 template <typename number>
2012 inline number &
2014 {
2015  Assert(cols != nullptr, ExcNotInitialized());
2016  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2017  ExcInvalidIndex(i, j));
2018  return val[cols->operator()(i, j)];
2019 }
2020 
2021 
2022 
2023 template <typename number>
2024 inline number
2025 SparseMatrix<number>::el(const size_type i, const size_type j) const
2026 {
2027  Assert(cols != nullptr, ExcNotInitialized());
2028  const size_type index = cols->operator()(i, j);
2029 
2030  if (index != SparsityPattern::invalid_entry)
2031  return val[index];
2032  else
2033  return 0;
2034 }
2035 
2036 
2037 
2038 template <typename number>
2039 inline number
2041 {
2042  Assert(cols != nullptr, ExcNotInitialized());
2043  Assert(m() == n(), ExcNotQuadratic());
2044  AssertIndexRange(i, m());
2045 
2046  // Use that the first element in each row of a quadratic matrix is the main
2047  // diagonal
2048  return val[cols->rowstart[i]];
2049 }
2050 
2051 
2052 
2053 template <typename number>
2054 inline number &
2056 {
2057  Assert(cols != nullptr, ExcNotInitialized());
2058  Assert(m() == n(), ExcNotQuadratic());
2059  AssertIndexRange(i, m());
2060 
2061  // Use that the first element in each row of a quadratic matrix is the main
2062  // diagonal
2063  return val[cols->rowstart[i]];
2064 }
2065 
2066 
2067 
2068 template <typename number>
2069 template <typename ForwardIterator>
2070 void
2071 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2072  const ForwardIterator end)
2073 {
2074  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2075  ExcIteratorRange(std::distance(begin, end), m()));
2076 
2077  // for use in the inner loop, we define an alias to the type of the inner
2078  // iterators
2079  using inner_iterator =
2080  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2081  size_type row = 0;
2082  for (ForwardIterator i = begin; i != end; ++i, ++row)
2083  {
2084  const inner_iterator end_of_row = i->end();
2085  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2086  // write entries
2087  set(row, j->first, j->second);
2088  };
2089 }
2090 
2091 
2092 //---------------------------------------------------------------------------
2093 
2094 
2095 namespace SparseMatrixIterators
2096 {
2097  template <typename number>
2098  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2099  const std::size_t index_within_matrix)
2100  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2101  index_within_matrix)
2102  , matrix(matrix)
2103  {}
2104 
2105 
2106 
2107  template <typename number>
2108  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2109  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2110  , matrix(matrix)
2111  {}
2112 
2113 
2114 
2115  template <typename number>
2116  inline Accessor<number, true>::Accessor(
2118  : SparsityPatternIterators::Accessor(a)
2119  , matrix(&a.get_matrix())
2120  {}
2121 
2122 
2123 
2124  template <typename number>
2125  inline number
2126  Accessor<number, true>::value() const
2127  {
2128  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2129  return matrix->val[linear_index];
2130  }
2131 
2132 
2133 
2134  template <typename number>
2135  inline const typename Accessor<number, true>::MatrixType &
2136  Accessor<number, true>::get_matrix() const
2137  {
2138  return *matrix;
2139  }
2140 
2141 
2142 
2143  template <typename number>
2144  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2145  const bool)
2146  : accessor(accessor)
2147  {}
2148 
2149 
2150  template <typename number>
2151  inline Accessor<number, false>::Reference::operator number() const
2152  {
2153  AssertIndexRange(accessor->linear_index,
2154  accessor->matrix->n_nonzero_elements());
2155  return accessor->matrix->val[accessor->linear_index];
2156  }
2157 
2158 
2159 
2160  template <typename number>
2161  inline const typename Accessor<number, false>::Reference &
2162  Accessor<number, false>::Reference::operator=(const number n) const
2163  {
2164  AssertIndexRange(accessor->linear_index,
2165  accessor->matrix->n_nonzero_elements());
2166  accessor->matrix->val[accessor->linear_index] = n;
2167  return *this;
2168  }
2169 
2170 
2171 
2172  template <typename number>
2173  inline const typename Accessor<number, false>::Reference &
2174  Accessor<number, false>::Reference::operator+=(const number n) const
2175  {
2176  AssertIndexRange(accessor->linear_index,
2177  accessor->matrix->n_nonzero_elements());
2178  accessor->matrix->val[accessor->linear_index] += n;
2179  return *this;
2180  }
2181 
2182 
2183 
2184  template <typename number>
2185  inline const typename Accessor<number, false>::Reference &
2186  Accessor<number, false>::Reference::operator-=(const number n) const
2187  {
2188  AssertIndexRange(accessor->linear_index,
2189  accessor->matrix->n_nonzero_elements());
2190  accessor->matrix->val[accessor->linear_index] -= n;
2191  return *this;
2192  }
2193 
2194 
2195 
2196  template <typename number>
2197  inline const typename Accessor<number, false>::Reference &
2198  Accessor<number, false>::Reference::operator*=(const number n) const
2199  {
2200  AssertIndexRange(accessor->linear_index,
2201  accessor->matrix->n_nonzero_elements());
2202  accessor->matrix->val[accessor->linear_index] *= n;
2203  return *this;
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 Accessor<number, false>::Accessor(MatrixType * matrix,
2222  const std::size_t index)
2223  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2224  , matrix(matrix)
2225  {}
2226 
2227 
2228 
2229  template <typename number>
2230  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2231  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2232  , matrix(matrix)
2233  {}
2234 
2235 
2236 
2237  template <typename number>
2238  inline typename Accessor<number, false>::Reference
2239  Accessor<number, false>::value() const
2240  {
2241  return Reference(this, true);
2242  }
2243 
2244 
2245 
2246  template <typename number>
2247  inline typename Accessor<number, false>::MatrixType &
2248  Accessor<number, false>::get_matrix() const
2249  {
2250  return *matrix;
2251  }
2252 
2253 
2254 
2255  template <typename number, bool Constness>
2256  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2257  const std::size_t index)
2258  : accessor(matrix, index)
2259  {}
2260 
2261 
2262 
2263  template <typename number, bool Constness>
2264  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2265  : accessor(matrix)
2266  {}
2267 
2268 
2269 
2270  template <typename number, bool Constness>
2271  inline Iterator<number, Constness>::Iterator(
2273  : accessor(*i)
2274  {}
2275 
2276 
2277 
2278  template <typename number, bool Constness>
2279  inline const Iterator<number, Constness> &
2280  Iterator<number, Constness>::
2282  {
2283  accessor = *i;
2284  return *this;
2285  }
2286 
2287 
2288 
2289  template <typename number, bool Constness>
2290  inline Iterator<number, Constness> &
2292  {
2293  accessor.advance();
2294  return *this;
2295  }
2296 
2297 
2298  template <typename number, bool Constness>
2299  inline Iterator<number, Constness>
2301  {
2302  const Iterator iter = *this;
2303  accessor.advance();
2304  return iter;
2305  }
2306 
2307 
2308  template <typename number, bool Constness>
2309  inline const Accessor<number, Constness> &Iterator<number, Constness>::
2310  operator*() const
2311  {
2312  return accessor;
2313  }
2314 
2315 
2316  template <typename number, bool Constness>
2317  inline const Accessor<number, Constness> *Iterator<number, Constness>::
2318  operator->() const
2319  {
2320  return &accessor;
2321  }
2322 
2323 
2324  template <typename number, bool Constness>
2325  inline bool
2326  Iterator<number, Constness>::operator==(const Iterator &other) const
2327  {
2328  return (accessor == other.accessor);
2329  }
2330 
2331 
2332  template <typename number, bool Constness>
2333  inline bool
2334  Iterator<number, Constness>::operator!=(const Iterator &other) const
2335  {
2336  return !(*this == other);
2337  }
2338 
2339 
2340  template <typename number, bool Constness>
2341  inline bool
2342  Iterator<number, Constness>::operator<(const Iterator &other) const
2343  {
2344  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2345  ExcInternalError());
2346 
2347  return (accessor < other.accessor);
2348  }
2349 
2350 
2351  template <typename number, bool Constness>
2352  inline bool
2353  Iterator<number, Constness>::operator>(const Iterator &other) const
2354  {
2355  return (other < *this);
2356  }
2357 
2358 
2359  template <typename number, bool Constness>
2360  inline int
2361  Iterator<number, Constness>::operator-(const Iterator &other) const
2362  {
2363  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2364  ExcInternalError());
2365 
2366  return (*this)->linear_index - other->linear_index;
2367  }
2368 
2369 
2370 
2371  template <typename number, bool Constness>
2372  inline Iterator<number, Constness>
2374  {
2375  Iterator x = *this;
2376  for (size_type i = 0; i < n; ++i)
2377  ++x;
2378 
2379  return x;
2380  }
2381 
2382 } // namespace SparseMatrixIterators
2383 
2384 
2385 
2386 template <typename number>
2389 {
2390  return const_iterator(this, 0);
2391 }
2392 
2393 
2394 template <typename number>
2397 {
2398  return const_iterator(this);
2399 }
2400 
2401 
2402 template <typename number>
2403 inline typename SparseMatrix<number>::iterator
2405 {
2406  return iterator(this, 0);
2407 }
2408 
2409 
2410 template <typename number>
2411 inline typename SparseMatrix<number>::iterator
2413 {
2414  return iterator(this, cols->rowstart[cols->rows]);
2415 }
2416 
2417 
2418 template <typename number>
2421 {
2422  AssertIndexRange(r, m());
2423 
2424  return const_iterator(this, cols->rowstart[r]);
2425 }
2426 
2427 
2428 
2429 template <typename number>
2431 SparseMatrix<number>::end(const size_type r) const
2432 {
2433  AssertIndexRange(r, m());
2434 
2435  return const_iterator(this, cols->rowstart[r + 1]);
2436 }
2437 
2438 
2439 
2440 template <typename number>
2441 inline typename SparseMatrix<number>::iterator
2443 {
2444  AssertIndexRange(r, m());
2445 
2446  return iterator(this, cols->rowstart[r]);
2447 }
2448 
2449 
2450 
2451 template <typename number>
2452 inline typename SparseMatrix<number>::iterator
2454 {
2455  AssertIndexRange(r, m());
2456 
2457  return iterator(this, cols->rowstart[r + 1]);
2458 }
2459 
2460 
2461 
2462 template <typename number>
2463 template <class StreamType>
2464 inline void
2465 SparseMatrix<number>::print(StreamType &out,
2466  const bool across,
2467  const bool diagonal_first) const
2468 {
2469  Assert(cols != nullptr, ExcNotInitialized());
2470  Assert(val != nullptr, ExcNotInitialized());
2471 
2472  bool hanging_diagonal = false;
2473  number diagonal = number();
2474 
2475  for (size_type i = 0; i < cols->rows; ++i)
2476  {
2477  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2478  {
2479  if (!diagonal_first && i == cols->colnums[j])
2480  {
2481  diagonal = val[j];
2482  hanging_diagonal = true;
2483  }
2484  else
2485  {
2486  if (hanging_diagonal && cols->colnums[j] > i)
2487  {
2488  if (across)
2489  out << ' ' << i << ',' << i << ':' << diagonal;
2490  else
2491  out << '(' << i << ',' << i << ") " << diagonal
2492  << std::endl;
2493  hanging_diagonal = false;
2494  }
2495  if (across)
2496  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2497  else
2498  out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2499  << std::endl;
2500  }
2501  }
2502  if (hanging_diagonal)
2503  {
2504  if (across)
2505  out << ' ' << i << ',' << i << ':' << diagonal;
2506  else
2507  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2508  hanging_diagonal = false;
2509  }
2510  }
2511  if (across)
2512  out << std::endl;
2513 }
2514 
2515 
2516 template <typename number>
2517 inline void
2519 {
2520  // nothing to do here
2521 }
2522 
2523 
2524 
2525 template <typename number>
2526 inline void
2528 {
2529  // nothing to do here
2530 }
2531 
2532 # endif // DOXYGEN
2533 
2534 
2535 /*---------------------------- sparse_matrix.h ---------------------------*/
2536 
2538 
2539 #endif
2540 /*---------------------------- sparse_matrix.h ---------------------------*/
typename Accessor< number, Constness >::MatrixType MatrixType
size_type m() const
SparseMatrix & operator/=(const number factor)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
typename numbers::NumberTraits< number >::real_type real_type
void prepare_add()
const_iterator end() const
Contents is actually a matrix.
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
std::unique_ptr< number[]> val
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void set(const size_type i, const size_type j, const number value)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
number value_type
static const char V
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
static ::ExceptionBase & ExcNotInitialized()
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
static ::ExceptionBase & ExcDivideByZero()
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
std::unique_ptr< std::size_t[]> rowstart
Matrix is diagonal.
size_type n() const
std::string compress(const std::string &input)
Definition: utilities.cc:392
size_type n() const
T sum(const T &t, const MPI_Comm &mpi_communicator)
number diag_element(const size_type i) const
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
SparseMatrixIterators::Iterator< number, false > iterator
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2818
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
Accessor< number, Constness > accessor
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
VectorType::value_type * end(VectorType &V)
void prepare_set()
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
Expression operator>(const Expression &lhs, const Expression &rhs)
void add(const size_type i, const size_type j, const number value)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SparseMatrix< double > SparseMatrix
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
types::global_dof_index size_type
std::unique_ptr< size_type[]> colnums
static ::ExceptionBase & ExcNotQuadratic()
unsigned int global_dof_index
Definition: types.h:76
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
number el(const size_type i, const size_type j) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
VectorType::value_type * begin(VectorType &V)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2792
size_type m() const
static const size_type invalid_entry
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
const SparsityPattern & get_sparsity_pattern() const
const_iterator begin() const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
const number & operator()(const size_type i, const size_type j) const
SparseMatrixIterators::Iterator< number, true > const_iterator
#define AssertIsFinite(number)
Definition: exceptions.h:1722
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::size_t max_len
const Accessor< number, Constness > & value_type
SparseMatrix & operator*=(const number factor)
static ::ExceptionBase & ExcInternalError()