Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
chunk_sparse_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_chunk_sparse_matrix_h
16#define dealii_chunk_sparse_matrix_h
17
18
19#include <deal.II/base/config.h>
20
23
27
28#include <iterator>
29#include <memory>
30
31
33
34// Forward declarations
35#ifndef DOXYGEN
36template <typename number>
37class Vector;
38template <typename number>
39class FullMatrix;
40#endif
41
52{
53 // forward declaration
54 template <typename number, bool Constness>
55 class Iterator;
56
67 template <typename number, bool Constness>
69 {
70 public:
74 number
75 value() const;
76
80 number &
82
88 get_matrix() const;
89 };
90
91
92
99 template <typename number>
101 {
102 public:
108
112 Accessor(MatrixType *matrix, const unsigned int row);
113
118
123
127 number
128 value() const;
129
134 const MatrixType &
135 get_matrix() const;
136
137 private:
142
147
148 // Make iterator class a friend.
149 template <typename, bool>
150 friend class Iterator;
151 };
152
153
160 template <typename number>
162 {
163 private:
187 {
188 public:
193 Reference(const Accessor *accessor, const bool dummy);
194
198 operator number() const;
199
203 const Reference &
204 operator=(const number n) const;
205
209 const Reference &
210 operator+=(const number n) const;
211
215 const Reference &
216 operator-=(const number n) const;
217
221 const Reference &
222 operator*=(const number n) const;
223
227 const Reference &
228 operator/=(const number n) const;
229
230 private:
236 };
237
238 public:
244
248 Accessor(MatrixType *matrix, const unsigned int row);
249
254
259 value() const;
260
265 MatrixType &
266 get_matrix() const;
267
268 private:
273
278
279 // Make iterator class a friend.
280 template <typename, bool>
281 friend class Iterator;
282 };
283
284
285
296 template <typename number, bool Constness>
298 {
299 public:
304
310
316
321 Iterator(MatrixType *matrix, const unsigned int row);
322
327
333
337 Iterator &
339
345
350 operator*() const;
351
356 operator->() const;
357
361 bool
362 operator==(const Iterator &) const;
363
367 bool
368 operator!=(const Iterator &) const;
369
377 bool
378 operator<(const Iterator &) const;
379
384 bool
385 operator>(const Iterator &) const;
386
393 int
394 operator-(const Iterator &p) const;
395
400 operator+(const unsigned int n) const;
401
402 private:
407 };
408
409} // namespace ChunkSparseMatrixIterators
410
412
413namespace std
414{
415 template <typename number, bool Constness>
417 ::ChunkSparseMatrixIterators::Iterator<number, Constness>>
418 {
419 using iterator_category = forward_iterator_tag;
420 using value_type = typename ::ChunkSparseMatrixIterators::
421 Iterator<number, Constness>::value_type;
422 using difference_type = typename ::ChunkSparseMatrixIterators::
423 Iterator<number, Constness>::difference_type;
424 };
425} // namespace std
426
428
429
430
447template <typename number>
448class ChunkSparseMatrix : public virtual Subscriptor
449{
450public:
455
460 using value_type = number;
461
472
478
486
493 struct Traits
494 {
499 static const bool zero_addition_can_be_elided = true;
500 };
501
517
527
541 explicit ChunkSparseMatrix(const ChunkSparsityPattern &sparsity);
542
550 const IdentityMatrix &id);
551
556 virtual ~ChunkSparseMatrix() override;
557
569
578
589 operator=(const double d);
590
604 virtual void
605 reinit(const ChunkSparsityPattern &sparsity);
606
612 virtual void
623 bool
624 empty() const;
625
631 m() const;
632
638 n() const;
639
647
657
668
673 std::size_t
675
686 void
687 set(const size_type i, const size_type j, const number value);
688
694 void
695 add(const size_type i, const size_type j, const number value);
696
706 template <typename number2>
707 void
708 add(const size_type row,
709 const size_type n_cols,
710 const size_type *col_indices,
711 const number2 *values,
712 const bool elide_zero_values = true,
713 const bool col_indices_are_sorted = false);
714
719 operator*=(const number factor);
720
725 operator/=(const number factor);
726
739 void
741
758 template <typename somenumber>
761
779 template <typename ForwardIterator>
780 void
781 copy_from(const ForwardIterator begin, const ForwardIterator end);
782
788 template <typename somenumber>
789 void
791
803 template <typename somenumber>
804 void
805 add(const number factor, const ChunkSparseMatrix<somenumber> &matrix);
806
826 number
827 operator()(const size_type i, const size_type j) const;
828
841 number
842 el(const size_type i, const size_type j) const;
843
853 number
854 diag_element(const size_type i) const;
855
865 void
867 const size_type array_length,
868 size_type &row_length,
869 size_type *column_indices,
870 number *values) const;
871
891 template <class OutVector, class InVector>
892 void
893 vmult(OutVector &dst, const InVector &src) const;
894
910 template <class OutVector, class InVector>
911 void
912 Tvmult(OutVector &dst, const InVector &src) const;
913
928 template <class OutVector, class InVector>
929 void
930 vmult_add(OutVector &dst, const InVector &src) const;
931
947 template <class OutVector, class InVector>
948 void
949 Tvmult_add(OutVector &dst, const InVector &src) const;
950
966 template <typename somenumber>
967 somenumber
969
973 template <typename somenumber>
974 somenumber
976 const Vector<somenumber> &v) const;
984 template <typename somenumber>
985 somenumber
987 const Vector<somenumber> &x,
988 const Vector<somenumber> &b) const;
989
1003 real_type
1004 l1_norm() const;
1005
1013 real_type
1015
1020 real_type
1033 template <typename somenumber>
1034 void
1036 const Vector<somenumber> &src,
1037 const number omega = 1.) const;
1038
1042 template <typename somenumber>
1043 void
1045 const Vector<somenumber> &src,
1046 const number om = 1.) const;
1047
1051 template <typename somenumber>
1052 void
1054 const Vector<somenumber> &src,
1055 const number om = 1.) const;
1056
1060 template <typename somenumber>
1061 void
1063 const Vector<somenumber> &src,
1064 const number om = 1.) const;
1065
1071 template <typename somenumber>
1072 void
1073 SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1074
1079 template <typename somenumber>
1080 void
1081 SOR(Vector<somenumber> &v, const number om = 1.) const;
1082
1087 template <typename somenumber>
1088 void
1089 TSOR(Vector<somenumber> &v, const number om = 1.) const;
1090
1101 template <typename somenumber>
1102 void
1104 const std::vector<size_type> &permutation,
1105 const std::vector<size_type> &inverse_permutation,
1106 const number om = 1.) const;
1107
1118 template <typename somenumber>
1119 void
1121 const std::vector<size_type> &permutation,
1122 const std::vector<size_type> &inverse_permutation,
1123 const number om = 1.) const;
1124
1129 template <typename somenumber>
1130 void
1132 const Vector<somenumber> &b,
1133 const number om = 1.) const;
1134
1139 template <typename somenumber>
1140 void
1142 const Vector<somenumber> &b,
1143 const number om = 1.) const;
1144
1149 template <typename somenumber>
1150 void
1152 const Vector<somenumber> &b,
1153 const number om = 1.) const;
1170 begin() const;
1171
1181 end() const;
1182
1192 iterator
1194
1203 iterator
1205
1221 begin(const unsigned int r) const;
1222
1238 end(const unsigned int r) const;
1239
1254 iterator
1255 begin(const unsigned int r);
1256
1271 iterator
1272 end(const unsigned int r);
1283 void
1284 print(std::ostream &out) const;
1285
1308 void
1309 print_formatted(std::ostream &out,
1310 const unsigned int precision = 3,
1311 const bool scientific = true,
1312 const unsigned int width = 0,
1313 const char *zero_string = " ",
1314 const double denominator = 1.,
1315 const char *separator = " ") const;
1316
1322 void
1323 print_pattern(std::ostream &out, const double threshold = 0.) const;
1324
1335 void
1336 block_write(std::ostream &out) const;
1337
1354 void
1355 block_read(std::istream &in);
1366 int,
1367 int,
1368 << "You are trying to access the matrix entry with index <"
1369 << arg1 << ',' << arg2
1370 << ">, but this entry does not exist in the sparsity pattern "
1371 "of this matrix."
1372 "\n\n"
1373 "The most common cause for this problem is that you used "
1374 "a method to build the sparsity pattern that did not "
1375 "(completely) take into account all of the entries you "
1376 "will later try to write into. An example would be "
1377 "building a sparsity pattern that does not include "
1378 "the entries you will write into due to constraints "
1379 "on degrees of freedom such as hanging nodes or periodic "
1380 "boundary conditions. In such cases, building the "
1381 "sparsity pattern will succeed, but you will get errors "
1382 "such as the current one at one point or other when "
1383 "trying to write into the entries of the matrix.");
1392 int,
1393 int,
1394 << "The iterators denote a range of " << arg1
1395 << " elements, but the given number of rows was " << arg2);
1400 "You are attempting an operation on two vectors that "
1401 "are the same object, but the operation requires that the "
1402 "two objects are in fact different.");
1404private:
1411
1419 std::unique_ptr<number[]> val;
1420
1428
1432 size_type
1433 compute_location(const size_type i, const size_type j) const;
1434
1435 // make all other sparse matrices friends
1436 template <typename somenumber>
1437 friend class ChunkSparseMatrix;
1438
1439 // Also give access to internal details to the iterator/accessor classes.
1440 template <typename, bool>
1442 template <typename, bool>
1444};
1445
1448#ifndef DOXYGEN
1449/*---------------------- Inline functions -----------------------------------*/
1450
1451
1452
1453template <typename number>
1456{
1457 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1458 return cols->rows;
1459}
1460
1461
1462template <typename number>
1465{
1466 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1467 return cols->cols;
1468}
1469
1470
1471
1472template <typename number>
1473inline const ChunkSparsityPattern &
1475{
1476 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1477 return *cols;
1478}
1479
1480
1481
1482template <typename number>
1485 const size_type j) const
1486{
1487 const size_type chunk_size = cols->get_chunk_size();
1488 const size_type chunk_index =
1489 cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1490
1491 if (chunk_index == ChunkSparsityPattern::invalid_entry)
1493 else
1494 {
1495 return (chunk_index * chunk_size * chunk_size +
1496 (i % chunk_size) * chunk_size + (j % chunk_size));
1497 }
1498}
1499
1500
1501template <typename number>
1502inline void
1503ChunkSparseMatrix<number>::set(const size_type i,
1504 const size_type j,
1505 const number value)
1506{
1508
1509 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1510 // it is allowed to set elements of the matrix that are not part of the
1511 // sparsity pattern, if the value to which we set it is zero
1512 const size_type index = compute_location(i, j);
1513 Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1514 ExcInvalidIndex(i, j));
1515
1516 if (index != SparsityPattern::invalid_entry)
1517 val[index] = value;
1518}
1519
1520
1521
1522template <typename number>
1523inline void
1524ChunkSparseMatrix<number>::add(const size_type i,
1525 const size_type j,
1526 const number value)
1527{
1529
1530 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1531
1532 if (std::abs(value) != 0.)
1533 {
1534 const size_type index = compute_location(i, j);
1536 ExcInvalidIndex(i, j));
1537
1538 val[index] += value;
1539 }
1540}
1541
1542
1543
1544template <typename number>
1545template <typename number2>
1546inline void
1547ChunkSparseMatrix<number>::add(const size_type row,
1548 const size_type n_cols,
1549 const size_type *col_indices,
1550 const number2 *values,
1551 const bool /*elide_zero_values*/,
1552 const bool /*col_indices_are_sorted*/)
1553{
1554 // TODO: could be done more efficiently...
1555 for (size_type col = 0; col < n_cols; ++col)
1556 add(row, col_indices[col], static_cast<number>(values[col]));
1557}
1558
1559
1560
1561template <typename number>
1563ChunkSparseMatrix<number>::operator*=(const number factor)
1564{
1565 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1566 Assert(val != nullptr, ExcNotInitialized());
1567
1568 const size_type chunk_size = cols->get_chunk_size();
1569
1570 // multiply all elements of the matrix with the given factor. this includes
1571 // the padding elements in chunks that overlap the boundaries of the actual
1572 // matrix -- but since multiplication with a number does not violate the
1573 // invariant of keeping these elements at zero nothing can happen
1574 number *val_ptr = val.get();
1575 const number *const end_ptr =
1576 val.get() +
1577 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1578 while (val_ptr != end_ptr)
1579 *val_ptr++ *= factor;
1580
1581 return *this;
1582}
1583
1584
1585
1586template <typename number>
1588ChunkSparseMatrix<number>::operator/=(const number factor)
1589{
1590 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1591 Assert(val != nullptr, ExcNotInitialized());
1592 Assert(std::abs(factor) != 0, ExcDivideByZero());
1593
1594 const number factor_inv = 1. / factor;
1595
1596 const size_type chunk_size = cols->get_chunk_size();
1597
1598 // multiply all elements of the matrix with the given factor. this includes
1599 // the padding elements in chunks that overlap the boundaries of the actual
1600 // matrix -- but since multiplication with a number does not violate the
1601 // invariant of keeping these elements at zero nothing can happen
1602 number *val_ptr = val.get();
1603 const number *const end_ptr =
1604 val.get() +
1605 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1606
1607 while (val_ptr != end_ptr)
1608 *val_ptr++ *= factor_inv;
1609
1610 return *this;
1611}
1612
1613
1614
1615template <typename number>
1616inline number
1618 const size_type j) const
1619{
1620 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1621 AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1622 ExcInvalidIndex(i, j));
1623 return val[compute_location(i, j)];
1624}
1625
1626
1627
1628template <typename number>
1629inline number
1630ChunkSparseMatrix<number>::el(const size_type i, const size_type j) const
1631{
1632 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1633 const size_type index = compute_location(i, j);
1634
1636 return val[index];
1637 else
1638 return 0;
1639}
1640
1641
1642
1643template <typename number>
1644inline number
1645ChunkSparseMatrix<number>::diag_element(const size_type i) const
1646{
1647 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1648 Assert(m() == n(), ExcNotQuadratic());
1649 AssertIndexRange(i, m());
1650
1651 // Use that the first element in each row of a quadratic matrix is the main
1652 // diagonal of the chunk sparsity pattern
1653 const size_type chunk_size = cols->get_chunk_size();
1654 return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1655 chunk_size +
1656 (i % chunk_size) * chunk_size + (i % chunk_size)];
1657}
1658
1659
1660
1661template <typename number>
1662template <typename ForwardIterator>
1663inline void
1664ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1665 const ForwardIterator end)
1666{
1667 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1668 ExcIteratorRange(std::distance(begin, end), m()));
1669
1670 // for use in the inner loop, we define an alias to the type of the inner
1671 // iterators
1672 using inner_iterator =
1673 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1674 size_type row = 0;
1675 for (ForwardIterator i = begin; i != end; ++i, ++row)
1676 {
1677 const inner_iterator end_of_row = i->end();
1678 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1679 // write entries
1680 set(row, j->first, j->second);
1681 }
1682}
1683
1684
1685
1686//---------------------------------------------------------------------------
1687
1688
1690{
1691 template <typename number>
1692 inline Accessor<number, true>::Accessor(const MatrixType *matrix,
1693 const unsigned int row)
1694 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1695 row)
1696 , matrix(matrix)
1697 {}
1698
1699
1700
1701 template <typename number>
1702 inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1703 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1704 , matrix(matrix)
1705 {}
1706
1707
1708
1709 template <typename number>
1710 inline Accessor<number, true>::Accessor(
1712 : ChunkSparsityPatternIterators::Accessor(a)
1713 , matrix(&a.get_matrix())
1714 {}
1715
1716
1717
1718 template <typename number>
1719 inline number
1720 Accessor<number, true>::value() const
1721 {
1722 const unsigned int chunk_size =
1723 matrix->get_sparsity_pattern().get_chunk_size();
1724 return matrix->val[reduced_index() * chunk_size * chunk_size +
1725 chunk_row * chunk_size + chunk_col];
1726 }
1727
1728
1729
1730 template <typename number>
1731 inline const typename Accessor<number, true>::MatrixType &
1732 Accessor<number, true>::get_matrix() const
1733 {
1734 return *matrix;
1735 }
1736
1737
1738
1739 template <typename number>
1740 inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1741 const bool)
1742 : accessor(accessor)
1743 {}
1744
1745
1746 template <typename number>
1747 inline Accessor<number, false>::Reference::operator number() const
1748 {
1749 const unsigned int chunk_size =
1750 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1751 return accessor->matrix
1752 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1753 accessor->chunk_row * chunk_size + accessor->chunk_col];
1754 }
1755
1756
1757
1758 template <typename number>
1759 inline const typename Accessor<number, false>::Reference &
1760 Accessor<number, false>::Reference::operator=(const number n) const
1761 {
1762 const unsigned int chunk_size =
1763 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1764 accessor->matrix
1765 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1766 accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1767 return *this;
1768 }
1769
1770
1771
1772 template <typename number>
1773 inline const typename Accessor<number, false>::Reference &
1774 Accessor<number, false>::Reference::operator+=(const number n) const
1775 {
1776 const unsigned int chunk_size =
1777 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1778 accessor->matrix
1779 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1780 accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1781 return *this;
1782 }
1783
1784
1785
1786 template <typename number>
1787 inline const typename Accessor<number, false>::Reference &
1788 Accessor<number, false>::Reference::operator-=(const number n) const
1789 {
1790 const unsigned int chunk_size =
1791 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1792 accessor->matrix
1793 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1794 accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1795 return *this;
1796 }
1797
1798
1799
1800 template <typename number>
1801 inline const typename Accessor<number, false>::Reference &
1802 Accessor<number, false>::Reference::operator*=(const number n) const
1803 {
1804 const unsigned int chunk_size =
1805 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1806 accessor->matrix
1807 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1808 accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1809 return *this;
1810 }
1811
1812
1813
1814 template <typename number>
1815 inline const typename Accessor<number, false>::Reference &
1816 Accessor<number, false>::Reference::operator/=(const number n) const
1817 {
1818 const unsigned int chunk_size =
1819 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1820 accessor->matrix
1821 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1822 accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1823 return *this;
1824 }
1825
1826
1827
1828 template <typename number>
1829 inline Accessor<number, false>::Accessor(MatrixType *matrix,
1830 const unsigned int row)
1831 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1832 row)
1833 , matrix(matrix)
1834 {}
1835
1836
1837
1838 template <typename number>
1839 inline Accessor<number, false>::Accessor(MatrixType *matrix)
1840 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1841 , matrix(matrix)
1842 {}
1843
1844
1845
1846 template <typename number>
1847 inline typename Accessor<number, false>::Reference
1848 Accessor<number, false>::value() const
1849 {
1850 return Reference(this, true);
1851 }
1852
1853
1854
1855 template <typename number>
1856 inline typename Accessor<number, false>::MatrixType &
1857 Accessor<number, false>::get_matrix() const
1858 {
1859 return *matrix;
1860 }
1861
1862
1863
1864 template <typename number, bool Constness>
1865 inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
1866 const unsigned int row)
1867 : accessor(matrix, row)
1868 {}
1869
1870
1871
1872 template <typename number, bool Constness>
1873 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1874 : accessor(matrix)
1875 {}
1876
1877
1878
1879 template <typename number, bool Constness>
1880 inline Iterator<number, Constness>::Iterator(
1882 : accessor(*i)
1883 {}
1884
1885
1886
1887 template <typename number, bool Constness>
1888 inline Iterator<number, Constness> &
1889 Iterator<number, Constness>::operator++()
1890 {
1891 accessor.advance();
1892 return *this;
1893 }
1894
1895
1896 template <typename number, bool Constness>
1897 inline Iterator<number, Constness>
1898 Iterator<number, Constness>::operator++(int)
1899 {
1900 const Iterator iter = *this;
1901 accessor.advance();
1902 return iter;
1903 }
1904
1905
1906 template <typename number, bool Constness>
1907 inline const Accessor<number, Constness> &
1908 Iterator<number, Constness>::operator*() const
1909 {
1910 return accessor;
1911 }
1912
1913
1914 template <typename number, bool Constness>
1915 inline const Accessor<number, Constness> *
1916 Iterator<number, Constness>::operator->() const
1917 {
1918 return &accessor;
1919 }
1920
1921
1922 template <typename number, bool Constness>
1923 inline bool
1924 Iterator<number, Constness>::operator==(const Iterator &other) const
1925 {
1926 return (accessor == other.accessor);
1927 }
1928
1929
1930 template <typename number, bool Constness>
1931 inline bool
1932 Iterator<number, Constness>::operator!=(const Iterator &other) const
1933 {
1934 return !(*this == other);
1935 }
1936
1937
1938 template <typename number, bool Constness>
1939 inline bool
1940 Iterator<number, Constness>::operator<(const Iterator &other) const
1941 {
1942 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1944
1945 return (accessor < other.accessor);
1946 }
1947
1948
1949 template <typename number, bool Constness>
1950 inline bool
1951 Iterator<number, Constness>::operator>(const Iterator &other) const
1952 {
1953 return (other < *this);
1954 }
1955
1956
1957 template <typename number, bool Constness>
1958 inline int
1959 Iterator<number, Constness>::operator-(const Iterator &other) const
1960 {
1961 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1963
1964 // TODO: can be optimized
1965 int difference = 0;
1966 if (*this < other)
1967 {
1968 Iterator copy = *this;
1969 while (copy != other)
1970 {
1971 ++copy;
1972 --difference;
1973 }
1974 }
1975 else
1976 {
1977 Iterator copy = other;
1978 while (copy != *this)
1979 {
1980 ++copy;
1981 ++difference;
1982 }
1983 }
1984 return difference;
1985 }
1986
1987
1988
1989 template <typename number, bool Constness>
1990 inline Iterator<number, Constness>
1991 Iterator<number, Constness>::operator+(const unsigned int n) const
1992 {
1993 Iterator x = *this;
1994 for (unsigned int i = 0; i < n; ++i)
1995 ++x;
1996
1997 return x;
1998 }
1999
2000} // namespace ChunkSparseMatrixIterators
2001
2002
2003
2004template <typename number>
2007{
2008 return const_iterator(this, 0);
2009}
2010
2011
2012template <typename number>
2015{
2016 return const_iterator(this);
2017}
2018
2019
2020template <typename number>
2023{
2024 return iterator(this, 0);
2025}
2026
2027
2028template <typename number>
2031{
2032 return iterator(this);
2033}
2034
2035
2036template <typename number>
2038ChunkSparseMatrix<number>::begin(const unsigned int r) const
2039{
2040 AssertIndexRange(r, m());
2041 return const_iterator(this, r);
2042}
2043
2044
2045
2046template <typename number>
2048ChunkSparseMatrix<number>::end(const unsigned int r) const
2049{
2050 AssertIndexRange(r, m());
2051 return const_iterator(this, r + 1);
2052}
2053
2054
2055
2056template <typename number>
2058ChunkSparseMatrix<number>::begin(const unsigned int r)
2059{
2060 AssertIndexRange(r, m());
2061 return iterator(this, r);
2062}
2063
2064
2065
2066template <typename number>
2068ChunkSparseMatrix<number>::end(const unsigned int r)
2069{
2070 AssertIndexRange(r, m());
2071 return iterator(this, r + 1);
2072}
2073
2074
2075
2076#endif // DOXYGEN
2077
2079
2080#endif
Reference(const Accessor *accessor, const bool dummy)
Accessor(MatrixType *matrix, const unsigned int row)
Accessor(MatrixType *matrix, const unsigned int row)
Accessor(const ChunkSparseMatrixIterators::Accessor< number, false > &a)
const ChunkSparseMatrix< number > & get_matrix() const
bool operator>(const Iterator &) const
Iterator(const ChunkSparseMatrixIterators::Iterator< number, false > &i)
const Accessor< number, Constness > & value_type
bool operator<(const Iterator &) const
Iterator(MatrixType *matrix, const unsigned int row)
const Accessor< number, Constness > * operator->() const
bool operator==(const Iterator &) const
Accessor< number, Constness > accessor
const Accessor< number, Constness > & operator*() const
int operator-(const Iterator &p) const
typename Accessor< number, Constness >::MatrixType MatrixType
bool operator!=(const Iterator &) const
Iterator operator+(const unsigned int n) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) 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 char *separator=" ") const
real_type linfty_norm() const
size_type compute_location(const size_type i, const size_type j) const
ChunkSparseMatrix(const ChunkSparsityPattern &sparsity)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
ChunkSparseMatrix< number > & operator=(const IdentityMatrix &id)
void Tvmult(OutVector &dst, const InVector &src) const
const_iterator begin() const
void add(const size_type i, const size_type j, const number value)
void set(const size_type i, const size_type j, const number value)
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const_iterator begin(const unsigned int r) const
void print(std::ostream &out) const
bool empty() const
ChunkSparseMatrix & operator/=(const number factor)
iterator begin()
ChunkSparseMatrix & operator=(const double d)
virtual void reinit(const ChunkSparsityPattern &sparsity)
void extract_row_copy(const size_type row, const size_type array_length, size_type &row_length, size_type *column_indices, number *values) const
iterator end(const unsigned int r)
void print_pattern(std::ostream &out, const double threshold=0.) 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 copy_from(const ForwardIterator begin, const ForwardIterator end)
real_type l1_norm() const
const_iterator end() const
void vmult_add(OutVector &dst, const InVector &src) const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
virtual ~ChunkSparseMatrix() override
std::unique_ptr< number[]> val
void block_write(std::ostream &out) const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
size_type m() const
size_type n_actually_nonzero_elements() const
ChunkSparseMatrix< number > & operator=(const ChunkSparseMatrix< number > &)
ChunkSparseMatrix(const ChunkSparsityPattern &sparsity, const IdentityMatrix &id)
void vmult(OutVector &dst, const InVector &src) const
void TSOR(Vector< somenumber > &v, const number om=1.) const
void add(const number factor, const ChunkSparseMatrix< somenumber > &matrix)
size_type n_nonzero_elements() const
void copy_from(const FullMatrix< somenumber > &matrix)
virtual void clear()
ChunkSparseMatrix & operator*=(const number factor)
void block_read(std::istream &in)
number operator()(const size_type i, const size_type j) const
void Tvmult_add(OutVector &dst, const InVector &src) const
number el(const size_type i, const size_type j) const
SmartPointer< const ChunkSparsityPattern, ChunkSparseMatrix< number > > cols
ChunkSparseMatrix< number > & copy_from(const ChunkSparseMatrix< somenumber > &source)
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
typename numbers::NumberTraits< number >::real_type real_type
ChunkSparseMatrix(const ChunkSparseMatrix &)
number diag_element(const size_type i) const
void SOR(Vector< somenumber > &v, const number om=1.) const
iterator begin(const unsigned int r)
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
std::size_t memory_consumption() const
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
real_type frobenius_norm() const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
size_type n() const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
const_iterator end(const unsigned int r) const
const ChunkSparsityPattern & get_sparsity_pattern() const
Accessor(const ChunkSparsityPattern *matrix, const size_type row)
static const size_type invalid_entry
static constexpr size_type invalid_entry
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcDifferentChunkSparsityPatterns()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
#define AssertThrow(cond, exc)
constexpr int chunk_size
Definition cuda_size.h:34
@ matrix
Contents is actually a matrix.
VectorType::value_type * end(VectorType &V)
void copy(const T *begin, const T *end, U *dest)
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:81
static const bool zero_addition_can_be_elided
typename ::ChunkSparseMatrixIterators:: Iterator< number, Constness >::value_type value_type
typename ::ChunkSparseMatrixIterators:: Iterator< number, Constness >::difference_type difference_type