Reference documentation for deal.II version Git 05e4468 2017-09-21 10:18:23 +0200
sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2017 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 at
12 // the top level of the deal.II distribution.
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 #include <deal.II/base/subscriptor.h>
22 #include <deal.II/base/smartpointer.h>
23 #include <deal.II/lac/sparsity_pattern.h>
24 #include <deal.II/lac/identity_matrix.h>
25 #include <deal.II/lac/exceptions.h>
26 #include <deal.II/lac/vector.h>
27 
28 #include <memory>
29 
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 template <typename number> class Vector;
34 template <typename number> class FullMatrix;
35 template <typename Matrix> class BlockMatrixBase;
36 template <typename number> class SparseILU;
37 
38 #ifdef DEAL_II_WITH_TRILINOS
39 namespace TrilinosWrappers
40 {
41  class SparseMatrix;
42 }
43 #endif
44 
55 {
60 
61  // forward declaration
62  template <typename number, bool Constness>
63  class Iterator;
64 
75  template <typename number, bool Constness>
77  {
78  public:
82  number value() const;
83 
87  number &value();
88 
93  const SparseMatrix<number> &get_matrix () const;
94  };
95 
96 
97 
104  template <typename number>
105  class Accessor<number,true> : public SparsityPatternIterators::Accessor
106  {
107  public:
113 
117  Accessor (MatrixType *matrix,
118  const std::size_t index_within_matrix);
119 
123  Accessor (MatrixType *matrix);
124 
129 
133  number value() const;
134 
139  const MatrixType &get_matrix () const;
140 
141  private:
145  MatrixType *matrix;
146 
151 
155  template <typename, bool>
156  friend class Iterator;
157  };
158 
159 
166  template <typename number>
167  class Accessor<number,false> : public SparsityPatternIterators::Accessor
168  {
169  private:
194  class Reference
195  {
196  public:
201  Reference (const Accessor *accessor,
202  const bool dummy);
203 
207  operator number () const;
208 
212  const Reference &operator = (const number n) const;
213 
217  const Reference &operator += (const number n) const;
218 
222  const Reference &operator -= (const number n) const;
223 
227  const Reference &operator *= (const number n) const;
228 
232  const Reference &operator /= (const number n) const;
233 
234  private:
240  };
241 
242  public:
248 
252  Accessor (MatrixType *matrix,
253  const std::size_t index);
254 
258  Accessor (MatrixType *matrix);
259 
263  Reference value() const;
264 
269  MatrixType &get_matrix () const;
270 
271  private:
275  MatrixType *matrix;
276 
281 
285  template <typename, bool>
286  friend class Iterator;
287  };
288 
289 
290 
320  template <typename number, bool Constness>
321  class Iterator
322  {
323  public:
327  typedef
330 
335  typedef
337 
342  Iterator (MatrixType *matrix,
343  const std::size_t index_within_matrix);
344 
348  Iterator (MatrixType *matrix);
349 
355 
360 
364  Iterator operator++ (int);
365 
369  const Accessor<number,Constness> &operator* () const;
370 
375 
379  bool operator == (const Iterator &) const;
380 
384  bool operator != (const Iterator &) const;
385 
393  bool operator < (const Iterator &) const;
394 
399  bool operator > (const Iterator &) const;
400 
407  int operator - (const Iterator &p) const;
408 
412  Iterator operator + (const size_type n) const;
413 
414  private:
419  };
420 
421 }
422 
428 //TODO: Add multithreading to the other vmult functions.
429 
458 template <typename number>
459 class SparseMatrix : public virtual Subscriptor
460 {
461 public:
466 
471  typedef number value_type;
472 
483 
488  typedef
491 
498  typedef
501 
508  struct Traits
509  {
514  static const bool zero_addition_can_be_elided = true;
515  };
516 
531  SparseMatrix ();
532 
541  SparseMatrix (const SparseMatrix &);
542 
551 
565  explicit SparseMatrix (const SparsityPattern &sparsity);
566 
573  SparseMatrix (const SparsityPattern &sparsity,
574  const IdentityMatrix &id);
575 
580  virtual ~SparseMatrix ();
581 
592 
593 #ifdef DEAL_II_WITH_CXX11
594 
599 #endif
600 
608  operator= (const IdentityMatrix &id);
609 
621  SparseMatrix &operator = (const double d);
622 
636  virtual void reinit (const SparsityPattern &sparsity);
637 
643  virtual void clear ();
645 
653  bool empty () const;
654 
659  size_type m () const;
660 
665  size_type n () const;
666 
670  size_type get_row_length (const size_type row) const;
671 
677  std::size_t n_nonzero_elements () const;
678 
688  std::size_t n_actually_nonzero_elements (const double threshold = 0.) const;
689 
698  const SparsityPattern &get_sparsity_pattern () const;
699 
704  std::size_t memory_consumption () const;
705 
710 
712 
721  void set (const size_type i,
722  const size_type j,
723  const number value);
724 
740  template <typename number2>
741  void set (const std::vector<size_type> &indices,
742  const FullMatrix<number2> &full_matrix,
743  const bool elide_zero_values = false);
744 
750  template <typename number2>
751  void set (const std::vector<size_type> &row_indices,
752  const std::vector<size_type> &col_indices,
753  const FullMatrix<number2> &full_matrix,
754  const bool elide_zero_values = false);
755 
766  template <typename number2>
767  void set (const size_type row,
768  const std::vector<size_type> &col_indices,
769  const std::vector<number2> &values,
770  const bool elide_zero_values = false);
771 
781  template <typename number2>
782  void set (const size_type row,
783  const size_type n_cols,
784  const size_type *col_indices,
785  const number2 *values,
786  const bool elide_zero_values = false);
787 
793  void add (const size_type i,
794  const size_type j,
795  const number value);
796 
811  template <typename number2>
812  void add (const std::vector<size_type> &indices,
813  const FullMatrix<number2> &full_matrix,
814  const bool elide_zero_values = true);
815 
821  template <typename number2>
822  void add (const std::vector<size_type> &row_indices,
823  const std::vector<size_type> &col_indices,
824  const FullMatrix<number2> &full_matrix,
825  const bool elide_zero_values = true);
826 
836  template <typename number2>
837  void add (const size_type row,
838  const std::vector<size_type> &col_indices,
839  const std::vector<number2> &values,
840  const bool elide_zero_values = true);
841 
851  template <typename number2>
852  void add (const size_type row,
853  const size_type n_cols,
854  const size_type *col_indices,
855  const number2 *values,
856  const bool elide_zero_values = true,
857  const bool col_indices_are_sorted = false);
858 
862  SparseMatrix &operator *= (const number factor);
863 
867  SparseMatrix &operator /= (const number factor);
868 
881  void symmetrize ();
882 
899  template <typename somenumber>
901  copy_from (const SparseMatrix<somenumber> &source);
902 
919  template <typename ForwardIterator>
920  void copy_from (const ForwardIterator begin,
921  const ForwardIterator end);
922 
928  template <typename somenumber>
929  void copy_from (const FullMatrix<somenumber> &matrix);
930 
931 #ifdef DEAL_II_WITH_TRILINOS
932 
943 #endif
944 
956  template <typename somenumber>
957  void add (const number factor,
958  const SparseMatrix<somenumber> &matrix);
959 
961 
965 
979  number operator () (const size_type i,
980  const size_type j) const;
981 
994  number el (const size_type i,
995  const size_type j) const;
996 
1006  number diag_element (const size_type i) const;
1007 
1012  number &diag_element (const size_type i);
1013 
1015 
1035  template <class OutVector, class InVector>
1036  void vmult (OutVector &dst,
1037  const InVector &src) const;
1038 
1054  template <class OutVector, class InVector>
1055  void Tvmult (OutVector &dst,
1056  const InVector &src) const;
1057 
1074  template <class OutVector, class InVector>
1075  void vmult_add (OutVector &dst,
1076  const InVector &src) const;
1077 
1093  template <class OutVector, class InVector>
1094  void Tvmult_add (OutVector &dst,
1095  const InVector &src) const;
1096 
1114  template <typename somenumber>
1115  somenumber matrix_norm_square (const Vector<somenumber> &v) const;
1116 
1122  template <typename somenumber>
1123  somenumber matrix_scalar_product (const Vector<somenumber> &u,
1124  const Vector<somenumber> &v) const;
1125 
1135  template <typename somenumber>
1136  somenumber residual (Vector<somenumber> &dst,
1137  const Vector<somenumber> &x,
1138  const Vector<somenumber> &b) const;
1139 
1175  template <typename numberB, typename numberC>
1176  void mmult (SparseMatrix<numberC> &C,
1177  const SparseMatrix<numberB> &B,
1178  const Vector<number> &V = Vector<number>(),
1179  const bool rebuild_sparsity_pattern = true) const;
1180 
1205  template <typename numberB, typename numberC>
1206  void Tmmult (SparseMatrix<numberC> &C,
1207  const SparseMatrix<numberB> &B,
1208  const Vector<number> &V = Vector<number>(),
1209  const bool rebuild_sparsity_pattern = true) const;
1210 
1212 
1216 
1224  real_type l1_norm () const;
1225 
1233  real_type linfty_norm () const;
1234 
1239  real_type frobenius_norm () const;
1241 
1245 
1251  template <typename somenumber>
1253  const Vector<somenumber> &src,
1254  const number omega = 1.) const;
1255 
1262  template <typename somenumber>
1264  const Vector<somenumber> &src,
1265  const number omega = 1.,
1266  const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>()) const;
1267 
1271  template <typename somenumber>
1273  const Vector<somenumber> &src,
1274  const number om = 1.) const;
1275 
1279  template <typename somenumber>
1281  const Vector<somenumber> &src,
1282  const number om = 1.) const;
1283 
1289  template <typename somenumber>
1290  void SSOR (Vector<somenumber> &v,
1291  const number omega = 1.) const;
1292 
1297  template <typename somenumber>
1298  void SOR (Vector<somenumber> &v,
1299  const number om = 1.) const;
1300 
1305  template <typename somenumber>
1306  void TSOR (Vector<somenumber> &v,
1307  const number om = 1.) const;
1308 
1319  template <typename somenumber>
1320  void PSOR (Vector<somenumber> &v,
1321  const std::vector<size_type> &permutation,
1322  const std::vector<size_type> &inverse_permutation,
1323  const number om = 1.) const;
1324 
1335  template <typename somenumber>
1336  void TPSOR (Vector<somenumber> &v,
1337  const std::vector<size_type> &permutation,
1338  const std::vector<size_type> &inverse_permutation,
1339  const number om = 1.) const;
1340 
1346  template <typename somenumber>
1348  const Vector<somenumber> &b,
1349  const number om = 1.) const;
1350 
1355  template <typename somenumber>
1356  void SOR_step (Vector<somenumber> &v,
1357  const Vector<somenumber> &b,
1358  const number om = 1.) const;
1359 
1364  template <typename somenumber>
1365  void TSOR_step (Vector<somenumber> &v,
1366  const Vector<somenumber> &b,
1367  const number om = 1.) const;
1368 
1373  template <typename somenumber>
1374  void SSOR_step (Vector<somenumber> &v,
1375  const Vector<somenumber> &b,
1376  const number om = 1.) const;
1378 
1382 
1389  const_iterator begin () const;
1390 
1394  iterator begin ();
1395 
1399  const_iterator end () const;
1400 
1404  iterator end ();
1405 
1415  const_iterator begin (const size_type r) const;
1416 
1420  iterator begin (const size_type r);
1421 
1431  const_iterator end (const size_type r) const;
1432 
1436  iterator end (const size_type r);
1438 
1442 
1454  template <class StreamType>
1455  void print (StreamType &out,
1456  const bool across = false,
1457  const bool diagonal_first = true) const;
1458 
1479  void print_formatted (std::ostream &out,
1480  const unsigned int precision = 3,
1481  const bool scientific = true,
1482  const unsigned int width = 0,
1483  const char *zero_string = " ",
1484  const double denominator = 1.) const;
1485 
1491  void print_pattern(std::ostream &out,
1492  const double threshold = 0.) const;
1493 
1504  void block_write (std::ostream &out) const;
1505 
1522  void block_read (std::istream &in);
1524 
1533  int, int,
1534  << "You are trying to access the matrix entry with index <"
1535  << arg1 << ',' << arg2
1536  << ">, but this entry does not exist in the sparsity pattern "
1537  "of this matrix."
1538  "\n\n"
1539  "The most common cause for this problem is that you used "
1540  "a method to build the sparsity pattern that did not "
1541  "(completely) take into account all of the entries you "
1542  "will later try to write into. An example would be "
1543  "building a sparsity pattern that does not include "
1544  "the entries you will write into due to constraints "
1545  "on degrees of freedom such as hanging nodes or periodic "
1546  "boundary conditions. In such cases, building the "
1547  "sparsity pattern will succeed, but you will get errors "
1548  "such as the current one at one point or other when "
1549  "trying to write into the entries of the matrix.");
1554  "When copying one sparse matrix into another, "
1555  "or when adding one sparse matrix to another, "
1556  "both matrices need to refer to the same "
1557  "sparsity pattern.");
1562  int, int,
1563  << "The iterators denote a range of " << arg1
1564  << " elements, but the given number of rows was " << arg2);
1570 
1571 protected:
1582  void prepare_add();
1583 
1588  void prepare_set();
1589 
1590 private:
1597 
1605  std::unique_ptr<number[]> val;
1606 
1613  std::size_t max_len;
1614 
1615  // make all other sparse matrices friends
1616  template <typename somenumber> friend class SparseMatrix;
1617  template <typename somenumber> friend class SparseLUDecomposition;
1618  template <typename> friend class SparseILU;
1619 
1623  template <typename> friend class BlockMatrixBase;
1624 
1628  template <typename,bool> friend class SparseMatrixIterators::Iterator;
1629  template <typename,bool> friend class SparseMatrixIterators::Accessor;
1630 };
1631 
1632 #ifndef DOXYGEN
1633 /*---------------------- Inline functions -----------------------------------*/
1634 
1635 
1636 
1637 template <typename number>
1638 inline
1640 {
1641  Assert (cols != nullptr, ExcNotInitialized());
1642  return cols->rows;
1643 }
1644 
1645 
1646 template <typename number>
1647 inline
1649 {
1650  Assert (cols != nullptr, ExcNotInitialized());
1651  return cols->cols;
1652 }
1653 
1654 
1655 // Inline the set() and add() functions, since they will be called frequently.
1656 template <typename number>
1657 inline
1658 void
1659 SparseMatrix<number>::set (const size_type i,
1660  const size_type j,
1661  const number value)
1662 {
1663  AssertIsFinite(value);
1664 
1665  const size_type index = cols->operator()(i, j);
1666 
1667  // it is allowed to set elements of the matrix that are not part of the
1668  // sparsity pattern, if the value to which we set it is zero
1669  if (index == SparsityPattern::invalid_entry)
1670  {
1671  Assert ((index != SparsityPattern::invalid_entry) ||
1672  (value == number()),
1673  ExcInvalidIndex(i, j));
1674  return;
1675  }
1676 
1677  val[index] = value;
1678 }
1679 
1680 
1681 
1682 template <typename number>
1683 template <typename number2>
1684 inline
1685 void
1686 SparseMatrix<number>::set (const std::vector<size_type> &indices,
1687  const FullMatrix<number2> &values,
1688  const bool elide_zero_values)
1689 {
1690  Assert (indices.size() == values.m(),
1691  ExcDimensionMismatch(indices.size(), values.m()));
1692  Assert (values.m() == values.n(), ExcNotQuadratic());
1693 
1694  for (size_type i=0; i<indices.size(); ++i)
1695  set (indices[i], indices.size(), &indices[0], &values(i,0),
1696  elide_zero_values);
1697 }
1698 
1699 
1700 
1701 template <typename number>
1702 template <typename number2>
1703 inline
1704 void
1705 SparseMatrix<number>::set (const std::vector<size_type> &row_indices,
1706  const std::vector<size_type> &col_indices,
1707  const FullMatrix<number2> &values,
1708  const bool elide_zero_values)
1709 {
1710  Assert (row_indices.size() == values.m(),
1711  ExcDimensionMismatch(row_indices.size(), values.m()));
1712  Assert (col_indices.size() == values.n(),
1713  ExcDimensionMismatch(col_indices.size(), values.n()));
1714 
1715  for (size_type i=0; i<row_indices.size(); ++i)
1716  set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1717  elide_zero_values);
1718 }
1719 
1720 
1721 
1722 template <typename number>
1723 template <typename number2>
1724 inline
1725 void
1726 SparseMatrix<number>::set (const size_type row,
1727  const std::vector<size_type> &col_indices,
1728  const std::vector<number2> &values,
1729  const bool elide_zero_values)
1730 {
1731  Assert (col_indices.size() == values.size(),
1732  ExcDimensionMismatch(col_indices.size(), values.size()));
1733 
1734  set (row, col_indices.size(), &col_indices[0], &values[0],
1735  elide_zero_values);
1736 }
1737 
1738 
1739 
1740 template <typename number>
1741 inline
1742 void
1743 SparseMatrix<number>::add (const size_type i,
1744  const size_type j,
1745  const number value)
1746 {
1747  AssertIsFinite(value);
1748 
1749  if (value == number())
1750  return;
1751 
1752  const size_type index = cols->operator()(i, j);
1753 
1754  // it is allowed to add elements to the matrix that are not part of the
1755  // sparsity pattern, if the value to which we set it is zero
1756  if (index == SparsityPattern::invalid_entry)
1757  {
1758  Assert ((index != SparsityPattern::invalid_entry) ||
1759  (value == number()),
1760  ExcInvalidIndex(i, j));
1761  return;
1762  }
1763 
1764  val[index] += value;
1765 }
1766 
1767 
1768 
1769 template <typename number>
1770 template <typename number2>
1771 inline
1772 void
1773 SparseMatrix<number>::add (const std::vector<size_type> &indices,
1774  const FullMatrix<number2> &values,
1775  const bool elide_zero_values)
1776 {
1777  Assert (indices.size() == values.m(),
1778  ExcDimensionMismatch(indices.size(), values.m()));
1779  Assert (values.m() == values.n(), ExcNotQuadratic());
1780 
1781  for (size_type i=0; i<indices.size(); ++i)
1782  add (indices[i], indices.size(), &indices[0], &values(i,0),
1783  elide_zero_values);
1784 }
1785 
1786 
1787 
1788 template <typename number>
1789 template <typename number2>
1790 inline
1791 void
1792 SparseMatrix<number>::add (const std::vector<size_type> &row_indices,
1793  const std::vector<size_type> &col_indices,
1794  const FullMatrix<number2> &values,
1795  const bool elide_zero_values)
1796 {
1797  Assert (row_indices.size() == values.m(),
1798  ExcDimensionMismatch(row_indices.size(), values.m()));
1799  Assert (col_indices.size() == values.n(),
1800  ExcDimensionMismatch(col_indices.size(), values.n()));
1801 
1802  for (size_type i=0; i<row_indices.size(); ++i)
1803  add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1804  elide_zero_values);
1805 }
1806 
1807 
1808 
1809 template <typename number>
1810 template <typename number2>
1811 inline
1812 void
1813 SparseMatrix<number>::add (const size_type row,
1814  const std::vector<size_type> &col_indices,
1815  const std::vector<number2> &values,
1816  const bool elide_zero_values)
1817 {
1818  Assert (col_indices.size() == values.size(),
1819  ExcDimensionMismatch(col_indices.size(), values.size()));
1820 
1821  add (row, col_indices.size(), &col_indices[0], &values[0],
1822  elide_zero_values);
1823 }
1824 
1825 
1826 
1827 template <typename number>
1828 inline
1830 SparseMatrix<number>::operator *= (const number factor)
1831 {
1832  Assert (cols != nullptr, ExcNotInitialized());
1833  Assert (val != nullptr, ExcNotInitialized());
1834 
1835  number *val_ptr = &val[0];
1836  const number *const end_ptr = &val[cols->n_nonzero_elements()];
1837 
1838  while (val_ptr != end_ptr)
1839  *val_ptr++ *= factor;
1840 
1841  return *this;
1842 }
1843 
1844 
1845 
1846 template <typename number>
1847 inline
1849 SparseMatrix<number>::operator /= (const number factor)
1850 {
1851  Assert (cols != nullptr, ExcNotInitialized());
1852  Assert (val != nullptr, ExcNotInitialized());
1853  Assert (factor != number(), ExcDivideByZero());
1854 
1855  const number factor_inv = number(1.) / factor;
1856 
1857  number *val_ptr = &val[0];
1858  const number *const end_ptr = &val[cols->n_nonzero_elements()];
1859 
1860  while (val_ptr != end_ptr)
1861  *val_ptr++ *= factor_inv;
1862 
1863  return *this;
1864 }
1865 
1866 
1867 
1868 template <typename number>
1869 inline
1870 number SparseMatrix<number>::operator () (const size_type i,
1871  const size_type j) const
1872 {
1873  Assert (cols != nullptr, ExcNotInitialized());
1874  Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
1875  ExcInvalidIndex(i,j));
1876  return val[cols->operator()(i,j)];
1877 }
1878 
1879 
1880 
1881 template <typename number>
1882 inline
1883 number SparseMatrix<number>::el (const size_type i,
1884  const size_type j) const
1885 {
1886  Assert (cols != nullptr, ExcNotInitialized());
1887  const size_type index = cols->operator()(i,j);
1888 
1889  if (index != SparsityPattern::invalid_entry)
1890  return val[index];
1891  else
1892  return 0;
1893 }
1894 
1895 
1896 
1897 template <typename number>
1898 inline
1899 number SparseMatrix<number>::diag_element (const size_type i) const
1900 {
1901  Assert (cols != nullptr, ExcNotInitialized());
1902  Assert (m() == n(), ExcNotQuadratic());
1903  AssertIndexRange(i, m());
1904 
1905  // Use that the first element in each row of a quadratic matrix is the main
1906  // diagonal
1907  return val[cols->rowstart[i]];
1908 }
1909 
1910 
1911 
1912 template <typename number>
1913 inline
1914 number &SparseMatrix<number>::diag_element (const size_type i)
1915 {
1916  Assert (cols != nullptr, ExcNotInitialized());
1917  Assert (m() == n(), ExcNotQuadratic());
1918  AssertIndexRange(i, m());
1919 
1920  // Use that the first element in each row of a quadratic matrix is the main
1921  // diagonal
1922  return val[cols->rowstart[i]];
1923 }
1924 
1925 
1926 
1927 template <typename number>
1928 template <typename ForwardIterator>
1929 void
1930 SparseMatrix<number>::copy_from (const ForwardIterator begin,
1931  const ForwardIterator end)
1932 {
1933  Assert (static_cast<size_type>(std::distance (begin, end)) == m(),
1934  ExcIteratorRange (std::distance (begin, end), m()));
1935 
1936  // for use in the inner loop, we define a typedef to the type of the inner
1937  // iterators
1938  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1939  size_type row=0;
1940  for (ForwardIterator i=begin; i!=end; ++i, ++row)
1941  {
1942  const inner_iterator end_of_row = i->end();
1943  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1944  // write entries
1945  set (row, j->first, j->second);
1946  };
1947 }
1948 
1949 
1950 //---------------------------------------------------------------------------
1951 
1952 
1953 namespace SparseMatrixIterators
1954 {
1955  template <typename number>
1956  inline
1958  Accessor (const MatrixType *matrix,
1959  const std::size_t index_within_matrix)
1960  :
1961  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1962  index_within_matrix),
1963  matrix (matrix)
1964  {}
1965 
1966 
1967 
1968  template <typename number>
1969  inline
1970  Accessor<number,true>::
1971  Accessor (const MatrixType *matrix)
1972  :
1973  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1974  matrix (matrix)
1975  {}
1976 
1977 
1978 
1979  template <typename number>
1980  inline
1981  Accessor<number,true>::
1983  :
1984  SparsityPatternIterators::Accessor (a),
1985  matrix (&a.get_matrix())
1986  {}
1987 
1988 
1989 
1990  template <typename number>
1991  inline
1992  number
1993  Accessor<number, true>::value () const
1994  {
1995  AssertIndexRange(index_within_sparsity, matrix->n_nonzero_elements());
1996  return matrix->val[index_within_sparsity];
1997  }
1998 
1999 
2000 
2001  template <typename number>
2002  inline
2003  const typename Accessor<number, true>::MatrixType &
2004  Accessor<number, true>::get_matrix () const
2005  {
2006  return *matrix;
2007  }
2008 
2009 
2010 
2011  template <typename number>
2012  inline
2013  Accessor<number, false>::Reference::Reference (
2014  const Accessor *accessor,
2015  const bool)
2016  :
2017  accessor (accessor)
2018  {}
2019 
2020 
2021  template <typename number>
2022  inline
2023  Accessor<number, false>::Reference::operator number() const
2024  {
2025  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2026  return accessor->matrix->val[accessor->index_within_sparsity];
2027  }
2028 
2029 
2030 
2031  template <typename number>
2032  inline
2033  const typename Accessor<number, false>::Reference &
2034  Accessor<number, false>::Reference::operator = (const number n) const
2035  {
2036  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2037  accessor->matrix->val[accessor->index_within_sparsity] = n;
2038  return *this;
2039  }
2040 
2041 
2042 
2043  template <typename number>
2044  inline
2045  const typename Accessor<number, false>::Reference &
2046  Accessor<number, false>::Reference::operator += (const number n) const
2047  {
2048  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2049  accessor->matrix->val[accessor->index_within_sparsity] += n;
2050  return *this;
2051  }
2052 
2053 
2054 
2055  template <typename number>
2056  inline
2057  const typename Accessor<number, false>::Reference &
2058  Accessor<number, false>::Reference::operator -= (const number n) const
2059  {
2060  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2061  accessor->matrix->val[accessor->index_within_sparsity] -= n;
2062  return *this;
2063  }
2064 
2065 
2066 
2067  template <typename number>
2068  inline
2069  const typename Accessor<number, false>::Reference &
2070  Accessor<number, false>::Reference::operator *= (const number n) const
2071  {
2072  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2073  accessor->matrix->val[accessor->index_within_sparsity] *= n;
2074  return *this;
2075  }
2076 
2077 
2078 
2079  template <typename number>
2080  inline
2081  const typename Accessor<number, false>::Reference &
2082  Accessor<number, false>::Reference::operator /= (const number n) const
2083  {
2084  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2085  accessor->matrix->val[accessor->index_within_sparsity] /= n;
2086  return *this;
2087  }
2088 
2089 
2090 
2091  template <typename number>
2092  inline
2093  Accessor<number,false>::
2094  Accessor (MatrixType *matrix,
2095  const std::size_t index)
2096  :
2097  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2098  index),
2099  matrix (matrix)
2100  {}
2101 
2102 
2103 
2104  template <typename number>
2105  inline
2106  Accessor<number,false>::
2107  Accessor (MatrixType *matrix)
2108  :
2109  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
2110  matrix (matrix)
2111  {}
2112 
2113 
2114 
2115  template <typename number>
2116  inline
2117  typename Accessor<number, false>::Reference
2118  Accessor<number, false>::value() const
2119  {
2120  return Reference(this,true);
2121  }
2122 
2123 
2124 
2125 
2126  template <typename number>
2127  inline
2128  typename Accessor<number, false>::MatrixType &
2129  Accessor<number, false>::get_matrix () const
2130  {
2131  return *matrix;
2132  }
2133 
2134 
2135 
2136  template <typename number, bool Constness>
2137  inline
2138  Iterator<number, Constness>::
2139  Iterator (MatrixType *matrix,
2140  const std::size_t index)
2141  :
2142  accessor(matrix, index)
2143  {}
2144 
2145 
2146 
2147  template <typename number, bool Constness>
2148  inline
2149  Iterator<number, Constness>::
2150  Iterator (MatrixType *matrix)
2151  :
2152  accessor(matrix)
2153  {}
2154 
2155 
2156 
2157  template <typename number, bool Constness>
2158  inline
2159  Iterator<number, Constness>::
2161  :
2162  accessor(*i)
2163  {}
2164 
2165 
2166 
2167  template <typename number, bool Constness>
2168  inline
2169  Iterator<number, Constness> &
2170  Iterator<number,Constness>::operator++ ()
2171  {
2172  accessor.advance ();
2173  return *this;
2174  }
2175 
2176 
2177  template <typename number, bool Constness>
2178  inline
2179  Iterator<number,Constness>
2180  Iterator<number,Constness>::operator++ (int)
2181  {
2182  const Iterator iter = *this;
2183  accessor.advance ();
2184  return iter;
2185  }
2186 
2187 
2188  template <typename number, bool Constness>
2189  inline
2190  const Accessor<number,Constness> &
2191  Iterator<number,Constness>::operator* () const
2192  {
2193  return accessor;
2194  }
2195 
2196 
2197  template <typename number, bool Constness>
2198  inline
2199  const Accessor<number,Constness> *
2200  Iterator<number,Constness>::operator-> () const
2201  {
2202  return &accessor;
2203  }
2204 
2205 
2206  template <typename number, bool Constness>
2207  inline
2208  bool
2209  Iterator<number,Constness>::
2210  operator == (const Iterator &other) const
2211  {
2212  return (accessor == other.accessor);
2213  }
2214 
2215 
2216  template <typename number, bool Constness>
2217  inline
2218  bool
2219  Iterator<number,Constness>::
2220  operator != (const Iterator &other) const
2221  {
2222  return ! (*this == other);
2223  }
2224 
2225 
2226  template <typename number, bool Constness>
2227  inline
2228  bool
2229  Iterator<number,Constness>::
2230  operator < (const Iterator &other) const
2231  {
2232  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2233  ExcInternalError());
2234 
2235  return (accessor < other.accessor);
2236  }
2237 
2238 
2239  template <typename number, bool Constness>
2240  inline
2241  bool
2242  Iterator<number,Constness>::
2243  operator > (const Iterator &other) const
2244  {
2245  return (other < *this);
2246  }
2247 
2248 
2249  template <typename number, bool Constness>
2250  inline
2251  int
2252  Iterator<number,Constness>::
2253  operator - (const Iterator &other) const
2254  {
2255  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2256  ExcInternalError());
2257 
2258  return (*this)->index_within_sparsity - other->index_within_sparsity;
2259  }
2260 
2261 
2262 
2263  template <typename number, bool Constness>
2264  inline
2265  Iterator<number,Constness>
2266  Iterator<number,Constness>::
2267  operator + (const size_type n) const
2268  {
2269  Iterator x = *this;
2270  for (size_type i=0; i<n; ++i)
2271  ++x;
2272 
2273  return x;
2274  }
2275 
2276 }
2277 
2278 
2279 
2280 template <typename number>
2281 inline
2284 {
2285  return const_iterator(this, 0);
2286 }
2287 
2288 
2289 template <typename number>
2290 inline
2293 {
2294  return const_iterator(this);
2295 }
2296 
2297 
2298 template <typename number>
2299 inline
2302 {
2303  return iterator (this, 0);
2304 }
2305 
2306 
2307 template <typename number>
2308 inline
2311 {
2312  return iterator(this, cols->rowstart[cols->rows]);
2313 }
2314 
2315 
2316 template <typename number>
2317 inline
2319 SparseMatrix<number>::begin (const size_type r) const
2320 {
2321  Assert (r<m(), ExcIndexRange(r,0,m()));
2322 
2323  return const_iterator(this, cols->rowstart[r]);
2324 }
2325 
2326 
2327 
2328 template <typename number>
2329 inline
2331 SparseMatrix<number>::end (const size_type r) const
2332 {
2333  Assert (r<m(), ExcIndexRange(r,0,m()));
2334 
2335  return const_iterator(this, cols->rowstart[r+1]);
2336 }
2337 
2338 
2339 
2340 template <typename number>
2341 inline
2343 SparseMatrix<number>::begin (const size_type r)
2344 {
2345  Assert (r<m(), ExcIndexRange(r,0,m()));
2346 
2347  return iterator(this, cols->rowstart[r]);
2348 }
2349 
2350 
2351 
2352 template <typename number>
2353 inline
2355 SparseMatrix<number>::end (const size_type r)
2356 {
2357  Assert (r<m(), ExcIndexRange(r,0,m()));
2358 
2359  return iterator(this, cols->rowstart[r+1]);
2360 }
2361 
2362 
2363 
2364 template <typename number>
2365 template <class StreamType>
2366 inline
2367 void SparseMatrix<number>::print (StreamType &out,
2368  const bool across,
2369  const bool diagonal_first) const
2370 {
2371  Assert (cols != nullptr, ExcNotInitialized());
2372  Assert (val != nullptr, ExcNotInitialized());
2373 
2374  bool hanging_diagonal = false;
2375  number diagonal = number();
2376 
2377  for (size_type i=0; i<cols->rows; ++i)
2378  {
2379  for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
2380  {
2381  if (!diagonal_first && i == cols->colnums[j])
2382  {
2383  diagonal = val[j];
2384  hanging_diagonal = true;
2385  }
2386  else
2387  {
2388  if (hanging_diagonal && cols->colnums[j]>i)
2389  {
2390  if (across)
2391  out << ' ' << i << ',' << i << ':' << diagonal;
2392  else
2393  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2394  hanging_diagonal = false;
2395  }
2396  if (across)
2397  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2398  else
2399  out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << std::endl;
2400  }
2401  }
2402  if (hanging_diagonal)
2403  {
2404  if (across)
2405  out << ' ' << i << ',' << i << ':' << diagonal;
2406  else
2407  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2408  hanging_diagonal = false;
2409  }
2410  }
2411  if (across)
2412  out << std::endl;
2413 }
2414 
2415 
2416 template <typename number>
2417 inline
2418 void
2420 prepare_add()
2421 {
2422  //nothing to do here
2423 }
2424 
2425 
2426 
2427 template <typename number>
2428 inline
2429 void
2431 prepare_set()
2432 {
2433  //nothing to do here
2434 }
2435 
2436 #endif // DOXYGEN
2437 
2438 
2439 /*---------------------------- sparse_matrix.h ---------------------------*/
2440 
2441 DEAL_II_NAMESPACE_CLOSE
2442 
2443 #endif
2444 /*---------------------------- 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:605
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
Accessor< number, Constness >::MatrixType MatrixType
number operator()(const size_type i, const size_type j) const
Iterator operator+(const size_type n) const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
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:1199
void block_read(std::istream &in)
SparseMatrixIterators::Iterator< number, false > iterator
real_type frobenius_norm() const
static::ExceptionBase & ExcNotInitialized()
std::size_t n_nonzero_elements() const
real_type linfty_norm() const
static::ExceptionBase & ExcDifferentSparsityPatterns()
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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
static::ExceptionBase & ExcDivideByZero()
void SOR(Vector< somenumber > &v, const number om=1.) const
const_iterator begin() const
Definition: mpi.h:57
number el(const size_type i, const size_type j) const
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
unsigned int global_dof_index
Definition: types.h:88
const SparseMatrix< number > & get_matrix() const
#define Assert(cond, exc)
Definition: exceptions.h:337
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:59
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:582
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
const Accessor< number, Constness > & operator*() const
number value_type
real_type l1_norm() const
#define DeclException0(Exception0)
Definition: exceptions.h:570
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
static::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
void prepare_set()
const Accessor< number, Constness > * operator->() const
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 Accessor< number, Constness > & value_type
const_iterator end() const
Accessor< number, Constness > accessor
bool operator<(const Iterator &) const
int operator-(const Iterator &p) const
SparseMatrixIterators::Iterator< number, true > const_iterator
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
static::ExceptionBase & ExcNotQuadratic()
void vmult_add(OutVector &dst, const InVector &src) const
types::global_dof_index size_type
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
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
virtual ~SparseMatrix()
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)
bool operator==(const Iterator &) const
std::size_t n_actually_nonzero_elements(const double threshold=0.) const
virtual void clear()
#define AssertIsFinite(number)
Definition: exceptions.h:1211
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