Reference documentation for deal.II version Git c0262b0158 2021-10-21 11:22:17 -0600
\(\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 - 2021 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
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 &
397  operator++();
398 
402  Iterator
403  operator++(int);
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 
591  SparseMatrix();
592 
601  SparseMatrix(const SparseMatrix &);
602 
610  SparseMatrix(SparseMatrix<number> &&m) noexcept;
611 
625  explicit SparseMatrix(const SparsityPattern &sparsity);
626 
633  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
634 
639  virtual ~SparseMatrix() override;
640 
651  operator=(const SparseMatrix<number> &);
652 
658  operator=(SparseMatrix<number> &&m) noexcept;
659 
667  operator=(const IdentityMatrix &id);
668 
680  SparseMatrix &
681  operator=(const double d);
682 
696  virtual void
697  reinit(const SparsityPattern &sparsity);
698 
704  virtual void
705  clear();
707 
715  bool
716  empty() const;
717 
722  size_type
723  m() const;
724 
729  size_type
730  n() const;
731 
735  size_type
736  get_row_length(const size_type row) const;
737 
743  std::size_t
744  n_nonzero_elements() const;
745 
755  std::size_t
756  n_actually_nonzero_elements(const double threshold = 0.) const;
757 
766  const SparsityPattern &
767  get_sparsity_pattern() const;
768 
773  std::size_t
774  memory_consumption() const;
775 
780 
782 
791  void
792  set(const size_type i, const size_type j, const number value);
793 
809  template <typename number2>
810  void
811  set(const std::vector<size_type> &indices,
812  const FullMatrix<number2> & full_matrix,
813  const bool elide_zero_values = false);
814 
820  template <typename number2>
821  void
822  set(const std::vector<size_type> &row_indices,
823  const std::vector<size_type> &col_indices,
824  const FullMatrix<number2> & full_matrix,
825  const bool elide_zero_values = false);
826 
837  template <typename number2>
838  void
839  set(const size_type row,
840  const std::vector<size_type> &col_indices,
841  const std::vector<number2> & values,
842  const bool elide_zero_values = false);
843 
853  template <typename number2>
854  void
855  set(const size_type row,
856  const size_type n_cols,
857  const size_type *col_indices,
858  const number2 * values,
859  const bool elide_zero_values = false);
860 
866  void
867  add(const size_type i, const size_type j, const number value);
868 
883  template <typename number2>
884  void
885  add(const std::vector<size_type> &indices,
886  const FullMatrix<number2> & full_matrix,
887  const bool elide_zero_values = true);
888 
894  template <typename number2>
895  void
896  add(const std::vector<size_type> &row_indices,
897  const std::vector<size_type> &col_indices,
898  const FullMatrix<number2> & full_matrix,
899  const bool elide_zero_values = true);
900 
910  template <typename number2>
911  void
912  add(const size_type row,
913  const std::vector<size_type> &col_indices,
914  const std::vector<number2> & values,
915  const bool elide_zero_values = true);
916 
926  template <typename number2>
927  void
928  add(const size_type row,
929  const size_type n_cols,
930  const size_type *col_indices,
931  const number2 * values,
932  const bool elide_zero_values = true,
933  const bool col_indices_are_sorted = false);
934 
938  SparseMatrix &
939  operator*=(const number factor);
940 
944  SparseMatrix &
945  operator/=(const number factor);
946 
959  void
960  symmetrize();
961 
978  template <typename somenumber>
980  copy_from(const SparseMatrix<somenumber> &source);
981 
998  template <typename ForwardIterator>
999  void
1000  copy_from(const ForwardIterator begin, const ForwardIterator end);
1001 
1011  template <typename somenumber>
1012  void
1013  copy_from(const FullMatrix<somenumber> &matrix);
1014 
1015 # ifdef DEAL_II_WITH_TRILINOS
1016 
1026  copy_from(const TrilinosWrappers::SparseMatrix &matrix);
1027 # endif
1028 
1040  template <typename somenumber>
1041  void
1042  add(const number factor, const SparseMatrix<somenumber> &matrix);
1043 
1045 
1049 
1063  const number &
1064  operator()(const size_type i, const size_type j) const;
1065 
1069  number &
1070  operator()(const size_type i, const size_type j);
1071 
1084  number
1085  el(const size_type i, const size_type j) const;
1086 
1096  number
1097  diag_element(const size_type i) const;
1098 
1103  number &
1104  diag_element(const size_type i);
1105 
1107 
1127  template <class OutVector, class InVector>
1128  void
1129  vmult(OutVector &dst, const InVector &src) const;
1130 
1146  template <class OutVector, class InVector>
1147  void
1148  Tvmult(OutVector &dst, const InVector &src) const;
1149 
1166  template <class OutVector, class InVector>
1167  void
1168  vmult_add(OutVector &dst, const InVector &src) const;
1169 
1185  template <class OutVector, class InVector>
1186  void
1187  Tvmult_add(OutVector &dst, const InVector &src) const;
1188 
1206  template <typename somenumber>
1207  somenumber
1208  matrix_norm_square(const Vector<somenumber> &v) const;
1209 
1215  template <typename somenumber>
1216  somenumber
1217  matrix_scalar_product(const Vector<somenumber> &u,
1218  const Vector<somenumber> &v) const;
1219 
1229  template <typename somenumber>
1230  somenumber
1231  residual(Vector<somenumber> & dst,
1232  const Vector<somenumber> &x,
1233  const Vector<somenumber> &b) const;
1234 
1270  template <typename numberB, typename numberC>
1271  void
1272  mmult(SparseMatrix<numberC> & C,
1273  const SparseMatrix<numberB> &B,
1274  const Vector<number> & V = Vector<number>(),
1275  const bool rebuild_sparsity_pattern = true) const;
1276 
1301  template <typename numberB, typename numberC>
1302  void
1303  Tmmult(SparseMatrix<numberC> & C,
1304  const SparseMatrix<numberB> &B,
1305  const Vector<number> & V = Vector<number>(),
1306  const bool rebuild_sparsity_pattern = true) const;
1307 
1309 
1313 
1321  real_type
1322  l1_norm() const;
1323 
1331  real_type
1332  linfty_norm() const;
1333 
1338  real_type
1339  frobenius_norm() const;
1341 
1345 
1351  template <typename somenumber>
1352  void
1353  precondition_Jacobi(Vector<somenumber> & dst,
1354  const Vector<somenumber> &src,
1355  const number omega = 1.) const;
1356 
1363  template <typename somenumber>
1364  void
1365  precondition_SSOR(Vector<somenumber> & dst,
1366  const Vector<somenumber> & src,
1367  const number omega = 1.,
1368  const std::vector<std::size_t> &pos_right_of_diagonal =
1369  std::vector<std::size_t>()) const;
1370 
1374  template <typename somenumber>
1375  void
1376  precondition_SOR(Vector<somenumber> & dst,
1377  const Vector<somenumber> &src,
1378  const number om = 1.) const;
1379 
1383  template <typename somenumber>
1384  void
1385  precondition_TSOR(Vector<somenumber> & dst,
1386  const Vector<somenumber> &src,
1387  const number om = 1.) const;
1388 
1394  template <typename somenumber>
1395  void
1396  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1397 
1402  template <typename somenumber>
1403  void
1404  SOR(Vector<somenumber> &v, const number om = 1.) const;
1405 
1410  template <typename somenumber>
1411  void
1412  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1413 
1424  template <typename somenumber>
1425  void
1426  PSOR(Vector<somenumber> & v,
1427  const std::vector<size_type> &permutation,
1428  const std::vector<size_type> &inverse_permutation,
1429  const number om = 1.) const;
1430 
1441  template <typename somenumber>
1442  void
1443  TPSOR(Vector<somenumber> & v,
1444  const std::vector<size_type> &permutation,
1445  const std::vector<size_type> &inverse_permutation,
1446  const number om = 1.) const;
1447 
1453  template <typename somenumber>
1454  void
1455  Jacobi_step(Vector<somenumber> & v,
1456  const Vector<somenumber> &b,
1457  const number om = 1.) const;
1458 
1463  template <typename somenumber>
1464  void
1465  SOR_step(Vector<somenumber> & v,
1466  const Vector<somenumber> &b,
1467  const number om = 1.) const;
1468 
1473  template <typename somenumber>
1474  void
1475  TSOR_step(Vector<somenumber> & v,
1476  const Vector<somenumber> &b,
1477  const number om = 1.) const;
1478 
1483  template <typename somenumber>
1484  void
1485  SSOR_step(Vector<somenumber> & v,
1486  const Vector<somenumber> &b,
1487  const number om = 1.) const;
1489 
1493 
1501  begin() const;
1502 
1506  iterator
1507  begin();
1508 
1513  end() const;
1514 
1518  iterator
1519  end();
1520 
1531  begin(const size_type r) const;
1532 
1536  iterator
1537  begin(const size_type r);
1538 
1549  end(const size_type r) const;
1550 
1554  iterator
1555  end(const size_type r);
1557 
1561 
1573  template <class StreamType>
1574  void
1575  print(StreamType &out,
1576  const bool across = false,
1577  const bool diagonal_first = true) const;
1578 
1599  void
1600  print_formatted(std::ostream & out,
1601  const unsigned int precision = 3,
1602  const bool scientific = true,
1603  const unsigned int width = 0,
1604  const char * zero_string = " ",
1605  const double denominator = 1.) const;
1606 
1612  void
1613  print_pattern(std::ostream &out, const double threshold = 0.) const;
1614 
1623  void
1624  print_as_numpy_arrays(std::ostream & out,
1625  const unsigned int precision = 9) const;
1626 
1637  void
1638  block_write(std::ostream &out) const;
1639 
1656  void
1657  block_read(std::istream &in);
1659 
1667  DeclException2(ExcInvalidIndex,
1668  int,
1669  int,
1670  << "You are trying to access the matrix entry with index <"
1671  << arg1 << ',' << arg2
1672  << ">, but this entry does not exist in the sparsity pattern "
1673  "of this matrix."
1674  "\n\n"
1675  "The most common cause for this problem is that you used "
1676  "a method to build the sparsity pattern that did not "
1677  "(completely) take into account all of the entries you "
1678  "will later try to write into. An example would be "
1679  "building a sparsity pattern that does not include "
1680  "the entries you will write into due to constraints "
1681  "on degrees of freedom such as hanging nodes or periodic "
1682  "boundary conditions. In such cases, building the "
1683  "sparsity pattern will succeed, but you will get errors "
1684  "such as the current one at one point or other when "
1685  "trying to write into the entries of the matrix.");
1689  DeclExceptionMsg(ExcDifferentSparsityPatterns,
1690  "When copying one sparse matrix into another, "
1691  "or when adding one sparse matrix to another, "
1692  "both matrices need to refer to the same "
1693  "sparsity pattern.");
1697  DeclException2(ExcIteratorRange,
1698  int,
1699  int,
1700  << "The iterators denote a range of " << arg1
1701  << " elements, but the given number of rows was " << arg2);
1705  DeclExceptionMsg(ExcSourceEqualsDestination,
1706  "You are attempting an operation on two matrices that "
1707  "are the same object, but the operation requires that the "
1708  "two objects are in fact different.");
1710 
1711 protected:
1722  void
1723  prepare_add();
1724 
1729  void
1730  prepare_set();
1731 
1732 private:
1739 
1747  std::unique_ptr<number[]> val;
1748 
1755  std::size_t max_len;
1756 
1757  // make all other sparse matrices friends
1758  template <typename somenumber>
1759  friend class SparseMatrix;
1760  template <typename somenumber>
1762  template <typename>
1763  friend class SparseILU;
1764 
1765  // To allow it calling private prepare_add() and prepare_set().
1766  template <typename>
1767  friend class BlockMatrixBase;
1768 
1769  // Also give access to internal details to the iterator/accessor classes.
1770  template <typename, bool>
1772  template <typename, bool>
1774 
1775 # ifdef DEAL_II_WITH_MPI
1776  // Give access to internal datastructures to perform MPI operations.
1777  template <typename Number>
1778  friend void
1780  const MPI_Comm &,
1782 # endif
1783 };
1784 
1785 # ifndef DOXYGEN
1786 /*---------------------- Inline functions -----------------------------------*/
1787 
1788 
1789 
1790 template <typename number>
1791 inline typename SparseMatrix<number>::size_type
1793 {
1794  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1795  return cols->rows;
1796 }
1797 
1798 
1799 template <typename number>
1800 inline typename SparseMatrix<number>::size_type
1802 {
1803  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1804  return cols->cols;
1805 }
1806 
1807 
1808 // Inline the set() and add() functions, since they will be called frequently.
1809 template <typename number>
1810 inline void
1812  const size_type j,
1813  const number value)
1814 {
1815  AssertIsFinite(value);
1816 
1817  const size_type index = cols->operator()(i, j);
1818 
1819  // it is allowed to set elements of the matrix that are not part of the
1820  // sparsity pattern, if the value to which we set it is zero
1821  if (index == SparsityPattern::invalid_entry)
1822  {
1823  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1824  ExcInvalidIndex(i, j));
1825  return;
1826  }
1827 
1828  val[index] = value;
1829 }
1830 
1831 
1832 
1833 template <typename number>
1834 template <typename number2>
1835 inline void
1836 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1837  const FullMatrix<number2> & values,
1838  const bool elide_zero_values)
1839 {
1840  Assert(indices.size() == values.m(),
1841  ExcDimensionMismatch(indices.size(), values.m()));
1842  Assert(values.m() == values.n(), ExcNotQuadratic());
1843 
1844  for (size_type i = 0; i < indices.size(); ++i)
1845  set(indices[i],
1846  indices.size(),
1847  indices.data(),
1848  &values(i, 0),
1849  elide_zero_values);
1850 }
1851 
1852 
1853 
1854 template <typename number>
1855 template <typename number2>
1856 inline void
1857 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1858  const std::vector<size_type> &col_indices,
1859  const FullMatrix<number2> & values,
1860  const bool elide_zero_values)
1861 {
1862  Assert(row_indices.size() == values.m(),
1863  ExcDimensionMismatch(row_indices.size(), values.m()));
1864  Assert(col_indices.size() == values.n(),
1865  ExcDimensionMismatch(col_indices.size(), values.n()));
1866 
1867  for (size_type i = 0; i < row_indices.size(); ++i)
1868  set(row_indices[i],
1869  col_indices.size(),
1870  col_indices.data(),
1871  &values(i, 0),
1872  elide_zero_values);
1873 }
1874 
1875 
1876 
1877 template <typename number>
1878 template <typename number2>
1879 inline void
1881  const std::vector<size_type> &col_indices,
1882  const std::vector<number2> & values,
1883  const bool elide_zero_values)
1884 {
1885  Assert(col_indices.size() == values.size(),
1886  ExcDimensionMismatch(col_indices.size(), values.size()));
1887 
1888  set(row,
1889  col_indices.size(),
1890  col_indices.data(),
1891  values.data(),
1892  elide_zero_values);
1893 }
1894 
1895 
1896 
1897 template <typename number>
1898 inline void
1900  const size_type j,
1901  const number value)
1902 {
1903  AssertIsFinite(value);
1904 
1905  if (value == number())
1906  return;
1907 
1908  const size_type index = cols->operator()(i, j);
1909 
1910  // it is allowed to add elements to the matrix that are not part of the
1911  // sparsity pattern, if the value to which we set it is zero
1912  if (index == SparsityPattern::invalid_entry)
1913  {
1914  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1915  ExcInvalidIndex(i, j));
1916  return;
1917  }
1918 
1919  val[index] += value;
1920 }
1921 
1922 
1923 
1924 template <typename number>
1925 template <typename number2>
1926 inline void
1927 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1928  const FullMatrix<number2> & values,
1929  const bool elide_zero_values)
1930 {
1931  Assert(indices.size() == values.m(),
1932  ExcDimensionMismatch(indices.size(), values.m()));
1933  Assert(values.m() == values.n(), ExcNotQuadratic());
1934 
1935  for (size_type i = 0; i < indices.size(); ++i)
1936  add(indices[i],
1937  indices.size(),
1938  indices.data(),
1939  &values(i, 0),
1940  elide_zero_values);
1941 }
1942 
1943 
1944 
1945 template <typename number>
1946 template <typename number2>
1947 inline void
1948 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1949  const std::vector<size_type> &col_indices,
1950  const FullMatrix<number2> & values,
1951  const bool elide_zero_values)
1952 {
1953  Assert(row_indices.size() == values.m(),
1954  ExcDimensionMismatch(row_indices.size(), values.m()));
1955  Assert(col_indices.size() == values.n(),
1956  ExcDimensionMismatch(col_indices.size(), values.n()));
1957 
1958  for (size_type i = 0; i < row_indices.size(); ++i)
1959  add(row_indices[i],
1960  col_indices.size(),
1961  col_indices.data(),
1962  &values(i, 0),
1963  elide_zero_values);
1964 }
1965 
1966 
1967 
1968 template <typename number>
1969 template <typename number2>
1970 inline void
1972  const std::vector<size_type> &col_indices,
1973  const std::vector<number2> & values,
1974  const bool elide_zero_values)
1975 {
1976  Assert(col_indices.size() == values.size(),
1977  ExcDimensionMismatch(col_indices.size(), values.size()));
1978 
1979  add(row,
1980  col_indices.size(),
1981  col_indices.data(),
1982  values.data(),
1983  elide_zero_values);
1984 }
1985 
1986 
1987 
1988 template <typename number>
1989 inline SparseMatrix<number> &
1990 SparseMatrix<number>::operator*=(const number factor)
1991 {
1992  Assert(cols != nullptr, ExcNeedsSparsityPattern());
1993  Assert(val != nullptr, ExcNotInitialized());
1994 
1995  number * val_ptr = val.get();
1996  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1997 
1998  while (val_ptr != end_ptr)
1999  *val_ptr++ *= factor;
2000 
2001  return *this;
2002 }
2003 
2004 
2005 
2006 template <typename number>
2007 inline SparseMatrix<number> &
2008 SparseMatrix<number>::operator/=(const number factor)
2009 {
2010  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2011  Assert(val != nullptr, ExcNotInitialized());
2012  Assert(factor != number(), ExcDivideByZero());
2013 
2014  const number factor_inv = number(1.) / factor;
2015 
2016  number * val_ptr = val.get();
2017  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2018 
2019  while (val_ptr != end_ptr)
2020  *val_ptr++ *= factor_inv;
2021 
2022  return *this;
2023 }
2024 
2025 
2026 
2027 template <typename number>
2028 inline const number &
2030 {
2031  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2032  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2033  ExcInvalidIndex(i, j));
2034  return val[cols->operator()(i, j)];
2035 }
2036 
2037 
2038 
2039 template <typename number>
2040 inline number &
2042 {
2043  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2044  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2045  ExcInvalidIndex(i, j));
2046  return val[cols->operator()(i, j)];
2047 }
2048 
2049 
2050 
2051 template <typename number>
2052 inline number
2053 SparseMatrix<number>::el(const size_type i, const size_type j) const
2054 {
2055  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2056  const size_type index = cols->operator()(i, j);
2057 
2058  if (index != SparsityPattern::invalid_entry)
2059  return val[index];
2060  else
2061  return 0;
2062 }
2063 
2064 
2065 
2066 template <typename number>
2067 inline number
2069 {
2070  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2071  Assert(m() == n(), ExcNotQuadratic());
2072  AssertIndexRange(i, m());
2073 
2074  // Use that the first element in each row of a quadratic matrix is the main
2075  // diagonal
2076  return val[cols->rowstart[i]];
2077 }
2078 
2079 
2080 
2081 template <typename number>
2082 inline number &
2084 {
2085  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2086  Assert(m() == n(), ExcNotQuadratic());
2087  AssertIndexRange(i, m());
2088 
2089  // Use that the first element in each row of a quadratic matrix is the main
2090  // diagonal
2091  return val[cols->rowstart[i]];
2092 }
2093 
2094 
2095 
2096 template <typename number>
2097 template <typename ForwardIterator>
2098 void
2099 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2100  const ForwardIterator end)
2101 {
2102  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2103  ExcIteratorRange(std::distance(begin, end), m()));
2104 
2105  // for use in the inner loop, we define an alias to the type of the inner
2106  // iterators
2107  using inner_iterator =
2108  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2109  size_type row = 0;
2110  for (ForwardIterator i = begin; i != end; ++i, ++row)
2111  {
2112  const inner_iterator end_of_row = i->end();
2113  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2114  // write entries
2115  set(row, j->first, j->second);
2116  };
2117 }
2118 
2119 
2120 //---------------------------------------------------------------------------
2121 
2122 
2123 namespace SparseMatrixIterators
2124 {
2125  template <typename number>
2126  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2127  const std::size_t index_within_matrix)
2128  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2129  index_within_matrix)
2130  , matrix(matrix)
2131  {}
2132 
2133 
2134 
2135  template <typename number>
2136  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2137  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2138  , matrix(matrix)
2139  {}
2140 
2141 
2142 
2143  template <typename number>
2144  inline Accessor<number, true>::Accessor(
2146  : SparsityPatternIterators::Accessor(a)
2147  , matrix(&a.get_matrix())
2148  {}
2149 
2150 
2151 
2152  template <typename number>
2153  inline number
2154  Accessor<number, true>::value() const
2155  {
2156  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2157  return matrix->val[linear_index];
2158  }
2159 
2160 
2161 
2162  template <typename number>
2163  inline const typename Accessor<number, true>::MatrixType &
2164  Accessor<number, true>::get_matrix() const
2165  {
2166  return *matrix;
2167  }
2168 
2169 
2170 
2171  template <typename number>
2172  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2173  const bool)
2174  : accessor(accessor)
2175  {}
2176 
2177 
2178  template <typename number>
2179  inline Accessor<number, false>::Reference::operator number() const
2180  {
2181  AssertIndexRange(accessor->linear_index,
2182  accessor->matrix->n_nonzero_elements());
2183  return accessor->matrix->val[accessor->linear_index];
2184  }
2185 
2186 
2187 
2188  template <typename number>
2189  inline const typename Accessor<number, false>::Reference &
2190  Accessor<number, false>::Reference::operator=(const number n) const
2191  {
2192  AssertIndexRange(accessor->linear_index,
2193  accessor->matrix->n_nonzero_elements());
2194  accessor->matrix->val[accessor->linear_index] = n;
2195  return *this;
2196  }
2197 
2198 
2199 
2200  template <typename number>
2201  inline const typename Accessor<number, false>::Reference &
2202  Accessor<number, false>::Reference::operator+=(const number n) const
2203  {
2204  AssertIndexRange(accessor->linear_index,
2205  accessor->matrix->n_nonzero_elements());
2206  accessor->matrix->val[accessor->linear_index] += n;
2207  return *this;
2208  }
2209 
2210 
2211 
2212  template <typename number>
2213  inline const typename Accessor<number, false>::Reference &
2214  Accessor<number, false>::Reference::operator-=(const number n) const
2215  {
2216  AssertIndexRange(accessor->linear_index,
2217  accessor->matrix->n_nonzero_elements());
2218  accessor->matrix->val[accessor->linear_index] -= n;
2219  return *this;
2220  }
2221 
2222 
2223 
2224  template <typename number>
2225  inline const typename Accessor<number, false>::Reference &
2226  Accessor<number, false>::Reference::operator*=(const number n) const
2227  {
2228  AssertIndexRange(accessor->linear_index,
2229  accessor->matrix->n_nonzero_elements());
2230  accessor->matrix->val[accessor->linear_index] *= n;
2231  return *this;
2232  }
2233 
2234 
2235 
2236  template <typename number>
2237  inline const typename Accessor<number, false>::Reference &
2238  Accessor<number, false>::Reference::operator/=(const number n) const
2239  {
2240  AssertIndexRange(accessor->linear_index,
2241  accessor->matrix->n_nonzero_elements());
2242  accessor->matrix->val[accessor->linear_index] /= n;
2243  return *this;
2244  }
2245 
2246 
2247 
2248  template <typename number>
2249  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2250  const std::size_t index)
2251  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2252  , matrix(matrix)
2253  {}
2254 
2255 
2256 
2257  template <typename number>
2258  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2259  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2260  , matrix(matrix)
2261  {}
2262 
2263 
2264 
2265  template <typename number>
2266  inline typename Accessor<number, false>::Reference
2267  Accessor<number, false>::value() const
2268  {
2269  return Reference(this, true);
2270  }
2271 
2272 
2273 
2274  template <typename number>
2275  inline typename Accessor<number, false>::MatrixType &
2276  Accessor<number, false>::get_matrix() const
2277  {
2278  return *matrix;
2279  }
2280 
2281 
2282 
2283  template <typename number, bool Constness>
2284  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2285  const std::size_t index)
2286  : accessor(matrix, index)
2287  {}
2288 
2289 
2290 
2291  template <typename number, bool Constness>
2292  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2293  : accessor(matrix)
2294  {}
2295 
2296 
2297 
2298  template <typename number, bool Constness>
2299  inline Iterator<number, Constness>::Iterator(
2301  : accessor(*i)
2302  {}
2303 
2304 
2305 
2306  template <typename number, bool Constness>
2307  inline const Iterator<number, Constness> &
2308  Iterator<number, Constness>::operator=(
2310  {
2311  accessor = *i;
2312  return *this;
2313  }
2314 
2315 
2316 
2317  template <typename number, bool Constness>
2318  inline Iterator<number, Constness> &
2320  {
2321  accessor.advance();
2322  return *this;
2323  }
2324 
2325 
2326  template <typename number, bool Constness>
2327  inline Iterator<number, Constness>
2329  {
2330  const Iterator iter = *this;
2331  accessor.advance();
2332  return iter;
2333  }
2334 
2335 
2336  template <typename number, bool Constness>
2337  inline const Accessor<number, Constness> &
2339  {
2340  return accessor;
2341  }
2342 
2343 
2344  template <typename number, bool Constness>
2345  inline const Accessor<number, Constness> *
2346  Iterator<number, Constness>::operator->() const
2347  {
2348  return &accessor;
2349  }
2350 
2351 
2352  template <typename number, bool Constness>
2353  inline bool
2354  Iterator<number, Constness>::operator==(const Iterator &other) const
2355  {
2356  return (accessor == other.accessor);
2357  }
2358 
2359 
2360  template <typename number, bool Constness>
2361  inline bool
2362  Iterator<number, Constness>::operator!=(const Iterator &other) const
2363  {
2364  return !(*this == other);
2365  }
2366 
2367 
2368  template <typename number, bool Constness>
2369  inline bool
2370  Iterator<number, Constness>::operator<(const Iterator &other) const
2371  {
2372  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2373  ExcInternalError());
2374 
2375  return (accessor < other.accessor);
2376  }
2377 
2378 
2379  template <typename number, bool Constness>
2380  inline bool
2381  Iterator<number, Constness>::operator>(const Iterator &other) const
2382  {
2383  return (other < *this);
2384  }
2385 
2386 
2387  template <typename number, bool Constness>
2388  inline int
2389  Iterator<number, Constness>::operator-(const Iterator &other) const
2390  {
2391  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2392  ExcInternalError());
2393 
2394  return (*this)->linear_index - other->linear_index;
2395  }
2396 
2397 
2398 
2399  template <typename number, bool Constness>
2400  inline Iterator<number, Constness>
2402  {
2403  Iterator x = *this;
2404  for (size_type i = 0; i < n; ++i)
2405  ++x;
2406 
2407  return x;
2408  }
2409 
2410 } // namespace SparseMatrixIterators
2411 
2412 
2413 
2414 template <typename number>
2417 {
2418  return const_iterator(this, 0);
2419 }
2420 
2421 
2422 template <typename number>
2425 {
2426  return const_iterator(this);
2427 }
2428 
2429 
2430 template <typename number>
2431 inline typename SparseMatrix<number>::iterator
2433 {
2434  return iterator(this, 0);
2435 }
2436 
2437 
2438 template <typename number>
2439 inline typename SparseMatrix<number>::iterator
2441 {
2442  return iterator(this, cols->rowstart[cols->rows]);
2443 }
2444 
2445 
2446 template <typename number>
2449 {
2450  AssertIndexRange(r, m());
2451 
2452  return const_iterator(this, cols->rowstart[r]);
2453 }
2454 
2455 
2456 
2457 template <typename number>
2459 SparseMatrix<number>::end(const size_type r) const
2460 {
2461  AssertIndexRange(r, m());
2462 
2463  return const_iterator(this, cols->rowstart[r + 1]);
2464 }
2465 
2466 
2467 
2468 template <typename number>
2469 inline typename SparseMatrix<number>::iterator
2471 {
2472  AssertIndexRange(r, m());
2473 
2474  return iterator(this, cols->rowstart[r]);
2475 }
2476 
2477 
2478 
2479 template <typename number>
2480 inline typename SparseMatrix<number>::iterator
2482 {
2483  AssertIndexRange(r, m());
2484 
2485  return iterator(this, cols->rowstart[r + 1]);
2486 }
2487 
2488 
2489 
2490 template <typename number>
2491 template <class StreamType>
2492 inline void
2493 SparseMatrix<number>::print(StreamType &out,
2494  const bool across,
2495  const bool diagonal_first) const
2496 {
2497  Assert(cols != nullptr, ExcNeedsSparsityPattern());
2498  Assert(val != nullptr, ExcNotInitialized());
2499 
2500  bool hanging_diagonal = false;
2501  number diagonal = number();
2502 
2503  for (size_type i = 0; i < cols->rows; ++i)
2504  {
2505  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2506  {
2507  if (!diagonal_first && i == cols->colnums[j])
2508  {
2509  diagonal = val[j];
2510  hanging_diagonal = true;
2511  }
2512  else
2513  {
2514  if (hanging_diagonal && cols->colnums[j] > i)
2515  {
2516  if (across)
2517  out << ' ' << i << ',' << i << ':' << diagonal;
2518  else
2519  out << '(' << i << ',' << i << ") " << diagonal
2520  << std::endl;
2521  hanging_diagonal = false;
2522  }
2523  if (across)
2524  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2525  else
2526  out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2527  << std::endl;
2528  }
2529  }
2530  if (hanging_diagonal)
2531  {
2532  if (across)
2533  out << ' ' << i << ',' << i << ':' << diagonal;
2534  else
2535  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2536  hanging_diagonal = false;
2537  }
2538  }
2539  if (across)
2540  out << std::endl;
2541 }
2542 
2543 
2544 template <typename number>
2545 inline void
2547 {
2548  // nothing to do here
2549 }
2550 
2551 
2552 
2553 template <typename number>
2554 inline void
2556 {
2557  // nothing to do here
2558 }
2559 
2560 # endif // DOXYGEN
2561 
2562 
2563 /*---------------------------- sparse_matrix.h ---------------------------*/
2564 
2566 
2567 #endif
2568 /*---------------------------- 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:532
types::global_dof_index size_type
Definition: sparse_matrix.h:83
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:1720
number value_type
static const char V
STL namespace.
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:1461
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:487
SparseMatrixIterators::Iterator< number, false > iterator
Number linfty_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:3025
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
Accessor< number, Constness > accessor
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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
typename ::SparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
Number l1_norm(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:2999
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)
typename ::SparseMatrixIterators::Iterator< number, Constness >::value_type value_type
const number & operator()(const size_type i, const size_type j) const
static ::ExceptionBase & ExcNeedsSparsityPattern()
SparseMatrixIterators::Iterator< number, true > const_iterator
#define AssertIsFinite(number)
Definition: exceptions.h:1746
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()