deal.II version GIT relicensing-2238-gc05b561aad 2024-12-10 20:50:00+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
sparse_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2002 - 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_sparse_matrix_h
16#define dealii_sparse_matrix_h
17
18#include <deal.II/base/config.h>
19
23
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;
41template <typename Matrix>
42class BlockMatrixBase;
43template <typename number>
44class SparseILU;
45# ifdef DEAL_II_WITH_MPI
46namespace Utilities
47{
48 namespace MPI
49 {
50 template <typename Number>
51 void
53 }
54} // namespace Utilities
55# endif
56
57# ifdef DEAL_II_WITH_TRILINOS
58namespace TrilinosWrappers
59{
60 class SparseMatrix;
61}
62# endif
63#endif
64
75{
80
81 // forward declaration
82 template <typename number, bool Constness>
83 class Iterator;
84
95 template <typename number, bool Constness>
97 {
98 public:
102 number
103 value() const;
104
108 number &
110
116 get_matrix() const;
117 };
118
119
120
127 template <typename number>
128 class Accessor<number, true> : public SparsityPatternIterators::Accessor
129 {
130 public:
136
140 Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
141
146
151
155 number
156 value() const;
157
162 const MatrixType &
163 get_matrix() const;
164
165 private:
170
175
176 // Make iterator class a friend.
177 template <typename, bool>
178 friend class Iterator;
179 };
180
181
188 template <typename number>
189 class Accessor<number, false> : public SparsityPatternIterators::Accessor
190 {
191 private:
216 class Reference
217 {
218 public:
223 Reference(const Accessor *accessor, const bool dummy);
224
228 operator number() const;
229
233 const Reference &
234 operator=(const number n) const;
235
239 const Reference &
240 operator+=(const number n) const;
241
245 const Reference &
246 operator-=(const number n) const;
247
251 const Reference &
252 operator*=(const number n) const;
253
257 const Reference &
258 operator/=(const number n) const;
259
260 private:
266 };
267
268 public:
274
278 Accessor(MatrixType *matrix, const std::size_t index);
279
284
288 Reference
289 value() const;
290
295 MatrixType &
296 get_matrix() const;
297
298 private:
303
308
309 // Make iterator class a friend.
310 template <typename, bool>
311 friend class Iterator;
312 };
313
314
315
345 template <typename number, bool Constness>
347 {
348 public:
353
359
365
370 Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
371
376
382
388
392 Iterator &
394
400
405 operator*() const;
406
411 operator->() const;
412
416 bool
417 operator==(const Iterator &) const;
418
422 bool
423 operator!=(const Iterator &) const;
424
432 bool
433 operator<(const Iterator &) const;
434
439 bool
440 operator>(const Iterator &) const;
441
448 int
449 operator-(const Iterator &p) const;
450
455 operator+(const size_type n) const;
456
457 private:
462 };
463
464} // namespace SparseMatrixIterators
465
467
468namespace std
469{
470 template <typename number, bool Constness>
472 ::SparseMatrixIterators::Iterator<number, Constness>>
473 {
474 using iterator_category = forward_iterator_tag;
476 typename ::SparseMatrixIterators::Iterator<number,
477 Constness>::value_type;
478 using difference_type = typename ::SparseMatrixIterators::
479 Iterator<number, Constness>::difference_type;
480 };
481} // namespace std
482
484
490// TODO: Add multithreading to the other vmult functions.
491
518template <typename number>
520{
521public:
526
531 using value_type = number;
532
543
549
557
564 struct Traits
565 {
570 static const bool zero_addition_can_be_elided = true;
571 };
572
588
598
607
621 explicit SparseMatrix(const SparsityPattern &sparsity);
622
629 SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
630
635 virtual ~SparseMatrix() override;
636
648
655
664
677 operator=(const double d);
678
692 virtual void
693 reinit(const SparsityPattern &sparsity);
694
701 template <typename number2>
702 void
703 reinit(const SparseMatrix<number2> &sparse_matrix);
704
710 virtual void
721 bool
722 empty() const;
723
729 m() const;
730
736 n() const;
737
742 get_row_length(const size_type row) const;
743
749 std::size_t
751
761 std::size_t
762 n_actually_nonzero_elements(const double threshold = 0.) const;
763
772 const SparsityPattern &
774
779 std::size_t
781
786
797 void
798 set(const size_type i, const size_type j, const number value);
799
815 template <typename number2>
816 void
817 set(const std::vector<size_type> &indices,
818 const FullMatrix<number2> &full_matrix,
819 const bool elide_zero_values = false);
820
826 template <typename number2>
827 void
828 set(const std::vector<size_type> &row_indices,
829 const std::vector<size_type> &col_indices,
830 const FullMatrix<number2> &full_matrix,
831 const bool elide_zero_values = false);
832
843 template <typename number2>
844 void
845 set(const size_type row,
846 const std::vector<size_type> &col_indices,
847 const std::vector<number2> &values,
848 const bool elide_zero_values = false);
849
859 template <typename number2>
860 void
861 set(const size_type row,
862 const size_type n_cols,
863 const size_type *col_indices,
864 const number2 *values,
865 const bool elide_zero_values = false);
866
872 void
873 add(const size_type i, const size_type j, const number value);
874
889 template <typename number2>
890 void
891 add(const std::vector<size_type> &indices,
892 const FullMatrix<number2> &full_matrix,
893 const bool elide_zero_values = true);
894
900 template <typename number2>
901 void
902 add(const std::vector<size_type> &row_indices,
903 const std::vector<size_type> &col_indices,
904 const FullMatrix<number2> &full_matrix,
905 const bool elide_zero_values = true);
906
916 template <typename number2>
917 void
918 add(const size_type row,
919 const std::vector<size_type> &col_indices,
920 const std::vector<number2> &values,
921 const bool elide_zero_values = true);
922
932 template <typename number2>
933 void
934 add(const size_type row,
935 const size_type n_cols,
936 const size_type *col_indices,
937 const number2 *values,
938 const bool elide_zero_values = true,
939 const bool col_indices_are_sorted = false);
940
945 operator*=(const number factor);
946
951 operator/=(const number factor);
952
965 void
967
984 template <typename somenumber>
987
1004 template <typename ForwardIterator>
1005 void
1006 copy_from(const ForwardIterator begin, const ForwardIterator end);
1007
1017 template <typename somenumber>
1018 void
1020
1021#ifdef DEAL_II_WITH_TRILINOS
1033#endif
1034
1046 template <typename somenumber>
1047 void
1048 add(const number factor, const SparseMatrix<somenumber> &matrix);
1049
1069 const number &
1070 operator()(const size_type i, const size_type j) const;
1071
1075 number &
1076 operator()(const size_type i, const size_type j);
1077
1090 number
1091 el(const size_type i, const size_type j) const;
1092
1102 number
1103 diag_element(const size_type i) const;
1104
1109 number &
1111
1133 template <class OutVector, class InVector>
1134 void
1135 vmult(OutVector &dst, const InVector &src) const;
1136
1152 template <class OutVector, class InVector>
1153 void
1154 Tvmult(OutVector &dst, const InVector &src) const;
1155
1172 template <class OutVector, class InVector>
1173 void
1174 vmult_add(OutVector &dst, const InVector &src) const;
1175
1191 template <class OutVector, class InVector>
1192 void
1193 Tvmult_add(OutVector &dst, const InVector &src) const;
1194
1212 template <typename somenumber>
1213 somenumber
1215
1221 template <typename somenumber>
1222 somenumber
1224 const Vector<somenumber> &v) const;
1225
1235 template <typename somenumber>
1236 somenumber
1238 const Vector<somenumber> &x,
1239 const Vector<somenumber> &b) const;
1240
1276 template <typename numberB, typename numberC>
1277 void
1279 const SparseMatrix<numberB> &B,
1280 const Vector<number> &V = Vector<number>(),
1281 const bool rebuild_sparsity_pattern = true) const;
1282
1307 template <typename numberB, typename numberC>
1308 void
1310 const SparseMatrix<numberB> &B,
1311 const Vector<number> &V = Vector<number>(),
1312 const bool rebuild_sparsity_pattern = true) const;
1313
1327 real_type
1328 l1_norm() const;
1329
1337 real_type
1339
1344 real_type
1357 template <typename somenumber>
1358 void
1360 const Vector<somenumber> &src,
1361 const number omega = 1.) const;
1362
1369 template <typename somenumber>
1370 void
1372 const Vector<somenumber> &src,
1373 const number omega = 1.,
1374 const std::vector<std::size_t> &pos_right_of_diagonal =
1375 std::vector<std::size_t>()) const;
1376
1380 template <typename somenumber>
1381 void
1383 const Vector<somenumber> &src,
1384 const number omega = 1.) const;
1385
1389 template <typename somenumber>
1390 void
1392 const Vector<somenumber> &src,
1393 const number omega = 1.) const;
1394
1400 template <typename somenumber>
1401 void
1402 SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1403
1408 template <typename somenumber>
1409 void
1410 SOR(Vector<somenumber> &v, const number omega = 1.) const;
1411
1416 template <typename somenumber>
1417 void
1418 TSOR(Vector<somenumber> &v, const number omega = 1.) const;
1419
1430 template <typename somenumber>
1431 void
1433 const std::vector<size_type> &permutation,
1434 const std::vector<size_type> &inverse_permutation,
1435 const number omega = 1.) const;
1436
1447 template <typename somenumber>
1448 void
1450 const std::vector<size_type> &permutation,
1451 const std::vector<size_type> &inverse_permutation,
1452 const number omega = 1.) const;
1453
1459 template <typename somenumber>
1460 void
1462 const Vector<somenumber> &b,
1463 const number omega = 1.) const;
1464
1469 template <typename somenumber>
1470 void
1472 const Vector<somenumber> &b,
1473 const number omega = 1.) const;
1474
1479 template <typename somenumber>
1480 void
1482 const Vector<somenumber> &b,
1483 const number omega = 1.) const;
1484
1489 template <typename somenumber>
1490 void
1492 const Vector<somenumber> &b,
1493 const number omega = 1.) const;
1507 begin() const;
1508
1512 iterator
1514
1519 end() const;
1520
1524 iterator
1526
1537 begin(const size_type r) const;
1538
1542 iterator
1544
1555 end(const size_type r) const;
1556
1560 iterator
1561 end(const size_type r);
1579 template <typename StreamType>
1580 void
1581 print(StreamType &out,
1582 const bool across = false,
1583 const bool diagonal_first = true) const;
1584
1607 void
1608 print_formatted(std::ostream &out,
1609 const unsigned int precision = 3,
1610 const bool scientific = true,
1611 const unsigned int width = 0,
1612 const char *zero_string = " ",
1613 const double denominator = 1.,
1614 const char *separator = " ") const;
1615
1621 void
1622 print_pattern(std::ostream &out, const double threshold = 0.) const;
1623
1632 void
1633 print_as_numpy_arrays(std::ostream &out,
1634 const unsigned int precision = 9) const;
1635
1646 void
1647 block_write(std::ostream &out) const;
1648
1665 void
1666 block_read(std::istream &in);
1677 int,
1678 int,
1679 << "You are trying to access the matrix entry with index <"
1680 << arg1 << ',' << arg2
1681 << ">, but this entry does not exist in the sparsity pattern "
1682 "of this matrix."
1683 "\n\n"
1684 "The most common cause for this problem is that you used "
1685 "a method to build the sparsity pattern that did not "
1686 "(completely) take into account all of the entries you "
1687 "will later try to write into. An example would be "
1688 "building a sparsity pattern that does not include "
1689 "the entries you will write into due to constraints "
1690 "on degrees of freedom such as hanging nodes or periodic "
1691 "boundary conditions. In such cases, building the "
1692 "sparsity pattern will succeed, but you will get errors "
1693 "such as the current one at one point or other when "
1694 "trying to write into the entries of the matrix.");
1699 "When copying one sparse matrix into another, "
1700 "or when adding one sparse matrix to another, "
1701 "both matrices need to refer to the same "
1702 "sparsity pattern.");
1707 int,
1708 int,
1709 << "The iterators denote a range of " << arg1
1710 << " elements, but the given number of rows was " << arg2);
1715 "You are attempting an operation on two vectors that "
1716 "are the same object, but the operation requires that the "
1717 "two objects are in fact different.");
1720protected:
1731 void
1733
1738 void
1740
1741private:
1748
1756 std::unique_ptr<number[]> val;
1757
1764 std::size_t max_len;
1765
1766 // make all other sparse matrices friends
1767 template <typename somenumber>
1768 friend class SparseMatrix;
1769 template <typename somenumber>
1771 template <typename>
1772 friend class SparseILU;
1773
1774 // To allow it calling private prepare_add() and prepare_set().
1775 template <typename>
1776 friend class BlockMatrixBase;
1777
1778 // Also give access to internal details to the iterator/accessor classes.
1779 template <typename, bool>
1781 template <typename, bool>
1783
1784#ifdef DEAL_II_WITH_MPI
1785 // Give access to internal datastructures to perform MPI operations.
1786 template <typename Number>
1787 friend void
1789 const MPI_Comm,
1791#endif
1792};
1793
1794#ifndef DOXYGEN
1795/*---------------------- Inline functions -----------------------------------*/
1796
1797
1798
1799template <typename number>
1800template <typename number2>
1801void
1803{
1804 this->reinit(sparse_matrix.get_sparsity_pattern());
1805}
1806
1807
1808
1809template <typename number>
1810inline typename SparseMatrix<number>::size_type
1812{
1813 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1814 return cols->rows;
1815}
1816
1817
1818template <typename number>
1819inline typename SparseMatrix<number>::size_type
1821{
1822 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1823 return cols->cols;
1824}
1825
1826
1827// Inline the set() and add() functions, since they will be called frequently.
1828template <typename number>
1829inline void
1831 const size_type j,
1832 const number value)
1833{
1835
1836 const size_type index = cols->operator()(i, j);
1837
1838 // it is allowed to set elements of the matrix that are not part of the
1839 // sparsity pattern, if the value to which we set it is zero
1841 {
1842 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1843 ExcInvalidIndex(i, j));
1844 return;
1845 }
1846
1847 val[index] = value;
1848}
1849
1850
1851
1852template <typename number>
1853template <typename number2>
1854inline void
1855SparseMatrix<number>::set(const std::vector<size_type> &indices,
1856 const FullMatrix<number2> &values,
1857 const bool elide_zero_values)
1858{
1859 Assert(indices.size() == values.m(),
1860 ExcDimensionMismatch(indices.size(), values.m()));
1861 Assert(values.m() == values.n(), ExcNotQuadratic());
1862
1863 for (size_type i = 0; i < indices.size(); ++i)
1864 set(indices[i],
1865 indices.size(),
1866 indices.data(),
1867 &values(i, 0),
1868 elide_zero_values);
1869}
1870
1871
1872
1873template <typename number>
1874template <typename number2>
1875inline void
1876SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1877 const std::vector<size_type> &col_indices,
1878 const FullMatrix<number2> &values,
1879 const bool elide_zero_values)
1880{
1881 Assert(row_indices.size() == values.m(),
1882 ExcDimensionMismatch(row_indices.size(), values.m()));
1883 Assert(col_indices.size() == values.n(),
1884 ExcDimensionMismatch(col_indices.size(), values.n()));
1885
1886 for (size_type i = 0; i < row_indices.size(); ++i)
1887 set(row_indices[i],
1888 col_indices.size(),
1889 col_indices.data(),
1890 &values(i, 0),
1891 elide_zero_values);
1892}
1893
1894
1895
1896template <typename number>
1897template <typename number2>
1898inline void
1900 const std::vector<size_type> &col_indices,
1901 const std::vector<number2> &values,
1902 const bool elide_zero_values)
1903{
1904 Assert(col_indices.size() == values.size(),
1905 ExcDimensionMismatch(col_indices.size(), values.size()));
1906
1907 set(row,
1908 col_indices.size(),
1909 col_indices.data(),
1910 values.data(),
1911 elide_zero_values);
1912}
1913
1914
1915
1916template <typename number>
1917inline void
1919 const size_type j,
1920 const number value)
1921{
1923
1924 if (value == number())
1925 return;
1926
1927 const size_type index = cols->operator()(i, j);
1928
1929 // it is allowed to add elements to the matrix that are not part of the
1930 // sparsity pattern, if the value to which we set it is zero
1932 {
1933 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1934 ExcInvalidIndex(i, j));
1935 return;
1936 }
1937
1938 val[index] += value;
1939}
1940
1941
1942
1943template <typename number>
1944template <typename number2>
1945inline void
1946SparseMatrix<number>::add(const std::vector<size_type> &indices,
1947 const FullMatrix<number2> &values,
1948 const bool elide_zero_values)
1949{
1950 Assert(indices.size() == values.m(),
1951 ExcDimensionMismatch(indices.size(), values.m()));
1952 Assert(values.m() == values.n(), ExcNotQuadratic());
1953
1954 for (size_type i = 0; i < indices.size(); ++i)
1955 add(indices[i],
1956 indices.size(),
1957 indices.data(),
1958 &values(i, 0),
1959 elide_zero_values);
1960}
1961
1962
1963
1964template <typename number>
1965template <typename number2>
1966inline void
1967SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1968 const std::vector<size_type> &col_indices,
1969 const FullMatrix<number2> &values,
1970 const bool elide_zero_values)
1971{
1972 Assert(row_indices.size() == values.m(),
1973 ExcDimensionMismatch(row_indices.size(), values.m()));
1974 Assert(col_indices.size() == values.n(),
1975 ExcDimensionMismatch(col_indices.size(), values.n()));
1976
1977 for (size_type i = 0; i < row_indices.size(); ++i)
1978 add(row_indices[i],
1979 col_indices.size(),
1980 col_indices.data(),
1981 &values(i, 0),
1982 elide_zero_values);
1983}
1984
1985
1986
1987template <typename number>
1988template <typename number2>
1989inline void
1991 const std::vector<size_type> &col_indices,
1992 const std::vector<number2> &values,
1993 const bool elide_zero_values)
1994{
1995 Assert(col_indices.size() == values.size(),
1996 ExcDimensionMismatch(col_indices.size(), values.size()));
1997
1998 add(row,
1999 col_indices.size(),
2000 col_indices.data(),
2001 values.data(),
2002 elide_zero_values,
2003 std::is_sorted(col_indices.begin(), col_indices.end()));
2004}
2005
2006
2007
2008template <typename number>
2009inline SparseMatrix<number> &
2010SparseMatrix<number>::operator*=(const number factor)
2011{
2012 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2013 Assert(val != nullptr, ExcNotInitialized());
2014
2015 number *val_ptr = val.get();
2016 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2017
2018 while (val_ptr != end_ptr)
2019 *val_ptr++ *= factor;
2020
2021 return *this;
2022}
2023
2024
2025
2026template <typename number>
2027inline SparseMatrix<number> &
2028SparseMatrix<number>::operator/=(const number factor)
2029{
2030 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2031 Assert(val != nullptr, ExcNotInitialized());
2032 Assert(factor != number(), ExcDivideByZero());
2033
2034 const number factor_inv = number(1.) / factor;
2035
2036 number *val_ptr = val.get();
2037 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2038
2039 while (val_ptr != end_ptr)
2040 *val_ptr++ *= factor_inv;
2041
2042 return *this;
2043}
2044
2045
2046
2047template <typename number>
2048inline const number &
2050{
2051 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2052 Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2053 ExcInvalidIndex(i, j));
2054 return val[cols->operator()(i, j)];
2055}
2056
2057
2058
2059template <typename number>
2060inline number &
2062{
2063 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2064 Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2065 ExcInvalidIndex(i, j));
2066 return val[cols->operator()(i, j)];
2067}
2068
2069
2070
2071template <typename number>
2072inline number
2073SparseMatrix<number>::el(const size_type i, const size_type j) const
2074{
2075 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2076 const size_type index = cols->operator()(i, j);
2077
2079 return val[index];
2080 else
2081 return 0;
2082}
2083
2084
2085
2086template <typename number>
2087inline number
2089{
2090 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2091 Assert(m() == n(), ExcNotQuadratic());
2092 AssertIndexRange(i, m());
2093
2094 // Use that the first element in each row of a quadratic matrix is the main
2095 // diagonal
2096 return val[cols->rowstart[i]];
2097}
2098
2099
2100
2101template <typename number>
2102inline number &
2104{
2105 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2106 Assert(m() == n(), ExcNotQuadratic());
2107 AssertIndexRange(i, m());
2108
2109 // Use that the first element in each row of a quadratic matrix is the main
2110 // diagonal
2111 return val[cols->rowstart[i]];
2112}
2113
2114
2115
2116template <typename number>
2117template <typename ForwardIterator>
2118void
2119SparseMatrix<number>::copy_from(const ForwardIterator begin,
2120 const ForwardIterator end)
2121{
2122 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2123 ExcIteratorRange(std::distance(begin, end), m()));
2124
2125 // for use in the inner loop, we define an alias to the type of the inner
2126 // iterators
2127 using inner_iterator =
2128 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2129 size_type row = 0;
2130 for (ForwardIterator i = begin; i != end; ++i, ++row)
2131 {
2132 const inner_iterator end_of_row = i->end();
2133 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2134 // write entries
2135 set(row, j->first, j->second);
2136 };
2137}
2138
2139
2140//---------------------------------------------------------------------------
2141
2142
2143namespace SparseMatrixIterators
2144{
2145 template <typename number>
2146 inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2147 const std::size_t index_within_matrix)
2148 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2149 index_within_matrix)
2150 , matrix(matrix)
2151 {}
2152
2153
2154
2155 template <typename number>
2156 inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2157 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2158 , matrix(matrix)
2159 {}
2160
2161
2162
2163 template <typename number>
2164 inline Accessor<number, true>::Accessor(
2166 : SparsityPatternIterators::Accessor(a)
2167 , matrix(&a.get_matrix())
2168 {}
2169
2170
2171
2172 template <typename number>
2173 inline number
2174 Accessor<number, true>::value() const
2175 {
2176 AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2177 return matrix->val[linear_index];
2178 }
2179
2180
2181
2182 template <typename number>
2183 inline const typename Accessor<number, true>::MatrixType &
2184 Accessor<number, true>::get_matrix() const
2185 {
2186 return *matrix;
2187 }
2188
2189
2190
2191 template <typename number>
2192 inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2193 const bool)
2194 : accessor(accessor)
2195 {}
2196
2197
2198 template <typename number>
2199 inline Accessor<number, false>::Reference::operator number() const
2200 {
2201 AssertIndexRange(accessor->linear_index,
2202 accessor->matrix->n_nonzero_elements());
2203 return accessor->matrix->val[accessor->linear_index];
2204 }
2205
2206
2207
2208 template <typename number>
2209 inline const typename Accessor<number, false>::Reference &
2210 Accessor<number, false>::Reference::operator=(const number n) const
2211 {
2212 AssertIndexRange(accessor->linear_index,
2213 accessor->matrix->n_nonzero_elements());
2214 accessor->matrix->val[accessor->linear_index] = n;
2215 return *this;
2216 }
2217
2218
2219
2220 template <typename number>
2221 inline const typename Accessor<number, false>::Reference &
2222 Accessor<number, false>::Reference::operator+=(const number n) const
2223 {
2224 AssertIndexRange(accessor->linear_index,
2225 accessor->matrix->n_nonzero_elements());
2226 accessor->matrix->val[accessor->linear_index] += n;
2227 return *this;
2228 }
2229
2230
2231
2232 template <typename number>
2233 inline const typename Accessor<number, false>::Reference &
2234 Accessor<number, false>::Reference::operator-=(const number n) const
2235 {
2236 AssertIndexRange(accessor->linear_index,
2237 accessor->matrix->n_nonzero_elements());
2238 accessor->matrix->val[accessor->linear_index] -= n;
2239 return *this;
2240 }
2241
2242
2243
2244 template <typename number>
2245 inline const typename Accessor<number, false>::Reference &
2246 Accessor<number, false>::Reference::operator*=(const number n) const
2247 {
2248 AssertIndexRange(accessor->linear_index,
2249 accessor->matrix->n_nonzero_elements());
2250 accessor->matrix->val[accessor->linear_index] *= n;
2251 return *this;
2252 }
2253
2254
2255
2256 template <typename number>
2257 inline const typename Accessor<number, false>::Reference &
2258 Accessor<number, false>::Reference::operator/=(const number n) const
2259 {
2260 AssertIndexRange(accessor->linear_index,
2261 accessor->matrix->n_nonzero_elements());
2262 accessor->matrix->val[accessor->linear_index] /= n;
2263 return *this;
2264 }
2265
2266
2267
2268 template <typename number>
2269 inline Accessor<number, false>::Accessor(MatrixType *matrix,
2270 const std::size_t index)
2271 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2272 , matrix(matrix)
2273 {}
2274
2275
2276
2277 template <typename number>
2278 inline Accessor<number, false>::Accessor(MatrixType *matrix)
2279 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2280 , matrix(matrix)
2281 {}
2282
2283
2284
2285 template <typename number>
2286 inline typename Accessor<number, false>::Reference
2287 Accessor<number, false>::value() const
2288 {
2289 return Reference(this, true);
2290 }
2291
2292
2293
2294 template <typename number>
2295 inline typename Accessor<number, false>::MatrixType &
2296 Accessor<number, false>::get_matrix() const
2297 {
2298 return *matrix;
2299 }
2300
2301
2302
2303 template <typename number, bool Constness>
2304 inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
2305 const std::size_t index)
2306 : accessor(matrix, index)
2307 {}
2308
2309
2310
2311 template <typename number, bool Constness>
2312 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2313 : accessor(matrix)
2314 {}
2315
2316
2317
2318 template <typename number, bool Constness>
2319 inline Iterator<number, Constness>::Iterator(
2321 : accessor(*i)
2322 {}
2323
2324
2325
2326 template <typename number, bool Constness>
2327 inline const Iterator<number, Constness> &
2328 Iterator<number, Constness>::operator=(
2330 {
2331 accessor = *i;
2332 return *this;
2333 }
2334
2335
2336
2337 template <typename number, bool Constness>
2338 inline Iterator<number, Constness> &
2339 Iterator<number, Constness>::operator++()
2340 {
2341 accessor.advance();
2342 return *this;
2343 }
2344
2345
2346 template <typename number, bool Constness>
2347 inline Iterator<number, Constness>
2348 Iterator<number, Constness>::operator++(int)
2349 {
2350 const Iterator iter = *this;
2351 accessor.advance();
2352 return iter;
2353 }
2354
2355
2356 template <typename number, bool Constness>
2357 inline const Accessor<number, Constness> &
2358 Iterator<number, Constness>::operator*() const
2359 {
2360 return accessor;
2361 }
2362
2363
2364 template <typename number, bool Constness>
2365 inline const Accessor<number, Constness> *
2366 Iterator<number, Constness>::operator->() const
2367 {
2368 return &accessor;
2369 }
2370
2371
2372 template <typename number, bool Constness>
2373 inline bool
2374 Iterator<number, Constness>::operator==(const Iterator &other) const
2375 {
2376 return (accessor == other.accessor);
2377 }
2378
2379
2380 template <typename number, bool Constness>
2381 inline bool
2382 Iterator<number, Constness>::operator!=(const Iterator &other) const
2383 {
2384 return !(*this == other);
2385 }
2386
2387
2388 template <typename number, bool Constness>
2389 inline bool
2390 Iterator<number, Constness>::operator<(const Iterator &other) const
2391 {
2392 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2394
2395 return (accessor < other.accessor);
2396 }
2397
2398
2399 template <typename number, bool Constness>
2400 inline bool
2401 Iterator<number, Constness>::operator>(const Iterator &other) const
2402 {
2403 return (other < *this);
2404 }
2405
2406
2407 template <typename number, bool Constness>
2408 inline int
2409 Iterator<number, Constness>::operator-(const Iterator &other) const
2410 {
2411 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2413
2414 return (*this)->linear_index - other->linear_index;
2415 }
2416
2417
2418
2419 template <typename number, bool Constness>
2420 inline Iterator<number, Constness>
2421 Iterator<number, Constness>::operator+(const size_type n) const
2422 {
2423 Iterator x = *this;
2424 for (size_type i = 0; i < n; ++i)
2425 ++x;
2426
2427 return x;
2428 }
2429
2430} // namespace SparseMatrixIterators
2431
2432
2433
2434template <typename number>
2437{
2438 return const_iterator(this, 0);
2439}
2440
2441
2442template <typename number>
2445{
2446 return const_iterator(this);
2447}
2448
2449
2450template <typename number>
2451inline typename SparseMatrix<number>::iterator
2453{
2454 return iterator(this, 0);
2455}
2456
2457
2458template <typename number>
2459inline typename SparseMatrix<number>::iterator
2461{
2462 return iterator(this, cols->rowstart[cols->rows]);
2463}
2464
2465
2466template <typename number>
2468SparseMatrix<number>::begin(const size_type r) const
2469{
2470 AssertIndexRange(r, m());
2471
2472 return const_iterator(this, cols->rowstart[r]);
2473}
2474
2475
2476
2477template <typename number>
2479SparseMatrix<number>::end(const size_type r) const
2480{
2481 AssertIndexRange(r, m());
2482
2483 return const_iterator(this, cols->rowstart[r + 1]);
2484}
2485
2486
2487
2488template <typename number>
2489inline typename SparseMatrix<number>::iterator
2490SparseMatrix<number>::begin(const size_type r)
2491{
2492 AssertIndexRange(r, m());
2493
2494 return iterator(this, cols->rowstart[r]);
2495}
2496
2497
2498
2499template <typename number>
2500inline typename SparseMatrix<number>::iterator
2501SparseMatrix<number>::end(const size_type r)
2502{
2503 AssertIndexRange(r, m());
2504
2505 return iterator(this, cols->rowstart[r + 1]);
2506}
2507
2508
2509
2510template <typename number>
2511template <typename StreamType>
2512inline void
2513SparseMatrix<number>::print(StreamType &out,
2514 const bool across,
2515 const bool diagonal_first) const
2516{
2517 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2518 Assert(val != nullptr, ExcNotInitialized());
2519
2520 bool hanging_diagonal = false;
2521 number diagonal = number();
2522
2523 for (size_type i = 0; i < cols->rows; ++i)
2524 {
2525 for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2526 {
2527 if (!diagonal_first && i == cols->colnums[j])
2528 {
2529 diagonal = val[j];
2530 hanging_diagonal = true;
2531 }
2532 else
2533 {
2534 if (hanging_diagonal && cols->colnums[j] > i)
2535 {
2536 if (across)
2537 out << ' ' << i << ',' << i << ':' << diagonal;
2538 else
2539 out << '(' << i << ',' << i << ") " << diagonal
2540 << std::endl;
2541 hanging_diagonal = false;
2542 }
2543 if (across)
2544 out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2545 else
2546 out << '(' << i << ',' << cols->colnums[j] << ") " << val[j]
2547 << std::endl;
2548 }
2549 }
2550 if (hanging_diagonal)
2551 {
2552 if (across)
2553 out << ' ' << i << ',' << i << ':' << diagonal;
2554 else
2555 out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2556 hanging_diagonal = false;
2557 }
2558 }
2559 if (across)
2560 out << std::endl;
2561}
2562
2563
2564template <typename number>
2565inline void
2567{
2568 // nothing to do here
2569}
2570
2571
2572
2573template <typename number>
2574inline void
2576{
2577 // nothing to do here
2578}
2579
2580#endif // DOXYGEN
2581
2582
2584
2585#endif
const Reference & operator-=(const number n) const
Reference(const Accessor *accessor, const bool dummy)
const Reference & operator/=(const number n) const
const Reference & operator+=(const number n) const
const Reference & operator*=(const number n) const
const Reference & operator=(const number n) const
Accessor(MatrixType *matrix, const std::size_t index)
Accessor(const SparseMatrixIterators::Accessor< number, false > &a)
Accessor(MatrixType *matrix, const std::size_t index_within_matrix)
const SparseMatrix< number > & get_matrix() const
const Iterator< number, Constness > & operator=(const SparseMatrixIterators::Iterator< number, false > &i)
bool operator>(const Iterator &) const
bool operator==(const Iterator &) const
int operator-(const Iterator &p) const
bool operator<(const Iterator &) const
Iterator(MatrixType *matrix)
const Accessor< number, Constness > & value_type
Iterator operator+(const size_type n) const
const Accessor< number, Constness > & operator*() const
Accessor< number, Constness > accessor
Iterator(const SparseMatrixIterators::Iterator< number, false > &i)
Iterator(MatrixType *matrix, const std::size_t index_within_matrix)
const Accessor< number, Constness > * operator->() const
typename Accessor< number, Constness >::MatrixType MatrixType
bool operator!=(const Iterator &) const
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number2 > &values, const bool elide_zero_values=false)
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
SparseMatrix(const SparsityPattern &sparsity)
void add(const number factor, const SparseMatrix< somenumber > &matrix)
std::size_t n_nonzero_elements() const
size_type get_row_length(const size_type row) const
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
iterator end()
number & diag_element(const size_type i)
void prepare_add()
void Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
void Tvmult(OutVector &dst, const InVector &src) const
const_iterator begin(const size_type r) const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
const_iterator end() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
void set(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
void print_as_numpy_arrays(std::ostream &out, const unsigned int precision=9) const
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
ObserverPointer< const SparsityPattern, SparseMatrix< number > > cols
const SparsityPattern & get_sparsity_pattern() const
void set(const size_type i, const size_type j, const number value)
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
const_iterator begin() const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
virtual void clear()
virtual ~SparseMatrix() override
void vmult_add(OutVector &dst, const InVector &src) const
number diag_element(const size_type i) const
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=true)
void SSOR(Vector< somenumber > &v, const number omega=1.) const
somenumber matrix_norm_square(const Vector< somenumber > &v) const
SparseMatrix(SparseMatrix< number > &&m) noexcept
void TSOR(Vector< somenumber > &v, const number omega=1.) const
SparseMatrix< number > & copy_from(const TrilinosWrappers::SparseMatrix &matrix)
void symmetrize()
SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id)
void print_pattern(std::ostream &out, const double threshold=0.) const
void vmult(OutVector &dst, const InVector &src) const
std::size_t memory_consumption() const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number omega=1.) const
void block_write(std::ostream &out) const
void copy_from(const ForwardIterator begin, const ForwardIterator end)
const number & operator()(const size_type i, const size_type j) const
number & operator()(const size_type i, const size_type j)
void block_read(std::istream &in)
void prepare_set()
number el(const size_type i, const size_type j) const
SparseMatrix< number > & operator=(SparseMatrix< number > &&m) noexcept
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
const_iterator end(const size_type r) const
iterator begin()
void SOR(Vector< somenumber > &v, const number omega=1.) const
SparseMatrix & operator/=(const number factor)
std::unique_ptr< number[]> val
iterator begin(const size_type r)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
typename numbers::NumberTraits< number >::real_type real_type
void reinit(const SparseMatrix< number2 > &sparse_matrix)
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)
void add(const size_type i, const size_type j, const number value)
size_type n() const
size_type m() const
void copy_from(const FullMatrix< somenumber > &matrix)
std::size_t max_len
iterator end(const size_type r)
SparseMatrix & operator*=(const number factor)
void compress(VectorOperation::values)
SparseMatrix(const SparseMatrix &)
SparseMatrix< number > & operator=(const IdentityMatrix &id)
SparseMatrix & operator=(const double d)
real_type frobenius_norm() const
void Tvmult_add(OutVector &dst, const InVector &src) const
real_type linfty_norm() const
void add(const size_type row, const std::vector< size_type > &col_indices, const std::vector< number2 > &values, const bool elide_zero_values=true)
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
real_type l1_norm() const
std::size_t n_actually_nonzero_elements(const double threshold=0.) const
bool empty() const
void add(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=true)
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number omega=1.) const
virtual void reinit(const SparsityPattern &sparsity)
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=false)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number omega=1.) const
friend class ChunkSparsityPatternIterators::Accessor
SparsityPatternIterators::size_type size_type
static constexpr size_type invalid_entry
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcDifferentSparsityPatterns()
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcSourceEqualsDestination()
#define AssertIsFinite(number)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:534
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:489
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDivideByZero()
static ::ExceptionBase & ExcIteratorRange(int arg1, int arg2)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
STL namespace.
unsigned int global_dof_index
Definition types.h:81
static const bool zero_addition_can_be_elided
typename ::SparseMatrixIterators::Iterator< number, Constness >::difference_type difference_type
typename ::SparseMatrixIterators::Iterator< number, Constness >::value_type value_type