Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 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 
22 # include <deal.II/base/smartpointer.h>
23 # include <deal.II/base/subscriptor.h>
24 
25 # include <deal.II/lac/exceptions.h>
26 # include <deal.II/lac/identity_matrix.h>
27 # include <deal.II/lac/sparsity_pattern.h>
28 # include <deal.II/lac/vector_operation.h>
29 # ifdef DEAL_II_WITH_MPI
30 # include <mpi.h>
31 # endif
32 
33 # include <memory>
34 
35 
36 DEAL_II_NAMESPACE_OPEN
37 
38 template <typename number>
39 class Vector;
40 template <typename number>
41 class FullMatrix;
42 template <typename Matrix>
43 class BlockMatrixBase;
44 template <typename number>
45 class SparseILU;
46 # ifdef DEAL_II_WITH_MPI
47 namespace Utilities
48 {
49  namespace MPI
50  {
51  template <typename Number>
52  void
53  sum(const SparseMatrix<Number> &, const MPI_Comm &, SparseMatrix<Number> &);
54  }
55 } // namespace Utilities
56 # endif
57 
58 # ifdef DEAL_II_WITH_TRILINOS
59 namespace TrilinosWrappers
60 {
61  class SparseMatrix;
62 }
63 # endif
64 
75 {
80 
81  // forward declaration
82  template <typename number, bool Constness>
83  class Iterator;
84 
95  template <typename number, bool Constness>
97  {
98  public:
102  number
103  value() const;
104 
108  number &
109  value();
110 
115  const SparseMatrix<number> &
116  get_matrix() const;
117  };
118 
119 
120 
127  template <typename number>
128  class Accessor<number, true> : public SparsityPatternIterators::Accessor
129  {
130  public:
136 
140  Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
141 
145  Accessor(MatrixType *matrix);
146 
151 
155  number
156  value() const;
157 
162  const MatrixType &
163  get_matrix() const;
164 
165  private:
170 
175 
179  template <typename, bool>
180  friend class Iterator;
181  };
182 
183 
190  template <typename number>
191  class Accessor<number, false> : public SparsityPatternIterators::Accessor
192  {
193  private:
218  class Reference
219  {
220  public:
225  Reference(const Accessor *accessor, const bool dummy);
226 
230  operator number() const;
231 
235  const Reference &
236  operator=(const number n) const;
237 
241  const Reference &
242  operator+=(const number n) const;
243 
247  const Reference &
248  operator-=(const number n) const;
249 
253  const Reference &
254  operator*=(const number n) const;
255 
259  const Reference &
260  operator/=(const number n) const;
261 
262  private:
268  };
269 
270  public:
276 
280  Accessor(MatrixType *matrix, const std::size_t index);
281 
285  Accessor(MatrixType *matrix);
286 
290  Reference
291  value() const;
292 
297  MatrixType &
298  get_matrix() const;
299 
300  private:
305 
310 
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 
368  Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
369 
373  Iterator(MatrixType *matrix);
374 
380 
386 
390  Iterator &
391  operator++();
392 
396  Iterator
397  operator++(int);
398 
402  const Accessor<number, Constness> &operator*() const;
403 
407  const Accessor<number, Constness> *operator->() const;
408 
412  bool
413  operator==(const Iterator &) const;
414 
418  bool
419  operator!=(const Iterator &) const;
420 
428  bool
429  operator<(const Iterator &) const;
430 
435  bool
436  operator>(const Iterator &) const;
437 
444  int
445  operator-(const Iterator &p) const;
446 
450  Iterator
451  operator+(const size_type n) const;
452 
453  private:
458  };
459 
460 } // namespace SparseMatrixIterators
461 
467 // TODO: Add multithreading to the other vmult functions.
468 
497 template <typename number>
498 class SparseMatrix : public virtual Subscriptor
499 {
500 public:
505 
510  using value_type = number;
511 
522 
528 
536 
543  struct Traits
544  {
549  static const bool zero_addition_can_be_elided = true;
550  };
551 
566  SparseMatrix();
567 
576  SparseMatrix(const SparseMatrix &);
577 
585  SparseMatrix(SparseMatrix<number> &&m) noexcept;
586 
600  explicit SparseMatrix(const SparsityPattern &sparsity);
601 
608  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
609 
614  virtual ~SparseMatrix() override;
615 
626  operator=(const SparseMatrix<number> &);
627 
633  operator=(SparseMatrix<number> &&m) noexcept;
634 
642  operator=(const IdentityMatrix &id);
643 
655  SparseMatrix &
656  operator=(const double d);
657 
671  virtual void
672  reinit(const SparsityPattern &sparsity);
673 
679  virtual void
680  clear();
682 
690  bool
691  empty() const;
692 
697  size_type
698  m() const;
699 
704  size_type
705  n() const;
706 
710  size_type
711  get_row_length(const size_type row) const;
712 
718  std::size_t
719  n_nonzero_elements() const;
720 
730  std::size_t
731  n_actually_nonzero_elements(const double threshold = 0.) const;
732 
741  const SparsityPattern &
742  get_sparsity_pattern() const;
743 
748  std::size_t
749  memory_consumption() const;
750 
754  void compress(::VectorOperation::values);
755 
757 
766  void
767  set(const size_type i, const size_type j, const number value);
768 
784  template <typename number2>
785  void
786  set(const std::vector<size_type> &indices,
787  const FullMatrix<number2> & full_matrix,
788  const bool elide_zero_values = false);
789 
795  template <typename number2>
796  void
797  set(const std::vector<size_type> &row_indices,
798  const std::vector<size_type> &col_indices,
799  const FullMatrix<number2> & full_matrix,
800  const bool elide_zero_values = false);
801 
812  template <typename number2>
813  void
814  set(const size_type row,
815  const std::vector<size_type> &col_indices,
816  const std::vector<number2> & values,
817  const bool elide_zero_values = false);
818 
828  template <typename number2>
829  void
830  set(const size_type row,
831  const size_type n_cols,
832  const size_type *col_indices,
833  const number2 * values,
834  const bool elide_zero_values = false);
835 
841  void
842  add(const size_type i, const size_type j, const number value);
843 
858  template <typename number2>
859  void
860  add(const std::vector<size_type> &indices,
861  const FullMatrix<number2> & full_matrix,
862  const bool elide_zero_values = true);
863 
869  template <typename number2>
870  void
871  add(const std::vector<size_type> &row_indices,
872  const std::vector<size_type> &col_indices,
873  const FullMatrix<number2> & full_matrix,
874  const bool elide_zero_values = true);
875 
885  template <typename number2>
886  void
887  add(const size_type row,
888  const std::vector<size_type> &col_indices,
889  const std::vector<number2> & values,
890  const bool elide_zero_values = true);
891 
901  template <typename number2>
902  void
903  add(const size_type row,
904  const size_type n_cols,
905  const size_type *col_indices,
906  const number2 * values,
907  const bool elide_zero_values = true,
908  const bool col_indices_are_sorted = false);
909 
913  SparseMatrix &
914  operator*=(const number factor);
915 
919  SparseMatrix &
920  operator/=(const number factor);
921 
934  void
935  symmetrize();
936 
953  template <typename somenumber>
955  copy_from(const SparseMatrix<somenumber> &source);
956 
973  template <typename ForwardIterator>
974  void
975  copy_from(const ForwardIterator begin, const ForwardIterator end);
976 
982  template <typename somenumber>
983  void
984  copy_from(const FullMatrix<somenumber> &matrix);
985 
986 # ifdef DEAL_II_WITH_TRILINOS
987 
997  copy_from(const TrilinosWrappers::SparseMatrix &matrix);
998 # endif
999 
1011  template <typename somenumber>
1012  void
1013  add(const number factor, const SparseMatrix<somenumber> &matrix);
1014 
1016 
1020 
1034  const number &
1035  operator()(const size_type i, const size_type j) const;
1036 
1040  number &
1041  operator()(const size_type i, const size_type j);
1042 
1055  number
1056  el(const size_type i, const size_type j) const;
1057 
1067  number
1068  diag_element(const size_type i) const;
1069 
1074  number &
1075  diag_element(const size_type i);
1076 
1078 
1098  template <class OutVector, class InVector>
1099  void
1100  vmult(OutVector &dst, const InVector &src) const;
1101 
1117  template <class OutVector, class InVector>
1118  void
1119  Tvmult(OutVector &dst, const InVector &src) const;
1120 
1137  template <class OutVector, class InVector>
1138  void
1139  vmult_add(OutVector &dst, const InVector &src) const;
1140 
1156  template <class OutVector, class InVector>
1157  void
1158  Tvmult_add(OutVector &dst, const InVector &src) const;
1159 
1177  template <typename somenumber>
1178  somenumber
1179  matrix_norm_square(const Vector<somenumber> &v) const;
1180 
1186  template <typename somenumber>
1187  somenumber
1188  matrix_scalar_product(const Vector<somenumber> &u,
1189  const Vector<somenumber> &v) const;
1190 
1200  template <typename somenumber>
1201  somenumber
1202  residual(Vector<somenumber> & dst,
1203  const Vector<somenumber> &x,
1204  const Vector<somenumber> &b) const;
1205 
1241  template <typename numberB, typename numberC>
1242  void
1243  mmult(SparseMatrix<numberC> & C,
1244  const SparseMatrix<numberB> &B,
1245  const Vector<number> & V = Vector<number>(),
1246  const bool rebuild_sparsity_pattern = true) const;
1247 
1272  template <typename numberB, typename numberC>
1273  void
1274  Tmmult(SparseMatrix<numberC> & C,
1275  const SparseMatrix<numberB> &B,
1276  const Vector<number> & V = Vector<number>(),
1277  const bool rebuild_sparsity_pattern = true) const;
1278 
1280 
1284 
1292  real_type
1293  l1_norm() const;
1294 
1302  real_type
1303  linfty_norm() const;
1304 
1309  real_type
1310  frobenius_norm() const;
1312 
1316 
1322  template <typename somenumber>
1323  void
1324  precondition_Jacobi(Vector<somenumber> & dst,
1325  const Vector<somenumber> &src,
1326  const number omega = 1.) const;
1327 
1334  template <typename somenumber>
1335  void
1336  precondition_SSOR(Vector<somenumber> & dst,
1337  const Vector<somenumber> & src,
1338  const number omega = 1.,
1339  const std::vector<std::size_t> &pos_right_of_diagonal =
1340  std::vector<std::size_t>()) const;
1341 
1345  template <typename somenumber>
1346  void
1347  precondition_SOR(Vector<somenumber> & dst,
1348  const Vector<somenumber> &src,
1349  const number om = 1.) const;
1350 
1354  template <typename somenumber>
1355  void
1356  precondition_TSOR(Vector<somenumber> & dst,
1357  const Vector<somenumber> &src,
1358  const number om = 1.) const;
1359 
1365  template <typename somenumber>
1366  void
1367  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1368 
1373  template <typename somenumber>
1374  void
1375  SOR(Vector<somenumber> &v, const number om = 1.) const;
1376 
1381  template <typename somenumber>
1382  void
1383  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1384 
1395  template <typename somenumber>
1396  void
1397  PSOR(Vector<somenumber> & v,
1398  const std::vector<size_type> &permutation,
1399  const std::vector<size_type> &inverse_permutation,
1400  const number om = 1.) const;
1401 
1412  template <typename somenumber>
1413  void
1414  TPSOR(Vector<somenumber> & v,
1415  const std::vector<size_type> &permutation,
1416  const std::vector<size_type> &inverse_permutation,
1417  const number om = 1.) const;
1418 
1424  template <typename somenumber>
1425  void
1426  Jacobi_step(Vector<somenumber> & v,
1427  const Vector<somenumber> &b,
1428  const number om = 1.) const;
1429 
1434  template <typename somenumber>
1435  void
1436  SOR_step(Vector<somenumber> & v,
1437  const Vector<somenumber> &b,
1438  const number om = 1.) const;
1439 
1444  template <typename somenumber>
1445  void
1446  TSOR_step(Vector<somenumber> & v,
1447  const Vector<somenumber> &b,
1448  const number om = 1.) const;
1449 
1454  template <typename somenumber>
1455  void
1456  SSOR_step(Vector<somenumber> & v,
1457  const Vector<somenumber> &b,
1458  const number om = 1.) const;
1460 
1464 
1472  begin() const;
1473 
1477  iterator
1478  begin();
1479 
1484  end() const;
1485 
1489  iterator
1490  end();
1491 
1502  begin(const size_type r) const;
1503 
1507  iterator
1508  begin(const size_type r);
1509 
1520  end(const size_type r) const;
1521 
1525  iterator
1526  end(const size_type r);
1528 
1532 
1544  template <class StreamType>
1545  void
1546  print(StreamType &out,
1547  const bool across = false,
1548  const bool diagonal_first = true) const;
1549 
1570  void
1571  print_formatted(std::ostream & out,
1572  const unsigned int precision = 3,
1573  const bool scientific = true,
1574  const unsigned int width = 0,
1575  const char * zero_string = " ",
1576  const double denominator = 1.) const;
1577 
1583  void
1584  print_pattern(std::ostream &out, const double threshold = 0.) const;
1585 
1596  void
1597  block_write(std::ostream &out) const;
1598 
1615  void
1616  block_read(std::istream &in);
1618 
1626  DeclException2(ExcInvalidIndex,
1627  int,
1628  int,
1629  << "You are trying to access the matrix entry with index <"
1630  << arg1 << ',' << arg2
1631  << ">, but this entry does not exist in the sparsity pattern "
1632  "of this matrix."
1633  "\n\n"
1634  "The most common cause for this problem is that you used "
1635  "a method to build the sparsity pattern that did not "
1636  "(completely) take into account all of the entries you "
1637  "will later try to write into. An example would be "
1638  "building a sparsity pattern that does not include "
1639  "the entries you will write into due to constraints "
1640  "on degrees of freedom such as hanging nodes or periodic "
1641  "boundary conditions. In such cases, building the "
1642  "sparsity pattern will succeed, but you will get errors "
1643  "such as the current one at one point or other when "
1644  "trying to write into the entries of the matrix.");
1648  DeclExceptionMsg(ExcDifferentSparsityPatterns,
1649  "When copying one sparse matrix into another, "
1650  "or when adding one sparse matrix to another, "
1651  "both matrices need to refer to the same "
1652  "sparsity pattern.");
1656  DeclException2(ExcIteratorRange,
1657  int,
1658  int,
1659  << "The iterators denote a range of " << arg1
1660  << " elements, but the given number of rows was " << arg2);
1664  DeclExceptionMsg(ExcSourceEqualsDestination,
1665  "You are attempting an operation on two matrices that "
1666  "are the same object, but the operation requires that the "
1667  "two objects are in fact different.");
1669 
1670 protected:
1681  void
1682  prepare_add();
1683 
1688  void
1689  prepare_set();
1690 
1691 private:
1698 
1706  std::unique_ptr<number[]> val;
1707 
1714  std::size_t max_len;
1715 
1716  // make all other sparse matrices friends
1717  template <typename somenumber>
1718  friend class SparseMatrix;
1719  template <typename somenumber>
1720  friend class SparseLUDecomposition;
1721  template <typename>
1722  friend class SparseILU;
1723 
1727  template <typename>
1728  friend class BlockMatrixBase;
1729 
1733  template <typename, bool>
1735  template <typename, bool>
1736  friend class SparseMatrixIterators::Accessor;
1737 
1738 # ifdef DEAL_II_WITH_MPI
1739 
1742  template <typename Number>
1743  friend void
1745  const MPI_Comm &,
1747 # endif
1748 };
1749 
1750 # ifndef DOXYGEN
1751 /*---------------------- Inline functions -----------------------------------*/
1752 
1753 
1754 
1755 template <typename number>
1756 inline typename SparseMatrix<number>::size_type
1758 {
1759  Assert(cols != nullptr, ExcNotInitialized());
1760  return cols->rows;
1761 }
1762 
1763 
1764 template <typename number>
1765 inline typename SparseMatrix<number>::size_type
1767 {
1768  Assert(cols != nullptr, ExcNotInitialized());
1769  return cols->cols;
1770 }
1771 
1772 
1773 // Inline the set() and add() functions, since they will be called frequently.
1774 template <typename number>
1775 inline void
1777  const size_type j,
1778  const number value)
1779 {
1780  AssertIsFinite(value);
1781 
1782  const size_type index = cols->operator()(i, j);
1783 
1784  // it is allowed to set elements of the matrix that are not part of the
1785  // sparsity pattern, if the value to which we set it is zero
1786  if (index == SparsityPattern::invalid_entry)
1787  {
1788  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1789  ExcInvalidIndex(i, j));
1790  return;
1791  }
1792 
1793  val[index] = value;
1794 }
1795 
1796 
1797 
1798 template <typename number>
1799 template <typename number2>
1800 inline void
1801 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1802  const FullMatrix<number2> & values,
1803  const bool elide_zero_values)
1804 {
1805  Assert(indices.size() == values.m(),
1806  ExcDimensionMismatch(indices.size(), values.m()));
1807  Assert(values.m() == values.n(), ExcNotQuadratic());
1808 
1809  for (size_type i = 0; i < indices.size(); ++i)
1810  set(indices[i],
1811  indices.size(),
1812  indices.data(),
1813  &values(i, 0),
1814  elide_zero_values);
1815 }
1816 
1817 
1818 
1819 template <typename number>
1820 template <typename number2>
1821 inline void
1822 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1823  const std::vector<size_type> &col_indices,
1824  const FullMatrix<number2> & values,
1825  const bool elide_zero_values)
1826 {
1827  Assert(row_indices.size() == values.m(),
1828  ExcDimensionMismatch(row_indices.size(), values.m()));
1829  Assert(col_indices.size() == values.n(),
1830  ExcDimensionMismatch(col_indices.size(), values.n()));
1831 
1832  for (size_type i = 0; i < row_indices.size(); ++i)
1833  set(row_indices[i],
1834  col_indices.size(),
1835  col_indices.data(),
1836  &values(i, 0),
1837  elide_zero_values);
1838 }
1839 
1840 
1841 
1842 template <typename number>
1843 template <typename number2>
1844 inline void
1846  const std::vector<size_type> &col_indices,
1847  const std::vector<number2> & values,
1848  const bool elide_zero_values)
1849 {
1850  Assert(col_indices.size() == values.size(),
1851  ExcDimensionMismatch(col_indices.size(), values.size()));
1852 
1853  set(row,
1854  col_indices.size(),
1855  col_indices.data(),
1856  values.data(),
1857  elide_zero_values);
1858 }
1859 
1860 
1861 
1862 template <typename number>
1863 inline void
1865  const size_type j,
1866  const number value)
1867 {
1868  AssertIsFinite(value);
1869 
1870  if (value == number())
1871  return;
1872 
1873  const size_type index = cols->operator()(i, j);
1874 
1875  // it is allowed to add elements to the matrix that are not part of the
1876  // sparsity pattern, if the value to which we set it is zero
1877  if (index == SparsityPattern::invalid_entry)
1878  {
1879  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1880  ExcInvalidIndex(i, j));
1881  return;
1882  }
1883 
1884  val[index] += value;
1885 }
1886 
1887 
1888 
1889 template <typename number>
1890 template <typename number2>
1891 inline void
1892 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1893  const FullMatrix<number2> & values,
1894  const bool elide_zero_values)
1895 {
1896  Assert(indices.size() == values.m(),
1897  ExcDimensionMismatch(indices.size(), values.m()));
1898  Assert(values.m() == values.n(), ExcNotQuadratic());
1899 
1900  for (size_type i = 0; i < indices.size(); ++i)
1901  add(indices[i],
1902  indices.size(),
1903  indices.data(),
1904  &values(i, 0),
1905  elide_zero_values);
1906 }
1907 
1908 
1909 
1910 template <typename number>
1911 template <typename number2>
1912 inline void
1913 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1914  const std::vector<size_type> &col_indices,
1915  const FullMatrix<number2> & values,
1916  const bool elide_zero_values)
1917 {
1918  Assert(row_indices.size() == values.m(),
1919  ExcDimensionMismatch(row_indices.size(), values.m()));
1920  Assert(col_indices.size() == values.n(),
1921  ExcDimensionMismatch(col_indices.size(), values.n()));
1922 
1923  for (size_type i = 0; i < row_indices.size(); ++i)
1924  add(row_indices[i],
1925  col_indices.size(),
1926  col_indices.data(),
1927  &values(i, 0),
1928  elide_zero_values);
1929 }
1930 
1931 
1932 
1933 template <typename number>
1934 template <typename number2>
1935 inline void
1937  const std::vector<size_type> &col_indices,
1938  const std::vector<number2> & values,
1939  const bool elide_zero_values)
1940 {
1941  Assert(col_indices.size() == values.size(),
1942  ExcDimensionMismatch(col_indices.size(), values.size()));
1943 
1944  add(row,
1945  col_indices.size(),
1946  col_indices.data(),
1947  values.data(),
1948  elide_zero_values);
1949 }
1950 
1951 
1952 
1953 template <typename number>
1954 inline SparseMatrix<number> &
1955 SparseMatrix<number>::operator*=(const number factor)
1956 {
1957  Assert(cols != nullptr, ExcNotInitialized());
1958  Assert(val != nullptr, ExcNotInitialized());
1959 
1960  number * val_ptr = val.get();
1961  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1962 
1963  while (val_ptr != end_ptr)
1964  *val_ptr++ *= factor;
1965 
1966  return *this;
1967 }
1968 
1969 
1970 
1971 template <typename number>
1972 inline SparseMatrix<number> &
1973 SparseMatrix<number>::operator/=(const number factor)
1974 {
1975  Assert(cols != nullptr, ExcNotInitialized());
1976  Assert(val != nullptr, ExcNotInitialized());
1977  Assert(factor != number(), ExcDivideByZero());
1978 
1979  const number factor_inv = number(1.) / factor;
1980 
1981  number * val_ptr = val.get();
1982  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1983 
1984  while (val_ptr != end_ptr)
1985  *val_ptr++ *= factor_inv;
1986 
1987  return *this;
1988 }
1989 
1990 
1991 
1992 template <typename number>
1993 inline const number &
1995 {
1996  Assert(cols != nullptr, ExcNotInitialized());
1997  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
1998  ExcInvalidIndex(i, j));
1999  return val[cols->operator()(i, j)];
2000 }
2001 
2002 
2003 
2004 template <typename number>
2005 inline number &
2007 {
2008  Assert(cols != nullptr, ExcNotInitialized());
2009  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2010  ExcInvalidIndex(i, j));
2011  return val[cols->operator()(i, j)];
2012 }
2013 
2014 
2015 
2016 template <typename number>
2017 inline number
2018 SparseMatrix<number>::el(const size_type i, const size_type j) const
2019 {
2020  Assert(cols != nullptr, ExcNotInitialized());
2021  const size_type index = cols->operator()(i, j);
2022 
2023  if (index != SparsityPattern::invalid_entry)
2024  return val[index];
2025  else
2026  return 0;
2027 }
2028 
2029 
2030 
2031 template <typename number>
2032 inline number
2034 {
2035  Assert(cols != nullptr, ExcNotInitialized());
2036  Assert(m() == n(), ExcNotQuadratic());
2037  AssertIndexRange(i, m());
2038 
2039  // Use that the first element in each row of a quadratic matrix is the main
2040  // diagonal
2041  return val[cols->rowstart[i]];
2042 }
2043 
2044 
2045 
2046 template <typename number>
2047 inline number &
2049 {
2050  Assert(cols != nullptr, ExcNotInitialized());
2051  Assert(m() == n(), ExcNotQuadratic());
2052  AssertIndexRange(i, m());
2053 
2054  // Use that the first element in each row of a quadratic matrix is the main
2055  // diagonal
2056  return val[cols->rowstart[i]];
2057 }
2058 
2059 
2060 
2061 template <typename number>
2062 template <typename ForwardIterator>
2063 void
2064 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2065  const ForwardIterator end)
2066 {
2067  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2068  ExcIteratorRange(std::distance(begin, end), m()));
2069 
2070  // for use in the inner loop, we define an alias to the type of the inner
2071  // iterators
2072  using inner_iterator =
2073  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2074  size_type row = 0;
2075  for (ForwardIterator i = begin; i != end; ++i, ++row)
2076  {
2077  const inner_iterator end_of_row = i->end();
2078  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2079  // write entries
2080  set(row, j->first, j->second);
2081  };
2082 }
2083 
2084 
2085 //---------------------------------------------------------------------------
2086 
2087 
2088 namespace SparseMatrixIterators
2089 {
2090  template <typename number>
2091  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2092  const std::size_t index_within_matrix)
2093  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2094  index_within_matrix)
2095  , matrix(matrix)
2096  {}
2097 
2098 
2099 
2100  template <typename number>
2101  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2102  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2103  , matrix(matrix)
2104  {}
2105 
2106 
2107 
2108  template <typename number>
2109  inline Accessor<number, true>::Accessor(
2111  : SparsityPatternIterators::Accessor(a)
2112  , matrix(&a.get_matrix())
2113  {}
2114 
2115 
2116 
2117  template <typename number>
2118  inline number
2119  Accessor<number, true>::value() const
2120  {
2121  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2122  return matrix->val[linear_index];
2123  }
2124 
2125 
2126 
2127  template <typename number>
2128  inline const typename Accessor<number, true>::MatrixType &
2129  Accessor<number, true>::get_matrix() const
2130  {
2131  return *matrix;
2132  }
2133 
2134 
2135 
2136  template <typename number>
2137  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2138  const bool)
2139  : accessor(accessor)
2140  {}
2141 
2142 
2143  template <typename number>
2144  inline Accessor<number, false>::Reference::operator number() const
2145  {
2146  AssertIndexRange(accessor->linear_index,
2147  accessor->matrix->n_nonzero_elements());
2148  return accessor->matrix->val[accessor->linear_index];
2149  }
2150 
2151 
2152 
2153  template <typename number>
2154  inline const typename Accessor<number, false>::Reference &
2155  Accessor<number, false>::Reference::operator=(const number n) const
2156  {
2157  AssertIndexRange(accessor->linear_index,
2158  accessor->matrix->n_nonzero_elements());
2159  accessor->matrix->val[accessor->linear_index] = n;
2160  return *this;
2161  }
2162 
2163 
2164 
2165  template <typename number>
2166  inline const typename Accessor<number, false>::Reference &
2167  Accessor<number, false>::Reference::operator+=(const number n) const
2168  {
2169  AssertIndexRange(accessor->linear_index,
2170  accessor->matrix->n_nonzero_elements());
2171  accessor->matrix->val[accessor->linear_index] += n;
2172  return *this;
2173  }
2174 
2175 
2176 
2177  template <typename number>
2178  inline const typename Accessor<number, false>::Reference &
2179  Accessor<number, false>::Reference::operator-=(const number n) const
2180  {
2181  AssertIndexRange(accessor->linear_index,
2182  accessor->matrix->n_nonzero_elements());
2183  accessor->matrix->val[accessor->linear_index] -= n;
2184  return *this;
2185  }
2186 
2187 
2188 
2189  template <typename number>
2190  inline const typename Accessor<number, false>::Reference &
2191  Accessor<number, false>::Reference::operator*=(const number n) const
2192  {
2193  AssertIndexRange(accessor->linear_index,
2194  accessor->matrix->n_nonzero_elements());
2195  accessor->matrix->val[accessor->linear_index] *= n;
2196  return *this;
2197  }
2198 
2199 
2200 
2201  template <typename number>
2202  inline const typename Accessor<number, false>::Reference &
2203  Accessor<number, false>::Reference::operator/=(const number n) const
2204  {
2205  AssertIndexRange(accessor->linear_index,
2206  accessor->matrix->n_nonzero_elements());
2207  accessor->matrix->val[accessor->linear_index] /= n;
2208  return *this;
2209  }
2210 
2211 
2212 
2213  template <typename number>
2214  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2215  const std::size_t index)
2216  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2217  , matrix(matrix)
2218  {}
2219 
2220 
2221 
2222  template <typename number>
2223  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2224  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2225  , matrix(matrix)
2226  {}
2227 
2228 
2229 
2230  template <typename number>
2231  inline typename Accessor<number, false>::Reference
2232  Accessor<number, false>::value() const
2233  {
2234  return Reference(this, true);
2235  }
2236 
2237 
2238 
2239  template <typename number>
2240  inline typename Accessor<number, false>::MatrixType &
2241  Accessor<number, false>::get_matrix() const
2242  {
2243  return *matrix;
2244  }
2245 
2246 
2247 
2248  template <typename number, bool Constness>
2249  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2250  const std::size_t index)
2251  : accessor(matrix, index)
2252  {}
2253 
2254 
2255 
2256  template <typename number, bool Constness>
2257  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2258  : accessor(matrix)
2259  {}
2260 
2261 
2262 
2263  template <typename number, bool Constness>
2264  inline Iterator<number, Constness>::Iterator(
2266  : accessor(*i)
2267  {}
2268 
2269 
2270 
2271  template <typename number, bool Constness>
2272  inline const Iterator<number, Constness> &
2273  Iterator<number, Constness>::
2275  {
2276  accessor = *i;
2277  return *this;
2278  }
2279 
2280 
2281 
2282  template <typename number, bool Constness>
2283  inline Iterator<number, Constness> &
2284  Iterator<number, Constness>::operator++()
2285  {
2286  accessor.advance();
2287  return *this;
2288  }
2289 
2290 
2291  template <typename number, bool Constness>
2292  inline Iterator<number, Constness>
2293  Iterator<number, Constness>::operator++(int)
2294  {
2295  const Iterator iter = *this;
2296  accessor.advance();
2297  return iter;
2298  }
2299 
2300 
2301  template <typename number, bool Constness>
2302  inline const Accessor<number, Constness> &Iterator<number, Constness>::
2303  operator*() const
2304  {
2305  return accessor;
2306  }
2307 
2308 
2309  template <typename number, bool Constness>
2310  inline const Accessor<number, Constness> *Iterator<number, Constness>::
2311  operator->() const
2312  {
2313  return &accessor;
2314  }
2315 
2316 
2317  template <typename number, bool Constness>
2318  inline bool
2319  Iterator<number, Constness>::operator==(const Iterator &other) const
2320  {
2321  return (accessor == other.accessor);
2322  }
2323 
2324 
2325  template <typename number, bool Constness>
2326  inline bool
2327  Iterator<number, Constness>::operator!=(const Iterator &other) const
2328  {
2329  return !(*this == other);
2330  }
2331 
2332 
2333  template <typename number, bool Constness>
2334  inline bool
2335  Iterator<number, Constness>::operator<(const Iterator &other) const
2336  {
2337  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2338  ExcInternalError());
2339 
2340  return (accessor < other.accessor);
2341  }
2342 
2343 
2344  template <typename number, bool Constness>
2345  inline bool
2346  Iterator<number, Constness>::operator>(const Iterator &other) const
2347  {
2348  return (other < *this);
2349  }
2350 
2351 
2352  template <typename number, bool Constness>
2353  inline int
2354  Iterator<number, Constness>::operator-(const Iterator &other) const
2355  {
2356  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2357  ExcInternalError());
2358 
2359  return (*this)->linear_index - other->linear_index;
2360  }
2361 
2362 
2363 
2364  template <typename number, bool Constness>
2365  inline Iterator<number, Constness>
2367  {
2368  Iterator x = *this;
2369  for (size_type i = 0; i < n; ++i)
2370  ++x;
2371 
2372  return x;
2373  }
2374 
2375 } // namespace SparseMatrixIterators
2376 
2377 
2378 
2379 template <typename number>
2382 {
2383  return const_iterator(this, 0);
2384 }
2385 
2386 
2387 template <typename number>
2390 {
2391  return const_iterator(this);
2392 }
2393 
2394 
2395 template <typename number>
2396 inline typename SparseMatrix<number>::iterator
2398 {
2399  return iterator(this, 0);
2400 }
2401 
2402 
2403 template <typename number>
2404 inline typename SparseMatrix<number>::iterator
2406 {
2407  return iterator(this, cols->rowstart[cols->rows]);
2408 }
2409 
2410 
2411 template <typename number>
2414 {
2415  Assert(r < m(), ExcIndexRange(r, 0, m()));
2416 
2417  return const_iterator(this, cols->rowstart[r]);
2418 }
2419 
2420 
2421 
2422 template <typename number>
2424 SparseMatrix<number>::end(const size_type r) const
2425 {
2426  Assert(r < m(), ExcIndexRange(r, 0, m()));
2427 
2428  return const_iterator(this, cols->rowstart[r + 1]);
2429 }
2430 
2431 
2432 
2433 template <typename number>
2434 inline typename SparseMatrix<number>::iterator
2436 {
2437  Assert(r < m(), ExcIndexRange(r, 0, m()));
2438 
2439  return iterator(this, cols->rowstart[r]);
2440 }
2441 
2442 
2443 
2444 template <typename number>
2445 inline typename SparseMatrix<number>::iterator
2447 {
2448  Assert(r < m(), ExcIndexRange(r, 0, m()));
2449 
2450  return iterator(this, cols->rowstart[r + 1]);
2451 }
2452 
2453 
2454 
2455 template <typename number>
2456 template <class StreamType>
2457 inline void
2458 SparseMatrix<number>::print(StreamType &out,
2459  const bool across,
2460  const bool diagonal_first) const
2461 {
2462  Assert(cols != nullptr, ExcNotInitialized());
2463  Assert(val != nullptr, ExcNotInitialized());
2464 
2465  bool hanging_diagonal = false;
2466  number diagonal = number();
2467 
2468  for (size_type i = 0; i < cols->rows; ++i)
2469  {
2470  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2471  {
2472  if (!diagonal_first && i == cols->colnums[j])
2473  {
2474  diagonal = val[j];
2475  hanging_diagonal = true;
2476  }
2477  else
2478  {
2479  if (hanging_diagonal && cols->colnums[j] > i)
2480  {
2481  if (across)
2482  out << ' ' << i << ',' << i << ':' << diagonal;
2483  else
2484  out << '(' << i << ',' << i << ") " << diagonal
2485  << std::endl;
2486  hanging_diagonal = false;
2487  }
2488  if (across)
2489  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2490  else
2491  out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2492  << std::endl;
2493  }
2494  }
2495  if (hanging_diagonal)
2496  {
2497  if (across)
2498  out << ' ' << i << ',' << i << ':' << diagonal;
2499  else
2500  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2501  hanging_diagonal = false;
2502  }
2503  }
2504  if (across)
2505  out << std::endl;
2506 }
2507 
2508 
2509 template <typename number>
2510 inline void
2512 {
2513  // nothing to do here
2514 }
2515 
2516 
2517 
2518 template <typename number>
2519 inline void
2521 {
2522  // nothing to do here
2523 }
2524 
2525 # endif // DOXYGEN
2526 
2527 
2528 /*---------------------------- sparse_matrix.h ---------------------------*/
2529 
2530 DEAL_II_NAMESPACE_CLOSE
2531 
2532 #endif
2533 /*---------------------------- sparse_matrix.h ---------------------------*/
typename Accessor< number, Constness >::MatrixType MatrixType
size_type m() const
SparseMatrix & operator/=(const number factor)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
types::global_dof_index size_type
Definition: sparse_matrix.h:79
typename numbers::NumberTraits< number >::real_type real_type
void prepare_add()
const_iterator end() const
std::unique_ptr< number[]> val
void set(const size_type i, const size_type j, const number value)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
number value_type
static ::ExceptionBase & ExcNotInitialized()
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
static ::ExceptionBase & ExcDivideByZero()
LinearAlgebra::distributed::Vector< Number > Vector
std::unique_ptr< std::size_t[]> rowstart
size_type n() const
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:1407
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:496
SparseMatrixIterators::Iterator< number, false > iterator
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
Accessor< number, Constness > accessor
void prepare_set()
SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
void add(const size_type i, const size_type j, const number value)
types::global_dof_index size_type
std::unique_ptr< size_type[]> colnums
static ::ExceptionBase & ExcNotQuadratic()
unsigned int global_dof_index
Definition: types.h:89
Definition: cuda.h:31
number el(const size_type i, const size_type j) const
size_type m() const
static const size_type invalid_entry
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
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
#define AssertIsFinite(number)
Definition: exceptions.h:1669
std::size_t max_len
const Accessor< number, Constness > & value_type
SparseMatrix & operator*=(const number factor)
static ::ExceptionBase & ExcInternalError()