Reference documentation for deal.II version Git 437f94a 2018-06-23 09:16:19 -0400
sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2018 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:
169  MatrixType *matrix;
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:
304  MatrixType *matrix;
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 
384  Iterator &
385  operator++();
386 
390  Iterator
391  operator++(int);
392 
396  const Accessor<number, Constness> &operator*() const;
397 
402 
406  bool
407  operator==(const Iterator &) const;
408 
412  bool
413  operator!=(const Iterator &) const;
414 
422  bool
423  operator<(const Iterator &) const;
424 
429  bool
430  operator>(const Iterator &) const;
431 
438  int
439  operator-(const Iterator &p) const;
440 
444  Iterator
445  operator+(const size_type n) const;
446 
447  private:
452  };
453 
454 } // namespace SparseMatrixIterators
455 
461 // TODO: Add multithreading to the other vmult functions.
462 
491 template <typename number>
492 class SparseMatrix : public virtual Subscriptor
493 {
494 public:
499 
504  typedef number value_type;
505 
516 
522 
530 
537  struct Traits
538  {
543  static const bool zero_addition_can_be_elided = true;
544  };
545 
560  SparseMatrix();
561 
570  SparseMatrix(const SparseMatrix &);
571 
579  SparseMatrix(SparseMatrix<number> &&m) noexcept;
580 
594  explicit SparseMatrix(const SparsityPattern &sparsity);
595 
602  SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
603 
608  virtual ~SparseMatrix() override;
609 
619  SparseMatrix<number> &
620  operator=(const SparseMatrix<number> &);
621 
626  SparseMatrix<number> &
627  operator=(SparseMatrix<number> &&m) noexcept;
628 
635  SparseMatrix<number> &
636  operator=(const IdentityMatrix &id);
637 
649  SparseMatrix &
650  operator=(const double d);
651 
665  virtual void
666  reinit(const SparsityPattern &sparsity);
667 
673  virtual void
674  clear();
676 
684  bool
685  empty() const;
686 
691  size_type
692  m() const;
693 
698  size_type
699  n() const;
700 
704  size_type
705  get_row_length(const size_type row) const;
706 
712  std::size_t
713  n_nonzero_elements() const;
714 
724  std::size_t
725  n_actually_nonzero_elements(const double threshold = 0.) const;
726 
735  const SparsityPattern &
736  get_sparsity_pattern() const;
737 
742  std::size_t
743  memory_consumption() const;
744 
748  void compress(::VectorOperation::values);
749 
751 
760  void
761  set(const size_type i, const size_type j, const number value);
762 
778  template <typename number2>
779  void
780  set(const std::vector<size_type> &indices,
781  const FullMatrix<number2> & full_matrix,
782  const bool elide_zero_values = false);
783 
789  template <typename number2>
790  void
791  set(const std::vector<size_type> &row_indices,
792  const std::vector<size_type> &col_indices,
793  const FullMatrix<number2> & full_matrix,
794  const bool elide_zero_values = false);
795 
806  template <typename number2>
807  void
808  set(const size_type row,
809  const std::vector<size_type> &col_indices,
810  const std::vector<number2> & values,
811  const bool elide_zero_values = false);
812 
822  template <typename number2>
823  void
824  set(const size_type row,
825  const size_type n_cols,
826  const size_type *col_indices,
827  const number2 * values,
828  const bool elide_zero_values = false);
829 
835  void
836  add(const size_type i, const size_type j, const number value);
837 
852  template <typename number2>
853  void
854  add(const std::vector<size_type> &indices,
855  const FullMatrix<number2> & full_matrix,
856  const bool elide_zero_values = true);
857 
863  template <typename number2>
864  void
865  add(const std::vector<size_type> &row_indices,
866  const std::vector<size_type> &col_indices,
867  const FullMatrix<number2> & full_matrix,
868  const bool elide_zero_values = true);
869 
879  template <typename number2>
880  void
881  add(const size_type row,
882  const std::vector<size_type> &col_indices,
883  const std::vector<number2> & values,
884  const bool elide_zero_values = true);
885 
895  template <typename number2>
896  void
897  add(const size_type row,
898  const size_type n_cols,
899  const size_type *col_indices,
900  const number2 * values,
901  const bool elide_zero_values = true,
902  const bool col_indices_are_sorted = false);
903 
907  SparseMatrix &
908  operator*=(const number factor);
909 
913  SparseMatrix &
914  operator/=(const number factor);
915 
928  void
929  symmetrize();
930 
947  template <typename somenumber>
948  SparseMatrix<number> &
949  copy_from(const SparseMatrix<somenumber> &source);
950 
967  template <typename ForwardIterator>
968  void
969  copy_from(const ForwardIterator begin, const ForwardIterator end);
970 
976  template <typename somenumber>
977  void
978  copy_from(const FullMatrix<somenumber> &matrix);
979 
980 # ifdef DEAL_II_WITH_TRILINOS
981 
992 # endif
993 
1005  template <typename somenumber>
1006  void
1007  add(const number factor, const SparseMatrix<somenumber> &matrix);
1008 
1010 
1014 
1028  const number &
1029  operator()(const size_type i, const size_type j) const;
1030 
1034  number &
1035  operator()(const size_type i, const size_type j);
1036 
1049  number
1050  el(const size_type i, const size_type j) const;
1051 
1061  number
1062  diag_element(const size_type i) const;
1063 
1068  number &
1069  diag_element(const size_type i);
1070 
1072 
1092  template <class OutVector, class InVector>
1093  void
1094  vmult(OutVector &dst, const InVector &src) const;
1095 
1111  template <class OutVector, class InVector>
1112  void
1113  Tvmult(OutVector &dst, const InVector &src) const;
1114 
1131  template <class OutVector, class InVector>
1132  void
1133  vmult_add(OutVector &dst, const InVector &src) const;
1134 
1150  template <class OutVector, class InVector>
1151  void
1152  Tvmult_add(OutVector &dst, const InVector &src) const;
1153 
1171  template <typename somenumber>
1172  somenumber
1173  matrix_norm_square(const Vector<somenumber> &v) const;
1174 
1180  template <typename somenumber>
1181  somenumber
1182  matrix_scalar_product(const Vector<somenumber> &u,
1183  const Vector<somenumber> &v) const;
1184 
1194  template <typename somenumber>
1195  somenumber
1196  residual(Vector<somenumber> & dst,
1197  const Vector<somenumber> &x,
1198  const Vector<somenumber> &b) const;
1199 
1235  template <typename numberB, typename numberC>
1236  void
1238  const SparseMatrix<numberB> &B,
1239  const Vector<number> & V = Vector<number>(),
1240  const bool rebuild_sparsity_pattern = true) const;
1241 
1266  template <typename numberB, typename numberC>
1267  void
1269  const SparseMatrix<numberB> &B,
1270  const Vector<number> & V = Vector<number>(),
1271  const bool rebuild_sparsity_pattern = true) const;
1272 
1274 
1278 
1286  real_type
1287  l1_norm() const;
1288 
1296  real_type
1297  linfty_norm() const;
1298 
1303  real_type
1304  frobenius_norm() const;
1306 
1310 
1316  template <typename somenumber>
1317  void
1318  precondition_Jacobi(Vector<somenumber> & dst,
1319  const Vector<somenumber> &src,
1320  const number omega = 1.) const;
1321 
1328  template <typename somenumber>
1329  void
1330  precondition_SSOR(Vector<somenumber> & dst,
1331  const Vector<somenumber> & src,
1332  const number omega = 1.,
1333  const std::vector<std::size_t> &pos_right_of_diagonal =
1334  std::vector<std::size_t>()) const;
1335 
1339  template <typename somenumber>
1340  void
1341  precondition_SOR(Vector<somenumber> & dst,
1342  const Vector<somenumber> &src,
1343  const number om = 1.) const;
1344 
1348  template <typename somenumber>
1349  void
1350  precondition_TSOR(Vector<somenumber> & dst,
1351  const Vector<somenumber> &src,
1352  const number om = 1.) const;
1353 
1359  template <typename somenumber>
1360  void
1361  SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1362 
1367  template <typename somenumber>
1368  void
1369  SOR(Vector<somenumber> &v, const number om = 1.) const;
1370 
1375  template <typename somenumber>
1376  void
1377  TSOR(Vector<somenumber> &v, const number om = 1.) const;
1378 
1389  template <typename somenumber>
1390  void
1391  PSOR(Vector<somenumber> & v,
1392  const std::vector<size_type> &permutation,
1393  const std::vector<size_type> &inverse_permutation,
1394  const number om = 1.) const;
1395 
1406  template <typename somenumber>
1407  void
1408  TPSOR(Vector<somenumber> & v,
1409  const std::vector<size_type> &permutation,
1410  const std::vector<size_type> &inverse_permutation,
1411  const number om = 1.) const;
1412 
1418  template <typename somenumber>
1419  void
1420  Jacobi_step(Vector<somenumber> & v,
1421  const Vector<somenumber> &b,
1422  const number om = 1.) const;
1423 
1428  template <typename somenumber>
1429  void
1430  SOR_step(Vector<somenumber> & v,
1431  const Vector<somenumber> &b,
1432  const number om = 1.) const;
1433 
1438  template <typename somenumber>
1439  void
1440  TSOR_step(Vector<somenumber> & v,
1441  const Vector<somenumber> &b,
1442  const number om = 1.) const;
1443 
1448  template <typename somenumber>
1449  void
1450  SSOR_step(Vector<somenumber> & v,
1451  const Vector<somenumber> &b,
1452  const number om = 1.) const;
1454 
1458 
1465  const_iterator
1466  begin() const;
1467 
1471  iterator
1472  begin();
1473 
1477  const_iterator
1478  end() const;
1479 
1483  iterator
1484  end();
1485 
1495  const_iterator
1496  begin(const size_type r) const;
1497 
1501  iterator
1502  begin(const size_type r);
1503 
1513  const_iterator
1514  end(const size_type r) const;
1515 
1519  iterator
1520  end(const size_type r);
1522 
1526 
1538  template <class StreamType>
1539  void
1540  print(StreamType &out,
1541  const bool across = false,
1542  const bool diagonal_first = true) const;
1543 
1564  void
1565  print_formatted(std::ostream & out,
1566  const unsigned int precision = 3,
1567  const bool scientific = true,
1568  const unsigned int width = 0,
1569  const char * zero_string = " ",
1570  const double denominator = 1.) const;
1571 
1577  void
1578  print_pattern(std::ostream &out, const double threshold = 0.) const;
1579 
1590  void
1591  block_write(std::ostream &out) const;
1592 
1609  void
1610  block_read(std::istream &in);
1612 
1621  int,
1622  int,
1623  << "You are trying to access the matrix entry with index <"
1624  << arg1 << ',' << arg2
1625  << ">, but this entry does not exist in the sparsity pattern "
1626  "of this matrix."
1627  "\n\n"
1628  "The most common cause for this problem is that you used "
1629  "a method to build the sparsity pattern that did not "
1630  "(completely) take into account all of the entries you "
1631  "will later try to write into. An example would be "
1632  "building a sparsity pattern that does not include "
1633  "the entries you will write into due to constraints "
1634  "on degrees of freedom such as hanging nodes or periodic "
1635  "boundary conditions. In such cases, building the "
1636  "sparsity pattern will succeed, but you will get errors "
1637  "such as the current one at one point or other when "
1638  "trying to write into the entries of the matrix.");
1643  "When copying one sparse matrix into another, "
1644  "or when adding one sparse matrix to another, "
1645  "both matrices need to refer to the same "
1646  "sparsity pattern.");
1651  int,
1652  int,
1653  << "The iterators denote a range of " << arg1
1654  << " elements, but the given number of rows was " << arg2);
1659  "You are attempting an operation on two matrices that "
1660  "are the same object, but the operation requires that the "
1661  "two objects are in fact different.");
1663 
1664 protected:
1675  void
1676  prepare_add();
1677 
1682  void
1683  prepare_set();
1684 
1685 private:
1692 
1700  std::unique_ptr<number[]> val;
1701 
1708  std::size_t max_len;
1709 
1710  // make all other sparse matrices friends
1711  template <typename somenumber>
1712  friend class SparseMatrix;
1713  template <typename somenumber>
1714  friend class SparseLUDecomposition;
1715  template <typename>
1716  friend class SparseILU;
1717 
1721  template <typename>
1722  friend class BlockMatrixBase;
1723 
1727  template <typename, bool>
1729  template <typename, bool>
1730  friend class SparseMatrixIterators::Accessor;
1731 
1732 # ifdef DEAL_II_WITH_MPI
1733 
1736  template <typename Number>
1737  friend void
1739  const MPI_Comm &,
1741 # endif
1742 };
1743 
1744 # ifndef DOXYGEN
1745 /*---------------------- Inline functions -----------------------------------*/
1746 
1747 
1748 
1749 template <typename number>
1750 inline typename SparseMatrix<number>::size_type
1752 {
1753  Assert(cols != nullptr, ExcNotInitialized());
1754  return cols->rows;
1755 }
1756 
1757 
1758 template <typename number>
1759 inline typename SparseMatrix<number>::size_type
1761 {
1762  Assert(cols != nullptr, ExcNotInitialized());
1763  return cols->cols;
1764 }
1765 
1766 
1767 // Inline the set() and add() functions, since they will be called frequently.
1768 template <typename number>
1769 inline void
1770 SparseMatrix<number>::set(const size_type i,
1771  const size_type j,
1772  const number value)
1773 {
1774  AssertIsFinite(value);
1775 
1776  const size_type index = cols->operator()(i, j);
1777 
1778  // it is allowed to set elements of the matrix that are not part of the
1779  // sparsity pattern, if the value to which we set it is zero
1780  if (index == SparsityPattern::invalid_entry)
1781  {
1782  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1783  ExcInvalidIndex(i, j));
1784  return;
1785  }
1786 
1787  val[index] = value;
1788 }
1789 
1790 
1791 
1792 template <typename number>
1793 template <typename number2>
1794 inline void
1795 SparseMatrix<number>::set(const std::vector<size_type> &indices,
1796  const FullMatrix<number2> & values,
1797  const bool elide_zero_values)
1798 {
1799  Assert(indices.size() == values.m(),
1800  ExcDimensionMismatch(indices.size(), values.m()));
1801  Assert(values.m() == values.n(), ExcNotQuadratic());
1802 
1803  for (size_type i = 0; i < indices.size(); ++i)
1804  set(indices[i],
1805  indices.size(),
1806  indices.data(),
1807  &values(i, 0),
1808  elide_zero_values);
1809 }
1810 
1811 
1812 
1813 template <typename number>
1814 template <typename number2>
1815 inline void
1816 SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1817  const std::vector<size_type> &col_indices,
1818  const FullMatrix<number2> & values,
1819  const bool elide_zero_values)
1820 {
1821  Assert(row_indices.size() == values.m(),
1822  ExcDimensionMismatch(row_indices.size(), values.m()));
1823  Assert(col_indices.size() == values.n(),
1824  ExcDimensionMismatch(col_indices.size(), values.n()));
1825 
1826  for (size_type i = 0; i < row_indices.size(); ++i)
1827  set(row_indices[i],
1828  col_indices.size(),
1829  col_indices.data(),
1830  &values(i, 0),
1831  elide_zero_values);
1832 }
1833 
1834 
1835 
1836 template <typename number>
1837 template <typename number2>
1838 inline void
1839 SparseMatrix<number>::set(const size_type row,
1840  const std::vector<size_type> &col_indices,
1841  const std::vector<number2> & values,
1842  const bool elide_zero_values)
1843 {
1844  Assert(col_indices.size() == values.size(),
1845  ExcDimensionMismatch(col_indices.size(), values.size()));
1846 
1847  set(row,
1848  col_indices.size(),
1849  col_indices.data(),
1850  values.data(),
1851  elide_zero_values);
1852 }
1853 
1854 
1855 
1856 template <typename number>
1857 inline void
1858 SparseMatrix<number>::add(const size_type i,
1859  const size_type j,
1860  const number value)
1861 {
1862  AssertIsFinite(value);
1863 
1864  if (value == number())
1865  return;
1866 
1867  const size_type index = cols->operator()(i, j);
1868 
1869  // it is allowed to add elements to the matrix that are not part of the
1870  // sparsity pattern, if the value to which we set it is zero
1871  if (index == SparsityPattern::invalid_entry)
1872  {
1873  Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1874  ExcInvalidIndex(i, j));
1875  return;
1876  }
1877 
1878  val[index] += value;
1879 }
1880 
1881 
1882 
1883 template <typename number>
1884 template <typename number2>
1885 inline void
1886 SparseMatrix<number>::add(const std::vector<size_type> &indices,
1887  const FullMatrix<number2> & values,
1888  const bool elide_zero_values)
1889 {
1890  Assert(indices.size() == values.m(),
1891  ExcDimensionMismatch(indices.size(), values.m()));
1892  Assert(values.m() == values.n(), ExcNotQuadratic());
1893 
1894  for (size_type i = 0; i < indices.size(); ++i)
1895  add(indices[i],
1896  indices.size(),
1897  indices.data(),
1898  &values(i, 0),
1899  elide_zero_values);
1900 }
1901 
1902 
1903 
1904 template <typename number>
1905 template <typename number2>
1906 inline void
1907 SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1908  const std::vector<size_type> &col_indices,
1909  const FullMatrix<number2> & values,
1910  const bool elide_zero_values)
1911 {
1912  Assert(row_indices.size() == values.m(),
1913  ExcDimensionMismatch(row_indices.size(), values.m()));
1914  Assert(col_indices.size() == values.n(),
1915  ExcDimensionMismatch(col_indices.size(), values.n()));
1916 
1917  for (size_type i = 0; i < row_indices.size(); ++i)
1918  add(row_indices[i],
1919  col_indices.size(),
1920  col_indices.data(),
1921  &values(i, 0),
1922  elide_zero_values);
1923 }
1924 
1925 
1926 
1927 template <typename number>
1928 template <typename number2>
1929 inline void
1930 SparseMatrix<number>::add(const size_type row,
1931  const std::vector<size_type> &col_indices,
1932  const std::vector<number2> & values,
1933  const bool elide_zero_values)
1934 {
1935  Assert(col_indices.size() == values.size(),
1936  ExcDimensionMismatch(col_indices.size(), values.size()));
1937 
1938  add(row,
1939  col_indices.size(),
1940  col_indices.data(),
1941  values.data(),
1942  elide_zero_values);
1943 }
1944 
1945 
1946 
1947 template <typename number>
1948 inline SparseMatrix<number> &
1949 SparseMatrix<number>::operator*=(const number factor)
1950 {
1951  Assert(cols != nullptr, ExcNotInitialized());
1952  Assert(val != nullptr, ExcNotInitialized());
1953 
1954  number * val_ptr = val.get();
1955  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1956 
1957  while (val_ptr != end_ptr)
1958  *val_ptr++ *= factor;
1959 
1960  return *this;
1961 }
1962 
1963 
1964 
1965 template <typename number>
1966 inline SparseMatrix<number> &
1967 SparseMatrix<number>::operator/=(const number factor)
1968 {
1969  Assert(cols != nullptr, ExcNotInitialized());
1970  Assert(val != nullptr, ExcNotInitialized());
1971  Assert(factor != number(), ExcDivideByZero());
1972 
1973  const number factor_inv = number(1.) / factor;
1974 
1975  number * val_ptr = val.get();
1976  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1977 
1978  while (val_ptr != end_ptr)
1979  *val_ptr++ *= factor_inv;
1980 
1981  return *this;
1982 }
1983 
1984 
1985 
1986 template <typename number>
1987 inline const number &
1988 SparseMatrix<number>::operator()(const size_type i, const size_type j) const
1989 {
1990  Assert(cols != nullptr, ExcNotInitialized());
1991  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
1992  ExcInvalidIndex(i, j));
1993  return val[cols->operator()(i, j)];
1994 }
1995 
1996 
1997 
1998 template <typename number>
1999 inline number &
2000 SparseMatrix<number>::operator()(const size_type i, const size_type j)
2001 {
2002  Assert(cols != nullptr, ExcNotInitialized());
2003  Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2004  ExcInvalidIndex(i, j));
2005  return val[cols->operator()(i, j)];
2006 }
2007 
2008 
2009 
2010 template <typename number>
2011 inline number
2012 SparseMatrix<number>::el(const size_type i, const size_type j) const
2013 {
2014  Assert(cols != nullptr, ExcNotInitialized());
2015  const size_type index = cols->operator()(i, j);
2016 
2017  if (index != SparsityPattern::invalid_entry)
2018  return val[index];
2019  else
2020  return 0;
2021 }
2022 
2023 
2024 
2025 template <typename number>
2026 inline number
2027 SparseMatrix<number>::diag_element(const size_type i) const
2028 {
2029  Assert(cols != nullptr, ExcNotInitialized());
2030  Assert(m() == n(), ExcNotQuadratic());
2031  AssertIndexRange(i, m());
2032 
2033  // Use that the first element in each row of a quadratic matrix is the main
2034  // diagonal
2035  return val[cols->rowstart[i]];
2036 }
2037 
2038 
2039 
2040 template <typename number>
2041 inline number &
2042 SparseMatrix<number>::diag_element(const size_type i)
2043 {
2044  Assert(cols != nullptr, ExcNotInitialized());
2045  Assert(m() == n(), ExcNotQuadratic());
2046  AssertIndexRange(i, m());
2047 
2048  // Use that the first element in each row of a quadratic matrix is the main
2049  // diagonal
2050  return val[cols->rowstart[i]];
2051 }
2052 
2053 
2054 
2055 template <typename number>
2056 template <typename ForwardIterator>
2057 void
2058 SparseMatrix<number>::copy_from(const ForwardIterator begin,
2059  const ForwardIterator end)
2060 {
2061  Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2062  ExcIteratorRange(std::distance(begin, end), m()));
2063 
2064  // for use in the inner loop, we define a typedef to the type of the inner
2065  // iterators
2066  typedef
2067  typename std::iterator_traits<ForwardIterator>::value_type::const_iterator
2068  inner_iterator;
2069  size_type row = 0;
2070  for (ForwardIterator i = begin; i != end; ++i, ++row)
2071  {
2072  const inner_iterator end_of_row = i->end();
2073  for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2074  // write entries
2075  set(row, j->first, j->second);
2076  };
2077 }
2078 
2079 
2080 //---------------------------------------------------------------------------
2081 
2082 
2083 namespace SparseMatrixIterators
2084 {
2085  template <typename number>
2086  inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2087  const std::size_t index_within_matrix)
2088  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2089  index_within_matrix)
2090  , matrix(matrix)
2091  {}
2092 
2093 
2094 
2095  template <typename number>
2096  inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2097  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2098  , matrix(matrix)
2099  {}
2100 
2101 
2102 
2103  template <typename number>
2104  inline Accessor<number, true>::Accessor(
2106  : SparsityPatternIterators::Accessor(a)
2107  , matrix(&a.get_matrix())
2108  {}
2109 
2110 
2111 
2112  template <typename number>
2113  inline number
2114  Accessor<number, true>::value() const
2115  {
2116  AssertIndexRange(index_within_sparsity, matrix->n_nonzero_elements());
2117  return matrix->val[index_within_sparsity];
2118  }
2119 
2120 
2121 
2122  template <typename number>
2123  inline const typename Accessor<number, true>::MatrixType &
2124  Accessor<number, true>::get_matrix() const
2125  {
2126  return *matrix;
2127  }
2128 
2129 
2130 
2131  template <typename number>
2132  inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2133  const bool)
2134  : accessor(accessor)
2135  {}
2136 
2137 
2138  template <typename number>
2139  inline Accessor<number, false>::Reference::operator number() const
2140  {
2141  AssertIndexRange(accessor->index_within_sparsity,
2142  accessor->matrix->n_nonzero_elements());
2143  return accessor->matrix->val[accessor->index_within_sparsity];
2144  }
2145 
2146 
2147 
2148  template <typename number>
2149  inline const typename Accessor<number, false>::Reference &
2150  Accessor<number, false>::Reference::operator=(const number n) const
2151  {
2152  AssertIndexRange(accessor->index_within_sparsity,
2153  accessor->matrix->n_nonzero_elements());
2154  accessor->matrix->val[accessor->index_within_sparsity] = n;
2155  return *this;
2156  }
2157 
2158 
2159 
2160  template <typename number>
2161  inline const typename Accessor<number, false>::Reference &
2162  Accessor<number, false>::Reference::operator+=(const number n) const
2163  {
2164  AssertIndexRange(accessor->index_within_sparsity,
2165  accessor->matrix->n_nonzero_elements());
2166  accessor->matrix->val[accessor->index_within_sparsity] += n;
2167  return *this;
2168  }
2169 
2170 
2171 
2172  template <typename number>
2173  inline const typename Accessor<number, false>::Reference &
2174  Accessor<number, false>::Reference::operator-=(const number n) const
2175  {
2176  AssertIndexRange(accessor->index_within_sparsity,
2177  accessor->matrix->n_nonzero_elements());
2178  accessor->matrix->val[accessor->index_within_sparsity] -= n;
2179  return *this;
2180  }
2181 
2182 
2183 
2184  template <typename number>
2185  inline const typename Accessor<number, false>::Reference &
2186  Accessor<number, false>::Reference::operator*=(const number n) const
2187  {
2188  AssertIndexRange(accessor->index_within_sparsity,
2189  accessor->matrix->n_nonzero_elements());
2190  accessor->matrix->val[accessor->index_within_sparsity] *= n;
2191  return *this;
2192  }
2193 
2194 
2195 
2196  template <typename number>
2197  inline const typename Accessor<number, false>::Reference &
2198  Accessor<number, false>::Reference::operator/=(const number n) const
2199  {
2200  AssertIndexRange(accessor->index_within_sparsity,
2201  accessor->matrix->n_nonzero_elements());
2202  accessor->matrix->val[accessor->index_within_sparsity] /= n;
2203  return *this;
2204  }
2205 
2206 
2207 
2208  template <typename number>
2209  inline Accessor<number, false>::Accessor(MatrixType * matrix,
2210  const std::size_t index)
2211  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2212  , matrix(matrix)
2213  {}
2214 
2215 
2216 
2217  template <typename number>
2218  inline Accessor<number, false>::Accessor(MatrixType *matrix)
2219  : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2220  , matrix(matrix)
2221  {}
2222 
2223 
2224 
2225  template <typename number>
2226  inline typename Accessor<number, false>::Reference
2227  Accessor<number, false>::value() const
2228  {
2229  return Reference(this, true);
2230  }
2231 
2232 
2233 
2234  template <typename number>
2235  inline typename Accessor<number, false>::MatrixType &
2236  Accessor<number, false>::get_matrix() const
2237  {
2238  return *matrix;
2239  }
2240 
2241 
2242 
2243  template <typename number, bool Constness>
2244  inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2245  const std::size_t index)
2246  : accessor(matrix, index)
2247  {}
2248 
2249 
2250 
2251  template <typename number, bool Constness>
2252  inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2253  : accessor(matrix)
2254  {}
2255 
2256 
2257 
2258  template <typename number, bool Constness>
2259  inline Iterator<number, Constness>::Iterator(
2261  : accessor(*i)
2262  {}
2263 
2264 
2265 
2266  template <typename number, bool Constness>
2267  inline Iterator<number, Constness> &
2268  Iterator<number, Constness>::operator++()
2269  {
2270  accessor.advance();
2271  return *this;
2272  }
2273 
2274 
2275  template <typename number, bool Constness>
2276  inline Iterator<number, Constness>
2277  Iterator<number, Constness>::operator++(int)
2278  {
2279  const Iterator iter = *this;
2280  accessor.advance();
2281  return iter;
2282  }
2283 
2284 
2285  template <typename number, bool Constness>
2286  inline const Accessor<number, Constness> &Iterator<number, Constness>::
2287  operator*() const
2288  {
2289  return accessor;
2290  }
2291 
2292 
2293  template <typename number, bool Constness>
2294  inline const Accessor<number, Constness> *Iterator<number, Constness>::
2295  operator->() const
2296  {
2297  return &accessor;
2298  }
2299 
2300 
2301  template <typename number, bool Constness>
2302  inline bool
2303  Iterator<number, Constness>::operator==(const Iterator &other) const
2304  {
2305  return (accessor == other.accessor);
2306  }
2307 
2308 
2309  template <typename number, bool Constness>
2310  inline bool
2311  Iterator<number, Constness>::operator!=(const Iterator &other) const
2312  {
2313  return !(*this == other);
2314  }
2315 
2316 
2317  template <typename number, bool Constness>
2318  inline bool
2319  Iterator<number, Constness>::operator<(const Iterator &other) const
2320  {
2321  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2322  ExcInternalError());
2323 
2324  return (accessor < other.accessor);
2325  }
2326 
2327 
2328  template <typename number, bool Constness>
2329  inline bool
2330  Iterator<number, Constness>::operator>(const Iterator &other) const
2331  {
2332  return (other < *this);
2333  }
2334 
2335 
2336  template <typename number, bool Constness>
2337  inline int
2338  Iterator<number, Constness>::operator-(const Iterator &other) const
2339  {
2340  Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2341  ExcInternalError());
2342 
2343  return (*this)->index_within_sparsity - other->index_within_sparsity;
2344  }
2345 
2346 
2347 
2348  template <typename number, bool Constness>
2349  inline Iterator<number, Constness>
2350  Iterator<number, Constness>::operator+(const size_type n) const
2351  {
2352  Iterator x = *this;
2353  for (size_type i = 0; i < n; ++i)
2354  ++x;
2355 
2356  return x;
2357  }
2358 
2359 } // namespace SparseMatrixIterators
2360 
2361 
2362 
2363 template <typename number>
2366 {
2367  return const_iterator(this, 0);
2368 }
2369 
2370 
2371 template <typename number>
2374 {
2375  return const_iterator(this);
2376 }
2377 
2378 
2379 template <typename number>
2380 inline typename SparseMatrix<number>::iterator
2382 {
2383  return iterator(this, 0);
2384 }
2385 
2386 
2387 template <typename number>
2388 inline typename SparseMatrix<number>::iterator
2390 {
2391  return iterator(this, cols->rowstart[cols->rows]);
2392 }
2393 
2394 
2395 template <typename number>
2397 SparseMatrix<number>::begin(const size_type r) const
2398 {
2399  Assert(r < m(), ExcIndexRange(r, 0, m()));
2400 
2401  return const_iterator(this, cols->rowstart[r]);
2402 }
2403 
2404 
2405 
2406 template <typename number>
2408 SparseMatrix<number>::end(const size_type r) const
2409 {
2410  Assert(r < m(), ExcIndexRange(r, 0, m()));
2411 
2412  return const_iterator(this, cols->rowstart[r + 1]);
2413 }
2414 
2415 
2416 
2417 template <typename number>
2418 inline typename SparseMatrix<number>::iterator
2419 SparseMatrix<number>::begin(const size_type r)
2420 {
2421  Assert(r < m(), ExcIndexRange(r, 0, m()));
2422 
2423  return iterator(this, cols->rowstart[r]);
2424 }
2425 
2426 
2427 
2428 template <typename number>
2429 inline typename SparseMatrix<number>::iterator
2430 SparseMatrix<number>::end(const size_type r)
2431 {
2432  Assert(r < m(), ExcIndexRange(r, 0, m()));
2433 
2434  return iterator(this, cols->rowstart[r + 1]);
2435 }
2436 
2437 
2438 
2439 template <typename number>
2440 template <class StreamType>
2441 inline void
2442 SparseMatrix<number>::print(StreamType &out,
2443  const bool across,
2444  const bool diagonal_first) const
2445 {
2446  Assert(cols != nullptr, ExcNotInitialized());
2447  Assert(val != nullptr, ExcNotInitialized());
2448 
2449  bool hanging_diagonal = false;
2450  number diagonal = number();
2451 
2452  for (size_type i = 0; i < cols->rows; ++i)
2453  {
2454  for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2455  {
2456  if (!diagonal_first && i == cols->colnums[j])
2457  {
2458  diagonal = val[j];
2459  hanging_diagonal = true;
2460  }
2461  else
2462  {
2463  if (hanging_diagonal && cols->colnums[j] > i)
2464  {
2465  if (across)
2466  out << ' ' << i << ',' << i << ':' << diagonal;
2467  else
2468  out << '(' << i << ',' << i << ") " << diagonal
2469  << std::endl;
2470  hanging_diagonal = false;
2471  }
2472  if (across)
2473  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2474  else
2475  out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2476  << std::endl;
2477  }
2478  }
2479  if (hanging_diagonal)
2480  {
2481  if (across)
2482  out << ' ' << i << ',' << i << ':' << diagonal;
2483  else
2484  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2485  hanging_diagonal = false;
2486  }
2487  }
2488  if (across)
2489  out << std::endl;
2490 }
2491 
2492 
2493 template <typename number>
2494 inline void
2496 {
2497  // nothing to do here
2498 }
2499 
2500 
2501 
2502 template <typename number>
2503 inline void
2505 {
2506  // nothing to do here
2507 }
2508 
2509 # endif // DOXYGEN
2510 
2511 
2512 /*---------------------------- sparse_matrix.h ---------------------------*/
2513 
2514 DEAL_II_NAMESPACE_CLOSE
2515 
2516 #endif
2517 /*---------------------------- sparse_matrix.h ---------------------------*/
bool operator>(const Iterator &) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
SparseMatrix & operator/=(const number factor)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
bool operator!=(const Iterator &) const
void Tvmult(OutVector &dst, const InVector &src) const
void vmult(OutVector &dst, const InVector &src) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:420
static::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void Tvmult_add(OutVector &dst, const InVector &src) const
numbers::NumberTraits< number >::real_type real_type
void prepare_add()
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
Iterator operator+(const size_type n) const
static const size_type invalid_entry
void SSOR(Vector< somenumber > &v, const number omega=1.) const
size_type m() const
std::unique_ptr< number[]> val
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
bool empty() const
Iterator(MatrixType *matrix, const std::size_t index_within_matrix)
void set(const size_type i, const size_type j, const number value)
size_type n() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1407
void block_read(std::istream &in)
const Accessor< number, Constness > & value_type
real_type frobenius_norm() const
STL namespace.
static::ExceptionBase & ExcNotInitialized()
std::size_t n_nonzero_elements() const
real_type linfty_norm() const
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static::ExceptionBase & ExcDifferentSparsityPatterns()
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
virtual void reinit(const SparsityPattern &sparsity)
void Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
const Accessor< number, Constness > & operator*() const
static::ExceptionBase & ExcDivideByZero()
void SOR(Vector< somenumber > &v, const number om=1.) const
LinearAlgebra::distributed::Vector< Number > Vector
const_iterator begin() const
Matrix is diagonal.
number el(const size_type i, const size_type j) const
Accessor< number, Constness >::MatrixType MatrixType
void print_pattern(std::ostream &out, const double threshold=0.) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type n() const
void symmetrize()
number diag_element(const size_type i) const
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:89
const SparseMatrix< number > & get_matrix() const
#define Assert(cond, exc)
Definition: exceptions.h:1227
std::size_t memory_consumption() const
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
size_type get_row_length(const size_type row) const
types::global_dof_index size_type
Definition: sparse_matrix.h:79
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:397
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
number value_type
real_type l1_norm() const
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
Accessor< number, Constness > accessor
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const Accessor< number, Constness > * operator->() const
static::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
void prepare_set()
size_type m() const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void add(const size_type i, const size_type j, const number value)
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
const_iterator end() const
bool operator<(const Iterator &) const
int operator-(const Iterator &p) const
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
static::ExceptionBase & ExcNotQuadratic()
Definition: cuda.h:32
void vmult_add(OutVector &dst, const InVector &src) const
types::global_dof_index size_type
SparseMatrixIterators::Iterator< number, false > iterator
SparseMatrixIterators::Iterator< number, true > const_iterator
void TSOR(Vector< somenumber > &v, const number om=1.) const
const SparsityPattern & get_sparsity_pattern() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void block_write(std::ostream &out) const
Accessor(const SparsityPattern *matrix, const std::size_t index_within_sparsity)
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const
static const bool zero_addition_can_be_elided
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
const number & operator()(const size_type i, const size_type j) const
bool operator==(const Iterator &) const
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
std::size_t n_actually_nonzero_elements(const double threshold=0.) const
virtual void clear()
#define AssertIsFinite(number)
Definition: exceptions.h:1423
std::size_t max_len
void compress(::VectorOperation::values)
static::ExceptionBase & ExcSourceEqualsDestination()
SparseMatrix & operator*=(const number factor)
static::ExceptionBase & ExcInternalError()
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const