Reference documentation for deal.II version Git d9b14c2d54 2019-09-14 18:13:29 -0600
\(\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 // 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 
148  Accessor(MatrixType *matrix);
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 
286  Accessor(MatrixType *matrix);
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 
372  Iterator(MatrixType *matrix);
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 
496 template <typename number>
497 class SparseMatrix : public virtual Subscriptor
498 {
499 public:
504 
509  using value_type = number;
510 
521 
527 
535 
542  struct Traits
543  {
548  static const bool zero_addition_can_be_elided = true;
549  };
550 
565  SparseMatrix();
566 
575  SparseMatrix(const SparseMatrix &);
576 
584  SparseMatrix(SparseMatrix<number> &&m) noexcept;
585 
599  explicit SparseMatrix(const SparsityPattern &sparsity);
600 
607  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
608 
613  virtual ~SparseMatrix() override;
614 
625  operator=(const SparseMatrix<number> &);
626 
632  operator=(SparseMatrix<number> &&m) noexcept;
633 
641  operator=(const IdentityMatrix &id);
642 
654  SparseMatrix &
655  operator=(const double d);
656 
670  virtual void
671  reinit(const SparsityPattern &sparsity);
672 
678  virtual void
679  clear();
681 
689  bool
690  empty() const;
691 
696  size_type
697  m() const;
698 
703  size_type
704  n() const;
705 
709  size_type
710  get_row_length(const size_type row) const;
711 
717  std::size_t
718  n_nonzero_elements() const;
719 
729  std::size_t
730  n_actually_nonzero_elements(const double threshold = 0.) const;
731 
740  const SparsityPattern &
741  get_sparsity_pattern() const;
742 
747  std::size_t
748  memory_consumption() const;
749 
753  void compress(::VectorOperation::values);
754 
756 
765  void
766  set(const size_type i, const size_type j, const number value);
767 
783  template <typename number2>
784  void
785  set(const std::vector<size_type> &indices,
786  const FullMatrix<number2> & full_matrix,
787  const bool elide_zero_values = false);
788 
794  template <typename number2>
795  void
796  set(const std::vector<size_type> &row_indices,
797  const std::vector<size_type> &col_indices,
798  const FullMatrix<number2> & full_matrix,
799  const bool elide_zero_values = false);
800 
811  template <typename number2>
812  void
813  set(const size_type row,
814  const std::vector<size_type> &col_indices,
815  const std::vector<number2> & values,
816  const bool elide_zero_values = false);
817 
827  template <typename number2>
828  void
829  set(const size_type row,
830  const size_type n_cols,
831  const size_type *col_indices,
832  const number2 * values,
833  const bool elide_zero_values = false);
834 
840  void
841  add(const size_type i, const size_type j, const number value);
842 
857  template <typename number2>
858  void
859  add(const std::vector<size_type> &indices,
860  const FullMatrix<number2> & full_matrix,
861  const bool elide_zero_values = true);
862 
868  template <typename number2>
869  void
870  add(const std::vector<size_type> &row_indices,
871  const std::vector<size_type> &col_indices,
872  const FullMatrix<number2> & full_matrix,
873  const bool elide_zero_values = true);
874 
884  template <typename number2>
885  void
886  add(const size_type row,
887  const std::vector<size_type> &col_indices,
888  const std::vector<number2> & values,
889  const bool elide_zero_values = true);
890 
900  template <typename number2>
901  void
902  add(const size_type row,
903  const size_type n_cols,
904  const size_type *col_indices,
905  const number2 * values,
906  const bool elide_zero_values = true,
907  const bool col_indices_are_sorted = false);
908 
912  SparseMatrix &
913  operator*=(const number factor);
914 
918  SparseMatrix &
919  operator/=(const number factor);
920 
933  void
934  symmetrize();
935 
952  template <typename somenumber>
954  copy_from(const SparseMatrix<somenumber> &source);
955 
972  template <typename ForwardIterator>
973  void
974  copy_from(const ForwardIterator begin, const ForwardIterator end);
975 
981  template <typename somenumber>
982  void
983  copy_from(const FullMatrix<somenumber> &matrix);
984 
985 # ifdef DEAL_II_WITH_TRILINOS
986 
996  copy_from(const TrilinosWrappers::SparseMatrix &matrix);
997 # endif
998 
1010  template <typename somenumber>
1011  void
1012  add(const number factor, const SparseMatrix<somenumber> &matrix);
1013 
1015 
1019 
1033  const number &
1034  operator()(const size_type i, const size_type j) const;
1035 
1039  number &
1040  operator()(const size_type i, const size_type j);
1041 
1054  number
1055  el(const size_type i, const size_type j) const;
1056 
1066  number
1067  diag_element(const size_type i) const;
1068 
1073  number &
1074  diag_element(const size_type i);
1075 
1077 
1097  template <class OutVector, class InVector>
1098  void
1099  vmult(OutVector &dst, const InVector &src) const;
1100 
1116  template <class OutVector, class InVector>
1117  void
1118  Tvmult(OutVector &dst, const InVector &src) const;
1119 
1136  template <class OutVector, class InVector>
1137  void
1138  vmult_add(OutVector &dst, const InVector &src) const;
1139 
1155  template <class OutVector, class InVector>
1156  void
1157  Tvmult_add(OutVector &dst, const InVector &src) const;
1158 
1176  template <typename somenumber>
1177  somenumber
1178  matrix_norm_square(const Vector<somenumber> &v) const;
1179 
1185  template <typename somenumber>
1186  somenumber
1187  matrix_scalar_product(const Vector<somenumber> &u,
1188  const Vector<somenumber> &v) const;
1189 
1199  template <typename somenumber>
1200  somenumber
1201  residual(Vector<somenumber> & dst,
1202  const Vector<somenumber> &x,
1203  const Vector<somenumber> &b) const;
1204 
1240  template <typename numberB, typename numberC>
1241  void
1242  mmult(SparseMatrix<numberC> & C,
1243  const SparseMatrix<numberB> &B,
1244  const Vector<number> & V = Vector<number>(),
1245  const bool rebuild_sparsity_pattern = true) const;
1246 
1271  template <typename numberB, typename numberC>
1272  void
1273  Tmmult(SparseMatrix<numberC> & C,
1274  const SparseMatrix<numberB> &B,
1275  const Vector<number> & V = Vector<number>(),
1276  const bool rebuild_sparsity_pattern = true) const;
1277 
1279 
1283 
1291  real_type
1292  l1_norm() const;
1293 
1301  real_type
1302  linfty_norm() const;
1303 
1308  real_type
1309  frobenius_norm() const;
1311 
1315 
1321  template <typename somenumber>
1322  void
1323  precondition_Jacobi(Vector<somenumber> & dst,
1324  const Vector<somenumber> &src,
1325  const number omega = 1.) const;
1326 
1333  template <typename somenumber>
1334  void
1335  precondition_SSOR(Vector<somenumber> & dst,
1336  const Vector<somenumber> & src,
1337  const number omega = 1.,
1338  const std::vector<std::size_t> &pos_right_of_diagonal =
1339  std::vector<std::size_t>()) const;
1340 
1344  template <typename somenumber>
1345  void
1346  precondition_SOR(Vector<somenumber> & dst,
1347  const Vector<somenumber> &src,
1348  const number om = 1.) const;
1349 
1353  template <typename somenumber>
1354  void
1355  precondition_TSOR(Vector<somenumber> & dst,
1356  const Vector<somenumber> &src,
1357  const number om = 1.) const;
1358 
1364  template <typename somenumber>
1365  void
1366  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1367 
1372  template <typename somenumber>
1373  void
1374  SOR(Vector<somenumber> &v, const number om = 1.) const;
1375 
1380  template <typename somenumber>
1381  void
1382  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1383 
1394  template <typename somenumber>
1395  void
1396  PSOR(Vector<somenumber> & v,
1397  const std::vector<size_type> &permutation,
1398  const std::vector<size_type> &inverse_permutation,
1399  const number om = 1.) const;
1400 
1411  template <typename somenumber>
1412  void
1413  TPSOR(Vector<somenumber> & v,
1414  const std::vector<size_type> &permutation,
1415  const std::vector<size_type> &inverse_permutation,
1416  const number om = 1.) const;
1417 
1423  template <typename somenumber>
1424  void
1425  Jacobi_step(Vector<somenumber> & v,
1426  const Vector<somenumber> &b,
1427  const number om = 1.) const;
1428 
1433  template <typename somenumber>
1434  void
1435  SOR_step(Vector<somenumber> & v,
1436  const Vector<somenumber> &b,
1437  const number om = 1.) const;
1438 
1443  template <typename somenumber>
1444  void
1445  TSOR_step(Vector<somenumber> & v,
1446  const Vector<somenumber> &b,
1447  const number om = 1.) const;
1448 
1453  template <typename somenumber>
1454  void
1455  SSOR_step(Vector<somenumber> & v,
1456  const Vector<somenumber> &b,
1457  const number om = 1.) const;
1459 
1463 
1471  begin() const;
1472 
1476  iterator
1477  begin();
1478 
1483  end() const;
1484 
1488  iterator
1489  end();
1490 
1501  begin(const size_type r) const;
1502 
1506  iterator
1507  begin(const size_type r);
1508 
1519  end(const size_type r) const;
1520 
1524  iterator
1525  end(const size_type r);
1527 
1531 
1543  template <class StreamType>
1544  void
1545  print(StreamType &out,
1546  const bool across = false,
1547  const bool diagonal_first = true) const;
1548 
1569  void
1570  print_formatted(std::ostream & out,
1571  const unsigned int precision = 3,
1572  const bool scientific = true,
1573  const unsigned int width = 0,
1574  const char * zero_string = " ",
1575  const double denominator = 1.) const;
1576 
1582  void
1583  print_pattern(std::ostream &out, const double threshold = 0.) const;
1584 
1593  void
1594  print_as_numpy_arrays(std::ostream & out,
1595  const unsigned int precision = 9) const;
1596 
1607  void
1608  block_write(std::ostream &out) const;
1609 
1626  void
1627  block_read(std::istream &in);
1629 
1637  DeclException2(ExcInvalidIndex,
1638  int,
1639  int,
1640  << "You are trying to access the matrix entry with index <"
1641  << arg1 << ',' << arg2
1642  << ">, but this entry does not exist in the sparsity pattern "
1643  "of this matrix."
1644  "\n\n"
1645  "The most common cause for this problem is that you used "
1646  "a method to build the sparsity pattern that did not "
1647  "(completely) take into account all of the entries you "
1648  "will later try to write into. An example would be "
1649  "building a sparsity pattern that does not include "
1650  "the entries you will write into due to constraints "
1651  "on degrees of freedom such as hanging nodes or periodic "
1652  "boundary conditions. In such cases, building the "
1653  "sparsity pattern will succeed, but you will get errors "
1654  "such as the current one at one point or other when "
1655  "trying to write into the entries of the matrix.");
1659  DeclExceptionMsg(ExcDifferentSparsityPatterns,
1660  "When copying one sparse matrix into another, "
1661  "or when adding one sparse matrix to another, "
1662  "both matrices need to refer to the same "
1663  "sparsity pattern.");
1667  DeclException2(ExcIteratorRange,
1668  int,
1669  int,
1670  << "The iterators denote a range of " << arg1
1671  << " elements, but the given number of rows was " << arg2);
1675  DeclExceptionMsg(ExcSourceEqualsDestination,
1676  "You are attempting an operation on two matrices that "
1677  "are the same object, but the operation requires that the "
1678  "two objects are in fact different.");
1680 
1681 protected:
1692  void
1693  prepare_add();
1694 
1699  void
1700  prepare_set();
1701 
1702 private:
1709 
1717  std::unique_ptr<number[]> val;
1718 
1725  std::size_t max_len;
1726 
1727  // make all other sparse matrices friends
1728  template <typename somenumber>
1729  friend class SparseMatrix;
1730  template <typename somenumber>
1731  friend class SparseLUDecomposition;
1732  template <typename>
1733  friend class SparseILU;
1734 
1735  // To allow it calling private prepare_add() and prepare_set().
1736  template <typename>
1737  friend class BlockMatrixBase;
1738 
1739  // Also give access to internal details to the iterator/accessor classes.
1740  template <typename, bool>
1741  friend class SparseMatrixIterators::Iterator;
1742  template <typename, bool>
1743  friend class SparseMatrixIterators::Accessor;
1744 
1745 # ifdef DEAL_II_WITH_MPI
1746  // Give access to internal datastructures to perform MPI operations.
1747  template <typename Number>
1748  friend void
1750  const MPI_Comm &,
1752 # endif
1753 };
1754 
1755 # ifndef DOXYGEN
1756 /*---------------------- Inline functions -----------------------------------*/
1757 
1758 
1759 
1760 template <typename number>
1761 inline typename SparseMatrix<number>::size_type
1763 {
1764  Assert(cols != nullptr, ExcNotInitialized());
1765  return cols->rows;
1766 }
1767 
1768 
1769 template <typename number>
1770 inline typename SparseMatrix<number>::size_type
1772 {
1773  Assert(cols != nullptr, ExcNotInitialized());
1774  return cols->cols;
1775 }
1776 
1777 
1778 // Inline the set() and add() functions, since they will be called frequently.
1779 template <typename number>
1780 inline void
1782  const size_type j,
1783  const number value)
1784 {
1785  AssertIsFinite(value);
1786 
1787  const size_type index = cols->operator()(i, j);
1788 
1789  // it is allowed to set elements of the matrix that are not part of the
1790  // sparsity pattern, if the value to which we set it is zero
1791  if (index == SparsityPattern::invalid_entry)
1792  {
1793  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1794  ExcInvalidIndex(i, j));
1795  return;
1796  }
1797 
1798  val[index] = value;
1799 }
1800 
1801 
1802 
1803 template <typename number>
1804 template <typename number2>
1805 inline void
1806 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1807  const FullMatrix<number2> & values,
1808  const bool elide_zero_values)
1809 {
1810  Assert(indices.size() == values.m(),
1811  ExcDimensionMismatch(indices.size(), values.m()));
1812  Assert(values.m() == values.n(), ExcNotQuadratic());
1813 
1814  for (size_type i = 0; i < indices.size(); ++i)
1815  set(indices[i],
1816  indices.size(),
1817  indices.data(),
1818  &values(i, 0),
1819  elide_zero_values);
1820 }
1821 
1822 
1823 
1824 template <typename number>
1825 template <typename number2>
1826 inline void
1827 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1828  const std::vector<size_type> &col_indices,
1829  const FullMatrix<number2> & values,
1830  const bool elide_zero_values)
1831 {
1832  Assert(row_indices.size() == values.m(),
1833  ExcDimensionMismatch(row_indices.size(), values.m()));
1834  Assert(col_indices.size() == values.n(),
1835  ExcDimensionMismatch(col_indices.size(), values.n()));
1836 
1837  for (size_type i = 0; i < row_indices.size(); ++i)
1838  set(row_indices[i],
1839  col_indices.size(),
1840  col_indices.data(),
1841  &values(i, 0),
1842  elide_zero_values);
1843 }
1844 
1845 
1846 
1847 template <typename number>
1848 template <typename number2>
1849 inline void
1851  const std::vector<size_type> &col_indices,
1852  const std::vector<number2> & values,
1853  const bool elide_zero_values)
1854 {
1855  Assert(col_indices.size() == values.size(),
1856  ExcDimensionMismatch(col_indices.size(), values.size()));
1857 
1858  set(row,
1859  col_indices.size(),
1860  col_indices.data(),
1861  values.data(),
1862  elide_zero_values);
1863 }
1864 
1865 
1866 
1867 template <typename number>
1868 inline void
1870  const size_type j,
1871  const number value)
1872 {
1873  AssertIsFinite(value);
1874 
1875  if (value == number())
1876  return;
1877 
1878  const size_type index = cols->operator()(i, j);
1879 
1880  // it is allowed to add elements to the matrix that are not part of the
1881  // sparsity pattern, if the value to which we set it is zero
1882  if (index == SparsityPattern::invalid_entry)
1883  {
1884  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1885  ExcInvalidIndex(i, j));
1886  return;
1887  }
1888 
1889  val[index] += value;
1890 }
1891 
1892 
1893 
1894 template <typename number>
1895 template <typename number2>
1896 inline void
1897 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1898  const FullMatrix<number2> & values,
1899  const bool elide_zero_values)
1900 {
1901  Assert(indices.size() == values.m(),
1902  ExcDimensionMismatch(indices.size(), values.m()));
1903  Assert(values.m() == values.n(), ExcNotQuadratic());
1904 
1905  for (size_type i = 0; i < indices.size(); ++i)
1906  add(indices[i],
1907  indices.size(),
1908  indices.data(),
1909  &values(i, 0),
1910  elide_zero_values);
1911 }
1912 
1913 
1914 
1915 template <typename number>
1916 template <typename number2>
1917 inline void
1918 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1919  const std::vector<size_type> &col_indices,
1920  const FullMatrix<number2> & values,
1921  const bool elide_zero_values)
1922 {
1923  Assert(row_indices.size() == values.m(),
1924  ExcDimensionMismatch(row_indices.size(), values.m()));
1925  Assert(col_indices.size() == values.n(),
1926  ExcDimensionMismatch(col_indices.size(), values.n()));
1927 
1928  for (size_type i = 0; i < row_indices.size(); ++i)
1929  add(row_indices[i],
1930  col_indices.size(),
1931  col_indices.data(),
1932  &values(i, 0),
1933  elide_zero_values);
1934 }
1935 
1936 
1937 
1938 template <typename number>
1939 template <typename number2>
1940 inline void
1942  const std::vector<size_type> &col_indices,
1943  const std::vector<number2> & values,
1944  const bool elide_zero_values)
1945 {
1946  Assert(col_indices.size() == values.size(),
1947  ExcDimensionMismatch(col_indices.size(), values.size()));
1948 
1949  add(row,
1950  col_indices.size(),
1951  col_indices.data(),
1952  values.data(),
1953  elide_zero_values);
1954 }
1955 
1956 
1957 
1958 template <typename number>
1959 inline SparseMatrix<number> &
1960 SparseMatrix<number>::operator*=(const number factor)
1961 {
1962  Assert(cols != nullptr, ExcNotInitialized());
1963  Assert(val != nullptr, ExcNotInitialized());
1964 
1965  number * val_ptr = val.get();
1966  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1967 
1968  while (val_ptr != end_ptr)
1969  *val_ptr++ *= factor;
1970 
1971  return *this;
1972 }
1973 
1974 
1975 
1976 template <typename number>
1977 inline SparseMatrix<number> &
1978 SparseMatrix<number>::operator/=(const number factor)
1979 {
1980  Assert(cols != nullptr, ExcNotInitialized());
1981  Assert(val != nullptr, ExcNotInitialized());
1982  Assert(factor != number(), ExcDivideByZero());
1983 
1984  const number factor_inv = number(1.) / factor;
1985 
1986  number * val_ptr = val.get();
1987  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1988 
1989  while (val_ptr != end_ptr)
1990  *val_ptr++ *= factor_inv;
1991 
1992  return *this;
1993 }
1994 
1995 
1996 
1997 template <typename number>
1998 inline const number &
2000 {
2001  Assert(cols != nullptr, ExcNotInitialized());
2002  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2003  ExcInvalidIndex(i, j));
2004  return val[cols->operator()(i, j)];
2005 }
2006 
2007 
2008 
2009 template <typename number>
2010 inline number &
2012 {
2013  Assert(cols != nullptr, ExcNotInitialized());
2014  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2015  ExcInvalidIndex(i, j));
2016  return val[cols->operator()(i, j)];
2017 }
2018 
2019 
2020 
2021 template <typename number>
2022 inline number
2023 SparseMatrix<number>::el(const size_type i, const size_type j) const
2024 {
2025  Assert(cols != nullptr, ExcNotInitialized());
2026  const size_type index = cols->operator()(i, j);
2027 
2028  if (index != SparsityPattern::invalid_entry)
2029  return val[index];
2030  else
2031  return 0;
2032 }
2033 
2034 
2035 
2036 template <typename number>
2037 inline number
2039 {
2040  Assert(cols != nullptr, ExcNotInitialized());
2041  Assert(m() == n(), ExcNotQuadratic());
2042  AssertIndexRange(i, m());
2043 
2044  // Use that the first element in each row of a quadratic matrix is the main
2045  // diagonal
2046  return val[cols->rowstart[i]];
2047 }
2048 
2049 
2050 
2051 template <typename number>
2052 inline number &
2054 {
2055  Assert(cols != nullptr, ExcNotInitialized());
2056  Assert(m() == n(), ExcNotQuadratic());
2057  AssertIndexRange(i, m());
2058 
2059  // Use that the first element in each row of a quadratic matrix is the main
2060  // diagonal
2061  return val[cols->rowstart[i]];
2062 }
2063 
2064 
2065 
2066 template <typename number>
2067 template <typename ForwardIterator>
2068 void
2069 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2070  const ForwardIterator end)
2071 {
2072  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2073  ExcIteratorRange(std::distance(begin, end), m()));
2074 
2075  // for use in the inner loop, we define an alias to the type of the inner
2076  // iterators
2077  using inner_iterator =
2078  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2079  size_type row = 0;
2080  for (ForwardIterator i = begin; i != end; ++i, ++row)
2081  {
2082  const inner_iterator end_of_row = i->end();
2083  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2084  // write entries
2085  set(row, j->first, j->second);
2086  };
2087 }
2088 
2089 
2090 //---------------------------------------------------------------------------
2091 
2092 
2093 namespace SparseMatrixIterators
2094 {
2095  template <typename number>
2096  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2097  const std::size_t index_within_matrix)
2098  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2099  index_within_matrix)
2100  , matrix(matrix)
2101  {}
2102 
2103 
2104 
2105  template <typename number>
2106  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2107  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2108  , matrix(matrix)
2109  {}
2110 
2111 
2112 
2113  template <typename number>
2114  inline Accessor<number, true>::Accessor(
2116  : SparsityPatternIterators::Accessor(a)
2117  , matrix(&a.get_matrix())
2118  {}
2119 
2120 
2121 
2122  template <typename number>
2123  inline number
2124  Accessor<number, true>::value() const
2125  {
2126  AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2127  return matrix->val[linear_index];
2128  }
2129 
2130 
2131 
2132  template <typename number>
2133  inline const typename Accessor<number, true>::MatrixType &
2134  Accessor<number, true>::get_matrix() const
2135  {
2136  return *matrix;
2137  }
2138 
2139 
2140 
2141  template <typename number>
2142  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2143  const bool)
2144  : accessor(accessor)
2145  {}
2146 
2147 
2148  template <typename number>
2149  inline Accessor<number, false>::Reference::operator number() const
2150  {
2151  AssertIndexRange(accessor->linear_index,
2152  accessor->matrix->n_nonzero_elements());
2153  return accessor->matrix->val[accessor->linear_index];
2154  }
2155 
2156 
2157 
2158  template <typename number>
2159  inline const typename Accessor<number, false>::Reference &
2160  Accessor<number, false>::Reference::operator=(const number n) const
2161  {
2162  AssertIndexRange(accessor->linear_index,
2163  accessor->matrix->n_nonzero_elements());
2164  accessor->matrix->val[accessor->linear_index] = n;
2165  return *this;
2166  }
2167 
2168 
2169 
2170  template <typename number>
2171  inline const typename Accessor<number, false>::Reference &
2172  Accessor<number, false>::Reference::operator+=(const number n) const
2173  {
2174  AssertIndexRange(accessor->linear_index,
2175  accessor->matrix->n_nonzero_elements());
2176  accessor->matrix->val[accessor->linear_index] += n;
2177  return *this;
2178  }
2179 
2180 
2181 
2182  template <typename number>
2183  inline const typename Accessor<number, false>::Reference &
2184  Accessor<number, false>::Reference::operator-=(const number n) const
2185  {
2186  AssertIndexRange(accessor->linear_index,
2187  accessor->matrix->n_nonzero_elements());
2188  accessor->matrix->val[accessor->linear_index] -= n;
2189  return *this;
2190  }
2191 
2192 
2193 
2194  template <typename number>
2195  inline const typename Accessor<number, false>::Reference &
2196  Accessor<number, false>::Reference::operator*=(const number n) const
2197  {
2198  AssertIndexRange(accessor->linear_index,
2199  accessor->matrix->n_nonzero_elements());
2200  accessor->matrix->val[accessor->linear_index] *= n;
2201  return *this;
2202  }
2203 
2204 
2205 
2206  template <typename number>
2207  inline const typename Accessor<number, false>::Reference &
2208  Accessor<number, false>::Reference::operator/=(const number n) const
2209  {
2210  AssertIndexRange(accessor->linear_index,
2211  accessor->matrix->n_nonzero_elements());
2212  accessor->matrix->val[accessor->linear_index] /= n;
2213  return *this;
2214  }
2215 
2216 
2217 
2218  template <typename number>
2219  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2220  const std::size_t index)
2221  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2222  , matrix(matrix)
2223  {}
2224 
2225 
2226 
2227  template <typename number>
2228  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2229  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2230  , matrix(matrix)
2231  {}
2232 
2233 
2234 
2235  template <typename number>
2236  inline typename Accessor<number, false>::Reference
2237  Accessor<number, false>::value() const
2238  {
2239  return Reference(this, true);
2240  }
2241 
2242 
2243 
2244  template <typename number>
2245  inline typename Accessor<number, false>::MatrixType &
2246  Accessor<number, false>::get_matrix() const
2247  {
2248  return *matrix;
2249  }
2250 
2251 
2252 
2253  template <typename number, bool Constness>
2254  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2255  const std::size_t index)
2256  : accessor(matrix, index)
2257  {}
2258 
2259 
2260 
2261  template <typename number, bool Constness>
2262  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2263  : accessor(matrix)
2264  {}
2265 
2266 
2267 
2268  template <typename number, bool Constness>
2269  inline Iterator<number, Constness>::Iterator(
2271  : accessor(*i)
2272  {}
2273 
2274 
2275 
2276  template <typename number, bool Constness>
2277  inline const Iterator<number, Constness> &
2278  Iterator<number, Constness>::
2280  {
2281  accessor = *i;
2282  return *this;
2283  }
2284 
2285 
2286 
2287  template <typename number, bool Constness>
2288  inline Iterator<number, Constness> &
2289  Iterator<number, Constness>::operator++()
2290  {
2291  accessor.advance();
2292  return *this;
2293  }
2294 
2295 
2296  template <typename number, bool Constness>
2297  inline Iterator<number, Constness>
2298  Iterator<number, Constness>::operator++(int)
2299  {
2300  const Iterator iter = *this;
2301  accessor.advance();
2302  return iter;
2303  }
2304 
2305 
2306  template <typename number, bool Constness>
2307  inline const Accessor<number, Constness> &Iterator<number, Constness>::
2308  operator*() const
2309  {
2310  return accessor;
2311  }
2312 
2313 
2314  template <typename number, bool Constness>
2315  inline const Accessor<number, Constness> *Iterator<number, Constness>::
2316  operator->() const
2317  {
2318  return &accessor;
2319  }
2320 
2321 
2322  template <typename number, bool Constness>
2323  inline bool
2324  Iterator<number, Constness>::operator==(const Iterator &other) const
2325  {
2326  return (accessor == other.accessor);
2327  }
2328 
2329 
2330  template <typename number, bool Constness>
2331  inline bool
2332  Iterator<number, Constness>::operator!=(const Iterator &other) const
2333  {
2334  return !(*this == other);
2335  }
2336 
2337 
2338  template <typename number, bool Constness>
2339  inline bool
2340  Iterator<number, Constness>::operator<(const Iterator &other) const
2341  {
2342  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2343  ExcInternalError());
2344 
2345  return (accessor < other.accessor);
2346  }
2347 
2348 
2349  template <typename number, bool Constness>
2350  inline bool
2351  Iterator<number, Constness>::operator>(const Iterator &other) const
2352  {
2353  return (other < *this);
2354  }
2355 
2356 
2357  template <typename number, bool Constness>
2358  inline int
2359  Iterator<number, Constness>::operator-(const Iterator &other) const
2360  {
2361  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2362  ExcInternalError());
2363 
2364  return (*this)->linear_index - other->linear_index;
2365  }
2366 
2367 
2368 
2369  template <typename number, bool Constness>
2370  inline Iterator<number, Constness>
2372  {
2373  Iterator x = *this;
2374  for (size_type i = 0; i < n; ++i)
2375  ++x;
2376 
2377  return x;
2378  }
2379 
2380 } // namespace SparseMatrixIterators
2381 
2382 
2383 
2384 template <typename number>
2387 {
2388  return const_iterator(this, 0);
2389 }
2390 
2391 
2392 template <typename number>
2395 {
2396  return const_iterator(this);
2397 }
2398 
2399 
2400 template <typename number>
2401 inline typename SparseMatrix<number>::iterator
2403 {
2404  return iterator(this, 0);
2405 }
2406 
2407 
2408 template <typename number>
2409 inline typename SparseMatrix<number>::iterator
2411 {
2412  return iterator(this, cols->rowstart[cols->rows]);
2413 }
2414 
2415 
2416 template <typename number>
2419 {
2420  Assert(r < m(), ExcIndexRange(r, 0, m()));
2421 
2422  return const_iterator(this, cols->rowstart[r]);
2423 }
2424 
2425 
2426 
2427 template <typename number>
2429 SparseMatrix<number>::end(const size_type r) const
2430 {
2431  Assert(r < m(), ExcIndexRange(r, 0, m()));
2432 
2433  return const_iterator(this, cols->rowstart[r + 1]);
2434 }
2435 
2436 
2437 
2438 template <typename number>
2439 inline typename SparseMatrix<number>::iterator
2441 {
2442  Assert(r < m(), ExcIndexRange(r, 0, m()));
2443 
2444  return iterator(this, cols->rowstart[r]);
2445 }
2446 
2447 
2448 
2449 template <typename number>
2450 inline typename SparseMatrix<number>::iterator
2452 {
2453  Assert(r < m(), ExcIndexRange(r, 0, m()));
2454 
2455  return iterator(this, cols->rowstart[r + 1]);
2456 }
2457 
2458 
2459 
2460 template <typename number>
2461 template <class StreamType>
2462 inline void
2463 SparseMatrix<number>::print(StreamType &out,
2464  const bool across,
2465  const bool diagonal_first) const
2466 {
2467  Assert(cols != nullptr, ExcNotInitialized());
2468  Assert(val != nullptr, ExcNotInitialized());
2469 
2470  bool hanging_diagonal = false;
2471  number diagonal = number();
2472 
2473  for (size_type i = 0; i < cols->rows; ++i)
2474  {
2475  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2476  {
2477  if (!diagonal_first && i == cols->colnums[j])
2478  {
2479  diagonal = val[j];
2480  hanging_diagonal = true;
2481  }
2482  else
2483  {
2484  if (hanging_diagonal && cols->colnums[j] > i)
2485  {
2486  if (across)
2487  out << ' ' << i << ',' << i << ':' << diagonal;
2488  else
2489  out << '(' << i << ',' << i << ") " << diagonal
2490  << std::endl;
2491  hanging_diagonal = false;
2492  }
2493  if (across)
2494  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2495  else
2496  out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2497  << std::endl;
2498  }
2499  }
2500  if (hanging_diagonal)
2501  {
2502  if (across)
2503  out << ' ' << i << ',' << i << ':' << diagonal;
2504  else
2505  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2506  hanging_diagonal = false;
2507  }
2508  }
2509  if (across)
2510  out << std::endl;
2511 }
2512 
2513 
2514 template <typename number>
2515 inline void
2517 {
2518  // nothing to do here
2519 }
2520 
2521 
2522 
2523 template <typename number>
2524 inline void
2526 {
2527  // nothing to do here
2528 }
2529 
2530 # endif // DOXYGEN
2531 
2532 
2533 /*---------------------------- sparse_matrix.h ---------------------------*/
2534 
2535 DEAL_II_NAMESPACE_CLOSE
2536 
2537 #endif
2538 /*---------------------------- sparse_matrix.h ---------------------------*/
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
typename Accessor< number, Constness >::MatrixType MatrixType
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
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:82
typename numbers::NumberTraits< number >::real_type real_type
void prepare_add()
const_iterator end() const
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
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()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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()
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
DEAL_II_CONSTEXPR SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
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:1669
std::size_t max_len
const Accessor< number, Constness > & value_type
SparseMatrix & operator*=(const number factor)
static ::ExceptionBase & ExcInternalError()