Reference documentation for deal.II version GIT acb749f4f5 2024-02-21 16:10:02+00:00
\(\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// Copyright (C) 2008 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_chunk_sparse_matrix_h
17#define dealii_chunk_sparse_matrix_h
18
19
20#include <deal.II/base/config.h>
21
24
28
29#include <iterator>
30#include <memory>
31
32
34
35// Forward declarations
36#ifndef DOXYGEN
37template <typename number>
38class Vector;
39template <typename number>
40class FullMatrix;
41#endif
42
53{
54 // forward declaration
55 template <typename number, bool Constness>
56 class Iterator;
57
68 template <typename number, bool Constness>
70 {
71 public:
75 number
76 value() const;
77
81 number &
83
89 get_matrix() const;
90 };
91
92
93
100 template <typename number>
102 {
103 public:
109
113 Accessor(MatrixType *matrix, const unsigned int row);
114
119
124
128 number
129 value() const;
130
135 const MatrixType &
136 get_matrix() const;
137
138 private:
143
148
149 // Make iterator class a friend.
150 template <typename, bool>
151 friend class Iterator;
152 };
153
154
161 template <typename number>
163 {
164 private:
187 class Reference
188 {
189 public:
194 Reference(const Accessor *accessor, const bool dummy);
195
199 operator number() const;
200
204 const Reference &
205 operator=(const number n) const;
206
210 const Reference &
211 operator+=(const number n) const;
212
216 const Reference &
217 operator-=(const number n) const;
218
222 const Reference &
223 operator*=(const number n) const;
224
228 const Reference &
229 operator/=(const number n) const;
230
231 private:
237 };
238
239 public:
245
249 Accessor(MatrixType *matrix, const unsigned int row);
250
255
259 Reference
260 value() const;
261
266 MatrixType &
267 get_matrix() const;
268
269 private:
274
279
280 // Make iterator class a friend.
281 template <typename, bool>
282 friend class Iterator;
283 };
284
285
286
297 template <typename number, bool Constness>
299 {
300 public:
305
311
317
322 Iterator(MatrixType *matrix, const unsigned int row);
323
328
334
338 Iterator &
340
346
351 operator*() const;
352
357 operator->() const;
358
362 bool
363 operator==(const Iterator &) const;
364
368 bool
369 operator!=(const Iterator &) const;
370
378 bool
379 operator<(const Iterator &) const;
380
385 bool
386 operator>(const Iterator &) const;
387
394 int
395 operator-(const Iterator &p) const;
396
401 operator+(const unsigned int n) const;
402
403 private:
408 };
409
410} // namespace ChunkSparseMatrixIterators
411
413
414namespace std
415{
416 template <typename number, bool Constness>
418 ::ChunkSparseMatrixIterators::Iterator<number, Constness>>
419 {
420 using iterator_category = forward_iterator_tag;
421 using value_type = typename ::ChunkSparseMatrixIterators::
422 Iterator<number, Constness>::value_type;
423 using difference_type = typename ::ChunkSparseMatrixIterators::
424 Iterator<number, Constness>::difference_type;
425 };
426} // namespace std
427
429
430
431
448template <typename number>
449class ChunkSparseMatrix : public virtual Subscriptor
450{
451public:
456
461 using value_type = number;
462
473
479
487
494 struct Traits
495 {
500 static const bool zero_addition_can_be_elided = true;
501 };
502
518
528
542 explicit ChunkSparseMatrix(const ChunkSparsityPattern &sparsity);
543
551 const IdentityMatrix &id);
552
557 virtual ~ChunkSparseMatrix() override;
558
570
579
590 operator=(const double d);
591
605 virtual void
606 reinit(const ChunkSparsityPattern &sparsity);
607
613 virtual void
624 bool
625 empty() const;
626
632 m() const;
633
639 n() const;
640
648
658
669
674 std::size_t
676
687 void
688 set(const size_type i, const size_type j, const number value);
689
695 void
696 add(const size_type i, const size_type j, const number value);
697
707 template <typename number2>
708 void
709 add(const size_type row,
710 const size_type n_cols,
711 const size_type *col_indices,
712 const number2 *values,
713 const bool elide_zero_values = true,
714 const bool col_indices_are_sorted = false);
715
720 operator*=(const number factor);
721
726 operator/=(const number factor);
727
740 void
742
759 template <typename somenumber>
762
780 template <typename ForwardIterator>
781 void
782 copy_from(const ForwardIterator begin, const ForwardIterator end);
783
789 template <typename somenumber>
790 void
792
804 template <typename somenumber>
805 void
806 add(const number factor, const ChunkSparseMatrix<somenumber> &matrix);
807
827 number
828 operator()(const size_type i, const size_type j) const;
829
842 number
843 el(const size_type i, const size_type j) const;
844
854 number
855 diag_element(const size_type i) const;
856
866 void
868 const size_type array_length,
869 size_type &row_length,
870 size_type *column_indices,
871 number *values) const;
872
892 template <class OutVector, class InVector>
893 void
894 vmult(OutVector &dst, const InVector &src) const;
895
911 template <class OutVector, class InVector>
912 void
913 Tvmult(OutVector &dst, const InVector &src) const;
914
929 template <class OutVector, class InVector>
930 void
931 vmult_add(OutVector &dst, const InVector &src) const;
932
948 template <class OutVector, class InVector>
949 void
950 Tvmult_add(OutVector &dst, const InVector &src) const;
951
967 template <typename somenumber>
968 somenumber
970
974 template <typename somenumber>
975 somenumber
977 const Vector<somenumber> &v) const;
985 template <typename somenumber>
986 somenumber
988 const Vector<somenumber> &x,
989 const Vector<somenumber> &b) const;
990
1004 real_type
1005 l1_norm() const;
1006
1014 real_type
1016
1021 real_type
1034 template <typename somenumber>
1035 void
1037 const Vector<somenumber> &src,
1038 const number omega = 1.) const;
1039
1043 template <typename somenumber>
1044 void
1046 const Vector<somenumber> &src,
1047 const number om = 1.) const;
1048
1052 template <typename somenumber>
1053 void
1055 const Vector<somenumber> &src,
1056 const number om = 1.) const;
1057
1061 template <typename somenumber>
1062 void
1064 const Vector<somenumber> &src,
1065 const number om = 1.) const;
1066
1072 template <typename somenumber>
1073 void
1074 SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1075
1080 template <typename somenumber>
1081 void
1082 SOR(Vector<somenumber> &v, const number om = 1.) const;
1083
1088 template <typename somenumber>
1089 void
1090 TSOR(Vector<somenumber> &v, const number om = 1.) const;
1091
1102 template <typename somenumber>
1103 void
1105 const std::vector<size_type> &permutation,
1106 const std::vector<size_type> &inverse_permutation,
1107 const number om = 1.) const;
1108
1119 template <typename somenumber>
1120 void
1122 const std::vector<size_type> &permutation,
1123 const std::vector<size_type> &inverse_permutation,
1124 const number om = 1.) const;
1125
1130 template <typename somenumber>
1131 void
1133 const Vector<somenumber> &b,
1134 const number om = 1.) const;
1135
1140 template <typename somenumber>
1141 void
1143 const Vector<somenumber> &b,
1144 const number om = 1.) const;
1145
1150 template <typename somenumber>
1151 void
1153 const Vector<somenumber> &b,
1154 const number om = 1.) const;
1171 begin() const;
1172
1182 end() const;
1183
1193 iterator
1195
1204 iterator
1206
1222 begin(const unsigned int r) const;
1223
1239 end(const unsigned int r) const;
1240
1255 iterator
1256 begin(const unsigned int r);
1257
1272 iterator
1273 end(const unsigned int r);
1284 void
1285 print(std::ostream &out) const;
1286
1309 void
1310 print_formatted(std::ostream &out,
1311 const unsigned int precision = 3,
1312 const bool scientific = true,
1313 const unsigned int width = 0,
1314 const char *zero_string = " ",
1315 const double denominator = 1.,
1316 const char *separator = " ") const;
1317
1323 void
1324 print_pattern(std::ostream &out, const double threshold = 0.) const;
1325
1336 void
1337 block_write(std::ostream &out) const;
1338
1355 void
1356 block_read(std::istream &in);
1367 int,
1368 int,
1369 << "You are trying to access the matrix entry with index <"
1370 << arg1 << ',' << arg2
1371 << ">, but this entry does not exist in the sparsity pattern "
1372 "of this matrix."
1373 "\n\n"
1374 "The most common cause for this problem is that you used "
1375 "a method to build the sparsity pattern that did not "
1376 "(completely) take into account all of the entries you "
1377 "will later try to write into. An example would be "
1378 "building a sparsity pattern that does not include "
1379 "the entries you will write into due to constraints "
1380 "on degrees of freedom such as hanging nodes or periodic "
1381 "boundary conditions. In such cases, building the "
1382 "sparsity pattern will succeed, but you will get errors "
1383 "such as the current one at one point or other when "
1384 "trying to write into the entries of the matrix.");
1393 int,
1394 int,
1395 << "The iterators denote a range of " << arg1
1396 << " elements, but the given number of rows was " << arg2);
1401 "You are attempting an operation on two vectors that "
1402 "are the same object, but the operation requires that the "
1403 "two objects are in fact different.");
1405private:
1412
1420 std::unique_ptr<number[]> val;
1421
1429
1433 size_type
1434 compute_location(const size_type i, const size_type j) const;
1435
1436 // make all other sparse matrices friends
1437 template <typename somenumber>
1438 friend class ChunkSparseMatrix;
1439
1440 // Also give access to internal details to the iterator/accessor classes.
1441 template <typename, bool>
1443 template <typename, bool>
1445};
1446
1449#ifndef DOXYGEN
1450/*---------------------- Inline functions -----------------------------------*/
1451
1452
1453
1454template <typename number>
1457{
1458 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1459 return cols->rows;
1460}
1461
1462
1463template <typename number>
1466{
1467 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1468 return cols->cols;
1469}
1470
1471
1472
1473template <typename number>
1474inline const ChunkSparsityPattern &
1476{
1477 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1478 return *cols;
1479}
1480
1481
1482
1483template <typename number>
1486 const size_type j) const
1487{
1488 const size_type chunk_size = cols->get_chunk_size();
1489 const size_type chunk_index =
1490 cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1491
1492 if (chunk_index == ChunkSparsityPattern::invalid_entry)
1494 else
1495 {
1496 return (chunk_index * chunk_size * chunk_size +
1497 (i % chunk_size) * chunk_size + (j % chunk_size));
1498 }
1499}
1500
1501
1502template <typename number>
1503inline void
1505 const size_type j,
1506 const number value)
1507{
1509
1510 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1511 // it is allowed to set elements of the matrix that are not part of the
1512 // sparsity pattern, if the value to which we set it is zero
1513 const size_type index = compute_location(i, j);
1514 Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1515 ExcInvalidIndex(i, j));
1516
1517 if (index != SparsityPattern::invalid_entry)
1518 val[index] = value;
1519}
1520
1521
1522
1523template <typename number>
1524inline void
1526 const size_type j,
1527 const number value)
1528{
1530
1531 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1532
1533 if (std::abs(value) != 0.)
1534 {
1535 const size_type index = compute_location(i, j);
1537 ExcInvalidIndex(i, j));
1538
1539 val[index] += value;
1540 }
1541}
1542
1543
1544
1545template <typename number>
1546template <typename number2>
1547inline void
1549 const size_type n_cols,
1550 const size_type *col_indices,
1551 const number2 *values,
1552 const bool /*elide_zero_values*/,
1553 const bool /*col_indices_are_sorted*/)
1554{
1555 // TODO: could be done more efficiently...
1556 for (size_type col = 0; col < n_cols; ++col)
1557 add(row, col_indices[col], static_cast<number>(values[col]));
1558}
1559
1560
1561
1562template <typename number>
1564ChunkSparseMatrix<number>::operator*=(const number factor)
1565{
1566 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1567 Assert(val != nullptr, ExcNotInitialized());
1568
1569 const size_type chunk_size = cols->get_chunk_size();
1570
1571 // multiply all elements of the matrix with the given factor. this includes
1572 // the padding elements in chunks that overlap the boundaries of the actual
1573 // matrix -- but since multiplication with a number does not violate the
1574 // invariant of keeping these elements at zero nothing can happen
1575 number *val_ptr = val.get();
1576 const number *const end_ptr =
1577 val.get() +
1578 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1579 while (val_ptr != end_ptr)
1580 *val_ptr++ *= factor;
1581
1582 return *this;
1583}
1584
1585
1586
1587template <typename number>
1589ChunkSparseMatrix<number>::operator/=(const number factor)
1590{
1591 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1592 Assert(val != nullptr, ExcNotInitialized());
1593 Assert(std::abs(factor) != 0, ExcDivideByZero());
1594
1595 const number factor_inv = 1. / factor;
1596
1597 const size_type chunk_size = cols->get_chunk_size();
1598
1599 // multiply all elements of the matrix with the given factor. this includes
1600 // the padding elements in chunks that overlap the boundaries of the actual
1601 // matrix -- but since multiplication with a number does not violate the
1602 // invariant of keeping these elements at zero nothing can happen
1603 number *val_ptr = val.get();
1604 const number *const end_ptr =
1605 val.get() +
1606 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1607
1608 while (val_ptr != end_ptr)
1609 *val_ptr++ *= factor_inv;
1610
1611 return *this;
1612}
1613
1614
1615
1616template <typename number>
1617inline number
1619 const size_type j) const
1620{
1621 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1622 AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1623 ExcInvalidIndex(i, j));
1624 return val[compute_location(i, j)];
1625}
1626
1627
1628
1629template <typename number>
1630inline number
1632{
1633 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1634 const size_type index = compute_location(i, j);
1635
1637 return val[index];
1638 else
1639 return 0;
1640}
1641
1642
1643
1644template <typename number>
1645inline number
1647{
1648 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1649 Assert(m() == n(), ExcNotQuadratic());
1650 AssertIndexRange(i, m());
1651
1652 // Use that the first element in each row of a quadratic matrix is the main
1653 // diagonal of the chunk sparsity pattern
1654 const size_type chunk_size = cols->get_chunk_size();
1655 return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1656 chunk_size +
1657 (i % chunk_size) * chunk_size + (i % chunk_size)];
1658}
1659
1660
1661
1662template <typename number>
1663template <typename ForwardIterator>
1664inline void
1665ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1666 const ForwardIterator end)
1667{
1668 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1669 ExcIteratorRange(std::distance(begin, end), m()));
1670
1671 // for use in the inner loop, we define an alias to the type of the inner
1672 // iterators
1673 using inner_iterator =
1674 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1675 size_type row = 0;
1676 for (ForwardIterator i = begin; i != end; ++i, ++row)
1677 {
1678 const inner_iterator end_of_row = i->end();
1679 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1680 // write entries
1681 set(row, j->first, j->second);
1682 }
1683}
1684
1685
1686
1687//---------------------------------------------------------------------------
1688
1689
1691{
1692 template <typename number>
1693 inline Accessor<number, true>::Accessor(const MatrixType *matrix,
1694 const unsigned int row)
1695 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1696 row)
1697 , matrix(matrix)
1698 {}
1699
1700
1701
1702 template <typename number>
1703 inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1704 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1705 , matrix(matrix)
1706 {}
1707
1708
1709
1710 template <typename number>
1711 inline Accessor<number, true>::Accessor(
1713 : ChunkSparsityPatternIterators::Accessor(a)
1714 , matrix(&a.get_matrix())
1715 {}
1716
1717
1718
1719 template <typename number>
1720 inline number
1721 Accessor<number, true>::value() const
1722 {
1723 const unsigned int chunk_size =
1724 matrix->get_sparsity_pattern().get_chunk_size();
1725 return matrix->val[reduced_index() * chunk_size * chunk_size +
1726 chunk_row * chunk_size + chunk_col];
1727 }
1728
1729
1730
1731 template <typename number>
1732 inline const typename Accessor<number, true>::MatrixType &
1733 Accessor<number, true>::get_matrix() const
1734 {
1735 return *matrix;
1736 }
1737
1738
1739
1740 template <typename number>
1741 inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1742 const bool)
1743 : accessor(accessor)
1744 {}
1745
1746
1747 template <typename number>
1748 inline Accessor<number, false>::Reference::operator number() const
1749 {
1750 const unsigned int chunk_size =
1751 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1752 return accessor->matrix
1753 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1754 accessor->chunk_row * chunk_size + accessor->chunk_col];
1755 }
1756
1757
1758
1759 template <typename number>
1760 inline const typename Accessor<number, false>::Reference &
1761 Accessor<number, false>::Reference::operator=(const number n) const
1762 {
1763 const unsigned int chunk_size =
1764 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1765 accessor->matrix
1766 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1767 accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1768 return *this;
1769 }
1770
1771
1772
1773 template <typename number>
1774 inline const typename Accessor<number, false>::Reference &
1775 Accessor<number, false>::Reference::operator+=(const number n) const
1776 {
1777 const unsigned int chunk_size =
1778 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1779 accessor->matrix
1780 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1781 accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1782 return *this;
1783 }
1784
1785
1786
1787 template <typename number>
1788 inline const typename Accessor<number, false>::Reference &
1789 Accessor<number, false>::Reference::operator-=(const number n) const
1790 {
1791 const unsigned int chunk_size =
1792 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1793 accessor->matrix
1794 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1795 accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1796 return *this;
1797 }
1798
1799
1800
1801 template <typename number>
1802 inline const typename Accessor<number, false>::Reference &
1803 Accessor<number, false>::Reference::operator*=(const number n) const
1804 {
1805 const unsigned int chunk_size =
1806 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1807 accessor->matrix
1808 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1809 accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1810 return *this;
1811 }
1812
1813
1814
1815 template <typename number>
1816 inline const typename Accessor<number, false>::Reference &
1817 Accessor<number, false>::Reference::operator/=(const number n) const
1818 {
1819 const unsigned int chunk_size =
1820 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1821 accessor->matrix
1822 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1823 accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1824 return *this;
1825 }
1826
1827
1828
1829 template <typename number>
1830 inline Accessor<number, false>::Accessor(MatrixType *matrix,
1831 const unsigned int row)
1832 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1833 row)
1834 , matrix(matrix)
1835 {}
1836
1837
1838
1839 template <typename number>
1840 inline Accessor<number, false>::Accessor(MatrixType *matrix)
1841 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1842 , matrix(matrix)
1843 {}
1844
1845
1846
1847 template <typename number>
1848 inline typename Accessor<number, false>::Reference
1849 Accessor<number, false>::value() const
1850 {
1851 return Reference(this, true);
1852 }
1853
1854
1855
1856 template <typename number>
1857 inline typename Accessor<number, false>::MatrixType &
1858 Accessor<number, false>::get_matrix() const
1859 {
1860 return *matrix;
1861 }
1862
1863
1864
1865 template <typename number, bool Constness>
1866 inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
1867 const unsigned int row)
1868 : accessor(matrix, row)
1869 {}
1870
1871
1872
1873 template <typename number, bool Constness>
1874 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1875 : accessor(matrix)
1876 {}
1877
1878
1879
1880 template <typename number, bool Constness>
1881 inline Iterator<number, Constness>::Iterator(
1883 : accessor(*i)
1884 {}
1885
1886
1887
1888 template <typename number, bool Constness>
1889 inline Iterator<number, Constness> &
1890 Iterator<number, Constness>::operator++()
1891 {
1892 accessor.advance();
1893 return *this;
1894 }
1895
1896
1897 template <typename number, bool Constness>
1898 inline Iterator<number, Constness>
1899 Iterator<number, Constness>::operator++(int)
1900 {
1901 const Iterator iter = *this;
1902 accessor.advance();
1903 return iter;
1904 }
1905
1906
1907 template <typename number, bool Constness>
1908 inline const Accessor<number, Constness> &
1909 Iterator<number, Constness>::operator*() const
1910 {
1911 return accessor;
1912 }
1913
1914
1915 template <typename number, bool Constness>
1916 inline const Accessor<number, Constness> *
1917 Iterator<number, Constness>::operator->() const
1918 {
1919 return &accessor;
1920 }
1921
1922
1923 template <typename number, bool Constness>
1924 inline bool
1925 Iterator<number, Constness>::operator==(const Iterator &other) const
1926 {
1927 return (accessor == other.accessor);
1928 }
1929
1930
1931 template <typename number, bool Constness>
1932 inline bool
1933 Iterator<number, Constness>::operator!=(const Iterator &other) const
1934 {
1935 return !(*this == other);
1936 }
1937
1938
1939 template <typename number, bool Constness>
1940 inline bool
1941 Iterator<number, Constness>::operator<(const Iterator &other) const
1942 {
1943 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1945
1946 return (accessor < other.accessor);
1947 }
1948
1949
1950 template <typename number, bool Constness>
1951 inline bool
1952 Iterator<number, Constness>::operator>(const Iterator &other) const
1953 {
1954 return (other < *this);
1955 }
1956
1957
1958 template <typename number, bool Constness>
1959 inline int
1960 Iterator<number, Constness>::operator-(const Iterator &other) const
1961 {
1962 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1964
1965 // TODO: can be optimized
1966 int difference = 0;
1967 if (*this < other)
1968 {
1969 Iterator copy = *this;
1970 while (copy != other)
1971 {
1972 ++copy;
1973 --difference;
1974 }
1975 }
1976 else
1977 {
1978 Iterator copy = other;
1979 while (copy != *this)
1980 {
1981 ++copy;
1982 ++difference;
1983 }
1984 }
1985 return difference;
1986 }
1987
1988
1989
1990 template <typename number, bool Constness>
1991 inline Iterator<number, Constness>
1992 Iterator<number, Constness>::operator+(const unsigned int n) const
1993 {
1994 Iterator x = *this;
1995 for (unsigned int i = 0; i < n; ++i)
1996 ++x;
1997
1998 return x;
1999 }
2000
2001} // namespace ChunkSparseMatrixIterators
2002
2003
2004
2005template <typename number>
2008{
2009 return const_iterator(this, 0);
2010}
2011
2012
2013template <typename number>
2016{
2017 return const_iterator(this);
2018}
2019
2020
2021template <typename number>
2024{
2025 return iterator(this, 0);
2026}
2027
2028
2029template <typename number>
2032{
2033 return iterator(this);
2034}
2035
2036
2037template <typename number>
2039ChunkSparseMatrix<number>::begin(const unsigned int r) const
2040{
2041 AssertIndexRange(r, m());
2042 return const_iterator(this, r);
2043}
2044
2045
2046
2047template <typename number>
2049ChunkSparseMatrix<number>::end(const unsigned int r) const
2050{
2051 AssertIndexRange(r, m());
2052 return const_iterator(this, r + 1);
2053}
2054
2055
2056
2057template <typename number>
2059ChunkSparseMatrix<number>::begin(const unsigned int r)
2060{
2061 AssertIndexRange(r, m());
2062 return iterator(this, r);
2063}
2064
2065
2066
2067template <typename number>
2069ChunkSparseMatrix<number>::end(const unsigned int r)
2070{
2071 AssertIndexRange(r, m());
2072 return iterator(this, r + 1);
2073}
2074
2075
2076
2077#endif // DOXYGEN
2078
2080
2081#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:472
static ::ExceptionBase & ExcDifferentChunkSparsityPatterns()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:540
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:495
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:35
@ matrix
Contents is actually a matrix.
void copy(const T *begin, const T *end, U *dest)
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:611
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:82
static const bool zero_addition_can_be_elided
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
typename ::ChunkSparseMatrixIterators::Iterator< number, Constness >::value_type value_type