Reference documentation for deal.II version 8.4.1
sparse_matrix.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2016 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 DEAL_II_NAMESPACE_OPEN
29 
30 template <typename number> class Vector;
31 template <typename number> class FullMatrix;
32 template <typename Matrix> class BlockMatrixBase;
33 template <typename number> class SparseILU;
34 
35 #ifdef DEAL_II_WITH_TRILINOS
36 namespace TrilinosWrappers
37 {
38  class SparseMatrix;
39 }
40 #endif
41 
52 {
57 
58  // forward declaration
59  template <typename number, bool Constness>
60  class Iterator;
61 
72  template <typename number, bool Constness>
74  {
75  public:
79  number value() const;
80 
84  number &value();
85 
90  const SparseMatrix<number> &get_matrix () const;
91  };
92 
93 
94 
101  template <typename number>
102  class Accessor<number,true> : public SparsityPatternIterators::Accessor
103  {
104  public:
110 
114  Accessor (MatrixType *matrix,
115  const std::size_t index_within_matrix);
116 
120  Accessor (MatrixType *matrix);
121 
126 
130  number value() const;
131 
136  const MatrixType &get_matrix () const;
137 
138  private:
142  MatrixType *matrix;
143 
148 
152  template <typename, bool>
153  friend class Iterator;
154  };
155 
156 
163  template <typename number>
164  class Accessor<number,false> : public SparsityPatternIterators::Accessor
165  {
166  private:
191  class Reference
192  {
193  public:
198  Reference (const Accessor *accessor,
199  const bool dummy);
200 
204  operator number () const;
205 
209  const Reference &operator = (const number n) const;
210 
214  const Reference &operator += (const number n) const;
215 
219  const Reference &operator -= (const number n) const;
220 
224  const Reference &operator *= (const number n) const;
225 
229  const Reference &operator /= (const number n) const;
230 
231  private:
237  };
238 
239  public:
245 
249  Accessor (MatrixType *matrix,
250  const std::size_t index);
251 
255  Accessor (MatrixType *matrix);
256 
260  Reference value() const;
261 
266  MatrixType &get_matrix () const;
267 
268  private:
272  MatrixType *matrix;
273 
278 
282  template <typename, bool>
283  friend class Iterator;
284 
289  template <typename> friend class ::SparseMatrix;
290  };
291 
292 
293 
323  template <typename number, bool Constness>
324  class Iterator
325  {
326  public:
330  typedef
333 
338  typedef
340 
345  Iterator (MatrixType *matrix,
346  const std::size_t index_within_matrix);
347 
351  Iterator (MatrixType *matrix);
352 
358 
363 
367  Iterator operator++ (int);
368 
372  const Accessor<number,Constness> &operator* () const;
373 
378 
382  bool operator == (const Iterator &) const;
383 
387  bool operator != (const Iterator &) const;
388 
396  bool operator < (const Iterator &) const;
397 
402  bool operator > (const Iterator &) const;
403 
410  int operator - (const Iterator &p) const;
411 
415  Iterator operator + (const size_type n) const;
416 
417  private:
422  };
423 
424 }
425 
431 //TODO: Add multithreading to the other vmult functions.
432 
461 template <typename number>
462 class SparseMatrix : public virtual Subscriptor
463 {
464 public:
469 
474  typedef number value_type;
475 
486 
491  typedef
494 
501  typedef
504 
511  struct Traits
512  {
517  static const bool zero_addition_can_be_elided = true;
518  };
519 
534  SparseMatrix ();
535 
544  SparseMatrix (const SparseMatrix &);
545 
559  explicit SparseMatrix (const SparsityPattern &sparsity);
560 
567  SparseMatrix (const SparsityPattern &sparsity,
568  const IdentityMatrix &id);
569 
574  virtual ~SparseMatrix ();
575 
586 
594  operator= (const IdentityMatrix &id);
595 
607  SparseMatrix &operator = (const double d);
608 
622  virtual void reinit (const SparsityPattern &sparsity);
623 
629  virtual void clear ();
631 
639  bool empty () const;
640 
645  size_type m () const;
646 
651  size_type n () const;
652 
656  size_type get_row_length (const size_type row) const;
657 
663  size_type n_nonzero_elements () const;
664 
674  size_type n_actually_nonzero_elements (const double threshold = 0.) const;
675 
684  const SparsityPattern &get_sparsity_pattern () const;
685 
690  std::size_t memory_consumption () const;
691 
695  void compress (::VectorOperation::values);
696 
698 
707  void set (const size_type i,
708  const size_type j,
709  const number value);
710 
726  template <typename number2>
727  void set (const std::vector<size_type> &indices,
728  const FullMatrix<number2> &full_matrix,
729  const bool elide_zero_values = false);
730 
736  template <typename number2>
737  void set (const std::vector<size_type> &row_indices,
738  const std::vector<size_type> &col_indices,
739  const FullMatrix<number2> &full_matrix,
740  const bool elide_zero_values = false);
741 
752  template <typename number2>
753  void set (const size_type row,
754  const std::vector<size_type> &col_indices,
755  const std::vector<number2> &values,
756  const bool elide_zero_values = false);
757 
767  template <typename number2>
768  void set (const size_type row,
769  const size_type n_cols,
770  const size_type *col_indices,
771  const number2 *values,
772  const bool elide_zero_values = false);
773 
779  void add (const size_type i,
780  const size_type j,
781  const number value);
782 
797  template <typename number2>
798  void add (const std::vector<size_type> &indices,
799  const FullMatrix<number2> &full_matrix,
800  const bool elide_zero_values = true);
801 
807  template <typename number2>
808  void add (const std::vector<size_type> &row_indices,
809  const std::vector<size_type> &col_indices,
810  const FullMatrix<number2> &full_matrix,
811  const bool elide_zero_values = true);
812 
822  template <typename number2>
823  void add (const size_type row,
824  const std::vector<size_type> &col_indices,
825  const std::vector<number2> &values,
826  const bool elide_zero_values = true);
827 
837  template <typename number2>
838  void add (const size_type row,
839  const size_type n_cols,
840  const size_type *col_indices,
841  const number2 *values,
842  const bool elide_zero_values = true,
843  const bool col_indices_are_sorted = false);
844 
848  SparseMatrix &operator *= (const number factor);
849 
853  SparseMatrix &operator /= (const number factor);
854 
867  void symmetrize ();
868 
885  template <typename somenumber>
887  copy_from (const SparseMatrix<somenumber> &source);
888 
905  template <typename ForwardIterator>
906  void copy_from (const ForwardIterator begin,
907  const ForwardIterator end);
908 
914  template <typename somenumber>
915  void copy_from (const FullMatrix<somenumber> &matrix);
916 
917 #ifdef DEAL_II_WITH_TRILINOS
918 
929 #endif
930 
942  template <typename somenumber>
943  void add (const number factor,
944  const SparseMatrix<somenumber> &matrix);
945 
947 
951 
965  number operator () (const size_type i,
966  const size_type j) const;
967 
980  number el (const size_type i,
981  const size_type j) const;
982 
993  number diag_element (const size_type i) const;
994 
999  number &diag_element (const size_type i);
1000 
1002 
1022  template <class OutVector, class InVector>
1023  void vmult (OutVector &dst,
1024  const InVector &src) const;
1025 
1041  template <class OutVector, class InVector>
1042  void Tvmult (OutVector &dst,
1043  const InVector &src) const;
1044 
1061  template <class OutVector, class InVector>
1062  void vmult_add (OutVector &dst,
1063  const InVector &src) const;
1064 
1080  template <class OutVector, class InVector>
1081  void Tvmult_add (OutVector &dst,
1082  const InVector &src) const;
1083 
1101  template <typename somenumber>
1102  somenumber matrix_norm_square (const Vector<somenumber> &v) const;
1103 
1109  template <typename somenumber>
1110  somenumber matrix_scalar_product (const Vector<somenumber> &u,
1111  const Vector<somenumber> &v) const;
1112 
1122  template <typename somenumber>
1123  somenumber residual (Vector<somenumber> &dst,
1124  const Vector<somenumber> &x,
1125  const Vector<somenumber> &b) const;
1126 
1150  template <typename numberB, typename numberC>
1151  void mmult (SparseMatrix<numberC> &C,
1152  const SparseMatrix<numberB> &B,
1153  const Vector<number> &V = Vector<number>(),
1154  const bool rebuild_sparsity_pattern = true) const;
1155 
1180  template <typename numberB, typename numberC>
1181  void Tmmult (SparseMatrix<numberC> &C,
1182  const SparseMatrix<numberB> &B,
1183  const Vector<number> &V = Vector<number>(),
1184  const bool rebuild_sparsity_pattern = true) const;
1185 
1187 
1191 
1199  real_type l1_norm () const;
1200 
1208  real_type linfty_norm () const;
1209 
1214  real_type frobenius_norm () const;
1216 
1220 
1226  template <typename somenumber>
1228  const Vector<somenumber> &src,
1229  const number omega = 1.) const;
1230 
1237  template <typename somenumber>
1239  const Vector<somenumber> &src,
1240  const number omega = 1.,
1241  const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>()) const;
1242 
1246  template <typename somenumber>
1248  const Vector<somenumber> &src,
1249  const number om = 1.) const;
1250 
1254  template <typename somenumber>
1256  const Vector<somenumber> &src,
1257  const number om = 1.) const;
1258 
1264  template <typename somenumber>
1265  void SSOR (Vector<somenumber> &v,
1266  const number omega = 1.) const;
1267 
1272  template <typename somenumber>
1273  void SOR (Vector<somenumber> &v,
1274  const number om = 1.) const;
1275 
1280  template <typename somenumber>
1281  void TSOR (Vector<somenumber> &v,
1282  const number om = 1.) const;
1283 
1294  template <typename somenumber>
1295  void PSOR (Vector<somenumber> &v,
1296  const std::vector<size_type> &permutation,
1297  const std::vector<size_type> &inverse_permutation,
1298  const number om = 1.) const;
1299 
1310  template <typename somenumber>
1311  void TPSOR (Vector<somenumber> &v,
1312  const std::vector<size_type> &permutation,
1313  const std::vector<size_type> &inverse_permutation,
1314  const number om = 1.) const;
1315 
1321  template <typename somenumber>
1323  const Vector<somenumber> &b,
1324  const number om = 1.) const;
1325 
1330  template <typename somenumber>
1331  void SOR_step (Vector<somenumber> &v,
1332  const Vector<somenumber> &b,
1333  const number om = 1.) const;
1334 
1339  template <typename somenumber>
1340  void TSOR_step (Vector<somenumber> &v,
1341  const Vector<somenumber> &b,
1342  const number om = 1.) const;
1343 
1348  template <typename somenumber>
1349  void SSOR_step (Vector<somenumber> &v,
1350  const Vector<somenumber> &b,
1351  const number om = 1.) const;
1353 
1357 
1364  const_iterator begin () const;
1365 
1369  iterator begin ();
1370 
1374  const_iterator end () const;
1375 
1379  iterator end ();
1380 
1390  const_iterator begin (const size_type r) const;
1391 
1395  iterator begin (const size_type r);
1396 
1406  const_iterator end (const size_type r) const;
1407 
1411  iterator end (const size_type r);
1413 
1417 
1429  template <class StreamType>
1430  void print (StreamType &out,
1431  const bool across = false,
1432  const bool diagonal_first = true) const;
1433 
1454  void print_formatted (std::ostream &out,
1455  const unsigned int precision = 3,
1456  const bool scientific = true,
1457  const unsigned int width = 0,
1458  const char *zero_string = " ",
1459  const double denominator = 1.) const;
1460 
1466  void print_pattern(std::ostream &out,
1467  const double threshold = 0.) const;
1468 
1479  void block_write (std::ostream &out) const;
1480 
1497  void block_read (std::istream &in);
1499 
1507  DeclException2 (ExcInvalidIndex,
1508  int, int,
1509  << "You are trying to access the matrix entry with index <"
1510  << arg1 << ',' << arg2
1511  << ">, but this entry does not exist in the sparsity pattern "
1512  "of this matrix."
1513  "\n\n"
1514  "The most common cause for this problem is that you used "
1515  "a method to build the sparsity pattern that did not "
1516  "(completely) take into account all of the entries you "
1517  "will later try to write into. An example would be "
1518  "building a sparsity pattern that does not include "
1519  "the entries you will write into due to constraints "
1520  "on degrees of freedom such as hanging nodes or periodic "
1521  "boundary conditions. In such cases, building the "
1522  "sparsity pattern will succeed, but you will get errors "
1523  "such as the current one at one point or other when "
1524  "trying to write into the entries of the matrix.");
1528  DeclExceptionMsg (ExcDifferentSparsityPatterns,
1529  "When copying one sparse matrix into another, "
1530  "or when adding one sparse matrix to another, "
1531  "both matrices need to refer to the same "
1532  "sparsity pattern.");
1536  DeclException2 (ExcIteratorRange,
1537  int, int,
1538  << "The iterators denote a range of " << arg1
1539  << " elements, but the given number of rows was " << arg2);
1543  DeclException0 (ExcSourceEqualsDestination);
1545 
1546 protected:
1557  void prepare_add();
1558 
1563  void prepare_set();
1564 
1565 private:
1572 
1579  number *val;
1580 
1587  std::size_t max_len;
1588 
1589  // make all other sparse matrices friends
1590  template <typename somenumber> friend class SparseMatrix;
1591  template <typename somenumber> friend class SparseLUDecomposition;
1592  template <typename> friend class SparseILU;
1593 
1597  template <typename> friend class BlockMatrixBase;
1598 
1602  template <typename,bool> friend class SparseMatrixIterators::Iterator;
1603  template <typename,bool> friend class SparseMatrixIterators::Accessor;
1604 
1605 #ifndef DEAL_II_MSVC
1606  // Visual studio is choking on the following friend declaration, probably
1607  // because Reference is only defined in a specialization. It looks like
1608  // the library is compiling without this line, though.
1609  template <typename number2> friend class SparseMatrixIterators::Accessor<number2, false>::Reference;
1610 #endif
1611 };
1612 
1613 #ifndef DOXYGEN
1614 /*---------------------- Inline functions -----------------------------------*/
1615 
1616 
1617 
1618 template <typename number>
1619 inline
1620 typename SparseMatrix<number>::size_type SparseMatrix<number>::m () const
1621 {
1622  Assert (cols != 0, ExcNotInitialized());
1623  return cols->rows;
1624 }
1625 
1626 
1627 template <typename number>
1628 inline
1630 {
1631  Assert (cols != 0, ExcNotInitialized());
1632  return cols->cols;
1633 }
1634 
1635 
1636 // Inline the set() and add() functions, since they will be called frequently.
1637 template <typename number>
1638 inline
1639 void
1640 SparseMatrix<number>::set (const size_type i,
1641  const size_type j,
1642  const number value)
1643 {
1644  AssertIsFinite(value);
1645 
1646  const size_type index = cols->operator()(i, j);
1647 
1648  // it is allowed to set elements of the matrix that are not part of the
1649  // sparsity pattern, if the value to which we set it is zero
1650  if (index == SparsityPattern::invalid_entry)
1651  {
1652  Assert ((index != SparsityPattern::invalid_entry) ||
1653  (value == number()),
1654  ExcInvalidIndex(i, j));
1655  return;
1656  }
1657 
1658  val[index] = value;
1659 }
1660 
1661 
1662 
1663 template <typename number>
1664 template <typename number2>
1665 inline
1666 void
1667 SparseMatrix<number>::set (const std::vector<size_type> &indices,
1668  const FullMatrix<number2> &values,
1669  const bool elide_zero_values)
1670 {
1671  Assert (indices.size() == values.m(),
1672  ExcDimensionMismatch(indices.size(), values.m()));
1673  Assert (values.m() == values.n(), ExcNotQuadratic());
1674 
1675  for (size_type i=0; i<indices.size(); ++i)
1676  set (indices[i], indices.size(), &indices[0], &values(i,0),
1677  elide_zero_values);
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> &row_indices,
1687  const std::vector<size_type> &col_indices,
1688  const FullMatrix<number2> &values,
1689  const bool elide_zero_values)
1690 {
1691  Assert (row_indices.size() == values.m(),
1692  ExcDimensionMismatch(row_indices.size(), values.m()));
1693  Assert (col_indices.size() == values.n(),
1694  ExcDimensionMismatch(col_indices.size(), values.n()));
1695 
1696  for (size_type i=0; i<row_indices.size(); ++i)
1697  set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1698  elide_zero_values);
1699 }
1700 
1701 
1702 
1703 template <typename number>
1704 template <typename number2>
1705 inline
1706 void
1707 SparseMatrix<number>::set (const size_type row,
1708  const std::vector<size_type> &col_indices,
1709  const std::vector<number2> &values,
1710  const bool elide_zero_values)
1711 {
1712  Assert (col_indices.size() == values.size(),
1713  ExcDimensionMismatch(col_indices.size(), values.size()));
1714 
1715  set (row, col_indices.size(), &col_indices[0], &values[0],
1716  elide_zero_values);
1717 }
1718 
1719 
1720 
1721 template <typename number>
1722 inline
1723 void
1724 SparseMatrix<number>::add (const size_type i,
1725  const size_type j,
1726  const number value)
1727 {
1728  AssertIsFinite(value);
1729 
1730  if (value == number())
1731  return;
1732 
1733  const size_type index = cols->operator()(i, j);
1734 
1735  // it is allowed to add elements to the matrix that are not part of the
1736  // sparsity pattern, if the value to which we set it is zero
1737  if (index == SparsityPattern::invalid_entry)
1738  {
1739  Assert ((index != SparsityPattern::invalid_entry) ||
1740  (value == number()),
1741  ExcInvalidIndex(i, j));
1742  return;
1743  }
1744 
1745  val[index] += value;
1746 }
1747 
1748 
1749 
1750 template <typename number>
1751 template <typename number2>
1752 inline
1753 void
1754 SparseMatrix<number>::add (const std::vector<size_type> &indices,
1755  const FullMatrix<number2> &values,
1756  const bool elide_zero_values)
1757 {
1758  Assert (indices.size() == values.m(),
1759  ExcDimensionMismatch(indices.size(), values.m()));
1760  Assert (values.m() == values.n(), ExcNotQuadratic());
1761 
1762  for (size_type i=0; i<indices.size(); ++i)
1763  add (indices[i], indices.size(), &indices[0], &values(i,0),
1764  elide_zero_values);
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> &row_indices,
1774  const std::vector<size_type> &col_indices,
1775  const FullMatrix<number2> &values,
1776  const bool elide_zero_values)
1777 {
1778  Assert (row_indices.size() == values.m(),
1779  ExcDimensionMismatch(row_indices.size(), values.m()));
1780  Assert (col_indices.size() == values.n(),
1781  ExcDimensionMismatch(col_indices.size(), values.n()));
1782 
1783  for (size_type i=0; i<row_indices.size(); ++i)
1784  add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1785  elide_zero_values);
1786 }
1787 
1788 
1789 
1790 template <typename number>
1791 template <typename number2>
1792 inline
1793 void
1794 SparseMatrix<number>::add (const size_type row,
1795  const std::vector<size_type> &col_indices,
1796  const std::vector<number2> &values,
1797  const bool elide_zero_values)
1798 {
1799  Assert (col_indices.size() == values.size(),
1800  ExcDimensionMismatch(col_indices.size(), values.size()));
1801 
1802  add (row, col_indices.size(), &col_indices[0], &values[0],
1803  elide_zero_values);
1804 }
1805 
1806 
1807 
1808 template <typename number>
1809 inline
1811 SparseMatrix<number>::operator *= (const number factor)
1812 {
1813  Assert (cols != 0, ExcNotInitialized());
1814  Assert (val != 0, ExcNotInitialized());
1815 
1816  number *val_ptr = &val[0];
1817  const number *const end_ptr = &val[cols->n_nonzero_elements()];
1818 
1819  while (val_ptr != end_ptr)
1820  *val_ptr++ *= factor;
1821 
1822  return *this;
1823 }
1824 
1825 
1826 
1827 template <typename number>
1828 inline
1830 SparseMatrix<number>::operator /= (const number factor)
1831 {
1832  Assert (cols != 0, ExcNotInitialized());
1833  Assert (val != 0, ExcNotInitialized());
1834  Assert (factor != number(), ExcDivideByZero());
1835 
1836  const number factor_inv = number(1.) / factor;
1837 
1838  number *val_ptr = &val[0];
1839  const number *const end_ptr = &val[cols->n_nonzero_elements()];
1840 
1841  while (val_ptr != end_ptr)
1842  *val_ptr++ *= factor_inv;
1843 
1844  return *this;
1845 }
1846 
1847 
1848 
1849 template <typename number>
1850 inline
1851 number SparseMatrix<number>::operator () (const size_type i,
1852  const size_type j) const
1853 {
1854  Assert (cols != 0, ExcNotInitialized());
1855  Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
1856  ExcInvalidIndex(i,j));
1857  return val[cols->operator()(i,j)];
1858 }
1859 
1860 
1861 
1862 template <typename number>
1863 inline
1864 number SparseMatrix<number>::el (const size_type i,
1865  const size_type j) const
1866 {
1867  Assert (cols != 0, ExcNotInitialized());
1868  const size_type index = cols->operator()(i,j);
1869 
1870  if (index != SparsityPattern::invalid_entry)
1871  return val[index];
1872  else
1873  return 0;
1874 }
1875 
1876 
1877 
1878 template <typename number>
1879 inline
1880 number SparseMatrix<number>::diag_element (const size_type i) const
1881 {
1882  Assert (cols != 0, ExcNotInitialized());
1883  Assert (m() == n(), ExcNotQuadratic());
1884  AssertIndexRange(i, m());
1885 
1886  // Use that the first element in each row of a quadratic matrix is the main
1887  // diagonal
1888  return val[cols->rowstart[i]];
1889 }
1890 
1891 
1892 
1893 template <typename number>
1894 inline
1895 number &SparseMatrix<number>::diag_element (const size_type i)
1896 {
1897  Assert (cols != 0, ExcNotInitialized());
1898  Assert (m() == n(), ExcNotQuadratic());
1899  AssertIndexRange(i, m());
1900 
1901  // Use that the first element in each row of a quadratic matrix is the main
1902  // diagonal
1903  return val[cols->rowstart[i]];
1904 }
1905 
1906 
1907 
1908 template <typename number>
1909 template <typename ForwardIterator>
1910 void
1911 SparseMatrix<number>::copy_from (const ForwardIterator begin,
1912  const ForwardIterator end)
1913 {
1914  Assert (static_cast<size_type>(std::distance (begin, end)) == m(),
1915  ExcIteratorRange (std::distance (begin, end), m()));
1916 
1917  // for use in the inner loop, we define a typedef to the type of the inner
1918  // iterators
1919  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1920  size_type row=0;
1921  for (ForwardIterator i=begin; i!=end; ++i, ++row)
1922  {
1923  const inner_iterator end_of_row = i->end();
1924  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1925  // write entries
1926  set (row, j->first, j->second);
1927  };
1928 }
1929 
1930 
1931 //---------------------------------------------------------------------------
1932 
1933 
1934 namespace SparseMatrixIterators
1935 {
1936  template <typename number>
1937  inline
1938  Accessor<number,true>::
1939  Accessor (const MatrixType *matrix,
1940  const std::size_t index_within_matrix)
1941  :
1942  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1943  index_within_matrix),
1944  matrix (matrix)
1945  {}
1946 
1947 
1948 
1949  template <typename number>
1950  inline
1951  Accessor<number,true>::
1952  Accessor (const MatrixType *matrix)
1953  :
1954  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1955  matrix (matrix)
1956  {}
1957 
1958 
1959 
1960  template <typename number>
1961  inline
1962  Accessor<number,true>::
1964  :
1965  SparsityPatternIterators::Accessor (a),
1966  matrix (&a.get_matrix())
1967  {}
1968 
1969 
1970 
1971  template <typename number>
1972  inline
1973  number
1974  Accessor<number, true>::value () const
1975  {
1976  AssertIndexRange(index_within_sparsity, matrix->n_nonzero_elements());
1977  return matrix->val[index_within_sparsity];
1978  }
1979 
1980 
1981 
1982  template <typename number>
1983  inline
1984  const typename Accessor<number, true>::MatrixType &
1985  Accessor<number, true>::get_matrix () const
1986  {
1987  return *matrix;
1988  }
1989 
1990 
1991 
1992  template <typename number>
1993  inline
1994  Accessor<number, false>::Reference::Reference (
1995  const Accessor *accessor,
1996  const bool)
1997  :
1998  accessor (accessor)
1999  {}
2000 
2001 
2002  template <typename number>
2003  inline
2004  Accessor<number, false>::Reference::operator number() const
2005  {
2006  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2007  return accessor->matrix->val[accessor->index_within_sparsity];
2008  }
2009 
2010 
2011 
2012  template <typename number>
2013  inline
2014  const typename Accessor<number, false>::Reference &
2015  Accessor<number, false>::Reference::operator = (const number n) const
2016  {
2017  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2018  accessor->matrix->val[accessor->index_within_sparsity] = n;
2019  return *this;
2020  }
2021 
2022 
2023 
2024  template <typename number>
2025  inline
2026  const typename Accessor<number, false>::Reference &
2027  Accessor<number, false>::Reference::operator += (const number n) const
2028  {
2029  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2030  accessor->matrix->val[accessor->index_within_sparsity] += n;
2031  return *this;
2032  }
2033 
2034 
2035 
2036  template <typename number>
2037  inline
2038  const typename Accessor<number, false>::Reference &
2039  Accessor<number, false>::Reference::operator -= (const number n) const
2040  {
2041  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2042  accessor->matrix->val[accessor->index_within_sparsity] -= n;
2043  return *this;
2044  }
2045 
2046 
2047 
2048  template <typename number>
2049  inline
2050  const typename Accessor<number, false>::Reference &
2051  Accessor<number, false>::Reference::operator *= (const number n) const
2052  {
2053  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2054  accessor->matrix->val[accessor->index_within_sparsity] *= n;
2055  return *this;
2056  }
2057 
2058 
2059 
2060  template <typename number>
2061  inline
2062  const typename Accessor<number, false>::Reference &
2063  Accessor<number, false>::Reference::operator /= (const number n) const
2064  {
2065  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2066  accessor->matrix->val[accessor->index_within_sparsity] /= n;
2067  return *this;
2068  }
2069 
2070 
2071 
2072  template <typename number>
2073  inline
2074  Accessor<number,false>::
2075  Accessor (MatrixType *matrix,
2076  const std::size_t index)
2077  :
2078  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2079  index),
2080  matrix (matrix)
2081  {}
2082 
2083 
2084 
2085  template <typename number>
2086  inline
2087  Accessor<number,false>::
2088  Accessor (MatrixType *matrix)
2089  :
2090  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
2091  matrix (matrix)
2092  {}
2093 
2094 
2095 
2096  template <typename number>
2097  inline
2098  typename Accessor<number, false>::Reference
2099  Accessor<number, false>::value() const
2100  {
2101  return Reference(this,true);
2102  }
2103 
2104 
2105 
2106 
2107  template <typename number>
2108  inline
2109  typename Accessor<number, false>::MatrixType &
2110  Accessor<number, false>::get_matrix () const
2111  {
2112  return *matrix;
2113  }
2114 
2115 
2116 
2117  template <typename number, bool Constness>
2118  inline
2119  Iterator<number, Constness>::
2120  Iterator (MatrixType *matrix,
2121  const std::size_t index)
2122  :
2123  accessor(matrix, index)
2124  {}
2125 
2126 
2127 
2128  template <typename number, bool Constness>
2129  inline
2130  Iterator<number, Constness>::
2131  Iterator (MatrixType *matrix)
2132  :
2133  accessor(matrix)
2134  {}
2135 
2136 
2137 
2138  template <typename number, bool Constness>
2139  inline
2140  Iterator<number, Constness>::
2142  :
2143  accessor(*i)
2144  {}
2145 
2146 
2147 
2148  template <typename number, bool Constness>
2149  inline
2150  Iterator<number, Constness> &
2151  Iterator<number,Constness>::operator++ ()
2152  {
2153  accessor.advance ();
2154  return *this;
2155  }
2156 
2157 
2158  template <typename number, bool Constness>
2159  inline
2160  Iterator<number,Constness>
2161  Iterator<number,Constness>::operator++ (int)
2162  {
2163  const Iterator iter = *this;
2164  accessor.advance ();
2165  return iter;
2166  }
2167 
2168 
2169  template <typename number, bool Constness>
2170  inline
2171  const Accessor<number,Constness> &
2172  Iterator<number,Constness>::operator* () const
2173  {
2174  return accessor;
2175  }
2176 
2177 
2178  template <typename number, bool Constness>
2179  inline
2180  const Accessor<number,Constness> *
2181  Iterator<number,Constness>::operator-> () const
2182  {
2183  return &accessor;
2184  }
2185 
2186 
2187  template <typename number, bool Constness>
2188  inline
2189  bool
2190  Iterator<number,Constness>::
2191  operator == (const Iterator &other) const
2192  {
2193  return (accessor == other.accessor);
2194  }
2195 
2196 
2197  template <typename number, bool Constness>
2198  inline
2199  bool
2200  Iterator<number,Constness>::
2201  operator != (const Iterator &other) const
2202  {
2203  return ! (*this == other);
2204  }
2205 
2206 
2207  template <typename number, bool Constness>
2208  inline
2209  bool
2210  Iterator<number,Constness>::
2211  operator < (const Iterator &other) const
2212  {
2213  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2214  ExcInternalError());
2215 
2216  return (accessor < other.accessor);
2217  }
2218 
2219 
2220  template <typename number, bool Constness>
2221  inline
2222  bool
2223  Iterator<number,Constness>::
2224  operator > (const Iterator &other) const
2225  {
2226  return (other < *this);
2227  }
2228 
2229 
2230  template <typename number, bool Constness>
2231  inline
2232  int
2233  Iterator<number,Constness>::
2234  operator - (const Iterator &other) const
2235  {
2236  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2237  ExcInternalError());
2238 
2239  return (*this)->index_within_sparsity - other->index_within_sparsity;
2240  }
2241 
2242 
2243 
2244  template <typename number, bool Constness>
2245  inline
2246  Iterator<number,Constness>
2247  Iterator<number,Constness>::
2248  operator + (const size_type n) const
2249  {
2250  Iterator x = *this;
2251  for (size_type i=0; i<n; ++i)
2252  ++x;
2253 
2254  return x;
2255  }
2256 
2257 }
2258 
2259 
2260 
2261 template <typename number>
2262 inline
2265 {
2266  return const_iterator(this, 0);
2267 }
2268 
2269 
2270 template <typename number>
2271 inline
2274 {
2275  return const_iterator(this);
2276 }
2277 
2278 
2279 template <typename number>
2280 inline
2283 {
2284  return iterator (this, 0);
2285 }
2286 
2287 
2288 template <typename number>
2289 inline
2292 {
2293  return iterator(this, cols->rowstart[cols->rows]);
2294 }
2295 
2296 
2297 template <typename number>
2298 inline
2300 SparseMatrix<number>::begin (const size_type r) const
2301 {
2302  Assert (r<m(), ExcIndexRange(r,0,m()));
2303 
2304  return const_iterator(this, cols->rowstart[r]);
2305 }
2306 
2307 
2308 
2309 template <typename number>
2310 inline
2312 SparseMatrix<number>::end (const size_type r) const
2313 {
2314  Assert (r<m(), ExcIndexRange(r,0,m()));
2315 
2316  return const_iterator(this, cols->rowstart[r+1]);
2317 }
2318 
2319 
2320 
2321 template <typename number>
2322 inline
2324 SparseMatrix<number>::begin (const size_type r)
2325 {
2326  Assert (r<m(), ExcIndexRange(r,0,m()));
2327 
2328  return iterator(this, cols->rowstart[r]);
2329 }
2330 
2331 
2332 
2333 template <typename number>
2334 inline
2336 SparseMatrix<number>::end (const size_type r)
2337 {
2338  Assert (r<m(), ExcIndexRange(r,0,m()));
2339 
2340  return iterator(this, cols->rowstart[r+1]);
2341 }
2342 
2343 
2344 
2345 template <typename number>
2346 template <class StreamType>
2347 inline
2348 void SparseMatrix<number>::print (StreamType &out,
2349  const bool across,
2350  const bool diagonal_first) const
2351 {
2352  Assert (cols != 0, ExcNotInitialized());
2353  Assert (val != 0, ExcNotInitialized());
2354 
2355  bool hanging_diagonal = false;
2356  number diagonal = number();
2357 
2358  for (size_type i=0; i<cols->rows; ++i)
2359  {
2360  for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
2361  {
2362  if (!diagonal_first && i == cols->colnums[j])
2363  {
2364  diagonal = val[j];
2365  hanging_diagonal = true;
2366  }
2367  else
2368  {
2369  if (hanging_diagonal && cols->colnums[j]>i)
2370  {
2371  if (across)
2372  out << ' ' << i << ',' << i << ':' << diagonal;
2373  else
2374  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2375  hanging_diagonal = false;
2376  }
2377  if (across)
2378  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2379  else
2380  out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << std::endl;
2381  }
2382  }
2383  if (hanging_diagonal)
2384  {
2385  if (across)
2386  out << ' ' << i << ',' << i << ':' << diagonal;
2387  else
2388  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2389  hanging_diagonal = false;
2390  }
2391  }
2392  if (across)
2393  out << std::endl;
2394 }
2395 
2396 
2397 template <typename number>
2398 inline
2399 void
2401 prepare_add()
2402 {
2403  //nothing to do here
2404 }
2405 
2406 
2407 
2408 template <typename number>
2409 inline
2410 void
2412 prepare_set()
2413 {
2414  //nothing to do here
2415 }
2416 
2417 #endif // DOXYGEN
2418 
2419 
2420 /*---------------------------- sparse_matrix.h ---------------------------*/
2421 
2422 DEAL_II_NAMESPACE_CLOSE
2423 
2424 #endif
2425 /*---------------------------- 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
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
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
DeclException0(ExcSourceEqualsDestination)
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:1081
void block_read(std::istream &in)
SparseMatrixIterators::Iterator< number, false > iterator
real_type frobenius_norm() const
size_type n_nonzero_elements() const
real_type linfty_norm() const
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
void SOR(Vector< somenumber > &v, const number om=1.) const
DeclExceptionMsg(ExcDifferentSparsityPatterns,"When copying one sparse matrix into another, ""or when adding one sparse matrix to another, ""both matrices need to refer to the same ""sparsity pattern.")
const_iterator begin() const
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:294
std::size_t memory_consumption() const
size_type get_row_length(const size_type row) const
types::global_dof_index size_type
Definition: sparse_matrix.h:56
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
size_type n_nonzero_elements() const
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
std::size_t * rowstart
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
virtual void clear()
#define AssertIsFinite(number)
Definition: exceptions.h:1096
DeclException2(ExcInvalidIndex, int, int,<< "You are trying to access the matrix entry with index <"<< arg1<< ','<< arg2<< ">, but this entry does not exist in the sparsity pattern ""of this matrix.""\n\n""The most common cause for this problem is that you used ""a method to build the sparsity pattern that did not ""(completely) take into account all of the entries you ""will later try to write into. An example would be ""building a sparsity pattern that does not include ""the entries you will write into due to constraints ""on degrees of freedom such as hanging nodes or periodic ""boundary conditions. In such cases, building the ""sparsity pattern will succeed, but you will get errors ""such as the current one at one point or other when ""trying to write into the entries of the matrix.")
size_type n_actually_nonzero_elements(const double threshold=0.) const
std::size_t max_len
void compress(::VectorOperation::values)
SparseMatrix & operator*=(const number factor)
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const