Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3065-g80d88f1806 2025-04-14 16: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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
1818
1819template <typename number>
1820inline typename SparseMatrix<number>::size_type
1822{
1823 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1824 return cols->cols;
1825}
1826
1827
1828
1829template <typename number>
1830inline const SparsityPattern &
1832{
1833 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1834 return *cols;
1835}
1836
1837
1838
1839// Inline the set() and add() functions, since they will be called frequently.
1840template <typename number>
1841inline void
1843 const size_type j,
1844 const number value)
1845{
1847
1848 const size_type index = cols->operator()(i, j);
1849
1850 // it is allowed to set elements of the matrix that are not part of the
1851 // sparsity pattern, if the value to which we set it is zero
1853 {
1854 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1855 ExcInvalidIndex(i, j));
1856 return;
1857 }
1858
1859 val[index] = value;
1860}
1861
1862
1863
1864template <typename number>
1865template <typename number2>
1866inline void
1867SparseMatrix<number>::set(const std::vector<size_type> &indices,
1868 const FullMatrix<number2> &values,
1869 const bool elide_zero_values)
1870{
1871 Assert(indices.size() == values.m(),
1872 ExcDimensionMismatch(indices.size(), values.m()));
1873 Assert(values.m() == values.n(), ExcNotQuadratic());
1874
1875 for (size_type i = 0; i < indices.size(); ++i)
1876 set(indices[i],
1877 indices.size(),
1878 indices.data(),
1879 &values(i, 0),
1880 elide_zero_values);
1881}
1882
1883
1884
1885template <typename number>
1886template <typename number2>
1887inline void
1888SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1889 const std::vector<size_type> &col_indices,
1890 const FullMatrix<number2> &values,
1891 const bool elide_zero_values)
1892{
1893 Assert(row_indices.size() == values.m(),
1894 ExcDimensionMismatch(row_indices.size(), values.m()));
1895 Assert(col_indices.size() == values.n(),
1896 ExcDimensionMismatch(col_indices.size(), values.n()));
1897
1898 for (size_type i = 0; i < row_indices.size(); ++i)
1899 set(row_indices[i],
1900 col_indices.size(),
1901 col_indices.data(),
1902 &values(i, 0),
1903 elide_zero_values);
1904}
1905
1906
1907
1908template <typename number>
1909template <typename number2>
1910inline void
1912 const std::vector<size_type> &col_indices,
1913 const std::vector<number2> &values,
1914 const bool elide_zero_values)
1915{
1916 Assert(col_indices.size() == values.size(),
1917 ExcDimensionMismatch(col_indices.size(), values.size()));
1918
1919 set(row,
1920 col_indices.size(),
1921 col_indices.data(),
1922 values.data(),
1923 elide_zero_values);
1924}
1925
1926
1927
1928template <typename number>
1929inline void
1931 const size_type j,
1932 const number value)
1933{
1935
1936 if (value == number())
1937 return;
1938
1939 const size_type index = cols->operator()(i, j);
1940
1941 // it is allowed to add elements to the matrix that are not part of the
1942 // sparsity pattern, if the value to which we set it is zero
1944 {
1945 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1946 ExcInvalidIndex(i, j));
1947 return;
1948 }
1949
1950 val[index] += value;
1951}
1952
1953
1954
1955template <typename number>
1956template <typename number2>
1957inline void
1958SparseMatrix<number>::add(const std::vector<size_type> &indices,
1959 const FullMatrix<number2> &values,
1960 const bool elide_zero_values)
1961{
1962 Assert(indices.size() == values.m(),
1963 ExcDimensionMismatch(indices.size(), values.m()));
1964 Assert(values.m() == values.n(), ExcNotQuadratic());
1965
1966 for (size_type i = 0; i < indices.size(); ++i)
1967 add(indices[i],
1968 indices.size(),
1969 indices.data(),
1970 &values(i, 0),
1971 elide_zero_values);
1972}
1973
1974
1975
1976template <typename number>
1977template <typename number2>
1978inline void
1979SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1980 const std::vector<size_type> &col_indices,
1981 const FullMatrix<number2> &values,
1982 const bool elide_zero_values)
1983{
1984 Assert(row_indices.size() == values.m(),
1985 ExcDimensionMismatch(row_indices.size(), values.m()));
1986 Assert(col_indices.size() == values.n(),
1987 ExcDimensionMismatch(col_indices.size(), values.n()));
1988
1989 for (size_type i = 0; i < row_indices.size(); ++i)
1990 add(row_indices[i],
1991 col_indices.size(),
1992 col_indices.data(),
1993 &values(i, 0),
1994 elide_zero_values);
1995}
1996
1997
1998
1999template <typename number>
2000template <typename number2>
2001inline void
2003 const std::vector<size_type> &col_indices,
2004 const std::vector<number2> &values,
2005 const bool elide_zero_values)
2006{
2007 Assert(col_indices.size() == values.size(),
2008 ExcDimensionMismatch(col_indices.size(), values.size()));
2009
2010 add(row,
2011 col_indices.size(),
2012 col_indices.data(),
2013 values.data(),
2014 elide_zero_values,
2015 std::is_sorted(col_indices.begin(), col_indices.end()));
2016}
2017
2018
2019
2020template <typename number>
2021inline SparseMatrix<number> &
2022SparseMatrix<number>::operator*=(const number factor)
2023{
2024 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2025 Assert(val != nullptr, ExcNotInitialized());
2026
2027 number *val_ptr = val.get();
2028 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2029
2030 while (val_ptr != end_ptr)
2031 *val_ptr++ *= factor;
2032
2033 return *this;
2034}
2035
2036
2037
2038template <typename number>
2039inline SparseMatrix<number> &
2040SparseMatrix<number>::operator/=(const number factor)
2041{
2042 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2043 Assert(val != nullptr, ExcNotInitialized());
2044 Assert(factor != number(), ExcDivideByZero());
2045
2046 const number factor_inv = number(1.) / factor;
2047
2048 number *val_ptr = val.get();
2049 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2050
2051 while (val_ptr != end_ptr)
2052 *val_ptr++ *= factor_inv;
2053
2054 return *this;
2055}
2056
2057
2058
2059template <typename number>
2060inline const 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 &
2074{
2075 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2076 Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2077 ExcInvalidIndex(i, j));
2078 return val[cols->operator()(i, j)];
2079}
2080
2081
2082
2083template <typename number>
2084inline number
2085SparseMatrix<number>::el(const size_type i, const size_type j) const
2086{
2087 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2088 const size_type index = cols->operator()(i, j);
2089
2091 return val[index];
2092 else
2093 return 0;
2094}
2095
2096
2097
2098template <typename number>
2099inline number
2101{
2102 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2103 Assert(m() == n(), ExcNotQuadratic());
2104 AssertIndexRange(i, m());
2105
2106 // Use that the first element in each row of a quadratic matrix is the main
2107 // diagonal
2108 return val[cols->rowstart[i]];
2109}
2110
2111
2112
2113template <typename number>
2114inline number &
2116{
2117 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2118 Assert(m() == n(), ExcNotQuadratic());
2119 AssertIndexRange(i, m());
2120
2121 // Use that the first element in each row of a quadratic matrix is the main
2122 // diagonal
2123 return val[cols->rowstart[i]];
2124}
2125
2126
2127
2128template <typename number>
2129template <typename ForwardIterator>
2130void
2131SparseMatrix<number>::copy_from(const ForwardIterator begin,
2132 const ForwardIterator end)
2133{
2134 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2135 ExcIteratorRange(std::distance(begin, end), m()));
2136
2137 // for use in the inner loop, we define an alias to the type of the inner
2138 // iterators
2139 using inner_iterator =
2140 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2141 size_type row = 0;
2142 for (ForwardIterator i = begin; i != end; ++i, ++row)
2143 {
2144 const inner_iterator end_of_row = i->end();
2145 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2146 // write entries
2147 set(row, j->first, j->second);
2148 };
2149}
2150
2151
2152//---------------------------------------------------------------------------
2153
2154
2155namespace SparseMatrixIterators
2156{
2157 template <typename number>
2158 inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2159 const std::size_t index_within_matrix)
2160 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2161 index_within_matrix)
2162 , matrix(matrix)
2163 {}
2164
2165
2166
2167 template <typename number>
2168 inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2169 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2170 , matrix(matrix)
2171 {}
2172
2173
2174
2175 template <typename number>
2176 inline Accessor<number, true>::Accessor(
2178 : SparsityPatternIterators::Accessor(a)
2179 , matrix(&a.get_matrix())
2180 {}
2181
2182
2183
2184 template <typename number>
2185 inline number
2186 Accessor<number, true>::value() const
2187 {
2188 AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2189 return matrix->val[linear_index];
2190 }
2191
2192
2193
2194 template <typename number>
2195 inline const typename Accessor<number, true>::MatrixType &
2196 Accessor<number, true>::get_matrix() const
2197 {
2198 return *matrix;
2199 }
2200
2201
2202
2203 template <typename number>
2204 inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2205 const bool)
2206 : accessor(accessor)
2207 {}
2208
2209
2210 template <typename number>
2211 inline Accessor<number, false>::Reference::operator number() const
2212 {
2213 AssertIndexRange(accessor->linear_index,
2214 accessor->matrix->n_nonzero_elements());
2215 return accessor->matrix->val[accessor->linear_index];
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 const typename Accessor<number, false>::Reference &
2270 Accessor<number, false>::Reference::operator/=(const number n) const
2271 {
2272 AssertIndexRange(accessor->linear_index,
2273 accessor->matrix->n_nonzero_elements());
2274 accessor->matrix->val[accessor->linear_index] /= n;
2275 return *this;
2276 }
2277
2278
2279
2280 template <typename number>
2281 inline Accessor<number, false>::Accessor(MatrixType *matrix,
2282 const std::size_t index)
2283 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2284 , matrix(matrix)
2285 {}
2286
2287
2288
2289 template <typename number>
2290 inline Accessor<number, false>::Accessor(MatrixType *matrix)
2291 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2292 , matrix(matrix)
2293 {}
2294
2295
2296
2297 template <typename number>
2298 inline typename Accessor<number, false>::Reference
2299 Accessor<number, false>::value() const
2300 {
2301 return Reference(this, true);
2302 }
2303
2304
2305
2306 template <typename number>
2307 inline typename Accessor<number, false>::MatrixType &
2308 Accessor<number, false>::get_matrix() const
2309 {
2310 return *matrix;
2311 }
2312
2313
2314
2315 template <typename number, bool Constness>
2316 inline Iterator<number, Constness>::Iterator(MatrixType *matrix,
2317 const std::size_t index)
2318 : accessor(matrix, index)
2319 {}
2320
2321
2322
2323 template <typename number, bool Constness>
2324 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2325 : accessor(matrix)
2326 {}
2327
2328
2329
2330 template <typename number, bool Constness>
2331 inline Iterator<number, Constness>::Iterator(
2333 : accessor(*i)
2334 {}
2335
2336
2337
2338 template <typename number, bool Constness>
2339 inline const Iterator<number, Constness> &
2340 Iterator<number, Constness>::operator=(
2342 {
2343 accessor = *i;
2344 return *this;
2345 }
2346
2347
2348
2349 template <typename number, bool Constness>
2350 inline Iterator<number, Constness> &
2351 Iterator<number, Constness>::operator++()
2352 {
2353 accessor.advance();
2354 return *this;
2355 }
2356
2357
2358 template <typename number, bool Constness>
2359 inline Iterator<number, Constness>
2360 Iterator<number, Constness>::operator++(int)
2361 {
2362 const Iterator iter = *this;
2363 accessor.advance();
2364 return iter;
2365 }
2366
2367
2368 template <typename number, bool Constness>
2369 inline const Accessor<number, Constness> &
2370 Iterator<number, Constness>::operator*() const
2371 {
2372 return accessor;
2373 }
2374
2375
2376 template <typename number, bool Constness>
2377 inline const Accessor<number, Constness> *
2378 Iterator<number, Constness>::operator->() const
2379 {
2380 return &accessor;
2381 }
2382
2383
2384 template <typename number, bool Constness>
2385 inline bool
2386 Iterator<number, Constness>::operator==(const Iterator &other) const
2387 {
2388 return (accessor == other.accessor);
2389 }
2390
2391
2392 template <typename number, bool Constness>
2393 inline bool
2394 Iterator<number, Constness>::operator!=(const Iterator &other) const
2395 {
2396 return !(*this == other);
2397 }
2398
2399
2400 template <typename number, bool Constness>
2401 inline bool
2402 Iterator<number, Constness>::operator<(const Iterator &other) const
2403 {
2404 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2406
2407 return (accessor < other.accessor);
2408 }
2409
2410
2411 template <typename number, bool Constness>
2412 inline bool
2413 Iterator<number, Constness>::operator>(const Iterator &other) const
2414 {
2415 return (other < *this);
2416 }
2417
2418
2419 template <typename number, bool Constness>
2420 inline int
2421 Iterator<number, Constness>::operator-(const Iterator &other) const
2422 {
2423 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2425
2426 return (*this)->linear_index - other->linear_index;
2427 }
2428
2429
2430
2431 template <typename number, bool Constness>
2432 inline Iterator<number, Constness>
2433 Iterator<number, Constness>::operator+(const size_type n) const
2434 {
2435 Iterator x = *this;
2436 for (size_type i = 0; i < n; ++i)
2437 ++x;
2438
2439 return x;
2440 }
2441
2442} // namespace SparseMatrixIterators
2443
2444
2445
2446template <typename number>
2449{
2450 return const_iterator(this, 0);
2451}
2452
2453
2454template <typename number>
2457{
2458 return const_iterator(this);
2459}
2460
2461
2462template <typename number>
2463inline typename SparseMatrix<number>::iterator
2465{
2466 return iterator(this, 0);
2467}
2468
2469
2470template <typename number>
2471inline typename SparseMatrix<number>::iterator
2473{
2474 return iterator(this, cols->rowstart[cols->rows]);
2475}
2476
2477
2478template <typename number>
2480SparseMatrix<number>::begin(const size_type r) const
2481{
2482 AssertIndexRange(r, m());
2483
2484 return const_iterator(this, cols->rowstart[r]);
2485}
2486
2487
2488
2489template <typename number>
2491SparseMatrix<number>::end(const size_type r) const
2492{
2493 AssertIndexRange(r, m());
2494
2495 return const_iterator(this, cols->rowstart[r + 1]);
2496}
2497
2498
2499
2500template <typename number>
2501inline typename SparseMatrix<number>::iterator
2502SparseMatrix<number>::begin(const size_type r)
2503{
2504 AssertIndexRange(r, m());
2505
2506 return iterator(this, cols->rowstart[r]);
2507}
2508
2509
2510
2511template <typename number>
2512inline typename SparseMatrix<number>::iterator
2513SparseMatrix<number>::end(const size_type r)
2514{
2515 AssertIndexRange(r, m());
2516
2517 return iterator(this, cols->rowstart[r + 1]);
2518}
2519
2520
2521
2522template <typename number>
2523template <typename StreamType>
2524inline void
2525SparseMatrix<number>::print(StreamType &out,
2526 const bool across,
2527 const bool diagonal_first) const
2528{
2529 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2530 Assert(val != nullptr, ExcNotInitialized());
2531
2532 bool hanging_diagonal = false;
2533 number diagonal = number();
2534
2535 for (size_type i = 0; i < cols->rows; ++i)
2536 {
2537 for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2538 {
2539 if (!diagonal_first && i == cols->colnums[j])
2540 {
2541 diagonal = val[j];
2542 hanging_diagonal = true;
2543 }
2544 else
2545 {
2546 if (hanging_diagonal && cols->colnums[j] > i)
2547 {
2548 if (across)
2549 out << ' ' << i << ',' << i << ':' << diagonal;
2550 else
2551 out << '(' << i << ',' << i << ") " << diagonal
2552 << std::endl;
2553 hanging_diagonal = false;
2554 }
2555 if (across)
2556 out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2557 else
2558 out << '(' << i << ',' << cols->colnums[j] << ") " << val[j]
2559 << std::endl;
2560 }
2561 }
2562 if (hanging_diagonal)
2563 {
2564 if (across)
2565 out << ' ' << i << ',' << i << ':' << diagonal;
2566 else
2567 out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2568 hanging_diagonal = false;
2569 }
2570 }
2571 if (across)
2572 out << std::endl;
2573}
2574
2575
2576template <typename number>
2577inline void
2579{
2580 // nothing to do here
2581}
2582
2583
2584
2585template <typename number>
2586inline void
2588{
2589 // nothing to do here
2590}
2591
2592#endif // DOXYGEN
2593
2594
2596
2597#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:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
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)
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
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:94
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