Reference documentation for deal.II version 8.1.0
sparse_matrix.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_matrix.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__sparse_matrix_h
18 #define __deal2__sparse_matrix_h
19 
20 
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/lac/sparsity_pattern.h>
25 #include <deal.II/lac/identity_matrix.h>
26 #include <deal.II/lac/exceptions.h>
27 // Included for VectorOperation
28 #include <deal.II/lac/vector.h>
29 
31 
32 template <typename number> class Vector;
33 template <typename number> class FullMatrix;
34 template <typename Matrix> class BlockMatrixBase;
35 template <typename number> class SparseILU;
36 
37 #ifdef DEAL_II_WITH_TRILINOS
38 namespace TrilinosWrappers
39 {
40  class SparseMatrix;
41 }
42 #endif
43 
54 {
59 
60  // forward declaration
61  template <typename number, bool Constness>
62  class Iterator;
63 
74  template <typename number, bool Constness>
76  {
77  public:
81  number value() const;
82 
86  number &value();
87 
92  const SparseMatrix<number> &get_matrix () const;
93  };
94 
95 
96 
103  template <typename number>
104  class Accessor<number,true> : public SparsityPatternIterators::Accessor
105  {
106  public:
112 
119  Accessor (MatrixType *matrix,
120  const size_type row,
121  const size_type index) DEAL_II_DEPRECATED;
122 
126  Accessor (MatrixType *matrix,
127  const std::size_t index_within_matrix);
128 
132  Accessor (MatrixType *matrix);
133 
138 
142  number value() const;
143 
148  MatrixType &get_matrix () const;
149 
150  private:
154  MatrixType *matrix;
155 
160 
164  template <typename, bool>
165  friend class Iterator;
166  };
167 
168 
175  template <typename number>
176  class Accessor<number,false> : public SparsityPatternIterators::Accessor
177  {
178  private:
203  class Reference
204  {
205  public:
210  Reference (const Accessor *accessor,
211  const bool dummy);
212 
216  operator number () const;
217 
221  const Reference &operator = (const number n) const;
222 
226  const Reference &operator += (const number n) const;
227 
231  const Reference &operator -= (const number n) const;
232 
236  const Reference &operator *= (const number n) const;
237 
241  const Reference &operator /= (const number n) const;
242 
243  private:
249  };
250 
251  public:
257 
261  Accessor (MatrixType *matrix,
262  const size_type row,
263  const size_type index);
264 
268  Accessor (MatrixType *matrix,
269  const std::size_t index);
270 
274  Accessor (MatrixType *matrix);
275 
279  Reference value() const;
280 
285  MatrixType &get_matrix () const;
286 
287  private:
291  MatrixType *matrix;
292 
297 
301  template <typename, bool>
302  friend class Iterator;
303 
308  };
309 
310 
311 
342  template <typename number, bool Constness>
343  class Iterator
344  {
345  public:
349  typedef
352 
357  typedef
359 
368  const size_type row,
369  const size_type index) DEAL_II_DEPRECATED;
370 
376  const std::size_t index_within_matrix);
377 
382 
388 
392  Iterator &operator++ ();
393 
397  Iterator operator++ (int);
398 
402  const Accessor<number,Constness> &operator* () const;
403 
407  const Accessor<number,Constness> *operator-> () const;
408 
412  bool operator == (const Iterator &) const;
413 
417  bool operator != (const Iterator &) const;
418 
426  bool operator < (const Iterator &) const;
427 
432  bool operator > (const Iterator &) const;
433 
440  int operator - (const Iterator &p) const;
441 
445  Iterator operator + (const size_type n) const;
446 
447  private:
452  };
453 
454 }
455 
461 //TODO: Add multithreading to the other vmult functions.
462 
492 template <typename number>
493 class SparseMatrix : public virtual Subscriptor
494 {
495 public:
500 
504  typedef number value_type;
505 
516 
522  typedef
525 
532  typedef
535 
542  struct Traits
543  {
548  static const bool zero_addition_can_be_elided = true;
549  };
550 
565  SparseMatrix ();
566 
575  SparseMatrix (const SparseMatrix &);
576 
590  explicit SparseMatrix (const SparsityPattern &sparsity);
591 
598  SparseMatrix (const SparsityPattern &sparsity,
599  const IdentityMatrix &id);
600 
605  virtual ~SparseMatrix ();
606 
616  SparseMatrix<number> &operator = (const SparseMatrix<number> &);
617 
625  operator= (const IdentityMatrix &id);
626 
636  SparseMatrix &operator = (const double d);
637 
651  virtual void reinit (const SparsityPattern &sparsity);
652 
658  virtual void clear ();
660 
668  bool empty () const;
669 
674  size_type m () const;
675 
680  size_type n () const;
681 
685  size_type get_row_length (const size_type row) const;
686 
692  size_type n_nonzero_elements () const;
693 
703  size_type n_actually_nonzero_elements (const double threshold = 0.) const;
704 
713  const SparsityPattern &get_sparsity_pattern () const;
714 
719  std::size_t memory_consumption () const;
720 
724  void compress (::VectorOperation::values);
725 
727 
736  void set (const size_type i,
737  const size_type j,
738  const number value);
739 
755  template <typename number2>
756  void set (const std::vector<size_type> &indices,
757  const FullMatrix<number2> &full_matrix,
758  const bool elide_zero_values = false);
759 
765  template <typename number2>
766  void set (const std::vector<size_type> &row_indices,
767  const std::vector<size_type> &col_indices,
768  const FullMatrix<number2> &full_matrix,
769  const bool elide_zero_values = false);
770 
781  template <typename number2>
782  void set (const size_type row,
783  const std::vector<size_type> &col_indices,
784  const std::vector<number2> &values,
785  const bool elide_zero_values = false);
786 
796  template <typename number2>
797  void set (const size_type row,
798  const size_type n_cols,
799  const size_type *col_indices,
800  const number2 *values,
801  const bool elide_zero_values = false);
802 
808  void add (const size_type i,
809  const size_type j,
810  const number value);
811 
826  template <typename number2>
827  void add (const std::vector<size_type> &indices,
828  const FullMatrix<number2> &full_matrix,
829  const bool elide_zero_values = true);
830 
836  template <typename number2>
837  void add (const std::vector<size_type> &row_indices,
838  const std::vector<size_type> &col_indices,
839  const FullMatrix<number2> &full_matrix,
840  const bool elide_zero_values = true);
841 
851  template <typename number2>
852  void add (const size_type row,
853  const std::vector<size_type> &col_indices,
854  const std::vector<number2> &values,
855  const bool elide_zero_values = true);
856 
866  template <typename number2>
867  void add (const size_type row,
868  const size_type n_cols,
869  const size_type *col_indices,
870  const number2 *values,
871  const bool elide_zero_values = true,
872  const bool col_indices_are_sorted = false);
873 
877  SparseMatrix &operator *= (const number factor);
878 
882  SparseMatrix &operator /= (const number factor);
883 
896  void symmetrize ();
897 
911  template <typename somenumber>
913  copy_from (const SparseMatrix<somenumber> &source);
914 
931  template <typename ForwardIterator>
932  void copy_from (const ForwardIterator begin,
933  const ForwardIterator end);
934 
940  template <typename somenumber>
941  void copy_from (const FullMatrix<somenumber> &matrix);
942 
943 #ifdef DEAL_II_WITH_TRILINOS
944 
954  copy_from (const TrilinosWrappers::SparseMatrix &matrix);
955 #endif
956 
968  template <typename somenumber>
969  void add (const number factor,
971 
973 
977 
991  number operator () (const size_type i,
992  const size_type j) const;
993 
1006  number el (const size_type i,
1007  const size_type j) const;
1008 
1019  number diag_element (const size_type i) const;
1020 
1025  number &diag_element (const size_type i);
1026 
1035  number raw_entry (const size_type row,
1036  const size_type index) const DEAL_II_DEPRECATED;
1037 
1051  number global_entry (const size_type i) const DEAL_II_DEPRECATED;
1052 
1058  number &global_entry (const size_type i) DEAL_II_DEPRECATED;
1059 
1061 
1078  template <class OutVector, class InVector>
1079  void vmult (OutVector &dst,
1080  const InVector &src) const;
1081 
1096  template <class OutVector, class InVector>
1097  void Tvmult (OutVector &dst,
1098  const InVector &src) const;
1099 
1113  template <class OutVector, class InVector>
1114  void vmult_add (OutVector &dst,
1115  const InVector &src) const;
1116 
1131  template <class OutVector, class InVector>
1132  void Tvmult_add (OutVector &dst,
1133  const InVector &src) const;
1134 
1150  template <typename somenumber>
1151  somenumber matrix_norm_square (const Vector<somenumber> &v) const;
1152 
1156  template <typename somenumber>
1157  somenumber matrix_scalar_product (const Vector<somenumber> &u,
1158  const Vector<somenumber> &v) const;
1159 
1167  template <typename somenumber>
1168  somenumber residual (Vector<somenumber> &dst,
1169  const Vector<somenumber> &x,
1170  const Vector<somenumber> &b) const;
1171 
1195  template <typename numberB, typename numberC>
1196  void mmult (SparseMatrix<numberC> &C,
1197  const SparseMatrix<numberB> &B,
1198  const Vector<number> &V = Vector<number>(),
1199  const bool rebuild_sparsity_pattern = true) const;
1200 
1225  template <typename numberB, typename numberC>
1226  void Tmmult (SparseMatrix<numberC> &C,
1227  const SparseMatrix<numberB> &B,
1228  const Vector<number> &V = Vector<number>(),
1229  const bool rebuild_sparsity_pattern = true) const;
1230 
1232 
1236 
1244  real_type l1_norm () const;
1245 
1254  real_type linfty_norm () const;
1255 
1260  real_type frobenius_norm () const;
1262 
1266 
1272  template <typename somenumber>
1273  void precondition_Jacobi (Vector<somenumber> &dst,
1274  const Vector<somenumber> &src,
1275  const number omega = 1.) const;
1276 
1283  template <typename somenumber>
1284  void precondition_SSOR (Vector<somenumber> &dst,
1285  const Vector<somenumber> &src,
1286  const number omega = 1.,
1287  const std::vector<std::size_t> &pos_right_of_diagonal=std::vector<std::size_t>()) const;
1288 
1292  template <typename somenumber>
1293  void precondition_SOR (Vector<somenumber> &dst,
1294  const Vector<somenumber> &src,
1295  const number om = 1.) const;
1296 
1300  template <typename somenumber>
1301  void precondition_TSOR (Vector<somenumber> &dst,
1302  const Vector<somenumber> &src,
1303  const number om = 1.) const;
1304 
1310  template <typename somenumber>
1311  void SSOR (Vector<somenumber> &v,
1312  const number omega = 1.) const;
1313 
1318  template <typename somenumber>
1319  void SOR (Vector<somenumber> &v,
1320  const number om = 1.) const;
1321 
1326  template <typename somenumber>
1327  void TSOR (Vector<somenumber> &v,
1328  const number om = 1.) const;
1329 
1340  template <typename somenumber>
1341  void PSOR (Vector<somenumber> &v,
1342  const std::vector<size_type> &permutation,
1343  const std::vector<size_type> &inverse_permutation,
1344  const number om = 1.) const;
1345 
1356  template <typename somenumber>
1357  void TPSOR (Vector<somenumber> &v,
1358  const std::vector<size_type> &permutation,
1359  const std::vector<size_type> &inverse_permutation,
1360  const number om = 1.) const;
1361 
1367  template <typename somenumber>
1368  void Jacobi_step (Vector<somenumber> &v,
1369  const Vector<somenumber> &b,
1370  const number om = 1.) const;
1371 
1376  template <typename somenumber>
1377  void SOR_step (Vector<somenumber> &v,
1378  const Vector<somenumber> &b,
1379  const number om = 1.) const;
1380 
1385  template <typename somenumber>
1386  void TSOR_step (Vector<somenumber> &v,
1387  const Vector<somenumber> &b,
1388  const number om = 1.) const;
1389 
1394  template <typename somenumber>
1395  void SSOR_step (Vector<somenumber> &v,
1396  const Vector<somenumber> &b,
1397  const number om = 1.) const;
1399 
1403 
1411  const_iterator begin () const;
1412 
1416  const_iterator end () const;
1417 
1425  iterator begin ();
1426 
1430  iterator end ();
1431 
1444  const_iterator begin (const size_type r) const;
1445 
1455  const_iterator end (const size_type r) const;
1456 
1469  iterator begin (const size_type r);
1470 
1480  iterator end (const size_type r);
1482 
1486 
1498  template <class STREAM>
1499  void print (STREAM &out,
1500  const bool across = false,
1501  const bool diagonal_first = true) const;
1502 
1523  void print_formatted (std::ostream &out,
1524  const unsigned int precision = 3,
1525  const bool scientific = true,
1526  const unsigned int width = 0,
1527  const char *zero_string = " ",
1528  const double denominator = 1.) const;
1529 
1535  void print_pattern(std::ostream &out,
1536  const double threshold = 0.) const;
1537 
1548  void block_write (std::ostream &out) const;
1549 
1566  void block_read (std::istream &in);
1568 
1574  DeclException2 (ExcInvalidIndex,
1575  int, int,
1576  << "The entry with index <" << arg1 << ',' << arg2
1577  << "> does not exist.");
1581  DeclException1 (ExcInvalidIndex1,
1582  int,
1583  << "The index " << arg1 << " is not in the allowed range.");
1587  DeclException0 (ExcDifferentSparsityPatterns);
1591  DeclException2 (ExcIteratorRange,
1592  int, int,
1593  << "The iterators denote a range of " << arg1
1594  << " elements, but the given number of rows was " << arg2);
1598  DeclException0 (ExcSourceEqualsDestination);
1600 
1601 protected:
1612  void prepare_add();
1613 
1618  void prepare_set();
1619 
1620 private:
1627 
1634  number *val;
1635 
1642  std::size_t max_len;
1643 
1644  // make all other sparse matrices friends
1645  template <typename somenumber> friend class SparseMatrix;
1646  template <typename somenumber> friend class SparseLUDecomposition;
1647  template <typename> friend class SparseILU;
1648 
1652  template <typename> friend class BlockMatrixBase;
1653 
1658  template <typename,bool> friend class SparseMatrixIterators::Iterator;
1659  template <typename,bool> friend class SparseMatrixIterators::Accessor;
1660 };
1661 
1662 #ifndef DOXYGEN
1663 /*---------------------- Inline functions -----------------------------------*/
1664 
1665 
1666 
1667 template <typename number>
1668 inline
1670 {
1671  Assert (cols != 0, ExcNotInitialized());
1672  return cols->rows;
1673 }
1674 
1675 
1676 template <typename number>
1677 inline
1679 {
1680  Assert (cols != 0, ExcNotInitialized());
1681  return cols->cols;
1682 }
1683 
1684 
1685 // Inline the set() and add() functions, since they will be called frequently.
1686 template <typename number>
1687 inline
1688 void
1689 SparseMatrix<number>::set (const size_type i,
1690  const size_type j,
1691  const number value)
1692 {
1694 
1695  const size_type index = cols->operator()(i, j);
1696 
1697  // it is allowed to set elements of the matrix that are not part of the
1698  // sparsity pattern, if the value to which we set it is zero
1699  if (index == SparsityPattern::invalid_entry)
1700  {
1701  Assert ((index != SparsityPattern::invalid_entry) ||
1702  (value == 0.),
1703  ExcInvalidIndex(i, j));
1704  return;
1705  }
1706 
1707  val[index] = value;
1708 }
1709 
1710 
1711 
1712 template <typename number>
1713 template <typename number2>
1714 inline
1715 void
1716 SparseMatrix<number>::set (const std::vector<size_type> &indices,
1717  const FullMatrix<number2> &values,
1718  const bool elide_zero_values)
1719 {
1720  Assert (indices.size() == values.m(),
1721  ExcDimensionMismatch(indices.size(), values.m()));
1722  Assert (values.m() == values.n(), ExcNotQuadratic());
1723 
1724  for (size_type i=0; i<indices.size(); ++i)
1725  set (indices[i], indices.size(), &indices[0], &values(i,0),
1726  elide_zero_values);
1727 }
1728 
1729 
1730 
1731 template <typename number>
1732 template <typename number2>
1733 inline
1734 void
1735 SparseMatrix<number>::set (const std::vector<size_type> &row_indices,
1736  const std::vector<size_type> &col_indices,
1737  const FullMatrix<number2> &values,
1738  const bool elide_zero_values)
1739 {
1740  Assert (row_indices.size() == values.m(),
1741  ExcDimensionMismatch(row_indices.size(), values.m()));
1742  Assert (col_indices.size() == values.n(),
1743  ExcDimensionMismatch(col_indices.size(), values.n()));
1744 
1745  for (size_type i=0; i<row_indices.size(); ++i)
1746  set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1747  elide_zero_values);
1748 }
1749 
1750 
1751 
1752 template <typename number>
1753 template <typename number2>
1754 inline
1755 void
1756 SparseMatrix<number>::set (const size_type row,
1757  const std::vector<size_type> &col_indices,
1758  const std::vector<number2> &values,
1759  const bool elide_zero_values)
1760 {
1761  Assert (col_indices.size() == values.size(),
1762  ExcDimensionMismatch(col_indices.size(), values.size()));
1763 
1764  set (row, col_indices.size(), &col_indices[0], &values[0],
1765  elide_zero_values);
1766 }
1767 
1768 
1769 
1770 template <typename number>
1771 inline
1772 void
1773 SparseMatrix<number>::add (const size_type i,
1774  const size_type j,
1775  const number value)
1776 {
1778 
1779  if (value == 0)
1780  return;
1781 
1782  const size_type index = cols->operator()(i, j);
1783 
1784  // it is allowed to add elements to the matrix that are not part of the
1785  // sparsity pattern, if the value to which we set it is zero
1786  if (index == SparsityPattern::invalid_entry)
1787  {
1788  Assert ((index != SparsityPattern::invalid_entry) ||
1789  (value == 0.),
1790  ExcInvalidIndex(i, j));
1791  return;
1792  }
1793 
1794  val[index] += value;
1795 }
1796 
1797 
1798 
1799 template <typename number>
1800 template <typename number2>
1801 inline
1802 void
1803 SparseMatrix<number>::add (const std::vector<size_type> &indices,
1804  const FullMatrix<number2> &values,
1805  const bool elide_zero_values)
1806 {
1807  Assert (indices.size() == values.m(),
1808  ExcDimensionMismatch(indices.size(), values.m()));
1809  Assert (values.m() == values.n(), ExcNotQuadratic());
1810 
1811  for (size_type i=0; i<indices.size(); ++i)
1812  add (indices[i], indices.size(), &indices[0], &values(i,0),
1813  elide_zero_values);
1814 }
1815 
1816 
1817 
1818 template <typename number>
1819 template <typename number2>
1820 inline
1821 void
1822 SparseMatrix<number>::add (const std::vector<size_type> &row_indices,
1823  const std::vector<size_type> &col_indices,
1824  const FullMatrix<number2> &values,
1825  const bool elide_zero_values)
1826 {
1827  Assert (row_indices.size() == values.m(),
1828  ExcDimensionMismatch(row_indices.size(), values.m()));
1829  Assert (col_indices.size() == values.n(),
1830  ExcDimensionMismatch(col_indices.size(), values.n()));
1831 
1832  for (size_type i=0; i<row_indices.size(); ++i)
1833  add (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
1834  elide_zero_values);
1835 }
1836 
1837 
1838 
1839 template <typename number>
1840 template <typename number2>
1841 inline
1842 void
1843 SparseMatrix<number>::add (const size_type row,
1844  const std::vector<size_type> &col_indices,
1845  const std::vector<number2> &values,
1846  const bool elide_zero_values)
1847 {
1848  Assert (col_indices.size() == values.size(),
1849  ExcDimensionMismatch(col_indices.size(), values.size()));
1850 
1851  add (row, col_indices.size(), &col_indices[0], &values[0],
1852  elide_zero_values);
1853 }
1854 
1855 
1856 
1857 template <typename number>
1858 inline
1860 SparseMatrix<number>::operator *= (const number factor)
1861 {
1862  Assert (cols != 0, ExcNotInitialized());
1863  Assert (val != 0, ExcNotInitialized());
1864 
1865  number *val_ptr = &val[0];
1866  const number *const end_ptr = &val[cols->n_nonzero_elements()];
1867 
1868  while (val_ptr != end_ptr)
1869  *val_ptr++ *= factor;
1870 
1871  return *this;
1872 }
1873 
1874 
1875 
1876 template <typename number>
1877 inline
1879 SparseMatrix<number>::operator /= (const number factor)
1880 {
1881  Assert (cols != 0, ExcNotInitialized());
1882  Assert (val != 0, ExcNotInitialized());
1883  Assert (factor !=0, ExcDivideByZero());
1884 
1885  const number factor_inv = 1. / factor;
1886 
1887  number *val_ptr = &val[0];
1888  const number *const end_ptr = &val[cols->n_nonzero_elements()];
1889 
1890  while (val_ptr != end_ptr)
1891  *val_ptr++ *= factor_inv;
1892 
1893  return *this;
1894 }
1895 
1896 
1897 
1898 template <typename number>
1899 inline
1900 number SparseMatrix<number>::operator () (const size_type i,
1901  const size_type j) const
1902 {
1903  Assert (cols != 0, ExcNotInitialized());
1904  Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
1905  ExcInvalidIndex(i,j));
1906  return val[cols->operator()(i,j)];
1907 }
1908 
1909 
1910 
1911 template <typename number>
1912 inline
1913 number SparseMatrix<number>::el (const size_type i,
1914  const size_type j) const
1915 {
1916  Assert (cols != 0, ExcNotInitialized());
1917  const size_type index = cols->operator()(i,j);
1918 
1919  if (index != SparsityPattern::invalid_entry)
1920  return val[index];
1921  else
1922  return 0;
1923 }
1924 
1925 
1926 
1927 template <typename number>
1928 inline
1929 number SparseMatrix<number>::diag_element (const size_type i) const
1930 {
1931  Assert (cols != 0, ExcNotInitialized());
1932  Assert (m() == n(), ExcNotQuadratic());
1933  Assert (i<m(), ExcInvalidIndex1(i));
1934 
1935  // Use that the first element in each row of a quadratic matrix is the main
1936  // diagonal
1937  return val[cols->rowstart[i]];
1938 }
1939 
1940 
1941 
1942 template <typename number>
1943 inline
1944 number &SparseMatrix<number>::diag_element (const size_type i)
1945 {
1946  Assert (cols != 0, ExcNotInitialized());
1947  Assert (m() == n(), ExcNotQuadratic());
1948  Assert (i<m(), ExcInvalidIndex1(i));
1949 
1950  // Use that the first element in each row of a quadratic matrix is the main
1951  // diagonal
1952  return val[cols->rowstart[i]];
1953 }
1954 
1955 
1956 
1957 template <typename number>
1958 inline
1959 number
1960 SparseMatrix<number>::raw_entry (const size_type row,
1961  const size_type index) const
1962 {
1963  Assert(row<cols->rows, ExcIndexRange(row,0,cols->rows));
1964  Assert(index<cols->row_length(row),
1965  ExcIndexRange(index,0,cols->row_length(row)));
1966 
1967  // this function will soon go away.
1968  return val[cols->rowstart[row]+index];
1969 }
1970 
1971 
1972 
1973 template <typename number>
1974 inline
1975 number SparseMatrix<number>::global_entry (const size_type j) const
1976 {
1977  Assert (cols != 0, ExcNotInitialized());
1978  Assert (j < cols->n_nonzero_elements(),
1979  ExcIndexRange (j, 0, cols->n_nonzero_elements()));
1980 
1981  return val[j];
1982 }
1983 
1984 
1985 
1986 template <typename number>
1987 inline
1988 number &SparseMatrix<number>::global_entry (const size_type j)
1989 {
1990  Assert (cols != 0, ExcNotInitialized());
1991  Assert (j < cols->n_nonzero_elements(),
1992  ExcIndexRange (j, 0, cols->n_nonzero_elements()));
1993 
1994  return val[j];
1995 }
1996 
1997 
1998 
1999 template <typename number>
2000 template <typename ForwardIterator>
2001 void
2002 SparseMatrix<number>::copy_from (const ForwardIterator begin,
2003  const ForwardIterator end)
2004 {
2005  Assert (static_cast<size_type>(std::distance (begin, end)) == m(),
2006  ExcIteratorRange (std::distance (begin, end), m()));
2007 
2008  // for use in the inner loop, we define a typedef to the type of the inner
2009  // iterators
2010  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
2011  size_type row=0;
2012  for (ForwardIterator i=begin; i!=end; ++i, ++row)
2013  {
2014  const inner_iterator end_of_row = i->end();
2015  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
2016  // write entries
2017  set (row, j->first, j->second);
2018  };
2019 }
2020 
2021 
2022 //---------------------------------------------------------------------------
2023 
2024 
2025 namespace SparseMatrixIterators
2026 {
2027  template <typename number>
2028  inline
2029  Accessor<number,true>::
2030  Accessor (const MatrixType *matrix,
2031  const size_type row,
2032  const size_type index)
2033  :
2034  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2035  row, index),
2036  matrix (matrix)
2037  {}
2038 
2039 
2040 
2041  template <typename number>
2042  inline
2043  Accessor<number,true>::
2044  Accessor (const MatrixType *matrix,
2045  const std::size_t index_within_matrix)
2046  :
2047  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2048  index_within_matrix),
2049  matrix (matrix)
2050  {}
2051 
2052 
2053 
2054  template <typename number>
2055  inline
2056  Accessor<number,true>::
2057  Accessor (const MatrixType *matrix)
2058  :
2059  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
2060  matrix (matrix)
2061  {}
2062 
2063 
2064 
2065  template <typename number>
2066  inline
2067  Accessor<number,true>::
2069  :
2070  SparsityPatternIterators::Accessor (a),
2071  matrix (&a.get_matrix())
2072  {}
2073 
2074 
2075 
2076  template <typename number>
2077  inline
2078  number
2079  Accessor<number, true>::value () const
2080  {
2081  AssertIndexRange(index_within_sparsity, matrix->n_nonzero_elements());
2082  return matrix->val[index_within_sparsity];
2083  }
2084 
2085 
2086 
2087  template <typename number>
2088  inline
2089  typename Accessor<number, true>::MatrixType &
2090  Accessor<number, true>::get_matrix () const
2091  {
2092  return *matrix;
2093  }
2094 
2095 
2096 
2097  template <typename number>
2098  inline
2099  Accessor<number, false>::Reference::Reference (
2100  const Accessor *accessor,
2101  const bool)
2102  :
2103  accessor (accessor)
2104  {}
2105 
2106 
2107  template <typename number>
2108  inline
2109  Accessor<number, false>::Reference::operator number() const
2110  {
2111  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2112  return accessor->matrix->val[accessor->index_within_sparsity];
2113  }
2114 
2115 
2116 
2117  template <typename number>
2118  inline
2119  const typename Accessor<number, false>::Reference &
2120  Accessor<number, false>::Reference::operator = (const number n) const
2121  {
2122  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2123  accessor->matrix->val[accessor->index_within_sparsity] = n;
2124  return *this;
2125  }
2126 
2127 
2128 
2129  template <typename number>
2130  inline
2131  const typename Accessor<number, false>::Reference &
2132  Accessor<number, false>::Reference::operator += (const number n) const
2133  {
2134  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2135  accessor->matrix->val[accessor->index_within_sparsity] += n;
2136  return *this;
2137  }
2138 
2139 
2140 
2141  template <typename number>
2142  inline
2143  const typename Accessor<number, false>::Reference &
2144  Accessor<number, false>::Reference::operator -= (const number n) const
2145  {
2146  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2147  accessor->matrix->val[accessor->index_within_sparsity] -= n;
2148  return *this;
2149  }
2150 
2151 
2152 
2153  template <typename number>
2154  inline
2155  const typename Accessor<number, false>::Reference &
2156  Accessor<number, false>::Reference::operator *= (const number n) const
2157  {
2158  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2159  accessor->matrix->val[accessor->index_within_sparsity] *= n;
2160  return *this;
2161  }
2162 
2163 
2164 
2165  template <typename number>
2166  inline
2167  const typename Accessor<number, false>::Reference &
2168  Accessor<number, false>::Reference::operator /= (const number n) const
2169  {
2170  AssertIndexRange(accessor->index_within_sparsity, accessor->matrix->n_nonzero_elements());
2171  accessor->matrix->val[accessor->index_within_sparsity] /= n;
2172  return *this;
2173  }
2174 
2175 
2176 
2177  template <typename number>
2178  inline
2179  Accessor<number,false>::
2180  Accessor (MatrixType *matrix,
2181  const size_type row,
2182  const size_type index)
2183  :
2184  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2185  row, index),
2186  matrix (matrix)
2187  {}
2188 
2189 
2190 
2191  template <typename number>
2192  inline
2193  Accessor<number,false>::
2194  Accessor (MatrixType *matrix,
2195  const std::size_t index)
2196  :
2197  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
2198  index),
2199  matrix (matrix)
2200  {}
2201 
2202 
2203 
2204  template <typename number>
2205  inline
2206  Accessor<number,false>::
2207  Accessor (MatrixType *matrix)
2208  :
2209  SparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
2210  matrix (matrix)
2211  {}
2212 
2213 
2214 
2215  template <typename number>
2216  inline
2217  typename Accessor<number, false>::Reference
2218  Accessor<number, false>::value() const
2219  {
2220  return Reference(this,true);
2221  }
2222 
2223 
2224 
2225 
2226  template <typename number>
2227  inline
2228  typename Accessor<number, false>::MatrixType &
2229  Accessor<number, false>::get_matrix () const
2230  {
2231  return *matrix;
2232  }
2233 
2234 
2235 
2236  template <typename number, bool Constness>
2237  inline
2238  Iterator<number, Constness>::
2239  Iterator (MatrixType *matrix,
2240  const size_type r,
2241  const size_type i)
2242  :
2243  accessor(matrix, r, i)
2244  {}
2245 
2246 
2247 
2248  template <typename number, bool Constness>
2249  inline
2250  Iterator<number, Constness>::
2251  Iterator (MatrixType *matrix,
2252  const std::size_t index)
2253  :
2254  accessor(matrix, index)
2255  {}
2256 
2257 
2258 
2259  template <typename number, bool Constness>
2260  inline
2261  Iterator<number, Constness>::
2262  Iterator (MatrixType *matrix)
2263  :
2264  accessor(matrix)
2265  {}
2266 
2267 
2268 
2269  template <typename number, bool Constness>
2270  inline
2271  Iterator<number, Constness>::
2273  :
2274  accessor(*i)
2275  {}
2276 
2277 
2278 
2279  template <typename number, bool Constness>
2280  inline
2281  Iterator<number, Constness> &
2282  Iterator<number,Constness>::operator++ ()
2283  {
2284  accessor.advance ();
2285  return *this;
2286  }
2287 
2288 
2289  template <typename number, bool Constness>
2290  inline
2291  Iterator<number,Constness>
2292  Iterator<number,Constness>::operator++ (int)
2293  {
2294  const Iterator iter = *this;
2295  accessor.advance ();
2296  return iter;
2297  }
2298 
2299 
2300  template <typename number, bool Constness>
2301  inline
2302  const Accessor<number,Constness> &
2303  Iterator<number,Constness>::operator* () const
2304  {
2305  return accessor;
2306  }
2307 
2308 
2309  template <typename number, bool Constness>
2310  inline
2311  const Accessor<number,Constness> *
2312  Iterator<number,Constness>::operator-> () const
2313  {
2314  return &accessor;
2315  }
2316 
2317 
2318  template <typename number, bool Constness>
2319  inline
2320  bool
2321  Iterator<number,Constness>::
2322  operator == (const Iterator &other) const
2323  {
2324  return (accessor == other.accessor);
2325  }
2326 
2327 
2328  template <typename number, bool Constness>
2329  inline
2330  bool
2331  Iterator<number,Constness>::
2332  operator != (const Iterator &other) const
2333  {
2334  return ! (*this == other);
2335  }
2336 
2337 
2338  template <typename number, bool Constness>
2339  inline
2340  bool
2341  Iterator<number,Constness>::
2342  operator < (const Iterator &other) const
2343  {
2344  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2345  ExcInternalError());
2346 
2347  return (accessor < other.accessor);
2348  }
2349 
2350 
2351  template <typename number, bool Constness>
2352  inline
2353  bool
2354  Iterator<number,Constness>::
2355  operator > (const Iterator &other) const
2356  {
2357  return (other < *this);
2358  }
2359 
2360 
2361  template <typename number, bool Constness>
2362  inline
2363  int
2364  Iterator<number,Constness>::
2365  operator - (const Iterator &other) const
2366  {
2367  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
2368  ExcInternalError());
2369 
2370  return (*this)->index_within_sparsity - other->index_within_sparsity;
2371  }
2372 
2373 
2374 
2375  template <typename number, bool Constness>
2376  inline
2377  Iterator<number,Constness>
2378  Iterator<number,Constness>::
2379  operator + (const size_type n) const
2380  {
2381  Iterator x = *this;
2382  for (size_type i=0; i<n; ++i)
2383  ++x;
2384 
2385  return x;
2386  }
2387 
2388 }
2389 
2390 
2391 
2392 template <typename number>
2393 inline
2396 {
2397  return const_iterator(this, 0);
2398 }
2399 
2400 
2401 template <typename number>
2402 inline
2405 {
2406  return const_iterator(this);
2407 }
2408 
2409 
2410 template <typename number>
2411 inline
2414 {
2415  return iterator (this, 0);
2416 }
2417 
2418 
2419 template <typename number>
2420 inline
2423 {
2424  return iterator(this, cols->rowstart[cols->rows]);
2425 }
2426 
2427 
2428 template <typename number>
2429 inline
2431 SparseMatrix<number>::begin (const size_type r) const
2432 {
2433  Assert (r<m(), ExcIndexRange(r,0,m()));
2434 
2435  return const_iterator(this, cols->rowstart[r]);
2436 }
2437 
2438 
2439 
2440 template <typename number>
2441 inline
2443 SparseMatrix<number>::end (const size_type r) const
2444 {
2445  Assert (r<m(), ExcIndexRange(r,0,m()));
2446 
2447  return const_iterator(this, cols->rowstart[r+1]);
2448 }
2449 
2450 
2451 
2452 template <typename number>
2453 inline
2455 SparseMatrix<number>::begin (const size_type r)
2456 {
2457  Assert (r<m(), ExcIndexRange(r,0,m()));
2458 
2459  return iterator(this, cols->rowstart[r]);
2460 }
2461 
2462 
2463 
2464 template <typename number>
2465 inline
2467 SparseMatrix<number>::end (const size_type r)
2468 {
2469  Assert (r<m(), ExcIndexRange(r,0,m()));
2470 
2471  return iterator(this, cols->rowstart[r+1]);
2472 }
2473 
2474 
2475 
2476 template <typename number>
2477 template <class STREAM>
2478 inline
2479 void SparseMatrix<number>::print (STREAM &out,
2480  const bool across,
2481  const bool diagonal_first) const
2482 {
2483  Assert (cols != 0, ExcNotInitialized());
2484  Assert (val != 0, ExcNotInitialized());
2485 
2486  bool hanging_diagonal = false;
2487  number diagonal = number();
2488 
2489  for (size_type i=0; i<cols->rows; ++i)
2490  {
2491  for (size_type j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
2492  {
2493  if (!diagonal_first && i == cols->colnums[j])
2494  {
2495  diagonal = val[j];
2496  hanging_diagonal = true;
2497  }
2498  else
2499  {
2500  if (hanging_diagonal && cols->colnums[j]>i)
2501  {
2502  if (across)
2503  out << ' ' << i << ',' << i << ':' << diagonal;
2504  else
2505  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2506  hanging_diagonal = false;
2507  }
2508  if (across)
2509  out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2510  else
2511  out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << std::endl;
2512  }
2513  }
2514  if (hanging_diagonal)
2515  {
2516  if (across)
2517  out << ' ' << i << ',' << i << ':' << diagonal;
2518  else
2519  out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2520  hanging_diagonal = false;
2521  }
2522  }
2523  if (across)
2524  out << std::endl;
2525 }
2526 
2527 
2528 template <typename number>
2529 inline
2530 void
2532 prepare_add()
2533 {
2534  //nothing to do here
2535 }
2536 
2537 
2538 
2539 template <typename number>
2540 inline
2541 void
2543 prepare_set()
2544 {
2545  //nothing to do here
2546 }
2547 
2548 #endif // DOXYGEN
2549 
2550 
2551 /*---------------------------- sparse_matrix.h ---------------------------*/
2552 
2553 DEAL_II_NAMESPACE_CLOSE
2554 
2555 #endif
2556 /*---------------------------- sparse_matrix.h ---------------------------*/
number raw_entry(const size_type row, const size_type index) const DEAL_II_DEPRECATED
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:525
numbers::NumberTraits< number >::real_type real_type
void prepare_add()
Accessor< number, Constness >::MatrixType MatrixType
size_type m() const
size_type n() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
SparseMatrixIterators::Iterator< number, false > iterator
bool is_finite(const double x)
const_iterator begin() const
void print(STREAM &out, const bool across=false, const bool diagonal_first=true) const
size_type n() const
TrilinosScalar operator()(const size_type i, const size_type j) const
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
unsigned int global_dof_index
Definition: types.h:100
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
Definition: exceptions.h:299
SparseMatrix & operator/=(const TrilinosScalar factor)
types::global_dof_index size_type
Definition: sparse_matrix.h:58
const SparsityPattern & get_sparsity_pattern() const
number value_type
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void prepare_set()
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
size_type m() const
const Accessor< number, Constness > & value_type
const_iterator end() const
Accessor< number, Constness > accessor
void copy_from(const SparseMatrix &source)
SparseMatrixIterators::Iterator< number, true > const_iterator
std::size_t * rowstart
SparseMatrix & operator*=(const TrilinosScalar factor)
TrilinosScalar el(const size_type i, const size_type j) const
types::global_dof_index size_type
::ExceptionBase & ExcNumberNotFinite()
number global_entry(const size_type i) const DEAL_II_DEPRECATED
TrilinosScalar diag_element(const size_type i) const
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
::ExceptionBase & ExcNotInitialized()
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
void set(const size_type i, const size_type j, const TrilinosScalar value)
std::size_t max_len
static const size_type invalid_entry