Reference documentation for deal.II version 9.3.3
\(\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\}}\)
sparse_matrix.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2021 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_sparse_matrix_h
17# define dealii_sparse_matrix_h
18
19
20# include <deal.II/base/config.h>
21
24
29# ifdef DEAL_II_WITH_MPI
30# include <mpi.h>
31# endif
32
33# include <memory>
34
35
37
38// Forward declarations
39# ifndef DOXYGEN
40template <typename number>
41class Vector;
42template <typename number>
43class FullMatrix;
44template <typename Matrix>
45class BlockMatrixBase;
46template <typename number>
47class SparseILU;
48# ifdef DEAL_II_WITH_MPI
49namespace Utilities
50{
51 namespace MPI
52 {
53 template <typename Number>
54 void
56 }
57} // namespace Utilities
58# endif
59
60# ifdef DEAL_II_WITH_TRILINOS
61namespace TrilinosWrappers
62{
63 class SparseMatrix;
64}
65# endif
66# endif
67
78{
83
84 // forward declaration
85 template <typename number, bool Constness>
86 class Iterator;
87
98 template <typename number, bool Constness>
100 {
101 public:
105 number
106 value() const;
107
111 number &
113
119 get_matrix() const;
120 };
121
122
123
130 template <typename number>
131 class Accessor<number, true> : public SparsityPatternIterators::Accessor
132 {
133 public:
139
143 Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
144
149
154
158 number
159 value() const;
160
165 const MatrixType &
166 get_matrix() const;
167
168 private:
173
178
179 // Make iterator class a friend.
180 template <typename, bool>
181 friend class Iterator;
182 };
183
184
191 template <typename number>
192 class Accessor<number, false> : public SparsityPatternIterators::Accessor
193 {
194 private:
219 class Reference
220 {
221 public:
226 Reference(const Accessor *accessor, const bool dummy);
227
231 operator number() const;
232
236 const Reference &
237 operator=(const number n) const;
238
242 const Reference &
243 operator+=(const number n) const;
244
248 const Reference &
249 operator-=(const number n) const;
250
254 const Reference &
255 operator*=(const number n) const;
256
260 const Reference &
261 operator/=(const number n) const;
262
263 private:
269 };
270
271 public:
277
281 Accessor(MatrixType *matrix, const std::size_t index);
282
287
291 Reference
292 value() const;
293
298 MatrixType &
299 get_matrix() const;
300
301 private:
306
311
312 // Make iterator class a friend.
313 template <typename, bool>
314 friend class Iterator;
315 };
316
317
318
348 template <typename number, bool Constness>
350 {
351 public:
356
362
367 Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
368
373
379
385
389 Iterator &
391
397
402
407
411 bool
412 operator==(const Iterator &) const;
413
417 bool
418 operator!=(const Iterator &) const;
419
427 bool
428 operator<(const Iterator &) const;
429
434 bool
435 operator>(const Iterator &) const;
436
443 int
444 operator-(const Iterator &p) const;
445
450 operator+(const size_type n) const;
451
452 private:
457 };
458
459} // namespace SparseMatrixIterators
460
466// TODO: Add multithreading to the other vmult functions.
467
494template <typename number>
495class SparseMatrix : public virtual Subscriptor
496{
497public:
502
507 using value_type = number;
508
519
525
533
540 struct Traits
541 {
546 static const bool zero_addition_can_be_elided = true;
547 };
548
564
574
583
597 explicit SparseMatrix(const SparsityPattern &sparsity);
598
605 SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
606
611 virtual ~SparseMatrix() override;
612
624
631
640
653 operator=(const double d);
654
668 virtual void
669 reinit(const SparsityPattern &sparsity);
670
676 virtual void
679
687 bool
688 empty() const;
689
695 m() const;
696
702 n() const;
703
708 get_row_length(const size_type row) const;
709
715 std::size_t
717
727 std::size_t
728 n_actually_nonzero_elements(const double threshold = 0.) const;
729
738 const SparsityPattern &
740
745 std::size_t
747
752
754
763 void
764 set(const size_type i, const size_type j, const number value);
765
781 template <typename number2>
782 void
783 set(const std::vector<size_type> &indices,
784 const FullMatrix<number2> & full_matrix,
785 const bool elide_zero_values = false);
786
792 template <typename number2>
793 void
794 set(const std::vector<size_type> &row_indices,
795 const std::vector<size_type> &col_indices,
796 const FullMatrix<number2> & full_matrix,
797 const bool elide_zero_values = false);
798
809 template <typename number2>
810 void
811 set(const size_type row,
812 const std::vector<size_type> &col_indices,
813 const std::vector<number2> & values,
814 const bool elide_zero_values = false);
815
825 template <typename number2>
826 void
827 set(const size_type row,
828 const size_type n_cols,
829 const size_type *col_indices,
830 const number2 * values,
831 const bool elide_zero_values = false);
832
838 void
839 add(const size_type i, const size_type j, const number value);
840
855 template <typename number2>
856 void
857 add(const std::vector<size_type> &indices,
858 const FullMatrix<number2> & full_matrix,
859 const bool elide_zero_values = true);
860
866 template <typename number2>
867 void
868 add(const std::vector<size_type> &row_indices,
869 const std::vector<size_type> &col_indices,
870 const FullMatrix<number2> & full_matrix,
871 const bool elide_zero_values = true);
872
882 template <typename number2>
883 void
884 add(const size_type row,
885 const std::vector<size_type> &col_indices,
886 const std::vector<number2> & values,
887 const bool elide_zero_values = true);
888
898 template <typename number2>
899 void
900 add(const size_type row,
901 const size_type n_cols,
902 const size_type *col_indices,
903 const number2 * values,
904 const bool elide_zero_values = true,
905 const bool col_indices_are_sorted = false);
906
911 operator*=(const number factor);
912
917 operator/=(const number factor);
918
931 void
933
950 template <typename somenumber>
953
970 template <typename ForwardIterator>
971 void
972 copy_from(const ForwardIterator begin, const ForwardIterator end);
973
983 template <typename somenumber>
984 void
986
987# ifdef DEAL_II_WITH_TRILINOS
999# endif
1000
1012 template <typename somenumber>
1013 void
1014 add(const number factor, const SparseMatrix<somenumber> &matrix);
1015
1017
1021
1035 const number &
1036 operator()(const size_type i, const size_type j) const;
1037
1041 number &
1042 operator()(const size_type i, const size_type j);
1043
1056 number
1057 el(const size_type i, const size_type j) const;
1058
1068 number
1069 diag_element(const size_type i) const;
1070
1075 number &
1077
1079
1099 template <class OutVector, class InVector>
1100 void
1101 vmult(OutVector &dst, const InVector &src) const;
1102
1118 template <class OutVector, class InVector>
1119 void
1120 Tvmult(OutVector &dst, const InVector &src) const;
1121
1138 template <class OutVector, class InVector>
1139 void
1140 vmult_add(OutVector &dst, const InVector &src) const;
1141
1157 template <class OutVector, class InVector>
1158 void
1159 Tvmult_add(OutVector &dst, const InVector &src) const;
1160
1178 template <typename somenumber>
1179 somenumber
1181
1187 template <typename somenumber>
1188 somenumber
1190 const Vector<somenumber> &v) const;
1191
1201 template <typename somenumber>
1202 somenumber
1204 const Vector<somenumber> &x,
1205 const Vector<somenumber> &b) const;
1206
1242 template <typename numberB, typename numberC>
1243 void
1245 const SparseMatrix<numberB> &B,
1246 const Vector<number> & V = Vector<number>(),
1247 const bool rebuild_sparsity_pattern = true) const;
1248
1273 template <typename numberB, typename numberC>
1274 void
1276 const SparseMatrix<numberB> &B,
1277 const Vector<number> & V = Vector<number>(),
1278 const bool rebuild_sparsity_pattern = true) const;
1279
1281
1285
1293 real_type
1294 l1_norm() const;
1295
1303 real_type
1305
1310 real_type
1313
1317
1323 template <typename somenumber>
1324 void
1326 const Vector<somenumber> &src,
1327 const number omega = 1.) const;
1328
1335 template <typename somenumber>
1336 void
1338 const Vector<somenumber> & src,
1339 const number omega = 1.,
1340 const std::vector<std::size_t> &pos_right_of_diagonal =
1341 std::vector<std::size_t>()) const;
1342
1346 template <typename somenumber>
1347 void
1349 const Vector<somenumber> &src,
1350 const number om = 1.) const;
1351
1355 template <typename somenumber>
1356 void
1358 const Vector<somenumber> &src,
1359 const number om = 1.) const;
1360
1366 template <typename somenumber>
1367 void
1368 SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1369
1374 template <typename somenumber>
1375 void
1376 SOR(Vector<somenumber> &v, const number om = 1.) const;
1377
1382 template <typename somenumber>
1383 void
1384 TSOR(Vector<somenumber> &v, const number om = 1.) const;
1385
1396 template <typename somenumber>
1397 void
1399 const std::vector<size_type> &permutation,
1400 const std::vector<size_type> &inverse_permutation,
1401 const number om = 1.) const;
1402
1413 template <typename somenumber>
1414 void
1416 const std::vector<size_type> &permutation,
1417 const std::vector<size_type> &inverse_permutation,
1418 const number om = 1.) const;
1419
1425 template <typename somenumber>
1426 void
1428 const Vector<somenumber> &b,
1429 const number om = 1.) const;
1430
1435 template <typename somenumber>
1436 void
1438 const Vector<somenumber> &b,
1439 const number om = 1.) const;
1440
1445 template <typename somenumber>
1446 void
1448 const Vector<somenumber> &b,
1449 const number om = 1.) const;
1450
1455 template <typename somenumber>
1456 void
1458 const Vector<somenumber> &b,
1459 const number om = 1.) const;
1461
1465
1473 begin() const;
1474
1478 iterator
1480
1485 end() const;
1486
1490 iterator
1492
1503 begin(const size_type r) const;
1504
1508 iterator
1510
1521 end(const size_type r) const;
1522
1526 iterator
1527 end(const size_type r);
1529
1533
1545 template <class StreamType>
1546 void
1547 print(StreamType &out,
1548 const bool across = false,
1549 const bool diagonal_first = true) const;
1550
1571 void
1572 print_formatted(std::ostream & out,
1573 const unsigned int precision = 3,
1574 const bool scientific = true,
1575 const unsigned int width = 0,
1576 const char * zero_string = " ",
1577 const double denominator = 1.) const;
1578
1584 void
1585 print_pattern(std::ostream &out, const double threshold = 0.) const;
1586
1595 void
1596 print_as_numpy_arrays(std::ostream & out,
1597 const unsigned int precision = 9) const;
1598
1609 void
1610 block_write(std::ostream &out) const;
1611
1628 void
1629 block_read(std::istream &in);
1631
1640 int,
1641 int,
1642 << "You are trying to access the matrix entry with index <"
1643 << arg1 << ',' << arg2
1644 << ">, but this entry does not exist in the sparsity pattern "
1645 "of this matrix."
1646 "\n\n"
1647 "The most common cause for this problem is that you used "
1648 "a method to build the sparsity pattern that did not "
1649 "(completely) take into account all of the entries you "
1650 "will later try to write into. An example would be "
1651 "building a sparsity pattern that does not include "
1652 "the entries you will write into due to constraints "
1653 "on degrees of freedom such as hanging nodes or periodic "
1654 "boundary conditions. In such cases, building the "
1655 "sparsity pattern will succeed, but you will get errors "
1656 "such as the current one at one point or other when "
1657 "trying to write into the entries of the matrix.");
1662 "When copying one sparse matrix into another, "
1663 "or when adding one sparse matrix to another, "
1664 "both matrices need to refer to the same "
1665 "sparsity pattern.");
1670 int,
1671 int,
1672 << "The iterators denote a range of " << arg1
1673 << " elements, but the given number of rows was " << arg2);
1678 "You are attempting an operation on two matrices that "
1679 "are the same object, but the operation requires that the "
1680 "two objects are in fact different.");
1682
1683protected:
1694 void
1696
1701 void
1703
1704private:
1711
1719 std::unique_ptr<number[]> val;
1720
1727 std::size_t max_len;
1728
1729 // make all other sparse matrices friends
1730 template <typename somenumber>
1731 friend class SparseMatrix;
1732 template <typename somenumber>
1734 template <typename>
1735 friend class SparseILU;
1736
1737 // To allow it calling private prepare_add() and prepare_set().
1738 template <typename>
1739 friend class BlockMatrixBase;
1740
1741 // Also give access to internal details to the iterator/accessor classes.
1742 template <typename, bool>
1744 template <typename, bool>
1746
1747# ifdef DEAL_II_WITH_MPI
1748 // Give access to internal datastructures to perform MPI operations.
1749 template <typename Number>
1750 friend void
1752 const MPI_Comm &,
1754# endif
1755};
1756
1757# ifndef DOXYGEN
1758/*---------------------- Inline functions -----------------------------------*/
1759
1760
1761
1762template <typename number>
1763inline typename SparseMatrix<number>::size_type
1765{
1766 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1767 return cols->rows;
1768}
1769
1770
1771template <typename number>
1772inline typename SparseMatrix<number>::size_type
1774{
1775 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1776 return cols->cols;
1777}
1778
1779
1780// Inline the set() and add() functions, since they will be called frequently.
1781template <typename number>
1782inline void
1784 const size_type j,
1785 const number value)
1786{
1788
1789 const size_type index = cols->operator()(i, j);
1790
1791 // it is allowed to set elements of the matrix that are not part of the
1792 // sparsity pattern, if the value to which we set it is zero
1794 {
1795 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1796 ExcInvalidIndex(i, j));
1797 return;
1798 }
1799
1800 val[index] = value;
1801}
1802
1803
1804
1805template <typename number>
1806template <typename number2>
1807inline void
1808SparseMatrix<number>::set(const std::vector<size_type> &indices,
1810 const bool elide_zero_values)
1811{
1812 Assert(indices.size() == values.m(),
1813 ExcDimensionMismatch(indices.size(), values.m()));
1814 Assert(values.m() == values.n(), ExcNotQuadratic());
1815
1816 for (size_type i = 0; i < indices.size(); ++i)
1817 set(indices[i],
1818 indices.size(),
1819 indices.data(),
1820 &values(i, 0),
1821 elide_zero_values);
1822}
1823
1824
1825
1826template <typename number>
1827template <typename number2>
1828inline void
1829SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1830 const std::vector<size_type> &col_indices,
1832 const bool elide_zero_values)
1833{
1834 Assert(row_indices.size() == values.m(),
1835 ExcDimensionMismatch(row_indices.size(), values.m()));
1836 Assert(col_indices.size() == values.n(),
1837 ExcDimensionMismatch(col_indices.size(), values.n()));
1838
1839 for (size_type i = 0; i < row_indices.size(); ++i)
1840 set(row_indices[i],
1841 col_indices.size(),
1842 col_indices.data(),
1843 &values(i, 0),
1844 elide_zero_values);
1845}
1846
1847
1848
1849template <typename number>
1850template <typename number2>
1851inline void
1853 const std::vector<size_type> &col_indices,
1854 const std::vector<number2> & values,
1855 const bool elide_zero_values)
1856{
1857 Assert(col_indices.size() == values.size(),
1858 ExcDimensionMismatch(col_indices.size(), values.size()));
1859
1860 set(row,
1861 col_indices.size(),
1862 col_indices.data(),
1863 values.data(),
1864 elide_zero_values);
1865}
1866
1867
1868
1869template <typename number>
1870inline void
1872 const size_type j,
1873 const number value)
1874{
1876
1877 if (value == number())
1878 return;
1879
1880 const size_type index = cols->operator()(i, j);
1881
1882 // it is allowed to add elements to the matrix that are not part of the
1883 // sparsity pattern, if the value to which we set it is zero
1885 {
1886 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1887 ExcInvalidIndex(i, j));
1888 return;
1889 }
1890
1891 val[index] += value;
1892}
1893
1894
1895
1896template <typename number>
1897template <typename number2>
1898inline void
1899SparseMatrix<number>::add(const std::vector<size_type> &indices,
1901 const bool elide_zero_values)
1902{
1903 Assert(indices.size() == values.m(),
1904 ExcDimensionMismatch(indices.size(), values.m()));
1905 Assert(values.m() == values.n(), ExcNotQuadratic());
1906
1907 for (size_type i = 0; i < indices.size(); ++i)
1908 add(indices[i],
1909 indices.size(),
1910 indices.data(),
1911 &values(i, 0),
1912 elide_zero_values);
1913}
1914
1915
1916
1917template <typename number>
1918template <typename number2>
1919inline void
1920SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1921 const std::vector<size_type> &col_indices,
1923 const bool elide_zero_values)
1924{
1925 Assert(row_indices.size() == values.m(),
1926 ExcDimensionMismatch(row_indices.size(), values.m()));
1927 Assert(col_indices.size() == values.n(),
1928 ExcDimensionMismatch(col_indices.size(), values.n()));
1929
1930 for (size_type i = 0; i < row_indices.size(); ++i)
1931 add(row_indices[i],
1932 col_indices.size(),
1933 col_indices.data(),
1934 &values(i, 0),
1935 elide_zero_values);
1936}
1937
1938
1939
1940template <typename number>
1941template <typename number2>
1942inline void
1944 const std::vector<size_type> &col_indices,
1945 const std::vector<number2> & values,
1946 const bool elide_zero_values)
1947{
1948 Assert(col_indices.size() == values.size(),
1949 ExcDimensionMismatch(col_indices.size(), values.size()));
1950
1951 add(row,
1952 col_indices.size(),
1953 col_indices.data(),
1954 values.data(),
1955 elide_zero_values);
1956}
1957
1958
1959
1960template <typename number>
1961inline SparseMatrix<number> &
1962SparseMatrix<number>::operator*=(const number factor)
1963{
1964 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1965 Assert(val != nullptr, ExcNotInitialized());
1966
1967 number * val_ptr = val.get();
1968 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1969
1970 while (val_ptr != end_ptr)
1971 *val_ptr++ *= factor;
1972
1973 return *this;
1974}
1975
1976
1977
1978template <typename number>
1979inline SparseMatrix<number> &
1980SparseMatrix<number>::operator/=(const number factor)
1981{
1982 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1983 Assert(val != nullptr, ExcNotInitialized());
1984 Assert(factor != number(), ExcDivideByZero());
1985
1986 const number factor_inv = number(1.) / factor;
1987
1988 number * val_ptr = val.get();
1989 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
1990
1991 while (val_ptr != end_ptr)
1992 *val_ptr++ *= factor_inv;
1993
1994 return *this;
1995}
1996
1997
1998
1999template <typename number>
2000inline const number &
2002{
2003 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2004 Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2005 ExcInvalidIndex(i, j));
2006 return val[cols->operator()(i, j)];
2007}
2008
2009
2010
2011template <typename number>
2012inline number &
2014{
2015 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2016 Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2017 ExcInvalidIndex(i, j));
2018 return val[cols->operator()(i, j)];
2019}
2020
2021
2022
2023template <typename number>
2024inline number
2025SparseMatrix<number>::el(const size_type i, const size_type j) const
2026{
2027 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2028 const size_type index = cols->operator()(i, j);
2029
2031 return val[index];
2032 else
2033 return 0;
2034}
2035
2036
2037
2038template <typename number>
2039inline number
2041{
2042 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2043 Assert(m() == n(), ExcNotQuadratic());
2044 AssertIndexRange(i, m());
2045
2046 // Use that the first element in each row of a quadratic matrix is the main
2047 // diagonal
2048 return val[cols->rowstart[i]];
2049}
2050
2051
2052
2053template <typename number>
2054inline number &
2056{
2057 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2058 Assert(m() == n(), ExcNotQuadratic());
2059 AssertIndexRange(i, m());
2060
2061 // Use that the first element in each row of a quadratic matrix is the main
2062 // diagonal
2063 return val[cols->rowstart[i]];
2064}
2065
2066
2067
2068template <typename number>
2069template <typename ForwardIterator>
2070void
2071SparseMatrix<number>::copy_from(const ForwardIterator begin,
2072 const ForwardIterator end)
2073{
2074 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2075 ExcIteratorRange(std::distance(begin, end), m()));
2076
2077 // for use in the inner loop, we define an alias to the type of the inner
2078 // iterators
2079 using inner_iterator =
2080 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2081 size_type row = 0;
2082 for (ForwardIterator i = begin; i != end; ++i, ++row)
2083 {
2084 const inner_iterator end_of_row = i->end();
2085 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2086 // write entries
2087 set(row, j->first, j->second);
2088 };
2089}
2090
2091
2092//---------------------------------------------------------------------------
2093
2094
2095namespace SparseMatrixIterators
2096{
2097 template <typename number>
2098 inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2099 const std::size_t index_within_matrix)
2100 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2101 index_within_matrix)
2102 , matrix(matrix)
2103 {}
2104
2105
2106
2107 template <typename number>
2108 inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2109 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2110 , matrix(matrix)
2111 {}
2112
2113
2114
2115 template <typename number>
2116 inline Accessor<number, true>::Accessor(
2118 : SparsityPatternIterators::Accessor(a)
2119 , matrix(&a.get_matrix())
2120 {}
2121
2122
2123
2124 template <typename number>
2125 inline number
2126 Accessor<number, true>::value() const
2127 {
2128 AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2129 return matrix->val[linear_index];
2130 }
2131
2132
2133
2134 template <typename number>
2135 inline const typename Accessor<number, true>::MatrixType &
2136 Accessor<number, true>::get_matrix() const
2137 {
2138 return *matrix;
2139 }
2140
2141
2142
2143 template <typename number>
2144 inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2145 const bool)
2146 : accessor(accessor)
2147 {}
2148
2149
2150 template <typename number>
2151 inline Accessor<number, false>::Reference::operator number() const
2152 {
2153 AssertIndexRange(accessor->linear_index,
2154 accessor->matrix->n_nonzero_elements());
2155 return accessor->matrix->val[accessor->linear_index];
2156 }
2157
2158
2159
2160 template <typename number>
2161 inline const typename Accessor<number, false>::Reference &
2162 Accessor<number, false>::Reference::operator=(const number n) const
2163 {
2164 AssertIndexRange(accessor->linear_index,
2165 accessor->matrix->n_nonzero_elements());
2166 accessor->matrix->val[accessor->linear_index] = n;
2167 return *this;
2168 }
2169
2170
2171
2172 template <typename number>
2173 inline const typename Accessor<number, false>::Reference &
2174 Accessor<number, false>::Reference::operator+=(const number n) const
2175 {
2176 AssertIndexRange(accessor->linear_index,
2177 accessor->matrix->n_nonzero_elements());
2178 accessor->matrix->val[accessor->linear_index] += n;
2179 return *this;
2180 }
2181
2182
2183
2184 template <typename number>
2185 inline const typename Accessor<number, false>::Reference &
2186 Accessor<number, false>::Reference::operator-=(const number n) const
2187 {
2188 AssertIndexRange(accessor->linear_index,
2189 accessor->matrix->n_nonzero_elements());
2190 accessor->matrix->val[accessor->linear_index] -= n;
2191 return *this;
2192 }
2193
2194
2195
2196 template <typename number>
2197 inline const typename Accessor<number, false>::Reference &
2198 Accessor<number, false>::Reference::operator*=(const number n) const
2199 {
2200 AssertIndexRange(accessor->linear_index,
2201 accessor->matrix->n_nonzero_elements());
2202 accessor->matrix->val[accessor->linear_index] *= n;
2203 return *this;
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 Accessor<number, false>::Accessor(MatrixType * matrix,
2222 const std::size_t index)
2223 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2224 , matrix(matrix)
2225 {}
2226
2227
2228
2229 template <typename number>
2230 inline Accessor<number, false>::Accessor(MatrixType *matrix)
2231 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2232 , matrix(matrix)
2233 {}
2234
2235
2236
2237 template <typename number>
2238 inline typename Accessor<number, false>::Reference
2239 Accessor<number, false>::value() const
2240 {
2241 return Reference(this, true);
2242 }
2243
2244
2245
2246 template <typename number>
2247 inline typename Accessor<number, false>::MatrixType &
2248 Accessor<number, false>::get_matrix() const
2249 {
2250 return *matrix;
2251 }
2252
2253
2254
2255 template <typename number, bool Constness>
2256 inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2257 const std::size_t index)
2258 : accessor(matrix, index)
2259 {}
2260
2261
2262
2263 template <typename number, bool Constness>
2264 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2265 : accessor(matrix)
2266 {}
2267
2268
2269
2270 template <typename number, bool Constness>
2271 inline Iterator<number, Constness>::Iterator(
2273 : accessor(*i)
2274 {}
2275
2276
2277
2278 template <typename number, bool Constness>
2279 inline const Iterator<number, Constness> &
2280 Iterator<number, Constness>::
2282 {
2283 accessor = *i;
2284 return *this;
2285 }
2286
2287
2288
2289 template <typename number, bool Constness>
2290 inline Iterator<number, Constness> &
2292 {
2293 accessor.advance();
2294 return *this;
2295 }
2296
2297
2298 template <typename number, bool Constness>
2299 inline Iterator<number, Constness>
2301 {
2302 const Iterator iter = *this;
2303 accessor.advance();
2304 return iter;
2305 }
2306
2307
2308 template <typename number, bool Constness>
2309 inline const Accessor<number, Constness> &Iterator<number, Constness>::
2310 operator*() const
2311 {
2312 return accessor;
2313 }
2314
2315
2316 template <typename number, bool Constness>
2317 inline const Accessor<number, Constness> *Iterator<number, Constness>::
2318 operator->() const
2319 {
2320 return &accessor;
2321 }
2322
2323
2324 template <typename number, bool Constness>
2325 inline bool
2326 Iterator<number, Constness>::operator==(const Iterator &other) const
2327 {
2328 return (accessor == other.accessor);
2329 }
2330
2331
2332 template <typename number, bool Constness>
2333 inline bool
2334 Iterator<number, Constness>::operator!=(const Iterator &other) const
2335 {
2336 return !(*this == other);
2337 }
2338
2339
2340 template <typename number, bool Constness>
2341 inline bool
2342 Iterator<number, Constness>::operator<(const Iterator &other) const
2343 {
2344 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2346
2347 return (accessor < other.accessor);
2348 }
2349
2350
2351 template <typename number, bool Constness>
2352 inline bool
2353 Iterator<number, Constness>::operator>(const Iterator &other) const
2354 {
2355 return (other < *this);
2356 }
2357
2358
2359 template <typename number, bool Constness>
2360 inline int
2361 Iterator<number, Constness>::operator-(const Iterator &other) const
2362 {
2363 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2365
2366 return (*this)->linear_index - other->linear_index;
2367 }
2368
2369
2370
2371 template <typename number, bool Constness>
2372 inline Iterator<number, Constness>
2374 {
2375 Iterator x = *this;
2376 for (size_type i = 0; i < n; ++i)
2377 ++x;
2378
2379 return x;
2380 }
2381
2382} // namespace SparseMatrixIterators
2383
2384
2385
2386template <typename number>
2389{
2390 return const_iterator(this, 0);
2391}
2392
2393
2394template <typename number>
2397{
2398 return const_iterator(this);
2399}
2400
2401
2402template <typename number>
2403inline typename SparseMatrix<number>::iterator
2405{
2406 return iterator(this, 0);
2407}
2408
2409
2410template <typename number>
2411inline typename SparseMatrix<number>::iterator
2413{
2414 return iterator(this, cols->rowstart[cols->rows]);
2415}
2416
2417
2418template <typename number>
2421{
2422 AssertIndexRange(r, m());
2423
2424 return const_iterator(this, cols->rowstart[r]);
2425}
2426
2427
2428
2429template <typename number>
2432{
2433 AssertIndexRange(r, m());
2434
2435 return const_iterator(this, cols->rowstart[r + 1]);
2436}
2437
2438
2439
2440template <typename number>
2441inline typename SparseMatrix<number>::iterator
2443{
2444 AssertIndexRange(r, m());
2445
2446 return iterator(this, cols->rowstart[r]);
2447}
2448
2449
2450
2451template <typename number>
2452inline typename SparseMatrix<number>::iterator
2454{
2455 AssertIndexRange(r, m());
2456
2457 return iterator(this, cols->rowstart[r + 1]);
2458}
2459
2460
2461
2462template <typename number>
2463template <class StreamType>
2464inline void
2465SparseMatrix<number>::print(StreamType &out,
2466 const bool across,
2467 const bool diagonal_first) const
2468{
2469 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2470 Assert(val != nullptr, ExcNotInitialized());
2471
2472 bool hanging_diagonal = false;
2473 number diagonal = number();
2474
2475 for (size_type i = 0; i < cols->rows; ++i)
2476 {
2477 for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2478 {
2479 if (!diagonal_first && i == cols->colnums[j])
2480 {
2481 diagonal = val[j];
2482 hanging_diagonal = true;
2483 }
2484 else
2485 {
2486 if (hanging_diagonal && cols->colnums[j] > i)
2487 {
2488 if (across)
2489 out << ' ' << i << ',' << i << ':' << diagonal;
2490 else
2491 out << '(' << i << ',' << i << ") " << diagonal
2492 << std::endl;
2493 hanging_diagonal = false;
2494 }
2495 if (across)
2496 out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2497 else
2498 out << "(" << i << "," << cols->colnums[j] << ") " << val[j]
2499 << std::endl;
2500 }
2501 }
2502 if (hanging_diagonal)
2503 {
2504 if (across)
2505 out << ' ' << i << ',' << i << ':' << diagonal;
2506 else
2507 out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2508 hanging_diagonal = false;
2509 }
2510 }
2511 if (across)
2512 out << std::endl;
2513}
2514
2515
2516template <typename number>
2517inline void
2519{
2520 // nothing to do here
2521}
2522
2523
2524
2525template <typename number>
2526inline void
2528{
2529 // nothing to do here
2530}
2531
2532# endif // DOXYGEN
2533
2534
2535/*---------------------------- sparse_matrix.h ---------------------------*/
2536
2538
2539#endif
2540/*---------------------------- sparse_matrix.h ---------------------------*/
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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 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
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
iterator end()
number & diag_element(const size_type i)
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
void compress(::VectorOperation::values)
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
const_iterator end() const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void set(const std::vector< size_type > &indices, const FullMatrix< number2 > &full_matrix, const bool elide_zero_values=false)
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
const SparsityPattern & get_sparsity_pattern() const
void set(const size_type i, const size_type j, const number value)
const_iterator begin() const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
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
SparseMatrix< number > & copy_from(const TrilinosWrappers::SparseMatrix &matrix)
void symmetrize()
SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id)
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void SOR(Vector< somenumber > &v, const number om=1.) const
void vmult(OutVector &dst, const InVector &src) const
std::size_t memory_consumption() const
void copy_from(const ForwardIterator begin, const ForwardIterator end)
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const number & operator()(const size_type i, const size_type j) const
number & operator()(const size_type i, const size_type j)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
number el(const size_type i, const size_type j) const
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) 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
const_iterator end(const size_type r) const
iterator begin()
SparseMatrix & operator/=(const number factor)
void TSOR(Vector< somenumber > &v, const number om=1.) const
iterator begin(const size_type r)
typename numbers::NumberTraits< number >::real_type real_type
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)
iterator end(const size_type r)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
SparseMatrix & operator*=(const number factor)
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)
number value_type
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)
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)
Definition: vector.h:110
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typenameProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
__global__ void set(Number *val, const Number s, const size_type N)
void prepare_add()
static ::ExceptionBase & ExcDifferentSparsityPatterns()
void print_as_numpy_arrays(std::ostream &out, const unsigned int precision=9) const
static ::ExceptionBase & ExcInvalidIndex(int arg1, int arg2)
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
static const bool zero_addition_can_be_elided
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
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcSourceEqualsDestination()
void print_pattern(std::ostream &out, const double threshold=0.) const
#define AssertIsFinite(number)
Definition: exceptions.h:1721
void block_write(std::ostream &out) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
void block_read(std::istream &in)
void prepare_set()
void print(StreamType &out, const bool across=false, const bool diagonal_first=true) const
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:493
Expression operator>(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcInternalError()
std::unique_ptr< number[]> val
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()
std::size_t max_len
static const size_type invalid_entry
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
static const char V
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
types::global_dof_index size_type
Definition: sparse_matrix.h:82
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
Definition: types.h:76
BarycentricPolynomial< dim, Number1 > operator-(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
BarycentricPolynomial< dim, Number1 > operator+(const Number2 &a, const BarycentricPolynomial< dim, Number1 > &bp)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)