Reference documentation for deal.II version 8.1.0
chunk_sparse_matrix.h
1 // ---------------------------------------------------------------------
2 // @f$Id: chunk_sparse_matrix.h 31414 2013-10-24 22:15:50Z bangerth @f$
3 //
4 // Copyright (C) 2008 - 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__chunk_sparse_matrix_h
18 #define __deal2__chunk_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/chunk_sparsity_pattern.h>
25 #include <deal.II/lac/identity_matrix.h>
26 #include <deal.II/lac/exceptions.h>
27 
29 
30 template<typename number> class Vector;
31 template<typename number> class FullMatrix;
32 
42 {
43  // forward declaration
44  template <typename number, bool Constness>
45  class Iterator;
46 
57  template <typename number, bool Constness>
59  {
60  public:
64  number value() const;
65 
69  number &value();
70 
75  const ChunkSparseMatrix<number> &get_matrix () const;
76  };
77 
78 
79 
86  template <typename number>
88  {
89  public:
95 
99  Accessor (MatrixType *matrix,
100  const unsigned int row);
101 
105  Accessor (MatrixType *matrix);
106 
111 
115  number value() const;
116 
121  MatrixType &get_matrix () const;
122 
123  private:
127  MatrixType *matrix;
128 
133 
137  template <typename, bool>
138  friend class Iterator;
139  };
140 
141 
148  template <typename number>
150  {
151  private:
174  class Reference
175  {
176  public:
181  Reference (const Accessor *accessor,
182  const bool dummy);
183 
187  operator number () const;
188 
192  const Reference &operator = (const number n) const;
193 
197  const Reference &operator += (const number n) const;
198 
202  const Reference &operator -= (const number n) const;
203 
207  const Reference &operator *= (const number n) const;
208 
212  const Reference &operator /= (const number n) const;
213 
214  private:
220  };
221 
222  public:
228 
232  Accessor (MatrixType *matrix,
233  const unsigned int row);
234 
238  Accessor (MatrixType *matrix);
239 
243  Reference value() const;
244 
249  MatrixType &get_matrix () const;
250 
251  private:
255  MatrixType *matrix;
256 
261 
265  template <typename, bool>
266  friend class Iterator;
267 
272  };
273 
274 
275 
286  template <typename number, bool Constness>
287  class Iterator
288  {
289  public:
293  typedef
296 
301  typedef
303 
308  Iterator (MatrixType *matrix,
309  const unsigned int row);
310 
314  Iterator (MatrixType *matrix);
315 
321 
326 
330  Iterator operator++ (int);
331 
335  const Accessor<number,Constness> &operator* () const;
336 
341 
345  bool operator == (const Iterator &) const;
346 
350  bool operator != (const Iterator &) const;
351 
359  bool operator < (const Iterator &) const;
360 
365  bool operator > (const Iterator &) const;
366 
373  int operator - (const Iterator &p) const;
374 
378  Iterator operator + (const unsigned int n) const;
379 
380  private:
385  };
386 
387 }
388 
389 
390 
408 template <typename number>
409 class ChunkSparseMatrix : public virtual Subscriptor
410 {
411 public:
416 
421  typedef number value_type;
422 
433 
439  typedef
442 
449  typedef
452 
459  struct Traits
460  {
465  static const bool zero_addition_can_be_elided = true;
466  };
467 
483 
493 
507  explicit ChunkSparseMatrix (const ChunkSparsityPattern &sparsity);
508 
515  ChunkSparseMatrix (const ChunkSparsityPattern &sparsity,
516  const IdentityMatrix &id);
517 
522  virtual ~ChunkSparseMatrix ();
523 
534 
542  operator= (const IdentityMatrix &id);
543 
553  ChunkSparseMatrix &operator = (const double d);
554 
568  virtual void reinit (const ChunkSparsityPattern &sparsity);
569 
575  virtual void clear ();
577 
585  bool empty () const;
586 
591  size_type m () const;
592 
597  size_type n () const;
598 
604  size_type n_nonzero_elements () const;
605 
613  size_type n_actually_nonzero_elements () const;
614 
624 
629  std::size_t memory_consumption () const;
630 
632 
641  void set (const size_type i,
642  const size_type j,
643  const number value);
644 
650  void add (const size_type i,
651  const size_type j,
652  const number value);
653 
663  template <typename number2>
664  void add (const size_type row,
665  const size_type n_cols,
666  const size_type *col_indices,
667  const number2 *values,
668  const bool elide_zero_values = true,
669  const bool col_indices_are_sorted = false);
670 
674  ChunkSparseMatrix &operator *= (const number factor);
675 
679  ChunkSparseMatrix &operator /= (const number factor);
680 
693  void symmetrize ();
694 
708  template <typename somenumber>
711 
729  template <typename ForwardIterator>
730  void copy_from (const ForwardIterator begin,
731  const ForwardIterator end);
732 
738  template <typename somenumber>
739  void copy_from (const FullMatrix<somenumber> &matrix);
740 
752  template <typename somenumber>
753  void add (const number factor,
754  const ChunkSparseMatrix<somenumber> &matrix);
755 
757 
761 
775  number operator () (const size_type i,
776  const size_type j) const;
777 
790  number el (const size_type i,
791  const size_type j) const;
792 
805  number diag_element (const size_type i) const;
806 
811  number &diag_element (const size_type i);
812 
814 
831  template <class OutVector, class InVector>
832  void vmult (OutVector &dst,
833  const InVector &src) const;
834 
849  template <class OutVector, class InVector>
850  void Tvmult (OutVector &dst,
851  const InVector &src) const;
852 
866  template <class OutVector, class InVector>
867  void vmult_add (OutVector &dst,
868  const InVector &src) const;
869 
884  template <class OutVector, class InVector>
885  void Tvmult_add (OutVector &dst,
886  const InVector &src) const;
887 
903  template <typename somenumber>
904  somenumber matrix_norm_square (const Vector<somenumber> &v) const;
905 
909  template <typename somenumber>
910  somenumber matrix_scalar_product (const Vector<somenumber> &u,
911  const Vector<somenumber> &v) const;
919  template <typename somenumber>
920  somenumber residual (Vector<somenumber> &dst,
921  const Vector<somenumber> &x,
922  const Vector<somenumber> &b) const;
923 
925 
929 
937  real_type l1_norm () const;
938 
946  real_type linfty_norm () const;
947 
952  real_type frobenius_norm () const;
954 
958 
964  template <typename somenumber>
966  const Vector<somenumber> &src,
967  const number omega = 1.) const;
968 
972  template <typename somenumber>
974  const Vector<somenumber> &src,
975  const number om = 1.) const;
976 
980  template <typename somenumber>
982  const Vector<somenumber> &src,
983  const number om = 1.) const;
984 
988  template <typename somenumber>
990  const Vector<somenumber> &src,
991  const number om = 1.) const;
992 
998  template <typename somenumber>
999  void SSOR (Vector<somenumber> &v,
1000  const number omega = 1.) const;
1001 
1006  template <typename somenumber>
1007  void SOR (Vector<somenumber> &v,
1008  const number om = 1.) const;
1009 
1014  template <typename somenumber>
1015  void TSOR (Vector<somenumber> &v,
1016  const number om = 1.) const;
1017 
1028  template <typename somenumber>
1029  void PSOR (Vector<somenumber> &v,
1030  const std::vector<size_type> &permutation,
1031  const std::vector<size_type> &inverse_permutation,
1032  const number om = 1.) const;
1033 
1044  template <typename somenumber>
1045  void TPSOR (Vector<somenumber> &v,
1046  const std::vector<size_type> &permutation,
1047  const std::vector<size_type> &inverse_permutation,
1048  const number om = 1.) const;
1049 
1054  template <typename somenumber>
1055  void SOR_step (Vector<somenumber> &v,
1056  const Vector<somenumber> &b,
1057  const number om = 1.) const;
1058 
1063  template <typename somenumber>
1064  void TSOR_step (Vector<somenumber> &v,
1065  const Vector<somenumber> &b,
1066  const number om = 1.) const;
1067 
1072  template <typename somenumber>
1073  void SSOR_step (Vector<somenumber> &v,
1074  const Vector<somenumber> &b,
1075  const number om = 1.) const;
1077 
1081 
1091  const_iterator begin () const;
1092 
1101  const_iterator end () const;
1102 
1112  iterator begin ();
1113 
1122  iterator end ();
1123 
1138  const_iterator begin (const unsigned int r) const;
1139 
1154  const_iterator end (const unsigned int r) const;
1155 
1170  iterator begin (const unsigned int r);
1171 
1186  iterator end (const unsigned int r);
1188 
1192 
1197  void print (std::ostream &out) const;
1198 
1219  void print_formatted (std::ostream &out,
1220  const unsigned int precision = 3,
1221  const bool scientific = true,
1222  const unsigned int width = 0,
1223  const char *zero_string = " ",
1224  const double denominator = 1.) const;
1225 
1231  void print_pattern(std::ostream &out,
1232  const double threshold = 0.) const;
1233 
1244  void block_write (std::ostream &out) const;
1245 
1262  void block_read (std::istream &in);
1264 
1270  DeclException2 (ExcInvalidIndex,
1271  int, int,
1272  << "The entry with index <" << arg1 << ',' << arg2
1273  << "> does not exist.");
1277  DeclException1 (ExcInvalidIndex1,
1278  int,
1279  << "The index " << arg1 << " is not in the allowed range.");
1283  DeclException0 (ExcDifferentChunkSparsityPatterns);
1287  DeclException2 (ExcIteratorRange,
1288  int, int,
1289  << "The iterators denote a range of " << arg1
1290  << " elements, but the given number of rows was " << arg2);
1294  DeclException0 (ExcSourceEqualsDestination);
1296 private:
1303 
1310  number *val;
1311 
1318  size_type max_len;
1319 
1323  size_type compute_location (const size_type i,
1324  const size_type j) const;
1325 
1326  // make all other sparse matrices friends
1327  template <typename somenumber> friend class ChunkSparseMatrix;
1328 
1333  template <typename,bool> friend class ChunkSparseMatrixIterators::Iterator;
1334  template <typename,bool> friend class ChunkSparseMatrixIterators::Accessor;
1335 };
1336 
1339 #ifndef DOXYGEN
1340 /*---------------------- Inline functions -----------------------------------*/
1341 
1342 
1343 
1344 template <typename number>
1345 inline
1348 {
1349  Assert (cols != 0, ExcNotInitialized());
1350  return cols->rows;
1351 }
1352 
1353 
1354 template <typename number>
1355 inline
1358 {
1359  Assert (cols != 0, ExcNotInitialized());
1360  return cols->cols;
1361 }
1362 
1363 
1364 
1365 template <typename number>
1366 inline
1367 const ChunkSparsityPattern &
1369 {
1370  Assert (cols != 0, ExcNotInitialized());
1371  return *cols;
1372 }
1373 
1374 
1375 
1376 template <typename number>
1377 inline
1380  const size_type j) const
1381 {
1382  const size_type chunk_size = cols->get_chunk_size();
1383  const size_type chunk_index
1384  = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
1385 
1386  if (chunk_index == ChunkSparsityPattern::invalid_entry)
1388  else
1389  {
1390  return (chunk_index * chunk_size * chunk_size
1391  +
1392  (i % chunk_size) * chunk_size
1393  +
1394  (j % chunk_size));
1395  }
1396 }
1397 
1398 
1399 template <typename number>
1400 inline
1402  const size_type j,
1403  const number value)
1404 {
1405 
1407 
1408  Assert (cols != 0, ExcNotInitialized());
1409  // it is allowed to set elements of the matrix that are not part of the
1410  // sparsity pattern, if the value to which we set it is zero
1411  const size_type index = compute_location(i,j);
1412  Assert ((index != SparsityPattern::invalid_entry) ||
1413  (value == 0.),
1414  ExcInvalidIndex(i,j));
1415 
1416  if (index != SparsityPattern::invalid_entry)
1417  val[index] = value;
1418 }
1419 
1420 
1421 
1422 template <typename number>
1423 inline
1425  const size_type j,
1426  const number value)
1427 {
1428 
1430 
1431  Assert (cols != 0, ExcNotInitialized());
1432 
1433  if (value != 0.)
1434  {
1435  const size_type index = compute_location(i,j);
1437  ExcInvalidIndex(i,j));
1438 
1439  val[index] += value;
1440  }
1441 }
1442 
1443 
1444 
1445 template <typename number>
1446 template <typename number2>
1447 inline
1449  const size_type n_cols,
1450  const size_type *col_indices,
1451  const number2 *values,
1452  const bool /*elide_zero_values*/,
1453  const bool /*col_indices_are_sorted*/)
1454 {
1455  // TODO: could be done more efficiently...
1456  for (size_type col=0; col<n_cols; ++col)
1457  add(row, col_indices[col], static_cast<number>(values[col]));
1458 }
1459 
1460 
1461 
1462 template <typename number>
1463 inline
1465 ChunkSparseMatrix<number>::operator *= (const number factor)
1466 {
1467  Assert (cols != 0, ExcNotInitialized());
1468  Assert (val != 0, ExcNotInitialized());
1469 
1470  const size_type chunk_size = cols->get_chunk_size();
1471 
1472  // multiply all elements of the matrix with the given factor. this includes
1473  // the padding elements in chunks that overlap the boundaries of the actual
1474  // matrix -- but since multiplication with a number does not violate the
1475  // invariant of keeping these elements at zero nothing can happen
1476  number *val_ptr = val;
1477  const number *const end_ptr = val +
1478  cols->sparsity_pattern.n_nonzero_elements()
1479  *
1480  chunk_size * chunk_size;
1481  while (val_ptr != end_ptr)
1482  *val_ptr++ *= factor;
1483 
1484  return *this;
1485 }
1486 
1487 
1488 
1489 template <typename number>
1490 inline
1492 ChunkSparseMatrix<number>::operator /= (const number factor)
1493 {
1494  Assert (cols != 0, ExcNotInitialized());
1495  Assert (val != 0, ExcNotInitialized());
1496  Assert (factor !=0, ExcDivideByZero());
1497 
1498  const number factor_inv = 1. / factor;
1499 
1500  const size_type chunk_size = cols->get_chunk_size();
1501 
1502  // multiply all elements of the matrix with the given factor. this includes
1503  // the padding elements in chunks that overlap the boundaries of the actual
1504  // matrix -- but since multiplication with a number does not violate the
1505  // invariant of keeping these elements at zero nothing can happen
1506  number *val_ptr = val;
1507  const number *const end_ptr = val +
1508  cols->sparsity_pattern.n_nonzero_elements()
1509  *
1510  chunk_size * chunk_size;
1511 
1512  while (val_ptr != end_ptr)
1513  *val_ptr++ *= factor_inv;
1514 
1515  return *this;
1516 }
1517 
1518 
1519 
1520 template <typename number>
1521 inline
1523  const size_type j) const
1524 {
1525  Assert (cols != 0, ExcNotInitialized());
1526  AssertThrow (compute_location(i,j) != SparsityPattern::invalid_entry,
1527  ExcInvalidIndex(i,j));
1528  return val[compute_location(i,j)];
1529 }
1530 
1531 
1532 
1533 template <typename number>
1534 inline
1536  const size_type j) const
1537 {
1538  Assert (cols != 0, ExcNotInitialized());
1539  const size_type index = compute_location(i,j);
1540 
1542  return val[index];
1543  else
1544  return 0;
1545 }
1546 
1547 
1548 
1549 template <typename number>
1550 inline
1552 {
1553  Assert (cols != 0, ExcNotInitialized());
1554  Assert (m() == n(), ExcNotQuadratic());
1555  Assert (i<m(), ExcInvalidIndex1(i));
1556 
1557  // Use that the first element in each row of a quadratic matrix is the main
1558  // diagonal of the chunk sparsity pattern
1559  const size_type chunk_size = cols->get_chunk_size();
1560  return val[cols->sparsity_pattern.rowstart[i/chunk_size]
1561  *
1562  chunk_size * chunk_size
1563  +
1564  (i % chunk_size) * chunk_size
1565  +
1566  (i % chunk_size)];
1567 }
1568 
1569 
1570 
1571 template <typename number>
1572 template <typename ForwardIterator>
1573 inline
1574 void
1575 ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
1576  const ForwardIterator end)
1577 {
1578  Assert (static_cast<size_type >(std::distance (begin, end)) == m(),
1579  ExcIteratorRange (std::distance (begin, end), m()));
1580 
1581  // for use in the inner loop, we define a typedef to the type of the inner
1582  // iterators
1583  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
1584  size_type row=0;
1585  for (ForwardIterator i=begin; i!=end; ++i, ++row)
1586  {
1587  const inner_iterator end_of_row = i->end();
1588  for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
1589  // write entries
1590  set (row, j->first, j->second);
1591  }
1592 }
1593 
1594 
1595 
1596 //---------------------------------------------------------------------------
1597 
1598 
1600 {
1601  template <typename number>
1602  inline
1604  Accessor (const MatrixType *matrix,
1605  const unsigned int row)
1606  :
1607  ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1608  row),
1609  matrix (matrix)
1610  {}
1611 
1612 
1613 
1614  template <typename number>
1615  inline
1616  Accessor<number,true>::
1617  Accessor (const MatrixType *matrix)
1618  :
1619  ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1620  matrix (matrix)
1621  {}
1622 
1623 
1624 
1625  template <typename number>
1626  inline
1627  Accessor<number,true>::
1629  :
1630  ChunkSparsityPatternIterators::Accessor (a),
1631  matrix (&a.get_matrix())
1632  {}
1633 
1634 
1635 
1636  template <typename number>
1637  inline
1638  number
1639  Accessor<number, true>::value () const
1640  {
1641  const unsigned int chunk_size = matrix->get_sparsity_pattern().get_chunk_size();
1642  return matrix->val[reduced_index() * chunk_size * chunk_size
1643  +
1644  chunk_row * chunk_size
1645  +
1646  chunk_col];
1647  }
1648 
1649 
1650 
1651  template <typename number>
1652  inline
1653  typename Accessor<number, true>::MatrixType &
1654  Accessor<number, true>::get_matrix () const
1655  {
1656  return *matrix;
1657  }
1658 
1659 
1660 
1661  template <typename number>
1662  inline
1663  Accessor<number, false>::Reference::Reference (
1664  const Accessor *accessor,
1665  const bool)
1666  :
1667  accessor (accessor)
1668  {}
1669 
1670 
1671  template <typename number>
1672  inline
1673  Accessor<number, false>::Reference::operator number() const
1674  {
1675  const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1676  return accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1677  +
1678  accessor->chunk_row * chunk_size
1679  +
1680  accessor->chunk_col];
1681  }
1682 
1683 
1684 
1685  template <typename number>
1686  inline
1687  const typename Accessor<number, false>::Reference &
1688  Accessor<number, false>::Reference::operator = (const number n) const
1689  {
1690  const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1691  accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1692  +
1693  accessor->chunk_row * chunk_size
1694  +
1695  accessor->chunk_col] = n;
1696  return *this;
1697  }
1698 
1699 
1700 
1701  template <typename number>
1702  inline
1703  const typename Accessor<number, false>::Reference &
1704  Accessor<number, false>::Reference::operator += (const number n) const
1705  {
1706  const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1707  accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1708  +
1709  accessor->chunk_row * chunk_size
1710  +
1711  accessor->chunk_col] += n;
1712  return *this;
1713  }
1714 
1715 
1716 
1717  template <typename number>
1718  inline
1719  const typename Accessor<number, false>::Reference &
1720  Accessor<number, false>::Reference::operator -= (const number n) const
1721  {
1722  const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1723  accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1724  +
1725  accessor->chunk_row * chunk_size
1726  +
1727  accessor->chunk_col] -= n;
1728  return *this;
1729  }
1730 
1731 
1732 
1733  template <typename number>
1734  inline
1735  const typename Accessor<number, false>::Reference &
1736  Accessor<number, false>::Reference::operator *= (const number n) const
1737  {
1738  const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1739  accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1740  +
1741  accessor->chunk_row * chunk_size
1742  +
1743  accessor->chunk_col] *= n;
1744  return *this;
1745  }
1746 
1747 
1748 
1749  template <typename number>
1750  inline
1751  const typename Accessor<number, false>::Reference &
1752  Accessor<number, false>::Reference::operator /= (const number n) const
1753  {
1754  const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
1755  accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
1756  +
1757  accessor->chunk_row * chunk_size
1758  +
1759  accessor->chunk_col] /= n;
1760  return *this;
1761  }
1762 
1763 
1764 
1765  template <typename number>
1766  inline
1767  Accessor<number,false>::
1768  Accessor (MatrixType *matrix,
1769  const unsigned int row)
1770  :
1771  ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
1772  row),
1773  matrix (matrix)
1774  {}
1775 
1776 
1777 
1778  template <typename number>
1779  inline
1780  Accessor<number,false>::
1781  Accessor (MatrixType *matrix)
1782  :
1783  ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
1784  matrix (matrix)
1785  {}
1786 
1787 
1788 
1789  template <typename number>
1790  inline
1791  typename Accessor<number, false>::Reference
1792  Accessor<number, false>::value() const
1793  {
1794  return Reference(this,true);
1795  }
1796 
1797 
1798 
1799 
1800  template <typename number>
1801  inline
1802  typename Accessor<number, false>::MatrixType &
1803  Accessor<number, false>::get_matrix () const
1804  {
1805  return *matrix;
1806  }
1807 
1808 
1809 
1810  template <typename number, bool Constness>
1811  inline
1812  Iterator<number, Constness>::
1813  Iterator (MatrixType *matrix,
1814  const unsigned int row)
1815  :
1816  accessor(matrix, row)
1817  {}
1818 
1819 
1820 
1821  template <typename number, bool Constness>
1822  inline
1823  Iterator<number, Constness>::
1824  Iterator (MatrixType *matrix)
1825  :
1826  accessor(matrix)
1827  {}
1828 
1829 
1830 
1831  template <typename number, bool Constness>
1832  inline
1833  Iterator<number, Constness>::
1835  :
1836  accessor(*i)
1837  {}
1838 
1839 
1840 
1841  template <typename number, bool Constness>
1842  inline
1843  Iterator<number, Constness> &
1844  Iterator<number,Constness>::operator++ ()
1845  {
1846  accessor.advance ();
1847  return *this;
1848  }
1849 
1850 
1851  template <typename number, bool Constness>
1852  inline
1853  Iterator<number,Constness>
1854  Iterator<number,Constness>::operator++ (int)
1855  {
1856  const Iterator iter = *this;
1857  accessor.advance ();
1858  return iter;
1859  }
1860 
1861 
1862  template <typename number, bool Constness>
1863  inline
1864  const Accessor<number,Constness> &
1865  Iterator<number,Constness>::operator* () const
1866  {
1867  return accessor;
1868  }
1869 
1870 
1871  template <typename number, bool Constness>
1872  inline
1873  const Accessor<number,Constness> *
1874  Iterator<number,Constness>::operator-> () const
1875  {
1876  return &accessor;
1877  }
1878 
1879 
1880  template <typename number, bool Constness>
1881  inline
1882  bool
1883  Iterator<number,Constness>::
1884  operator == (const Iterator &other) const
1885  {
1886  return (accessor == other.accessor);
1887  }
1888 
1889 
1890  template <typename number, bool Constness>
1891  inline
1892  bool
1893  Iterator<number,Constness>::
1894  operator != (const Iterator &other) const
1895  {
1896  return ! (*this == other);
1897  }
1898 
1899 
1900  template <typename number, bool Constness>
1901  inline
1902  bool
1903  Iterator<number,Constness>::
1904  operator < (const Iterator &other) const
1905  {
1906  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
1907  ExcInternalError());
1908 
1909  return (accessor < other.accessor);
1910  }
1911 
1912 
1913  template <typename number, bool Constness>
1914  inline
1915  bool
1916  Iterator<number,Constness>::
1917  operator > (const Iterator &other) const
1918  {
1919  return (other < *this);
1920  }
1921 
1922 
1923  template <typename number, bool Constness>
1924  inline
1925  int
1926  Iterator<number,Constness>::
1927  operator - (const Iterator &other) const
1928  {
1929  Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
1930  ExcInternalError());
1931 
1932  // TODO: can be optimized
1933  int difference = 0;
1934  if (*this < other)
1935  {
1936  Iterator copy = *this;
1937  while (copy != other)
1938  {
1939  ++copy;
1940  --difference;
1941  }
1942  }
1943  else
1944  {
1945  Iterator copy = other;
1946  while (copy != *this)
1947  {
1948  ++copy;
1949  ++difference;
1950  }
1951  }
1952  return difference;
1953  }
1954 
1955 
1956 
1957  template <typename number, bool Constness>
1958  inline
1959  Iterator<number,Constness>
1960  Iterator<number,Constness>::
1961  operator + (const unsigned int n) const
1962  {
1963  Iterator x = *this;
1964  for (unsigned int i=0; i<n; ++i)
1965  ++x;
1966 
1967  return x;
1968  }
1969 
1970 }
1971 
1972 
1973 
1974 template <typename number>
1975 inline
1978 {
1979  return const_iterator(this, 0);
1980 }
1981 
1982 
1983 template <typename number>
1984 inline
1987 {
1988  return const_iterator(this);
1989 }
1990 
1991 
1992 template <typename number>
1993 inline
1996 {
1997  return iterator(this, 0);
1998 }
1999 
2000 
2001 template <typename number>
2002 inline
2005 {
2006  return iterator(this);
2007 }
2008 
2009 
2010 template <typename number>
2011 inline
2013 ChunkSparseMatrix<number>::begin (const unsigned int r) const
2014 {
2015  Assert (r<m(), ExcIndexRange(r,0,m()));
2016  return const_iterator(this, r);
2017 }
2018 
2019 
2020 
2021 template <typename number>
2022 inline
2024 ChunkSparseMatrix<number>::end (const unsigned int r) const
2025 {
2026  Assert (r<m(), ExcIndexRange(r,0,m()));
2027  return const_iterator(this, r+1);
2028 }
2029 
2030 
2031 
2032 template <typename number>
2033 inline
2035 ChunkSparseMatrix<number>::begin (const unsigned int r)
2036 {
2037  Assert (r<m(), ExcIndexRange(r,0,m()));
2038  return iterator(this, r);
2039 }
2040 
2041 
2042 
2043 template <typename number>
2044 inline
2046 ChunkSparseMatrix<number>::end (const unsigned int r)
2047 {
2048  Assert (r<m(), ExcIndexRange(r,0,m()));
2049  return iterator(this, r+1);
2050 }
2051 
2052 
2053 
2054 
2055 #endif // DOXYGEN
2056 
2057 
2058 /*---------------------------- chunk_sparse_matrix.h ---------------------------*/
2059 
2060 DEAL_II_NAMESPACE_CLOSE
2061 
2062 #endif
2063 /*---------------------------- chunk_sparse_matrix.h ---------------------------*/
number el(const size_type i, const size_type j) const
Iterator operator+(const unsigned int n) const
bool operator>(const Iterator &) const
ChunkSparseMatrixIterators::Iterator< number, false > iterator
DeclException0(ExcDifferentChunkSparsityPatterns)
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void TSOR(Vector< somenumber > &v, const number om=1.) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
void print(std::ostream &out) const
ChunkSparseMatrixIterators::Iterator< number, true > const_iterator
void print_pattern(std::ostream &out, const double threshold=0.) const
const_iterator end() const
types::global_dof_index size_type
Iterator(MatrixType *matrix, const unsigned int row)
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
void SOR(Vector< somenumber > &v, const number om=1.) const
size_type n_actually_nonzero_elements() const
DeclException1(ExcInvalidIndex1, int,<< "The index "<< arg1<< " is not in the allowed range.")
bool is_finite(const double x)
void vmult_add(OutVector &dst, const InVector &src) const
size_type n_nonzero_elements() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void vmult(OutVector &dst, const InVector &src) const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
numbers::NumberTraits< number >::real_type real_type
bool operator!=(const Iterator &) const
const ChunkSparsityPattern & get_sparsity_pattern() const
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
bool operator==(const Iterator &) const
number operator()(const size_type i, const size_type j) const
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
void block_write(std::ostream &out) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
size_type n() const
ChunkSparseMatrix & operator*=(const number factor)
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
ChunkSparseMatrix< number > & operator=(const ChunkSparseMatrix< number > &)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type m() const
size_type compute_location(const size_type i, const size_type j) const
const Accessor< number, Constness > & value_type
Accessor(const ChunkSparsityPattern *matrix, const unsigned int row)
DeclException2(ExcInvalidIndex, int, int,<< "The entry with index <"<< arg1<< ','<< arg2<< "> does not exist.")
::ExceptionBase & ExcNumberNotFinite()
Accessor< number, Constness > accessor
static const bool zero_addition_can_be_elided
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const ChunkSparseMatrix< number > & get_matrix() const
Accessor< number, Constness >::MatrixType MatrixType
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
void block_read(std::istream &in)
ChunkSparseMatrix & operator/=(const number factor)
::ExceptionBase & ExcNotInitialized()
bool operator<(const Iterator &) const
number diag_element(const size_type i) const
void set(const size_type i, const size_type j, const number value)
static const size_type invalid_entry
std::size_t memory_consumption() const
virtual void reinit(const ChunkSparsityPattern &sparsity)
void add(const size_type i, const size_type j, const number value)
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
somenumber matrix_norm_square(const Vector< somenumber > &v) const
const_iterator begin() const
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
int operator-(const Iterator &p) const
const Accessor< number, Constness > & operator*() const
void Tvmult(OutVector &dst, const InVector &src) const
const Accessor< number, Constness > * operator->() const
void TPSOR(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
static const size_type invalid_entry