Reference documentation for deal.II version Git 2618e0f 2017-11-23 17:25:26 +0100
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_operation.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);
1569  "You are attempting an operation on two matrices that "
1570  "are the same object, but the operation requires that the "
1571  "two objects are in fact different.");
1573 
1574 protected:
1585  void prepare_add();
1586 
1591  void prepare_set();
1592 
1593 private:
1600 
1608  std::unique_ptr<number[]> val;
1609 
1616  std::size_t max_len;
1617 
1618  // make all other sparse matrices friends
1619  template <typename somenumber> friend class SparseMatrix;
1620  template <typename somenumber> friend class SparseLUDecomposition;
1621  template <typename> friend class SparseILU;
1622 
1626  template <typename> friend class BlockMatrixBase;
1627 
1631  template <typename,bool> friend class SparseMatrixIterators::Iterator;
1632  template <typename,bool> friend class SparseMatrixIterators::Accessor;
1633 };
1634 
1635 #ifndef DOXYGEN
1636 /*---------------------- Inline functions -----------------------------------*/
1637 
1638 
1639 
1640 template <typename number>
1641 inline
1643 {
1644  Assert (cols != nullptr, ExcNotInitialized());
1645  return cols->rows;
1646 }
1647 
1648 
1649 template <typename number>
1650 inline
1652 {
1653  Assert (cols != nullptr, ExcNotInitialized());
1654  return cols->cols;
1655 }
1656 
1657 
1658 // Inline the set() and add() functions, since they will be called frequently.
1659 template <typename number>
1660 inline
1661 void
1662 SparseMatrix<number>::set (const size_type i,
1663  const size_type j,
1664  const number value)
1665 {
1666  AssertIsFinite(value);
1667 
1668  const size_type index = cols->operator()(i, j);
1669 
1670  // it is allowed to set elements of the matrix that are not part of the
1671  // sparsity pattern, if the value to which we set it is zero
1672  if (index == SparsityPattern::invalid_entry)
1673  {
1674  Assert ((index != SparsityPattern::invalid_entry) ||
1675  (value == number()),
1676  ExcInvalidIndex(i, j));
1677  return;
1678  }
1679 
1680  val[index] = value;
1681 }
1682 
1683 
1684 
1685 template <typename number>
1686 template <typename number2>
1687 inline
1688 void
1689 SparseMatrix<number>::set (const std::vector<size_type> &indices,
1690  const FullMatrix<number2> &values,
1691  const bool elide_zero_values)
1692 {
1693  Assert (indices.size() == values.m(),
1694  ExcDimensionMismatch(indices.size(), values.m()));
1695  Assert (values.m() == values.n(), ExcNotQuadratic());
1696 
1697  for (size_type i=0; i<indices.size(); ++i)
1698  set (indices[i], indices.size(), indices.data(), &values(i,0),
1699  elide_zero_values);
1700 }
1701 
1702 
1703 
1704 template <typename number>
1705 template <typename number2>
1706 inline
1707 void
1708 SparseMatrix<number>::set (const std::vector<size_type> &row_indices,
1709  const std::vector<size_type> &col_indices,
1710  const FullMatrix<number2> &values,
1711  const bool elide_zero_values)
1712 {
1713  Assert (row_indices.size() == values.m(),
1714  ExcDimensionMismatch(row_indices.size(), values.m()));
1715  Assert (col_indices.size() == values.n(),
1716  ExcDimensionMismatch(col_indices.size(), values.n()));
1717 
1718  for (size_type i=0; i<row_indices.size(); ++i)
1719  set (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
1720  elide_zero_values);
1721 }
1722 
1723 
1724 
1725 template <typename number>
1726 template <typename number2>
1727 inline
1728 void
1729 SparseMatrix<number>::set (const size_type row,
1730  const std::vector<size_type> &col_indices,
1731  const std::vector<number2> &values,
1732  const bool elide_zero_values)
1733 {
1734  Assert (col_indices.size() == values.size(),
1735  ExcDimensionMismatch(col_indices.size(), values.size()));
1736 
1737  set (row, col_indices.size(), col_indices.data(), values.data(),
1738  elide_zero_values);
1739 }
1740 
1741 
1742 
1743 template <typename number>
1744 inline
1745 void
1746 SparseMatrix<number>::add (const size_type i,
1747  const size_type j,
1748  const number value)
1749 {
1750  AssertIsFinite(value);
1751 
1752  if (value == number())
1753  return;
1754 
1755  const size_type index = cols->operator()(i, j);
1756 
1757  // it is allowed to add elements to the matrix that are not part of the
1758  // sparsity pattern, if the value to which we set it is zero
1759  if (index == SparsityPattern::invalid_entry)
1760  {
1761  Assert ((index != SparsityPattern::invalid_entry) ||
1762  (value == number()),
1763  ExcInvalidIndex(i, j));
1764  return;
1765  }
1766 
1767  val[index] += value;
1768 }
1769 
1770 
1771 
1772 template <typename number>
1773 template <typename number2>
1774 inline
1775 void
1776 SparseMatrix<number>::add (const std::vector<size_type> &indices,
1777  const FullMatrix<number2> &values,
1778  const bool elide_zero_values)
1779 {
1780  Assert (indices.size() == values.m(),
1781  ExcDimensionMismatch(indices.size(), values.m()));
1782  Assert (values.m() == values.n(), ExcNotQuadratic());
1783 
1784  for (size_type i=0; i<indices.size(); ++i)
1785  add (indices[i], indices.size(), indices.data(), &values(i,0),
1786  elide_zero_values);
1787 }
1788 
1789 
1790 
1791 template <typename number>
1792 template <typename number2>
1793 inline
1794 void
1795 SparseMatrix<number>::add (const std::vector<size_type> &row_indices,
1796  const std::vector<size_type> &col_indices,
1797  const FullMatrix<number2> &values,
1798  const bool elide_zero_values)
1799 {
1800  Assert (row_indices.size() == values.m(),
1801  ExcDimensionMismatch(row_indices.size(), values.m()));
1802  Assert (col_indices.size() == values.n(),
1803  ExcDimensionMismatch(col_indices.size(), values.n()));
1804 
1805  for (size_type i=0; i<row_indices.size(); ++i)
1806  add (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
1807  elide_zero_values);
1808 }
1809 
1810 
1811 
1812 template <typename number>
1813 template <typename number2>
1814 inline
1815 void
1816 SparseMatrix<number>::add (const size_type row,
1817  const std::vector<size_type> &col_indices,
1818  const std::vector<number2> &values,
1819  const bool elide_zero_values)
1820 {
1821  Assert (col_indices.size() == values.size(),
1822  ExcDimensionMismatch(col_indices.size(), values.size()));
1823 
1824  add (row, col_indices.size(), col_indices.data(), values.data(),
1825  elide_zero_values);
1826 }
1827 
1828 
1829 
1830 template <typename number>
1831 inline
1833 SparseMatrix<number>::operator *= (const number factor)
1834 {
1835  Assert (cols != nullptr, ExcNotInitialized());
1836  Assert (val != nullptr, ExcNotInitialized());
1837 
1838  number *val_ptr = val.get();
1839  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1840 
1841  while (val_ptr != end_ptr)
1842  *val_ptr++ *= factor;
1843 
1844  return *this;
1845 }
1846 
1847 
1848 
1849 template <typename number>
1850 inline
1852 SparseMatrix<number>::operator /= (const number factor)
1853 {
1854  Assert (cols != nullptr, ExcNotInitialized());
1855  Assert (val != nullptr, ExcNotInitialized());
1856  Assert (factor != number(), ExcDivideByZero());
1857 
1858  const number factor_inv = number(1.) / factor;
1859 
1860  number *val_ptr = val.get();
1861  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1862 
1863  while (val_ptr != end_ptr)
1864  *val_ptr++ *= factor_inv;
1865 
1866  return *this;
1867 }
1868 
1869 
1870 
1871 template <typename number>
1872 inline
1873 number SparseMatrix<number>::operator () (const size_type i,
1874  const size_type j) const
1875 {
1876  Assert (cols != nullptr, ExcNotInitialized());
1877  Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
1878  ExcInvalidIndex(i,j));
1879  return val[cols->operator()(i,j)];
1880 }
1881 
1882 
1883 
1884 template <typename number>
1885 inline
1886 number SparseMatrix<number>::el (const size_type i,
1887  const size_type j) const
1888 {
1889  Assert (cols != nullptr, ExcNotInitialized());
1890  const size_type index = cols->operator()(i,j);
1891 
1892  if (index != SparsityPattern::invalid_entry)
1893  return val[index];
1894  else
1895  return 0;
1896 }
1897 
1898 
1899 
1900 template <typename number>
1901 inline
1902 number SparseMatrix<number>::diag_element (const size_type i) const
1903 {
1904  Assert (cols != nullptr, ExcNotInitialized());
1905  Assert (m() == n(), ExcNotQuadratic());
1906  AssertIndexRange(i, m());
1907 
1908  // Use that the first element in each row of a quadratic matrix is the main
1909  // diagonal
1910  return val[cols->rowstart[i]];
1911 }
1912 
1913 
1914 
1915 template <typename number>
1916 inline
1917 number &SparseMatrix<number>::diag_element (const size_type i)
1918 {
1919  Assert (cols != nullptr, ExcNotInitialized());
1920  Assert (m() == n(), ExcNotQuadratic());
1921  AssertIndexRange(i, m());
1922 
1923  // Use that the first element in each row of a quadratic matrix is the main
1924  // diagonal
1925  return val[cols->rowstart[i]];
1926 }
1927 
1928 
1929 
1930 template <typename number>
1931 template <typename ForwardIterator>
1932 void
1933 SparseMatrix<number>::copy_from (const ForwardIterator begin,
1934  const ForwardIterator end)
1935 {
1936  Assert (static_cast<size_type>(std::distance (begin, end)) == m(),
1937  ExcIteratorRange (std::distance (begin, end), m()));
1938 
1939  // for use in the inner loop, we define a typedef to the type of the inner
1940  // iterators
1941  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1942  size_type row=0;
1943  for (ForwardIterator i=begin; i!=end; ++i, ++row)
1944  {
1945  const inner_iterator end_of_row = i->end();
1946  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1947  // write entries
1948  set (row, j->first, j->second);
1949  };
1950 }
1951 
1952 
1953 //---------------------------------------------------------------------------
1954 
1955 
1956 namespace SparseMatrixIterators
1957 {
1958  template <typename number>
1959  inline
1961  Accessor (const MatrixType *matrix,
1962  const std::size_t index_within_matrix)
1963  :
1964  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1965  index_within_matrix),
1966  matrix (matrix)
1967  {}
1968 
1969 
1970 
1971  template <typename number>
1972  inline
1973  Accessor<number,true>::
1974  Accessor (const MatrixType *matrix)
1975  :
1976  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1977  matrix (matrix)
1978  {}
1979 
1980 
1981 
1982  template <typename number>
1983  inline
1984  Accessor<number,true>::
1986  :
1987  SparsityPatternIterators::Accessor (a),
1988  matrix (&a.get_matrix())
1989  {}
1990 
1991 
1992 
1993  template <typename number>
1994  inline
1995  number
1996  Accessor<number, true>::value () const
1997  {
1998  AssertIndexRange(index_within_sparsity, matrix->n_nonzero_elements());
1999  return matrix->val[index_within_sparsity];
2000  }
2001 
2002 
2003 
2004  template <typename number>
2005  inline
2006  const typename Accessor<number, true>::MatrixType &
2007  Accessor<number, true>::get_matrix () const
2008  {
2009  return *matrix;
2010  }
2011 
2012 
2013 
2014  template <typename number>
2015  inline
2016  Accessor<number, false>::Reference::Reference (
2017  const Accessor *accessor,
2018  const bool)
2019  :
2020  accessor (accessor)
2021  {}
2022 
2023 
2024  template <typename number>
2025  inline
2026  Accessor<number, false>::Reference::operator number() const
2027  {
2028  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2029  return accessor->matrix->val[accessor->index_within_sparsity];
2030  }
2031 
2032 
2033 
2034  template <typename number>
2035  inline
2036  const typename Accessor<number, false>::Reference &
2037  Accessor<number, false>::Reference::operator = (const number n) const
2038  {
2039  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2040  accessor->matrix->val[accessor->index_within_sparsity] = n;
2041  return *this;
2042  }
2043 
2044 
2045 
2046  template <typename number>
2047  inline
2048  const typename Accessor<number, false>::Reference &
2049  Accessor<number, false>::Reference::operator += (const number n) const
2050  {
2051  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2052  accessor->matrix->val[accessor->index_within_sparsity] += n;
2053  return *this;
2054  }
2055 
2056 
2057 
2058  template <typename number>
2059  inline
2060  const typename Accessor<number, false>::Reference &
2061  Accessor<number, false>::Reference::operator -= (const number n) const
2062  {
2063  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2064  accessor->matrix->val[accessor->index_within_sparsity] -= n;
2065  return *this;
2066  }
2067 
2068 
2069 
2070  template <typename number>
2071  inline
2072  const typename Accessor<number, false>::Reference &
2073  Accessor<number, false>::Reference::operator *= (const number n) const
2074  {
2075  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2076  accessor->matrix->val[accessor->index_within_sparsity] *= n;
2077  return *this;
2078  }
2079 
2080 
2081 
2082  template <typename number>
2083  inline
2084  const typename Accessor<number, false>::Reference &
2085  Accessor<number, false>::Reference::operator /= (const number n) const
2086  {
2087  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2088  accessor->matrix->val[accessor->index_within_sparsity] /= n;
2089  return *this;
2090  }
2091 
2092 
2093 
2094  template <typename number>
2095  inline
2096  Accessor<number,false>::
2097  Accessor (MatrixType *matrix,
2098  const std::size_t index)
2099  :
2100  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2101  index),
2102  matrix (matrix)
2103  {}
2104 
2105 
2106 
2107  template <typename number>
2108  inline
2109  Accessor<number,false>::
2110  Accessor (MatrixType *matrix)
2111  :
2112  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
2113  matrix (matrix)
2114  {}
2115 
2116 
2117 
2118  template <typename number>
2119  inline
2120  typename Accessor<number, false>::Reference
2121  Accessor<number, false>::value() const
2122  {
2123  return Reference(this,true);
2124  }
2125 
2126 
2127 
2128 
2129  template <typename number>
2130  inline
2131  typename Accessor<number, false>::MatrixType &
2132  Accessor<number, false>::get_matrix () const
2133  {
2134  return *matrix;
2135  }
2136 
2137 
2138 
2139  template <typename number, bool Constness>
2140  inline
2141  Iterator<number, Constness>::
2142  Iterator (MatrixType *matrix,
2143  const std::size_t index)
2144  :
2145  accessor(matrix, index)
2146  {}
2147 
2148 
2149 
2150  template <typename number, bool Constness>
2151  inline
2152  Iterator<number, Constness>::
2153  Iterator (MatrixType *matrix)
2154  :
2155  accessor(matrix)
2156  {}
2157 
2158 
2159 
2160  template <typename number, bool Constness>
2161  inline
2162  Iterator<number, Constness>::
2164  :
2165  accessor(*i)
2166  {}
2167 
2168 
2169 
2170  template <typename number, bool Constness>
2171  inline
2172  Iterator<number, Constness> &
2173  Iterator<number,Constness>::operator++ ()
2174  {
2175  accessor.advance ();
2176  return *this;
2177  }
2178 
2179 
2180  template <typename number, bool Constness>
2181  inline
2182  Iterator<number,Constness>
2183  Iterator<number,Constness>::operator++ (int)
2184  {
2185  const Iterator iter = *this;
2186  accessor.advance ();
2187  return iter;
2188  }
2189 
2190 
2191  template <typename number, bool Constness>
2192  inline
2193  const Accessor<number,Constness> &
2194  Iterator<number,Constness>::operator* () const
2195  {
2196  return accessor;
2197  }
2198 
2199 
2200  template <typename number, bool Constness>
2201  inline
2202  const Accessor<number,Constness> *
2203  Iterator<number,Constness>::operator-> () const
2204  {
2205  return &accessor;
2206  }
2207 
2208 
2209  template <typename number, bool Constness>
2210  inline
2211  bool
2212  Iterator<number,Constness>::
2213  operator == (const Iterator &other) const
2214  {
2215  return (accessor == other.accessor);
2216  }
2217 
2218 
2219  template <typename number, bool Constness>
2220  inline
2221  bool
2222  Iterator<number,Constness>::
2223  operator != (const Iterator &other) const
2224  {
2225  return ! (*this == other);
2226  }
2227 
2228 
2229  template <typename number, bool Constness>
2230  inline
2231  bool
2232  Iterator<number,Constness>::
2233  operator < (const Iterator &other) const
2234  {
2235  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2236  ExcInternalError());
2237 
2238  return (accessor < other.accessor);
2239  }
2240 
2241 
2242  template <typename number, bool Constness>
2243  inline
2244  bool
2245  Iterator<number,Constness>::
2246  operator > (const Iterator &other) const
2247  {
2248  return (other < *this);
2249  }
2250 
2251 
2252  template <typename number, bool Constness>
2253  inline
2254  int
2255  Iterator<number,Constness>::
2256  operator - (const Iterator &other) const
2257  {
2258  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2259  ExcInternalError());
2260 
2261  return (*this)->index_within_sparsity - other->index_within_sparsity;
2262  }
2263 
2264 
2265 
2266  template <typename number, bool Constness>
2267  inline
2268  Iterator<number,Constness>
2269  Iterator<number,Constness>::
2270  operator + (const size_type n) const
2271  {
2272  Iterator x = *this;
2273  for (size_type i=0; i<n; ++i)
2274  ++x;
2275 
2276  return x;
2277  }
2278 
2279 }
2280 
2281 
2282 
2283 template <typename number>
2284 inline
2287 {
2288  return const_iterator(this, 0);
2289 }
2290 
2291 
2292 template <typename number>
2293 inline
2296 {
2297  return const_iterator(this);
2298 }
2299 
2300 
2301 template <typename number>
2302 inline
2305 {
2306  return iterator (this, 0);
2307 }
2308 
2309 
2310 template <typename number>
2311 inline
2314 {
2315  return iterator(this, cols->rowstart[cols->rows]);
2316 }
2317 
2318 
2319 template <typename number>
2320 inline
2322 SparseMatrix<number>::begin (const size_type r) const
2323 {
2324  Assert (r<m(), ExcIndexRange(r,0,m()));
2325 
2326  return const_iterator(this, cols->rowstart[r]);
2327 }
2328 
2329 
2330 
2331 template <typename number>
2332 inline
2334 SparseMatrix<number>::end (const size_type r) const
2335 {
2336  Assert (r<m(), ExcIndexRange(r,0,m()));
2337 
2338  return const_iterator(this, cols->rowstart[r+1]);
2339 }
2340 
2341 
2342 
2343 template <typename number>
2344 inline
2346 SparseMatrix<number>::begin (const size_type r)
2347 {
2348  Assert (r<m(), ExcIndexRange(r,0,m()));
2349 
2350  return iterator(this, cols->rowstart[r]);
2351 }
2352 
2353 
2354 
2355 template <typename number>
2356 inline
2358 SparseMatrix<number>::end (const size_type r)
2359 {
2360  Assert (r<m(), ExcIndexRange(r,0,m()));
2361 
2362  return iterator(this, cols->rowstart[r+1]);
2363 }
2364 
2365 
2366 
2367 template <typename number>
2368 template <class StreamType>
2369 inline
2370 void SparseMatrix<number>::print (StreamType &out,
2371  const bool across,
2372  const bool diagonal_first) const
2373 {
2374  Assert (cols != nullptr, ExcNotInitialized());
2375  Assert (val != nullptr, ExcNotInitialized());
2376 
2377  bool hanging_diagonal = false;
2378  number diagonal = number();
2379 
2380  for (size_type i=0; i<cols->rows; ++i)
2381  {
2382  for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
2383  {
2384  if (!diagonal_first && i == cols->colnums[j])
2385  {
2386  diagonal = val[j];
2387  hanging_diagonal = true;
2388  }
2389  else
2390  {
2391  if (hanging_diagonal && cols->colnums[j]>i)
2392  {
2393  if (across)
2394  out << ' ' << i << ',' << i << ':' << diagonal;
2395  else
2396  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2397  hanging_diagonal = false;
2398  }
2399  if (across)
2400  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2401  else
2402  out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << std::endl;
2403  }
2404  }
2405  if (hanging_diagonal)
2406  {
2407  if (across)
2408  out << ' ' << i << ',' << i << ':' << diagonal;
2409  else
2410  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2411  hanging_diagonal = false;
2412  }
2413  }
2414  if (across)
2415  out << std::endl;
2416 }
2417 
2418 
2419 template <typename number>
2420 inline
2421 void
2423 prepare_add()
2424 {
2425  //nothing to do here
2426 }
2427 
2428 
2429 
2430 template <typename number>
2431 inline
2432 void
2434 prepare_set()
2435 {
2436  //nothing to do here
2437 }
2438 
2439 #endif // DOXYGEN
2440 
2441 
2442 /*---------------------------- sparse_matrix.h ---------------------------*/
2443 
2444 DEAL_II_NAMESPACE_CLOSE
2445 
2446 #endif
2447 /*---------------------------- 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
Contents is actually a matrix.
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:1207
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:59
Matrix is diagonal.
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
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:1219
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