Reference documentation for deal.II version Git 1ad260a 2018-02-22 16:12:58 +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 
598 
606  operator= (const IdentityMatrix &id);
607 
619  SparseMatrix &operator = (const double d);
620 
634  virtual void reinit (const SparsityPattern &sparsity);
635 
641  virtual void clear ();
643 
651  bool empty () const;
652 
657  size_type m () const;
658 
663  size_type n () const;
664 
668  size_type get_row_length (const size_type row) const;
669 
675  std::size_t n_nonzero_elements () const;
676 
686  std::size_t n_actually_nonzero_elements (const double threshold = 0.) const;
687 
696  const SparsityPattern &get_sparsity_pattern () const;
697 
702  std::size_t memory_consumption () const;
703 
708 
710 
719  void set (const size_type i,
720  const size_type j,
721  const number value);
722 
738  template <typename number2>
739  void set (const std::vector<size_type> &indices,
740  const FullMatrix<number2> &full_matrix,
741  const bool elide_zero_values = false);
742 
748  template <typename number2>
749  void set (const std::vector<size_type> &row_indices,
750  const std::vector<size_type> &col_indices,
751  const FullMatrix<number2> &full_matrix,
752  const bool elide_zero_values = false);
753 
764  template <typename number2>
765  void set (const size_type row,
766  const std::vector<size_type> &col_indices,
767  const std::vector<number2> &values,
768  const bool elide_zero_values = false);
769 
779  template <typename number2>
780  void set (const size_type row,
781  const size_type n_cols,
782  const size_type *col_indices,
783  const number2 *values,
784  const bool elide_zero_values = false);
785 
791  void add (const size_type i,
792  const size_type j,
793  const number value);
794 
809  template <typename number2>
810  void add (const std::vector<size_type> &indices,
811  const FullMatrix<number2> &full_matrix,
812  const bool elide_zero_values = true);
813 
819  template <typename number2>
820  void add (const std::vector<size_type> &row_indices,
821  const std::vector<size_type> &col_indices,
822  const FullMatrix<number2> &full_matrix,
823  const bool elide_zero_values = true);
824 
834  template <typename number2>
835  void add (const size_type row,
836  const std::vector<size_type> &col_indices,
837  const std::vector<number2> &values,
838  const bool elide_zero_values = true);
839 
849  template <typename number2>
850  void add (const size_type row,
851  const size_type n_cols,
852  const size_type *col_indices,
853  const number2 *values,
854  const bool elide_zero_values = true,
855  const bool col_indices_are_sorted = false);
856 
860  SparseMatrix &operator *= (const number factor);
861 
865  SparseMatrix &operator /= (const number factor);
866 
879  void symmetrize ();
880 
897  template <typename somenumber>
899  copy_from (const SparseMatrix<somenumber> &source);
900 
917  template <typename ForwardIterator>
918  void copy_from (const ForwardIterator begin,
919  const ForwardIterator end);
920 
926  template <typename somenumber>
927  void copy_from (const FullMatrix<somenumber> &matrix);
928 
929 #ifdef DEAL_II_WITH_TRILINOS
930 
941 #endif
942 
954  template <typename somenumber>
955  void add (const number factor,
956  const SparseMatrix<somenumber> &matrix);
957 
959 
963 
977  number operator () (const size_type i,
978  const size_type j) const;
979 
992  number el (const size_type i,
993  const size_type j) const;
994 
1004  number diag_element (const size_type i) const;
1005 
1010  number &diag_element (const size_type i);
1011 
1013 
1033  template <class OutVector, class InVector>
1034  void vmult (OutVector &dst,
1035  const InVector &src) const;
1036 
1052  template <class OutVector, class InVector>
1053  void Tvmult (OutVector &dst,
1054  const InVector &src) const;
1055 
1072  template <class OutVector, class InVector>
1073  void vmult_add (OutVector &dst,
1074  const InVector &src) const;
1075 
1091  template <class OutVector, class InVector>
1092  void Tvmult_add (OutVector &dst,
1093  const InVector &src) const;
1094 
1112  template <typename somenumber>
1113  somenumber matrix_norm_square (const Vector<somenumber> &v) const;
1114 
1120  template <typename somenumber>
1121  somenumber matrix_scalar_product (const Vector<somenumber> &u,
1122  const Vector<somenumber> &v) const;
1123 
1133  template <typename somenumber>
1134  somenumber residual (Vector<somenumber> &dst,
1135  const Vector<somenumber> &x,
1136  const Vector<somenumber> &b) const;
1137 
1173  template <typename numberB, typename numberC>
1174  void mmult (SparseMatrix<numberC> &C,
1175  const SparseMatrix<numberB> &B,
1176  const Vector<number> &V = Vector<number>(),
1177  const bool rebuild_sparsity_pattern = true) const;
1178 
1203  template <typename numberB, typename numberC>
1204  void Tmmult (SparseMatrix<numberC> &C,
1205  const SparseMatrix<numberB> &B,
1206  const Vector<number> &V = Vector<number>(),
1207  const bool rebuild_sparsity_pattern = true) const;
1208 
1210 
1214 
1222  real_type l1_norm () const;
1223 
1231  real_type linfty_norm () const;
1232 
1237  real_type frobenius_norm () const;
1239 
1243 
1249  template <typename somenumber>
1251  const Vector<somenumber> &src,
1252  const number omega = 1.) const;
1253 
1260  template <typename somenumber>
1262  const Vector<somenumber> &src,
1263  const number omega = 1.,
1264  const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>()) const;
1265 
1269  template <typename somenumber>
1271  const Vector<somenumber> &src,
1272  const number om = 1.) const;
1273 
1277  template <typename somenumber>
1279  const Vector<somenumber> &src,
1280  const number om = 1.) const;
1281 
1287  template <typename somenumber>
1288  void SSOR (Vector<somenumber> &v,
1289  const number omega = 1.) const;
1290 
1295  template <typename somenumber>
1296  void SOR (Vector<somenumber> &v,
1297  const number om = 1.) const;
1298 
1303  template <typename somenumber>
1304  void TSOR (Vector<somenumber> &v,
1305  const number om = 1.) const;
1306 
1317  template <typename somenumber>
1318  void PSOR (Vector<somenumber> &v,
1319  const std::vector<size_type> &permutation,
1320  const std::vector<size_type> &inverse_permutation,
1321  const number om = 1.) const;
1322 
1333  template <typename somenumber>
1334  void TPSOR (Vector<somenumber> &v,
1335  const std::vector<size_type> &permutation,
1336  const std::vector<size_type> &inverse_permutation,
1337  const number om = 1.) const;
1338 
1344  template <typename somenumber>
1346  const Vector<somenumber> &b,
1347  const number om = 1.) const;
1348 
1353  template <typename somenumber>
1354  void SOR_step (Vector<somenumber> &v,
1355  const Vector<somenumber> &b,
1356  const number om = 1.) const;
1357 
1362  template <typename somenumber>
1363  void TSOR_step (Vector<somenumber> &v,
1364  const Vector<somenumber> &b,
1365  const number om = 1.) const;
1366 
1371  template <typename somenumber>
1372  void SSOR_step (Vector<somenumber> &v,
1373  const Vector<somenumber> &b,
1374  const number om = 1.) const;
1376 
1380 
1387  const_iterator begin () const;
1388 
1392  iterator begin ();
1393 
1397  const_iterator end () const;
1398 
1402  iterator end ();
1403 
1413  const_iterator begin (const size_type r) const;
1414 
1418  iterator begin (const size_type r);
1419 
1429  const_iterator end (const size_type r) const;
1430 
1434  iterator end (const size_type r);
1436 
1440 
1452  template <class StreamType>
1453  void print (StreamType &out,
1454  const bool across = false,
1455  const bool diagonal_first = true) const;
1456 
1477  void print_formatted (std::ostream &out,
1478  const unsigned int precision = 3,
1479  const bool scientific = true,
1480  const unsigned int width = 0,
1481  const char *zero_string = " ",
1482  const double denominator = 1.) const;
1483 
1489  void print_pattern(std::ostream &out,
1490  const double threshold = 0.) const;
1491 
1502  void block_write (std::ostream &out) const;
1503 
1520  void block_read (std::istream &in);
1522 
1531  int, int,
1532  << "You are trying to access the matrix entry with index <"
1533  << arg1 << ',' << arg2
1534  << ">, but this entry does not exist in the sparsity pattern "
1535  "of this matrix."
1536  "\n\n"
1537  "The most common cause for this problem is that you used "
1538  "a method to build the sparsity pattern that did not "
1539  "(completely) take into account all of the entries you "
1540  "will later try to write into. An example would be "
1541  "building a sparsity pattern that does not include "
1542  "the entries you will write into due to constraints "
1543  "on degrees of freedom such as hanging nodes or periodic "
1544  "boundary conditions. In such cases, building the "
1545  "sparsity pattern will succeed, but you will get errors "
1546  "such as the current one at one point or other when "
1547  "trying to write into the entries of the matrix.");
1552  "When copying one sparse matrix into another, "
1553  "or when adding one sparse matrix to another, "
1554  "both matrices need to refer to the same "
1555  "sparsity pattern.");
1560  int, int,
1561  << "The iterators denote a range of " << arg1
1562  << " elements, but the given number of rows was " << arg2);
1567  "You are attempting an operation on two matrices that "
1568  "are the same object, but the operation requires that the "
1569  "two objects are in fact different.");
1571 
1572 protected:
1583  void prepare_add();
1584 
1589  void prepare_set();
1590 
1591 private:
1598 
1606  std::unique_ptr<number[]> val;
1607 
1614  std::size_t max_len;
1615 
1616  // make all other sparse matrices friends
1617  template <typename somenumber> friend class SparseMatrix;
1618  template <typename somenumber> friend class SparseLUDecomposition;
1619  template <typename> friend class SparseILU;
1620 
1624  template <typename> friend class BlockMatrixBase;
1625 
1629  template <typename,bool> friend class SparseMatrixIterators::Iterator;
1630  template <typename,bool> friend class SparseMatrixIterators::Accessor;
1631 };
1632 
1633 #ifndef DOXYGEN
1634 /*---------------------- Inline functions -----------------------------------*/
1635 
1636 
1637 
1638 template <typename number>
1639 inline
1641 {
1642  Assert (cols != nullptr, ExcNotInitialized());
1643  return cols->rows;
1644 }
1645 
1646 
1647 template <typename number>
1648 inline
1650 {
1651  Assert (cols != nullptr, ExcNotInitialized());
1652  return cols->cols;
1653 }
1654 
1655 
1656 // Inline the set() and add() functions, since they will be called frequently.
1657 template <typename number>
1658 inline
1659 void
1660 SparseMatrix<number>::set (const size_type i,
1661  const size_type j,
1662  const number value)
1663 {
1664  AssertIsFinite(value);
1665 
1666  const size_type index = cols->operator()(i, j);
1667 
1668  // it is allowed to set elements of the matrix that are not part of the
1669  // sparsity pattern, if the value to which we set it is zero
1670  if (index == SparsityPattern::invalid_entry)
1671  {
1672  Assert ((index != SparsityPattern::invalid_entry) ||
1673  (value == number()),
1674  ExcInvalidIndex(i, j));
1675  return;
1676  }
1677 
1678  val[index] = value;
1679 }
1680 
1681 
1682 
1683 template <typename number>
1684 template <typename number2>
1685 inline
1686 void
1687 SparseMatrix<number>::set (const std::vector<size_type> &indices,
1688  const FullMatrix<number2> &values,
1689  const bool elide_zero_values)
1690 {
1691  Assert (indices.size() == values.m(),
1692  ExcDimensionMismatch(indices.size(), values.m()));
1693  Assert (values.m() == values.n(), ExcNotQuadratic());
1694 
1695  for (size_type i=0; i<indices.size(); ++i)
1696  set (indices[i], indices.size(), indices.data(), &values(i,0),
1697  elide_zero_values);
1698 }
1699 
1700 
1701 
1702 template <typename number>
1703 template <typename number2>
1704 inline
1705 void
1706 SparseMatrix<number>::set (const std::vector<size_type> &row_indices,
1707  const std::vector<size_type> &col_indices,
1708  const FullMatrix<number2> &values,
1709  const bool elide_zero_values)
1710 {
1711  Assert (row_indices.size() == values.m(),
1712  ExcDimensionMismatch(row_indices.size(), values.m()));
1713  Assert (col_indices.size() == values.n(),
1714  ExcDimensionMismatch(col_indices.size(), values.n()));
1715 
1716  for (size_type i=0; i<row_indices.size(); ++i)
1717  set (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
1718  elide_zero_values);
1719 }
1720 
1721 
1722 
1723 template <typename number>
1724 template <typename number2>
1725 inline
1726 void
1727 SparseMatrix<number>::set (const size_type row,
1728  const std::vector<size_type> &col_indices,
1729  const std::vector<number2> &values,
1730  const bool elide_zero_values)
1731 {
1732  Assert (col_indices.size() == values.size(),
1733  ExcDimensionMismatch(col_indices.size(), values.size()));
1734 
1735  set (row, col_indices.size(), col_indices.data(), values.data(),
1736  elide_zero_values);
1737 }
1738 
1739 
1740 
1741 template <typename number>
1742 inline
1743 void
1744 SparseMatrix<number>::add (const size_type i,
1745  const size_type j,
1746  const number value)
1747 {
1748  AssertIsFinite(value);
1749 
1750  if (value == number())
1751  return;
1752 
1753  const size_type index = cols->operator()(i, j);
1754 
1755  // it is allowed to add elements to the matrix that are not part of the
1756  // sparsity pattern, if the value to which we set it is zero
1757  if (index == SparsityPattern::invalid_entry)
1758  {
1759  Assert ((index != SparsityPattern::invalid_entry) ||
1760  (value == number()),
1761  ExcInvalidIndex(i, j));
1762  return;
1763  }
1764 
1765  val[index] += value;
1766 }
1767 
1768 
1769 
1770 template <typename number>
1771 template <typename number2>
1772 inline
1773 void
1774 SparseMatrix<number>::add (const std::vector<size_type> &indices,
1775  const FullMatrix<number2> &values,
1776  const bool elide_zero_values)
1777 {
1778  Assert (indices.size() == values.m(),
1779  ExcDimensionMismatch(indices.size(), values.m()));
1780  Assert (values.m() == values.n(), ExcNotQuadratic());
1781 
1782  for (size_type i=0; i<indices.size(); ++i)
1783  add (indices[i], indices.size(), indices.data(), &values(i,0),
1784  elide_zero_values);
1785 }
1786 
1787 
1788 
1789 template <typename number>
1790 template <typename number2>
1791 inline
1792 void
1793 SparseMatrix<number>::add (const std::vector<size_type> &row_indices,
1794  const std::vector<size_type> &col_indices,
1795  const FullMatrix<number2> &values,
1796  const bool elide_zero_values)
1797 {
1798  Assert (row_indices.size() == values.m(),
1799  ExcDimensionMismatch(row_indices.size(), values.m()));
1800  Assert (col_indices.size() == values.n(),
1801  ExcDimensionMismatch(col_indices.size(), values.n()));
1802 
1803  for (size_type i=0; i<row_indices.size(); ++i)
1804  add (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
1805  elide_zero_values);
1806 }
1807 
1808 
1809 
1810 template <typename number>
1811 template <typename number2>
1812 inline
1813 void
1814 SparseMatrix<number>::add (const size_type row,
1815  const std::vector<size_type> &col_indices,
1816  const std::vector<number2> &values,
1817  const bool elide_zero_values)
1818 {
1819  Assert (col_indices.size() == values.size(),
1820  ExcDimensionMismatch(col_indices.size(), values.size()));
1821 
1822  add (row, col_indices.size(), col_indices.data(), values.data(),
1823  elide_zero_values);
1824 }
1825 
1826 
1827 
1828 template <typename number>
1829 inline
1831 SparseMatrix<number>::operator *= (const number factor)
1832 {
1833  Assert (cols != nullptr, ExcNotInitialized());
1834  Assert (val != nullptr, ExcNotInitialized());
1835 
1836  number *val_ptr = val.get();
1837  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1838 
1839  while (val_ptr != end_ptr)
1840  *val_ptr++ *= factor;
1841 
1842  return *this;
1843 }
1844 
1845 
1846 
1847 template <typename number>
1848 inline
1850 SparseMatrix<number>::operator /= (const number factor)
1851 {
1852  Assert (cols != nullptr, ExcNotInitialized());
1853  Assert (val != nullptr, ExcNotInitialized());
1854  Assert (factor != number(), ExcDivideByZero());
1855 
1856  const number factor_inv = number(1.) / factor;
1857 
1858  number *val_ptr = val.get();
1859  const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1860 
1861  while (val_ptr != end_ptr)
1862  *val_ptr++ *= factor_inv;
1863 
1864  return *this;
1865 }
1866 
1867 
1868 
1869 template <typename number>
1870 inline
1871 number SparseMatrix<number>::operator () (const size_type i,
1872  const size_type j) const
1873 {
1874  Assert (cols != nullptr, ExcNotInitialized());
1875  Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
1876  ExcInvalidIndex(i,j));
1877  return val[cols->operator()(i,j)];
1878 }
1879 
1880 
1881 
1882 template <typename number>
1883 inline
1884 number SparseMatrix<number>::el (const size_type i,
1885  const size_type j) const
1886 {
1887  Assert (cols != nullptr, ExcNotInitialized());
1888  const size_type index = cols->operator()(i,j);
1889 
1890  if (index != SparsityPattern::invalid_entry)
1891  return val[index];
1892  else
1893  return 0;
1894 }
1895 
1896 
1897 
1898 template <typename number>
1899 inline
1900 number SparseMatrix<number>::diag_element (const size_type i) const
1901 {
1902  Assert (cols != nullptr, ExcNotInitialized());
1903  Assert (m() == n(), ExcNotQuadratic());
1904  AssertIndexRange(i, m());
1905 
1906  // Use that the first element in each row of a quadratic matrix is the main
1907  // diagonal
1908  return val[cols->rowstart[i]];
1909 }
1910 
1911 
1912 
1913 template <typename number>
1914 inline
1915 number &SparseMatrix<number>::diag_element (const size_type i)
1916 {
1917  Assert (cols != nullptr, ExcNotInitialized());
1918  Assert (m() == n(), ExcNotQuadratic());
1919  AssertIndexRange(i, m());
1920 
1921  // Use that the first element in each row of a quadratic matrix is the main
1922  // diagonal
1923  return val[cols->rowstart[i]];
1924 }
1925 
1926 
1927 
1928 template <typename number>
1929 template <typename ForwardIterator>
1930 void
1931 SparseMatrix<number>::copy_from (const ForwardIterator begin,
1932  const ForwardIterator end)
1933 {
1934  Assert (static_cast<size_type>(std::distance (begin, end)) == m(),
1935  ExcIteratorRange (std::distance (begin, end), m()));
1936 
1937  // for use in the inner loop, we define a typedef to the type of the inner
1938  // iterators
1939  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1940  size_type row=0;
1941  for (ForwardIterator i=begin; i!=end; ++i, ++row)
1942  {
1943  const inner_iterator end_of_row = i->end();
1944  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1945  // write entries
1946  set (row, j->first, j->second);
1947  };
1948 }
1949 
1950 
1951 //---------------------------------------------------------------------------
1952 
1953 
1954 namespace SparseMatrixIterators
1955 {
1956  template <typename number>
1957  inline
1959  Accessor (const MatrixType *matrix,
1960  const std::size_t index_within_matrix)
1961  :
1962  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1963  index_within_matrix),
1964  matrix (matrix)
1965  {}
1966 
1967 
1968 
1969  template <typename number>
1970  inline
1971  Accessor<number,true>::
1972  Accessor (const MatrixType *matrix)
1973  :
1974  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1975  matrix (matrix)
1976  {}
1977 
1978 
1979 
1980  template <typename number>
1981  inline
1982  Accessor<number,true>::
1984  :
1985  SparsityPatternIterators::Accessor (a),
1986  matrix (&a.get_matrix())
1987  {}
1988 
1989 
1990 
1991  template <typename number>
1992  inline
1993  number
1994  Accessor<number, true>::value () const
1995  {
1996  AssertIndexRange(index_within_sparsity, matrix->n_nonzero_elements());
1997  return matrix->val[index_within_sparsity];
1998  }
1999 
2000 
2001 
2002  template <typename number>
2003  inline
2004  const typename Accessor<number, true>::MatrixType &
2005  Accessor<number, true>::get_matrix () const
2006  {
2007  return *matrix;
2008  }
2009 
2010 
2011 
2012  template <typename number>
2013  inline
2014  Accessor<number, false>::Reference::Reference (
2015  const Accessor *accessor,
2016  const bool)
2017  :
2018  accessor (accessor)
2019  {}
2020 
2021 
2022  template <typename number>
2023  inline
2024  Accessor<number, false>::Reference::operator number() const
2025  {
2026  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2027  return accessor->matrix->val[accessor->index_within_sparsity];
2028  }
2029 
2030 
2031 
2032  template <typename number>
2033  inline
2034  const typename Accessor<number, false>::Reference &
2035  Accessor<number, false>::Reference::operator = (const number n) const
2036  {
2037  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2038  accessor->matrix->val[accessor->index_within_sparsity] = n;
2039  return *this;
2040  }
2041 
2042 
2043 
2044  template <typename number>
2045  inline
2046  const typename Accessor<number, false>::Reference &
2047  Accessor<number, false>::Reference::operator += (const number n) const
2048  {
2049  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2050  accessor->matrix->val[accessor->index_within_sparsity] += n;
2051  return *this;
2052  }
2053 
2054 
2055 
2056  template <typename number>
2057  inline
2058  const typename Accessor<number, false>::Reference &
2059  Accessor<number, false>::Reference::operator -= (const number n) const
2060  {
2061  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2062  accessor->matrix->val[accessor->index_within_sparsity] -= n;
2063  return *this;
2064  }
2065 
2066 
2067 
2068  template <typename number>
2069  inline
2070  const typename Accessor<number, false>::Reference &
2071  Accessor<number, false>::Reference::operator *= (const number n) const
2072  {
2073  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2074  accessor->matrix->val[accessor->index_within_sparsity] *= n;
2075  return *this;
2076  }
2077 
2078 
2079 
2080  template <typename number>
2081  inline
2082  const typename Accessor<number, false>::Reference &
2083  Accessor<number, false>::Reference::operator /= (const number n) const
2084  {
2085  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2086  accessor->matrix->val[accessor->index_within_sparsity] /= n;
2087  return *this;
2088  }
2089 
2090 
2091 
2092  template <typename number>
2093  inline
2094  Accessor<number,false>::
2095  Accessor (MatrixType *matrix,
2096  const std::size_t index)
2097  :
2098  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2099  index),
2100  matrix (matrix)
2101  {}
2102 
2103 
2104 
2105  template <typename number>
2106  inline
2107  Accessor<number,false>::
2108  Accessor (MatrixType *matrix)
2109  :
2110  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
2111  matrix (matrix)
2112  {}
2113 
2114 
2115 
2116  template <typename number>
2117  inline
2118  typename Accessor<number, false>::Reference
2119  Accessor<number, false>::value() const
2120  {
2121  return Reference(this,true);
2122  }
2123 
2124 
2125 
2126 
2127  template <typename number>
2128  inline
2129  typename Accessor<number, false>::MatrixType &
2130  Accessor<number, false>::get_matrix () const
2131  {
2132  return *matrix;
2133  }
2134 
2135 
2136 
2137  template <typename number, bool Constness>
2138  inline
2139  Iterator<number, Constness>::
2140  Iterator (MatrixType *matrix,
2141  const std::size_t index)
2142  :
2143  accessor(matrix, index)
2144  {}
2145 
2146 
2147 
2148  template <typename number, bool Constness>
2149  inline
2150  Iterator<number, Constness>::
2151  Iterator (MatrixType *matrix)
2152  :
2153  accessor(matrix)
2154  {}
2155 
2156 
2157 
2158  template <typename number, bool Constness>
2159  inline
2160  Iterator<number, Constness>::
2162  :
2163  accessor(*i)
2164  {}
2165 
2166 
2167 
2168  template <typename number, bool Constness>
2169  inline
2170  Iterator<number, Constness> &
2171  Iterator<number,Constness>::operator++ ()
2172  {
2173  accessor.advance ();
2174  return *this;
2175  }
2176 
2177 
2178  template <typename number, bool Constness>
2179  inline
2180  Iterator<number,Constness>
2181  Iterator<number,Constness>::operator++ (int)
2182  {
2183  const Iterator iter = *this;
2184  accessor.advance ();
2185  return iter;
2186  }
2187 
2188 
2189  template <typename number, bool Constness>
2190  inline
2191  const Accessor<number,Constness> &
2192  Iterator<number,Constness>::operator* () const
2193  {
2194  return accessor;
2195  }
2196 
2197 
2198  template <typename number, bool Constness>
2199  inline
2200  const Accessor<number,Constness> *
2201  Iterator<number,Constness>::operator-> () const
2202  {
2203  return &accessor;
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  return (accessor == other.accessor);
2214  }
2215 
2216 
2217  template <typename number, bool Constness>
2218  inline
2219  bool
2220  Iterator<number,Constness>::
2221  operator != (const Iterator &other) const
2222  {
2223  return ! (*this == other);
2224  }
2225 
2226 
2227  template <typename number, bool Constness>
2228  inline
2229  bool
2230  Iterator<number,Constness>::
2231  operator < (const Iterator &other) const
2232  {
2233  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2234  ExcInternalError());
2235 
2236  return (accessor < other.accessor);
2237  }
2238 
2239 
2240  template <typename number, bool Constness>
2241  inline
2242  bool
2243  Iterator<number,Constness>::
2244  operator > (const Iterator &other) const
2245  {
2246  return (other < *this);
2247  }
2248 
2249 
2250  template <typename number, bool Constness>
2251  inline
2252  int
2253  Iterator<number,Constness>::
2254  operator - (const Iterator &other) const
2255  {
2256  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2257  ExcInternalError());
2258 
2259  return (*this)->index_within_sparsity - other->index_within_sparsity;
2260  }
2261 
2262 
2263 
2264  template <typename number, bool Constness>
2265  inline
2266  Iterator<number,Constness>
2267  Iterator<number,Constness>::
2268  operator + (const size_type n) const
2269  {
2270  Iterator x = *this;
2271  for (size_type i=0; i<n; ++i)
2272  ++x;
2273 
2274  return x;
2275  }
2276 
2277 }
2278 
2279 
2280 
2281 template <typename number>
2282 inline
2285 {
2286  return const_iterator(this, 0);
2287 }
2288 
2289 
2290 template <typename number>
2291 inline
2294 {
2295  return const_iterator(this);
2296 }
2297 
2298 
2299 template <typename number>
2300 inline
2303 {
2304  return iterator (this, 0);
2305 }
2306 
2307 
2308 template <typename number>
2309 inline
2312 {
2313  return iterator(this, cols->rowstart[cols->rows]);
2314 }
2315 
2316 
2317 template <typename number>
2318 inline
2320 SparseMatrix<number>::begin (const size_type r) const
2321 {
2322  Assert (r<m(), ExcIndexRange(r,0,m()));
2323 
2324  return const_iterator(this, cols->rowstart[r]);
2325 }
2326 
2327 
2328 
2329 template <typename number>
2330 inline
2332 SparseMatrix<number>::end (const size_type r) const
2333 {
2334  Assert (r<m(), ExcIndexRange(r,0,m()));
2335 
2336  return const_iterator(this, cols->rowstart[r+1]);
2337 }
2338 
2339 
2340 
2341 template <typename number>
2342 inline
2344 SparseMatrix<number>::begin (const size_type r)
2345 {
2346  Assert (r<m(), ExcIndexRange(r,0,m()));
2347 
2348  return iterator(this, cols->rowstart[r]);
2349 }
2350 
2351 
2352 
2353 template <typename number>
2354 inline
2356 SparseMatrix<number>::end (const size_type r)
2357 {
2358  Assert (r<m(), ExcIndexRange(r,0,m()));
2359 
2360  return iterator(this, cols->rowstart[r+1]);
2361 }
2362 
2363 
2364 
2365 template <typename number>
2366 template <class StreamType>
2367 inline
2368 void SparseMatrix<number>::print (StreamType &out,
2369  const bool across,
2370  const bool diagonal_first) const
2371 {
2372  Assert (cols != nullptr, ExcNotInitialized());
2373  Assert (val != nullptr, ExcNotInitialized());
2374 
2375  bool hanging_diagonal = false;
2376  number diagonal = number();
2377 
2378  for (size_type i=0; i<cols->rows; ++i)
2379  {
2380  for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
2381  {
2382  if (!diagonal_first && i == cols->colnums[j])
2383  {
2384  diagonal = val[j];
2385  hanging_diagonal = true;
2386  }
2387  else
2388  {
2389  if (hanging_diagonal && cols->colnums[j]>i)
2390  {
2391  if (across)
2392  out << ' ' << i << ',' << i << ':' << diagonal;
2393  else
2394  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2395  hanging_diagonal = false;
2396  }
2397  if (across)
2398  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2399  else
2400  out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << std::endl;
2401  }
2402  }
2403  if (hanging_diagonal)
2404  {
2405  if (across)
2406  out << ' ' << i << ',' << i << ':' << diagonal;
2407  else
2408  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2409  hanging_diagonal = false;
2410  }
2411  }
2412  if (across)
2413  out << std::endl;
2414 }
2415 
2416 
2417 template <typename number>
2418 inline
2419 void
2421 prepare_add()
2422 {
2423  //nothing to do here
2424 }
2425 
2426 
2427 
2428 template <typename number>
2429 inline
2430 void
2432 prepare_set()
2433 {
2434  //nothing to do here
2435 }
2436 
2437 #endif // DOXYGEN
2438 
2439 
2440 /*---------------------------- sparse_matrix.h ---------------------------*/
2441 
2442 DEAL_II_NAMESPACE_CLOSE
2443 
2444 #endif
2445 /*---------------------------- 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:617
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:1237
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
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:349
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:594
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:1253
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