Reference documentation for deal.II version Git bed997f895 2020-09-22 11:49:20 -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
55  sum(const SparseMatrix<Number> &, const MPI_Comm &, SparseMatrix<Number> &);
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 
979  template <typename somenumber>
980  void
981  copy_from(const FullMatrix<somenumber> &matrix);
982 
983 # ifdef DEAL_II_WITH_TRILINOS
984 
994  copy_from(const TrilinosWrappers::SparseMatrix &matrix);
995 # endif
996 
1008  template <typename somenumber>
1009  void
1010  add(const number factor, const SparseMatrix<somenumber> &matrix);
1011 
1013 
1017 
1031  const number &
1032  operator()(const size_type i, const size_type j) const;
1033 
1037  number &
1038  operator()(const size_type i, const size_type j);
1039 
1052  number
1053  el(const size_type i, const size_type j) const;
1054 
1064  number
1065  diag_element(const size_type i) const;
1066 
1071  number &
1072  diag_element(const size_type i);
1073 
1075 
1095  template <class OutVector, class InVector>
1096  void
1097  vmult(OutVector &dst, const InVector &src) const;
1098 
1114  template <class OutVector, class InVector>
1115  void
1116  Tvmult(OutVector &dst, const InVector &src) const;
1117 
1134  template <class OutVector, class InVector>
1135  void
1136  vmult_add(OutVector &dst, const InVector &src) const;
1137 
1153  template <class OutVector, class InVector>
1154  void
1155  Tvmult_add(OutVector &dst, const InVector &src) const;
1156 
1174  template <typename somenumber>
1175  somenumber
1176  matrix_norm_square(const Vector<somenumber> &v) const;
1177 
1183  template <typename somenumber>
1184  somenumber
1185  matrix_scalar_product(const Vector<somenumber> &u,
1186  const Vector<somenumber> &v) const;
1187 
1197  template <typename somenumber>
1198  somenumber
1199  residual(Vector<somenumber> & dst,
1200  const Vector<somenumber> &x,
1201  const Vector<somenumber> &b) const;
1202 
1238  template <typename numberB, typename numberC>
1239  void
1240  mmult(SparseMatrix<numberC> & C,
1241  const SparseMatrix<numberB> &B,
1242  const Vector<number> & V = Vector<number>(),
1243  const bool rebuild_sparsity_pattern = true) const;
1244 
1269  template <typename numberB, typename numberC>
1270  void
1271  Tmmult(SparseMatrix<numberC> & C,
1272  const SparseMatrix<numberB> &B,
1273  const Vector<number> & V = Vector<number>(),
1274  const bool rebuild_sparsity_pattern = true) const;
1275 
1277 
1281 
1289  real_type
1290  l1_norm() const;
1291 
1299  real_type
1300  linfty_norm() const;
1301 
1306  real_type
1307  frobenius_norm() const;
1309 
1313 
1319  template <typename somenumber>
1320  void
1321  precondition_Jacobi(Vector<somenumber> & dst,
1322  const Vector<somenumber> &src,
1323  const number omega = 1.) const;
1324 
1331  template <typename somenumber>
1332  void
1333  precondition_SSOR(Vector<somenumber> & dst,
1334  const Vector<somenumber> & src,
1335  const number omega = 1.,
1336  const std::vector<std::size_t> &pos_right_of_diagonal =
1337  std::vector<std::size_t>()) const;
1338 
1342  template <typename somenumber>
1343  void
1344  precondition_SOR(Vector<somenumber> & dst,
1345  const Vector<somenumber> &src,
1346  const number om = 1.) const;
1347 
1351  template <typename somenumber>
1352  void
1353  precondition_TSOR(Vector<somenumber> & dst,
1354  const Vector<somenumber> &src,
1355  const number om = 1.) const;
1356 
1362  template <typename somenumber>
1363  void
1364  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1365 
1370  template <typename somenumber>
1371  void
1372  SOR(Vector<somenumber> &v, const number om = 1.) const;
1373 
1378  template <typename somenumber>
1379  void
1380  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1381 
1392  template <typename somenumber>
1393  void
1394  PSOR(Vector<somenumber> & v,
1395  const std::vector<size_type> &permutation,
1396  const std::vector<size_type> &inverse_permutation,
1397  const number om = 1.) const;
1398 
1409  template <typename somenumber>
1410  void
1411  TPSOR(Vector<somenumber> & v,
1412  const std::vector<size_type> &permutation,
1413  const std::vector<size_type> &inverse_permutation,
1414  const number om = 1.) const;
1415 
1421  template <typename somenumber>
1422  void
1423  Jacobi_step(Vector<somenumber> & v,
1424  const Vector<somenumber> &b,
1425  const number om = 1.) const;
1426 
1431  template <typename somenumber>
1432  void
1433  SOR_step(Vector<somenumber> & v,
1434  const Vector<somenumber> &b,
1435  const number om = 1.) const;
1436 
1441  template <typename somenumber>
1442  void
1443  TSOR_step(Vector<somenumber> & v,
1444  const Vector<somenumber> &b,
1445  const number om = 1.) const;
1446 
1451  template <typename somenumber>
1452  void
1453  SSOR_step(Vector<somenumber> & v,
1454  const Vector<somenumber> &b,
1455  const number om = 1.) const;
1457 
1461 
1469  begin() const;
1470 
1474  iterator
1475  begin();
1476 
1481  end() const;
1482 
1486  iterator
1487  end();
1488 
1499  begin(const size_type r) const;
1500 
1504  iterator
1505  begin(const size_type r);
1506 
1517  end(const size_type r) const;
1518 
1522  iterator
1523  end(const size_type r);
1525 
1529 
1541  template <class StreamType>
1542  void
1543  print(StreamType &out,
1544  const bool across = false,
1545  const bool diagonal_first = true) const;
1546 
1567  void
1568  print_formatted(std::ostream & out,
1569  const unsigned int precision = 3,
1570  const bool scientific = true,
1571  const unsigned int width = 0,
1572  const char * zero_string = " ",
1573  const double denominator = 1.) const;
1574 
1580  void
1581  print_pattern(std::ostream &out, const double threshold = 0.) const;
1582 
1591  void
1592  print_as_numpy_arrays(std::ostream & out,
1593  const unsigned int precision = 9) const;
1594 
1605  void
1606  block_write(std::ostream &out) const;
1607 
1624  void
1625  block_read(std::istream &in);
1627 
1635  DeclException2(ExcInvalidIndex,
1636  int,
1637  int,
1638  << "You are trying to access the matrix entry with index <"
1639  << arg1 << ',' << arg2
1640  << ">, but this entry does not exist in the sparsity pattern "
1641  "of this matrix."
1642  "\n\n"
1643  "The most common cause for this problem is that you used "
1644  "a method to build the sparsity pattern that did not "
1645  "(completely) take into account all of the entries you "
1646  "will later try to write into. An example would be "
1647  "building a sparsity pattern that does not include "
1648  "the entries you will write into due to constraints "
1649  "on degrees of freedom such as hanging nodes or periodic "
1650  "boundary conditions. In such cases, building the "
1651  "sparsity pattern will succeed, but you will get errors "
1652  "such as the current one at one point or other when "
1653  "trying to write into the entries of the matrix.");
1657  DeclExceptionMsg(ExcDifferentSparsityPatterns,
1658  "When copying one sparse matrix into another, "
1659  "or when adding one sparse matrix to another, "
1660  "both matrices need to refer to the same "
1661  "sparsity pattern.");
1665  DeclException2(ExcIteratorRange,
1666  int,
1667  int,
1668  << "The iterators denote a range of " << arg1
1669  << " elements, but the given number of rows was " << arg2);
1673  DeclExceptionMsg(ExcSourceEqualsDestination,
1674  "You are attempting an operation on two matrices that "
1675  "are the same object, but the operation requires that the "
1676  "two objects are in fact different.");
1678 
1679 protected:
1690  void
1691  prepare_add();
1692 
1697  void
1698  prepare_set();
1699 
1700 private:
1707 
1715  std::unique_ptr<number[]> val;
1716 
1723  std::size_t max_len;
1724 
1725  // make all other sparse matrices friends
1726  template <typename somenumber>
1727  friend class SparseMatrix;
1728  template <typename somenumber>
1730  template <typename>
1731  friend class SparseILU;
1732 
1733  // To allow it calling private prepare_add() and prepare_set().
1734  template <typename>
1735  friend class BlockMatrixBase;
1736 
1737  // Also give access to internal details to the iterator/accessor classes.
1738  template <typename, bool>
1740  template <typename, bool>
1742 
1743 # ifdef DEAL_II_WITH_MPI
1744  // Give access to internal datastructures to perform MPI operations.
1745  template <typename Number>
1746  friend void
1748  const MPI_Comm &,
1750 # endif
1751 };
1752 
1753 # ifndef DOXYGEN
1754 /*---------------------- Inline functions -----------------------------------*/
1755 
1756 
1757 
1758 template <typename number>
1759 inline typename SparseMatrix<number>::size_type
1761 {
1762  Assert(cols != nullptr, ExcNotInitialized());
1763  return cols->rows;
1764 }
1765 
1766 
1767 template <typename number>
1768 inline typename SparseMatrix<number>::size_type
1770 {
1771  Assert(cols != nullptr, ExcNotInitialized());
1772  return cols->cols;
1773 }
1774 
1775 
1776 // Inline the set() and add() functions, since they will be called frequently.
1777 template <typename number>
1778 inline void
1780  const size_type j,
1781  const number value)
1782 {
1783  AssertIsFinite(value);
1784 
1785  const size_type index = cols->operator()(i, j);
1786 
1787  // it is allowed to set elements of the matrix that are not part of the
1788  // sparsity pattern, if the value to which we set it is zero
1789  if (index == SparsityPattern::invalid_entry)
1790  {
1791  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1792  ExcInvalidIndex(i, j));
1793  return;
1794  }
1795 
1796  val[index] = value;
1797 }
1798 
1799 
1800 
1801 template <typename number>
1802 template <typename number2>
1803 inline void
1804 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1805  const FullMatrix<number2> & values,
1806  const bool elide_zero_values)
1807 {
1808  Assert(indices.size() == values.m(),
1809  ExcDimensionMismatch(indices.size(), values.m()));
1810  Assert(values.m() == values.n(), ExcNotQuadratic());
1811 
1812  for (size_type i = 0; i < indices.size(); ++i)
1813  set(indices[i],
1814  indices.size(),
1815  indices.data(),
1816  &values(i, 0),
1817  elide_zero_values);
1818 }
1819 
1820 
1821 
1822 template <typename number>
1823 template <typename number2>
1824 inline void
1825 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1826  const std::vector<size_type> &col_indices,
1827  const FullMatrix<number2> & values,
1828  const bool elide_zero_values)
1829 {
1830  Assert(row_indices.size() == values.m(),
1831  ExcDimensionMismatch(row_indices.size(), values.m()));
1832  Assert(col_indices.size() == values.n(),
1833  ExcDimensionMismatch(col_indices.size(), values.n()));
1834 
1835  for (size_type i = 0; i < row_indices.size(); ++i)
1836  set(row_indices[i],
1837  col_indices.size(),
1838  col_indices.data(),
1839  &values(i, 0),
1840  elide_zero_values);
1841 }
1842 
1843 
1844 
1845 template <typename number>
1846 template <typename number2>
1847 inline void
1849  const std::vector<size_type> &col_indices,
1850  const std::vector<number2> & values,
1851  const bool elide_zero_values)
1852 {
1853  Assert(col_indices.size() == values.size(),
1854  ExcDimensionMismatch(col_indices.size(), values.size()));
1855 
1856  set(row,
1857  col_indices.size(),
1858  col_indices.data(),
1859  values.data(),
1860  elide_zero_values);
1861 }
1862 
1863 
1864 
1865 template <typename number>
1866 inline void
1868  const size_type j,
1869  const number value)
1870 {
1871  AssertIsFinite(value);
1872 
1873  if (value == number())
1874  return;
1875 
1876  const size_type index = cols->operator()(i, j);
1877 
1878  // it is allowed to add elements to the matrix that are not part of the
1879  // sparsity pattern, if the value to which we set it is zero
1880  if (index == SparsityPattern::invalid_entry)
1881  {
1882  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1883  ExcInvalidIndex(i, j));
1884  return;
1885  }
1886 
1887  val[index] += value;
1888 }
1889 
1890 
1891 
1892 template <typename number>
1893 template <typename number2>
1894 inline void
1895 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1896  const FullMatrix<number2> & values,
1897  const bool elide_zero_values)
1898 {
1899  Assert(indices.size() == values.m(),
1900  ExcDimensionMismatch(indices.size(), values.m()));
1901  Assert(values.m() == values.n(), ExcNotQuadratic());
1902 
1903  for (size_type i = 0; i < indices.size(); ++i)
1904  add(indices[i],
1905  indices.size(),
1906  indices.data(),
1907  &values(i, 0),
1908  elide_zero_values);
1909 }
1910 
1911 
1912 
1913 template <typename number>
1914 template <typename number2>
1915 inline void
1916 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1917  const std::vector<size_type> &col_indices,
1918  const FullMatrix<number2> & values,
1919  const bool elide_zero_values)
1920 {
1921  Assert(row_indices.size() == values.m(),
1922  ExcDimensionMismatch(row_indices.size(), values.m()));
1923  Assert(col_indices.size() == values.n(),
1924  ExcDimensionMismatch(col_indices.size(), values.n()));
1925 
1926  for (size_type i = 0; i < row_indices.size(); ++i)
1927  add(row_indices[i],
1928  col_indices.size(),
1929  col_indices.data(),
1930  &values(i, 0),
1931  elide_zero_values);
1932 }
1933 
1934 
1935 
1936 template <typename number>
1937 template <typename number2>
1938 inline void
1940  const std::vector<size_type> &col_indices,
1941  const std::vector<number2> & values,
1942  const bool elide_zero_values)
1943 {
1944  Assert(col_indices.size() == values.size(),
1945  ExcDimensionMismatch(col_indices.size(), values.size()));
1946 
1947  add(row,
1948  col_indices.size(),
1949  col_indices.data(),
1950  values.data(),
1951  elide_zero_values);
1952 }
1953 
1954 
1955 
1956 template <typename number>
1957 inline SparseMatrix<number> &
1958 SparseMatrix<number>::operator*=(const number factor)
1959 {
1960  Assert(cols != nullptr, ExcNotInitialized());
1961  Assert(val != nullptr, ExcNotInitialized());
1962 
1963  number * val_ptr = val.get();
1964  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1965 
1966  while (val_ptr != end_ptr)
1967  *val_ptr++ *= factor;
1968 
1969  return *this;
1970 }
1971 
1972 
1973 
1974 template <typename number>
1975 inline SparseMatrix<number> &
1976 SparseMatrix<number>::operator/=(const number factor)
1977 {
1978  Assert(cols != nullptr, ExcNotInitialized());
1979  Assert(val != nullptr, ExcNotInitialized());
1980  Assert(factor != number(), ExcDivideByZero());
1981 
1982  const number factor_inv = number(1.) / factor;
1983 
1984  number * val_ptr = val.get();
1985  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1986 
1987  while (val_ptr != end_ptr)
1988  *val_ptr++ *= factor_inv;
1989 
1990  return *this;
1991 }
1992 
1993 
1994 
1995 template <typename number>
1996 inline const number &
1998 {
1999  Assert(cols != nullptr, ExcNotInitialized());
2000  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2001  ExcInvalidIndex(i, j));
2002  return val[cols->operator()(i, j)];
2003 }
2004 
2005 
2006 
2007 template <typename number>
2008 inline number &
2010 {
2011  Assert(cols != nullptr, ExcNotInitialized());
2012  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2013  ExcInvalidIndex(i, j));
2014  return val[cols->operator()(i, j)];
2015 }
2016 
2017 
2018 
2019 template <typename number>
2020 inline number
2021 SparseMatrix<number>::el(const size_type i, const size_type j) const
2022 {
2023  Assert(cols != nullptr, ExcNotInitialized());
2024  const size_type index = cols->operator()(i, j);
2025 
2026  if (index != SparsityPattern::invalid_entry)
2027  return val[index];
2028  else
2029  return 0;
2030 }
2031 
2032 
2033 
2034 template <typename number>
2035 inline number
2037 {
2038  Assert(cols != nullptr, ExcNotInitialized());
2039  Assert(m() == n(), ExcNotQuadratic());
2040  AssertIndexRange(i, m());
2041 
2042  // Use that the first element in each row of a quadratic matrix is the main
2043  // diagonal
2044  return val[cols->rowstart[i]];
2045 }
2046 
2047 
2048 
2049 template <typename number>
2050 inline number &
2052 {
2053  Assert(cols != nullptr, ExcNotInitialized());
2054  Assert(m() == n(), ExcNotQuadratic());
2055  AssertIndexRange(i, m());
2056 
2057  // Use that the first element in each row of a quadratic matrix is the main
2058  // diagonal
2059  return val[cols->rowstart[i]];
2060 }
2061 
2062 
2063 
2064 template <typename number>
2065 template <typename ForwardIterator>
2066 void
2067 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2068  const ForwardIterator end)
2069 {
2070  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2071  ExcIteratorRange(std::distance(begin, end), m()));
2072 
2073  // for use in the inner loop, we define an alias to the type of the inner
2074  // iterators
2075  using inner_iterator =
2076  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2077  size_type row = 0;
2078  for (ForwardIterator i = begin; i != end; ++i, ++row)
2079  {
2080  const inner_iterator end_of_row = i->end();
2081  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2082  // write entries
2083  set(row, j->first, j->second);
2084  };
2085 }
2086 
2087 
2088 //---------------------------------------------------------------------------
2089 
2090 
2091 namespace SparseMatrixIterators
2092 {
2093  template <typename number>
2094  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2095  const std::size_t index_within_matrix)
2096  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2097  index_within_matrix)
2098  , matrix(matrix)
2099  {}
2100 
2101 
2102 
2103  template <typename number>
2104  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2105  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2106  , matrix(matrix)
2107  {}
2108 
2109 
2110 
2111  template <typename number>
2112  inline Accessor<number, true>::Accessor(
2114  : SparsityPatternIterators::Accessor(a)
2115  , matrix(&a.get_matrix())
2116  {}
2117 
2118 
2119 
2120  template <typename number>
2121  inline number
2122  Accessor<number, true>::value() const
2123  {
2124  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2125  return matrix->val[linear_index];
2126  }
2127 
2128 
2129 
2130  template <typename number>
2131  inline const typename Accessor<number, true>::MatrixType &
2132  Accessor<number, true>::get_matrix() const
2133  {
2134  return *matrix;
2135  }
2136 
2137 
2138 
2139  template <typename number>
2140  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2141  const bool)
2142  : accessor(accessor)
2143  {}
2144 
2145 
2146  template <typename number>
2147  inline Accessor<number, false>::Reference::operator number() const
2148  {
2149  AssertIndexRange(accessor->linear_index,
2150  accessor->matrix->n_nonzero_elements());
2151  return accessor->matrix->val[accessor->linear_index];
2152  }
2153 
2154 
2155 
2156  template <typename number>
2157  inline const typename Accessor<number, false>::Reference &
2158  Accessor<number, false>::Reference::operator=(const number n) const
2159  {
2160  AssertIndexRange(accessor->linear_index,
2161  accessor->matrix->n_nonzero_elements());
2162  accessor->matrix->val[accessor->linear_index] = n;
2163  return *this;
2164  }
2165 
2166 
2167 
2168  template <typename number>
2169  inline const typename Accessor<number, false>::Reference &
2170  Accessor<number, false>::Reference::operator+=(const number n) const
2171  {
2172  AssertIndexRange(accessor->linear_index,
2173  accessor->matrix->n_nonzero_elements());
2174  accessor->matrix->val[accessor->linear_index] += n;
2175  return *this;
2176  }
2177 
2178 
2179 
2180  template <typename number>
2181  inline const typename Accessor<number, false>::Reference &
2182  Accessor<number, false>::Reference::operator-=(const number n) const
2183  {
2184  AssertIndexRange(accessor->linear_index,
2185  accessor->matrix->n_nonzero_elements());
2186  accessor->matrix->val[accessor->linear_index] -= n;
2187  return *this;
2188  }
2189 
2190 
2191 
2192  template <typename number>
2193  inline const typename Accessor<number, false>::Reference &
2194  Accessor<number, false>::Reference::operator*=(const number n) const
2195  {
2196  AssertIndexRange(accessor->linear_index,
2197  accessor->matrix->n_nonzero_elements());
2198  accessor->matrix->val[accessor->linear_index] *= n;
2199  return *this;
2200  }
2201 
2202 
2203 
2204  template <typename number>
2205  inline const typename Accessor<number, false>::Reference &
2206  Accessor<number, false>::Reference::operator/=(const number n) const
2207  {
2208  AssertIndexRange(accessor->linear_index,
2209  accessor->matrix->n_nonzero_elements());
2210  accessor->matrix->val[accessor->linear_index] /= n;
2211  return *this;
2212  }
2213 
2214 
2215 
2216  template <typename number>
2217  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2218  const std::size_t index)
2219  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2220  , matrix(matrix)
2221  {}
2222 
2223 
2224 
2225  template <typename number>
2226  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2227  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2228  , matrix(matrix)
2229  {}
2230 
2231 
2232 
2233  template <typename number>
2234  inline typename Accessor<number, false>::Reference
2235  Accessor<number, false>::value() const
2236  {
2237  return Reference(this, true);
2238  }
2239 
2240 
2241 
2242  template <typename number>
2243  inline typename Accessor<number, false>::MatrixType &
2244  Accessor<number, false>::get_matrix() const
2245  {
2246  return *matrix;
2247  }
2248 
2249 
2250 
2251  template <typename number, bool Constness>
2252  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2253  const std::size_t index)
2254  : accessor(matrix, index)
2255  {}
2256 
2257 
2258 
2259  template <typename number, bool Constness>
2260  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2261  : accessor(matrix)
2262  {}
2263 
2264 
2265 
2266  template <typename number, bool Constness>
2267  inline Iterator<number, Constness>::Iterator(
2269  : accessor(*i)
2270  {}
2271 
2272 
2273 
2274  template <typename number, bool Constness>
2275  inline const Iterator<number, Constness> &
2276  Iterator<number, Constness>::
2278  {
2279  accessor = *i;
2280  return *this;
2281  }
2282 
2283 
2284 
2285  template <typename number, bool Constness>
2286  inline Iterator<number, Constness> &
2288  {
2289  accessor.advance();
2290  return *this;
2291  }
2292 
2293 
2294  template <typename number, bool Constness>
2295  inline Iterator<number, Constness>
2297  {
2298  const Iterator iter = *this;
2299  accessor.advance();
2300  return iter;
2301  }
2302 
2303 
2304  template <typename number, bool Constness>
2305  inline const Accessor<number, Constness> &Iterator<number, Constness>::
2306  operator*() const
2307  {
2308  return accessor;
2309  }
2310 
2311 
2312  template <typename number, bool Constness>
2313  inline const Accessor<number, Constness> *Iterator<number, Constness>::
2314  operator->() const
2315  {
2316  return &accessor;
2317  }
2318 
2319 
2320  template <typename number, bool Constness>
2321  inline bool
2322  Iterator<number, Constness>::operator==(const Iterator &other) const
2323  {
2324  return (accessor == other.accessor);
2325  }
2326 
2327 
2328  template <typename number, bool Constness>
2329  inline bool
2330  Iterator<number, Constness>::operator!=(const Iterator &other) const
2331  {
2332  return !(*this == other);
2333  }
2334 
2335 
2336  template <typename number, bool Constness>
2337  inline bool
2338  Iterator<number, Constness>::operator<(const Iterator &other) const
2339  {
2340  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2341  ExcInternalError());
2342 
2343  return (accessor < other.accessor);
2344  }
2345 
2346 
2347  template <typename number, bool Constness>
2348  inline bool
2349  Iterator<number, Constness>::operator>(const Iterator &other) const
2350  {
2351  return (other < *this);
2352  }
2353 
2354 
2355  template <typename number, bool Constness>
2356  inline int
2357  Iterator<number, Constness>::operator-(const Iterator &other) const
2358  {
2359  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2360  ExcInternalError());
2361 
2362  return (*this)->linear_index - other->linear_index;
2363  }
2364 
2365 
2366 
2367  template <typename number, bool Constness>
2368  inline Iterator<number, Constness>
2370  {
2371  Iterator x = *this;
2372  for (size_type i = 0; i < n; ++i)
2373  ++x;
2374 
2375  return x;
2376  }
2377 
2378 } // namespace SparseMatrixIterators
2379 
2380 
2381 
2382 template <typename number>
2385 {
2386  return const_iterator(this, 0);
2387 }
2388 
2389 
2390 template <typename number>
2393 {
2394  return const_iterator(this);
2395 }
2396 
2397 
2398 template <typename number>
2399 inline typename SparseMatrix<number>::iterator
2401 {
2402  return iterator(this, 0);
2403 }
2404 
2405 
2406 template <typename number>
2407 inline typename SparseMatrix<number>::iterator
2409 {
2410  return iterator(this, cols->rowstart[cols->rows]);
2411 }
2412 
2413 
2414 template <typename number>
2417 {
2418  AssertIndexRange(r, m());
2419 
2420  return const_iterator(this, cols->rowstart[r]);
2421 }
2422 
2423 
2424 
2425 template <typename number>
2427 SparseMatrix<number>::end(const size_type r) const
2428 {
2429  AssertIndexRange(r, m());
2430 
2431  return const_iterator(this, cols->rowstart[r + 1]);
2432 }
2433 
2434 
2435 
2436 template <typename number>
2437 inline typename SparseMatrix<number>::iterator
2439 {
2440  AssertIndexRange(r, m());
2441 
2442  return iterator(this, cols->rowstart[r]);
2443 }
2444 
2445 
2446 
2447 template <typename number>
2448 inline typename SparseMatrix<number>::iterator
2450 {
2451  AssertIndexRange(r, m());
2452 
2453  return iterator(this, cols->rowstart[r + 1]);
2454 }
2455 
2456 
2457 
2458 template <typename number>
2459 template <class StreamType>
2460 inline void
2461 SparseMatrix<number>::print(StreamType &out,
2462  const bool across,
2463  const bool diagonal_first) const
2464 {
2465  Assert(cols != nullptr, ExcNotInitialized());
2466  Assert(val != nullptr, ExcNotInitialized());
2467 
2468  bool hanging_diagonal = false;
2469  number diagonal = number();
2470 
2471  for (size_type i = 0; i < cols->rows; ++i)
2472  {
2473  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2474  {
2475  if (!diagonal_first && i == cols->colnums[j])
2476  {
2477  diagonal = val[j];
2478  hanging_diagonal = true;
2479  }
2480  else
2481  {
2482  if (hanging_diagonal && cols->colnums[j] > i)
2483  {
2484  if (across)
2485  out << ' ' << i << ',' << i << ':' << diagonal;
2486  else
2487  out << '(' << i << ',' << i << ") " << diagonal
2488  << std::endl;
2489  hanging_diagonal = false;
2490  }
2491  if (across)
2492  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2493  else
2494  out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2495  << std::endl;
2496  }
2497  }
2498  if (hanging_diagonal)
2499  {
2500  if (across)
2501  out << ' ' << i << ',' << i << ':' << diagonal;
2502  else
2503  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2504  hanging_diagonal = false;
2505  }
2506  }
2507  if (across)
2508  out << std::endl;
2509 }
2510 
2511 
2512 template <typename number>
2513 inline void
2515 {
2516  // nothing to do here
2517 }
2518 
2519 
2520 
2521 template <typename number>
2522 inline void
2524 {
2525  // nothing to do here
2526 }
2527 
2528 # endif // DOXYGEN
2529 
2530 
2531 /*---------------------------- sparse_matrix.h ---------------------------*/
2532 
2534 
2535 #endif
2536 /*---------------------------- 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
#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:1636
number value_type
static const char V
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
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:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
constexpr SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
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:2811
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
Accessor< number, Constness > accessor
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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
Definition: cuda.h:31
number el(const size_type i, const size_type j) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
VectorType::value_type * begin(VectorType &V)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2785
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:1667
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()