Reference documentation for deal.II version 9.5.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// 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
1307 void
1308 print_formatted(std::ostream & out,
1309 const unsigned int precision = 3,
1310 const bool scientific = true,
1311 const unsigned int width = 0,
1312 const char * zero_string = " ",
1313 const double denominator = 1.) const;
1314
1320 void
1321 print_pattern(std::ostream &out, const double threshold = 0.) const;
1322
1333 void
1334 block_write(std::ostream &out) const;
1335
1352 void
1353 block_read(std::istream &in);
1364 int,
1365 int,
1366 << "You are trying to access the matrix entry with index <"
1367 << arg1 << ',' << arg2
1368 << ">, but this entry does not exist in the sparsity pattern "
1369 "of this matrix."
1370 "\n\n"
1371 "The most common cause for this problem is that you used "
1372 "a method to build the sparsity pattern that did not "
1373 "(completely) take into account all of the entries you "
1374 "will later try to write into. An example would be "
1375 "building a sparsity pattern that does not include "
1376 "the entries you will write into due to constraints "
1377 "on degrees of freedom such as hanging nodes or periodic "
1378 "boundary conditions. In such cases, building the "
1379 "sparsity pattern will succeed, but you will get errors "
1380 "such as the current one at one point or other when "
1381 "trying to write into the entries of the matrix.");
1390 int,
1391 int,
1392 << "The iterators denote a range of " << arg1
1393 << " elements, but the given number of rows was " << arg2);
1398 "You are attempting an operation on two vectors that "
1399 "are the same object, but the operation requires that the "
1400 "two objects are in fact different.");
1402private:
1409
1417 std::unique_ptr<number[]> val;
1418
1426
1430 size_type
1431 compute_location(const size_type i, const size_type j) const;
1432
1433 // make all other sparse matrices friends
1434 template <typename somenumber>
1435 friend class ChunkSparseMatrix;
1436
1437 // Also give access to internal details to the iterator/accessor classes.
1438 template <typename, bool>
1440 template <typename, bool>
1442};
1443
1446#ifndef DOXYGEN
1447/*---------------------- Inline functions -----------------------------------*/
1448
1449
1450
1451template <typename number>
1454{
1455 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1456 return cols->rows;
1457}
1458
1459
1460template <typename number>
1463{
1464 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1465 return cols->cols;
1466}
1467
1468
1469
1470template <typename number>
1471inline const ChunkSparsityPattern &
1473{
1474 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1475 return *cols;
1476}
1477
1478
1479
1480template <typename number>
1483 const size_type j) const
1484{
1485 const size_type chunk_size = cols->get_chunk_size();
1486 const size_type chunk_index =
1487 cols->sparsity_pattern(i / chunk_size, j / chunk_size);
1488
1489 if (chunk_index == ChunkSparsityPattern::invalid_entry)
1491 else
1492 {
1493 return (chunk_index * chunk_size * chunk_size +
1494 (i % chunk_size) * chunk_size + (j % chunk_size));
1495 }
1496}
1497
1498
1499template <typename number>
1500inline void
1502 const size_type j,
1503 const number value)
1504{
1506
1507 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1508 // it is allowed to set elements of the matrix that are not part of the
1509 // sparsity pattern, if the value to which we set it is zero
1510 const size_type index = compute_location(i, j);
1511 Assert((index != SparsityPattern::invalid_entry) || (value == 0.),
1512 ExcInvalidIndex(i, j));
1513
1514 if (index != SparsityPattern::invalid_entry)
1515 val[index] = value;
1516}
1517
1518
1519
1520template <typename number>
1521inline void
1523 const size_type j,
1524 const number value)
1525{
1527
1528 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1529
1530 if (std::abs(value) != 0.)
1531 {
1532 const size_type index = compute_location(i, j);
1534 ExcInvalidIndex(i, j));
1535
1536 val[index] += value;
1537 }
1538}
1539
1540
1541
1542template <typename number>
1543template <typename number2>
1544inline void
1546 const size_type n_cols,
1547 const size_type *col_indices,
1548 const number2 * values,
1549 const bool /*elide_zero_values*/,
1550 const bool /*col_indices_are_sorted*/)
1551{
1552 // TODO: could be done more efficiently...
1553 for (size_type col = 0; col < n_cols; ++col)
1554 add(row, col_indices[col], static_cast<number>(values[col]));
1555}
1556
1557
1558
1559template <typename number>
1561ChunkSparseMatrix<number>::operator*=(const number factor)
1562{
1563 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1564 Assert(val != nullptr, ExcNotInitialized());
1565
1566 const size_type chunk_size = cols->get_chunk_size();
1567
1568 // multiply all elements of the matrix with the given factor. this includes
1569 // the padding elements in chunks that overlap the boundaries of the actual
1570 // matrix -- but since multiplication with a number does not violate the
1571 // invariant of keeping these elements at zero nothing can happen
1572 number * val_ptr = val.get();
1573 const number *const end_ptr =
1574 val.get() +
1575 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1576 while (val_ptr != end_ptr)
1577 *val_ptr++ *= factor;
1578
1579 return *this;
1580}
1581
1582
1583
1584template <typename number>
1586ChunkSparseMatrix<number>::operator/=(const number factor)
1587{
1588 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1589 Assert(val != nullptr, ExcNotInitialized());
1590 Assert(std::abs(factor) != 0, ExcDivideByZero());
1591
1592 const number factor_inv = 1. / factor;
1593
1594 const size_type chunk_size = cols->get_chunk_size();
1595
1596 // multiply all elements of the matrix with the given factor. this includes
1597 // the padding elements in chunks that overlap the boundaries of the actual
1598 // matrix -- but since multiplication with a number does not violate the
1599 // invariant of keeping these elements at zero nothing can happen
1600 number * val_ptr = val.get();
1601 const number *const end_ptr =
1602 val.get() +
1603 cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size;
1604
1605 while (val_ptr != end_ptr)
1606 *val_ptr++ *= factor_inv;
1607
1608 return *this;
1609}
1610
1611
1612
1613template <typename number>
1614inline number
1616 const size_type j) const
1617{
1618 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1619 AssertThrow(compute_location(i, j) != SparsityPattern::invalid_entry,
1620 ExcInvalidIndex(i, j));
1621 return val[compute_location(i, j)];
1622}
1623
1624
1625
1626template <typename number>
1627inline number
1629{
1630 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1631 const size_type index = compute_location(i, j);
1632
1634 return val[index];
1635 else
1636 return 0;
1637}
1638
1639
1640
1641template <typename number>
1642inline number
1644{
1645 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1646 Assert(m() == n(), ExcNotQuadratic());
1647 AssertIndexRange(i, m());
1648
1649 // Use that the first element in each row of a quadratic matrix is the main
1650 // diagonal of the chunk sparsity pattern
1651 const size_type chunk_size = cols->get_chunk_size();
1652 return val[cols->sparsity_pattern.rowstart[i / chunk_size] * chunk_size *
1653 chunk_size +
1654 (i % chunk_size) * chunk_size + (i % chunk_size)];
1655}
1656
1657
1658
1659template <typename number>
1660template <typename ForwardIterator>
1661inline void
1662ChunkSparseMatrix<number>::copy_from(const ForwardIterator begin,
1663 const ForwardIterator end)
1664{
1665 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
1666 ExcIteratorRange(std::distance(begin, end), m()));
1667
1668 // for use in the inner loop, we define an alias to the type of the inner
1669 // iterators
1670 using inner_iterator =
1671 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
1672 size_type row = 0;
1673 for (ForwardIterator i = begin; i != end; ++i, ++row)
1674 {
1675 const inner_iterator end_of_row = i->end();
1676 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
1677 // write entries
1678 set(row, j->first, j->second);
1679 }
1680}
1681
1682
1683
1684//---------------------------------------------------------------------------
1685
1686
1688{
1689 template <typename number>
1690 inline Accessor<number, true>::Accessor(const MatrixType * matrix,
1691 const unsigned int row)
1692 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1693 row)
1694 , matrix(matrix)
1695 {}
1696
1697
1698
1699 template <typename number>
1700 inline Accessor<number, true>::Accessor(const MatrixType *matrix)
1701 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1702 , matrix(matrix)
1703 {}
1704
1705
1706
1707 template <typename number>
1708 inline Accessor<number, true>::Accessor(
1710 : ChunkSparsityPatternIterators::Accessor(a)
1711 , matrix(&a.get_matrix())
1712 {}
1713
1714
1715
1716 template <typename number>
1717 inline number
1718 Accessor<number, true>::value() const
1719 {
1720 const unsigned int chunk_size =
1721 matrix->get_sparsity_pattern().get_chunk_size();
1722 return matrix->val[reduced_index() * chunk_size * chunk_size +
1723 chunk_row * chunk_size + chunk_col];
1724 }
1725
1726
1727
1728 template <typename number>
1729 inline const typename Accessor<number, true>::MatrixType &
1730 Accessor<number, true>::get_matrix() const
1731 {
1732 return *matrix;
1733 }
1734
1735
1736
1737 template <typename number>
1738 inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
1739 const bool)
1740 : accessor(accessor)
1741 {}
1742
1743
1744 template <typename number>
1745 inline Accessor<number, false>::Reference::operator number() const
1746 {
1747 const unsigned int chunk_size =
1748 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1749 return accessor->matrix
1750 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1751 accessor->chunk_row * chunk_size + accessor->chunk_col];
1752 }
1753
1754
1755
1756 template <typename number>
1757 inline const typename Accessor<number, false>::Reference &
1758 Accessor<number, false>::Reference::operator=(const number n) const
1759 {
1760 const unsigned int chunk_size =
1761 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1762 accessor->matrix
1763 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1764 accessor->chunk_row * chunk_size + accessor->chunk_col] = n;
1765 return *this;
1766 }
1767
1768
1769
1770 template <typename number>
1771 inline const typename Accessor<number, false>::Reference &
1772 Accessor<number, false>::Reference::operator+=(const number n) const
1773 {
1774 const unsigned int chunk_size =
1775 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1776 accessor->matrix
1777 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1778 accessor->chunk_row * chunk_size + accessor->chunk_col] += n;
1779 return *this;
1780 }
1781
1782
1783
1784 template <typename number>
1785 inline const typename Accessor<number, false>::Reference &
1786 Accessor<number, false>::Reference::operator-=(const number n) const
1787 {
1788 const unsigned int chunk_size =
1789 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1790 accessor->matrix
1791 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1792 accessor->chunk_row * chunk_size + accessor->chunk_col] -= n;
1793 return *this;
1794 }
1795
1796
1797
1798 template <typename number>
1799 inline const typename Accessor<number, false>::Reference &
1800 Accessor<number, false>::Reference::operator*=(const number n) const
1801 {
1802 const unsigned int chunk_size =
1803 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1804 accessor->matrix
1805 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1806 accessor->chunk_row * chunk_size + accessor->chunk_col] *= n;
1807 return *this;
1808 }
1809
1810
1811
1812 template <typename number>
1813 inline const typename Accessor<number, false>::Reference &
1814 Accessor<number, false>::Reference::operator/=(const number n) const
1815 {
1816 const unsigned int chunk_size =
1817 accessor->matrix->get_sparsity_pattern().get_chunk_size();
1818 accessor->matrix
1819 ->val[accessor->reduced_index() * chunk_size * chunk_size +
1820 accessor->chunk_row * chunk_size + accessor->chunk_col] /= n;
1821 return *this;
1822 }
1823
1824
1825
1826 template <typename number>
1827 inline Accessor<number, false>::Accessor(MatrixType * matrix,
1828 const unsigned int row)
1829 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
1830 row)
1831 , matrix(matrix)
1832 {}
1833
1834
1835
1836 template <typename number>
1837 inline Accessor<number, false>::Accessor(MatrixType *matrix)
1838 : ChunkSparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
1839 , matrix(matrix)
1840 {}
1841
1842
1843
1844 template <typename number>
1845 inline typename Accessor<number, false>::Reference
1846 Accessor<number, false>::value() const
1847 {
1848 return Reference(this, true);
1849 }
1850
1851
1852
1853 template <typename number>
1854 inline typename Accessor<number, false>::MatrixType &
1855 Accessor<number, false>::get_matrix() const
1856 {
1857 return *matrix;
1858 }
1859
1860
1861
1862 template <typename number, bool Constness>
1863 inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
1864 const unsigned int row)
1865 : accessor(matrix, row)
1866 {}
1867
1868
1869
1870 template <typename number, bool Constness>
1871 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
1872 : accessor(matrix)
1873 {}
1874
1875
1876
1877 template <typename number, bool Constness>
1878 inline Iterator<number, Constness>::Iterator(
1880 : accessor(*i)
1881 {}
1882
1883
1884
1885 template <typename number, bool Constness>
1886 inline Iterator<number, Constness> &
1887 Iterator<number, Constness>::operator++()
1888 {
1889 accessor.advance();
1890 return *this;
1891 }
1892
1893
1894 template <typename number, bool Constness>
1895 inline Iterator<number, Constness>
1896 Iterator<number, Constness>::operator++(int)
1897 {
1898 const Iterator iter = *this;
1899 accessor.advance();
1900 return iter;
1901 }
1902
1903
1904 template <typename number, bool Constness>
1905 inline const Accessor<number, Constness> &
1906 Iterator<number, Constness>::operator*() const
1907 {
1908 return accessor;
1909 }
1910
1911
1912 template <typename number, bool Constness>
1913 inline const Accessor<number, Constness> *
1914 Iterator<number, Constness>::operator->() const
1915 {
1916 return &accessor;
1917 }
1918
1919
1920 template <typename number, bool Constness>
1921 inline bool
1922 Iterator<number, Constness>::operator==(const Iterator &other) const
1923 {
1924 return (accessor == other.accessor);
1925 }
1926
1927
1928 template <typename number, bool Constness>
1929 inline bool
1930 Iterator<number, Constness>::operator!=(const Iterator &other) const
1931 {
1932 return !(*this == other);
1933 }
1934
1935
1936 template <typename number, bool Constness>
1937 inline bool
1938 Iterator<number, Constness>::operator<(const Iterator &other) const
1939 {
1940 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1942
1943 return (accessor < other.accessor);
1944 }
1945
1946
1947 template <typename number, bool Constness>
1948 inline bool
1949 Iterator<number, Constness>::operator>(const Iterator &other) const
1950 {
1951 return (other < *this);
1952 }
1953
1954
1955 template <typename number, bool Constness>
1956 inline int
1957 Iterator<number, Constness>::operator-(const Iterator &other) const
1958 {
1959 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
1961
1962 // TODO: can be optimized
1963 int difference = 0;
1964 if (*this < other)
1965 {
1966 Iterator copy = *this;
1967 while (copy != other)
1968 {
1969 ++copy;
1970 --difference;
1971 }
1972 }
1973 else
1974 {
1975 Iterator copy = other;
1976 while (copy != *this)
1977 {
1978 ++copy;
1979 ++difference;
1980 }
1981 }
1982 return difference;
1983 }
1984
1985
1986
1987 template <typename number, bool Constness>
1988 inline Iterator<number, Constness>
1989 Iterator<number, Constness>::operator+(const unsigned int n) const
1990 {
1991 Iterator x = *this;
1992 for (unsigned int i = 0; i < n; ++i)
1993 ++x;
1994
1995 return x;
1996 }
1997
1998} // namespace ChunkSparseMatrixIterators
1999
2000
2001
2002template <typename number>
2005{
2006 return const_iterator(this, 0);
2007}
2008
2009
2010template <typename number>
2013{
2014 return const_iterator(this);
2015}
2016
2017
2018template <typename number>
2021{
2022 return iterator(this, 0);
2023}
2024
2025
2026template <typename number>
2029{
2030 return iterator(this);
2031}
2032
2033
2034template <typename number>
2036ChunkSparseMatrix<number>::begin(const unsigned int r) const
2037{
2038 AssertIndexRange(r, m());
2039 return const_iterator(this, r);
2040}
2041
2042
2043
2044template <typename number>
2046ChunkSparseMatrix<number>::end(const unsigned int r) const
2047{
2048 AssertIndexRange(r, m());
2049 return const_iterator(this, r + 1);
2050}
2051
2052
2053
2054template <typename number>
2056ChunkSparseMatrix<number>::begin(const unsigned int r)
2057{
2058 AssertIndexRange(r, m());
2059 return iterator(this, r);
2060}
2061
2062
2063
2064template <typename number>
2066ChunkSparseMatrix<number>::end(const unsigned int r)
2067{
2068 AssertIndexRange(r, m());
2069 return iterator(this, r + 1);
2070}
2071
2072
2073
2074#endif // DOXYGEN
2075
2077
2078#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
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
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) 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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
__global__ void set(Number *val, const Number s, const size_type N)
#define DeclException0(Exception0)
Definition exceptions.h:465
static ::ExceptionBase & ExcDifferentChunkSparsityPatterns()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define Assert(cond, exc)
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
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.
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: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