Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
sparse_matrix.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_sparse_matrix_h
17#define dealii_sparse_matrix_h
18
19#include <deal.II/base/config.h>
20
24
29
30#include <iterator>
31#include <memory>
32
33
35
36// Forward declarations
37#ifndef DOXYGEN
38template <typename number>
39class Vector;
40template <typename number>
41class FullMatrix;
42template <typename Matrix>
43class BlockMatrixBase;
44template <typename number>
45class SparseILU;
46# ifdef DEAL_II_WITH_MPI
47namespace Utilities
48{
49 namespace MPI
50 {
51 template <typename Number>
52 void
54 }
55} // namespace Utilities
56# endif
57
58# ifdef DEAL_II_WITH_TRILINOS
59namespace TrilinosWrappers
60{
61 class SparseMatrix;
62}
63# endif
64#endif
65
76{
81
82 // forward declaration
83 template <typename number, bool Constness>
84 class Iterator;
85
96 template <typename number, bool Constness>
98 {
99 public:
103 number
104 value() const;
105
109 number &
111
117 get_matrix() const;
118 };
119
120
121
128 template <typename number>
129 class Accessor<number, true> : public SparsityPatternIterators::Accessor
130 {
131 public:
137
141 Accessor(MatrixType *matrix, const std::size_t index_within_matrix);
142
147
152
156 number
157 value() const;
158
163 const MatrixType &
164 get_matrix() const;
165
166 private:
171
176
177 // Make iterator class a friend.
178 template <typename, bool>
179 friend class Iterator;
180 };
181
182
189 template <typename number>
190 class Accessor<number, false> : public SparsityPatternIterators::Accessor
191 {
192 private:
217 class Reference
218 {
219 public:
224 Reference(const Accessor *accessor, const bool dummy);
225
229 operator number() const;
230
234 const Reference &
235 operator=(const number n) const;
236
240 const Reference &
241 operator+=(const number n) const;
242
246 const Reference &
247 operator-=(const number n) const;
248
252 const Reference &
253 operator*=(const number n) const;
254
258 const Reference &
259 operator/=(const number n) const;
260
261 private:
267 };
268
269 public:
275
279 Accessor(MatrixType *matrix, const std::size_t index);
280
285
289 Reference
290 value() const;
291
296 MatrixType &
297 get_matrix() const;
298
299 private:
304
309
310 // Make iterator class a friend.
311 template <typename, bool>
312 friend class Iterator;
313 };
314
315
316
346 template <typename number, bool Constness>
348 {
349 public:
354
360
366
371 Iterator(MatrixType *matrix, const std::size_t index_within_matrix);
372
377
383
389
393 Iterator &
395
401
406 operator*() const;
407
412 operator->() const;
413
417 bool
418 operator==(const Iterator &) const;
419
423 bool
424 operator!=(const Iterator &) const;
425
433 bool
434 operator<(const Iterator &) const;
435
440 bool
441 operator>(const Iterator &) const;
442
449 int
450 operator-(const Iterator &p) const;
451
456 operator+(const size_type n) const;
457
458 private:
463 };
464
465} // namespace SparseMatrixIterators
466
468
469namespace std
470{
471 template <typename number, bool Constness>
473 ::SparseMatrixIterators::Iterator<number, Constness>>
474 {
475 using iterator_category = forward_iterator_tag;
477 typename ::SparseMatrixIterators::Iterator<number,
478 Constness>::value_type;
479 using difference_type = typename ::SparseMatrixIterators::
480 Iterator<number, Constness>::difference_type;
481 };
482} // namespace std
483
485
491// TODO: Add multithreading to the other vmult functions.
492
519template <typename number>
520class SparseMatrix : public virtual Subscriptor
521{
522public:
527
532 using value_type = number;
533
544
550
558
565 struct Traits
566 {
571 static const bool zero_addition_can_be_elided = true;
572 };
573
589
599
608
622 explicit SparseMatrix(const SparsityPattern &sparsity);
623
630 SparseMatrix(const SparsityPattern &sparsity, const IdentityMatrix &id);
631
636 virtual ~SparseMatrix() override;
637
649
656
665
678 operator=(const double d);
679
693 virtual void
694 reinit(const SparsityPattern &sparsity);
695
702 template <typename number2>
703 void
704 reinit(const SparseMatrix<number2> &sparse_matrix);
705
711 virtual void
722 bool
723 empty() const;
724
730 m() const;
731
737 n() const;
738
743 get_row_length(const size_type row) const;
744
750 std::size_t
752
762 std::size_t
763 n_actually_nonzero_elements(const double threshold = 0.) const;
764
773 const SparsityPattern &
775
780 std::size_t
782
787
798 void
799 set(const size_type i, const size_type j, const number value);
800
816 template <typename number2>
817 void
818 set(const std::vector<size_type> &indices,
819 const FullMatrix<number2> & full_matrix,
820 const bool elide_zero_values = false);
821
827 template <typename number2>
828 void
829 set(const std::vector<size_type> &row_indices,
830 const std::vector<size_type> &col_indices,
831 const FullMatrix<number2> & full_matrix,
832 const bool elide_zero_values = false);
833
844 template <typename number2>
845 void
846 set(const size_type row,
847 const std::vector<size_type> &col_indices,
848 const std::vector<number2> & values,
849 const bool elide_zero_values = false);
850
860 template <typename number2>
861 void
862 set(const size_type row,
863 const size_type n_cols,
864 const size_type *col_indices,
865 const number2 * values,
866 const bool elide_zero_values = false);
867
873 void
874 add(const size_type i, const size_type j, const number value);
875
890 template <typename number2>
891 void
892 add(const std::vector<size_type> &indices,
893 const FullMatrix<number2> & full_matrix,
894 const bool elide_zero_values = true);
895
901 template <typename number2>
902 void
903 add(const std::vector<size_type> &row_indices,
904 const std::vector<size_type> &col_indices,
905 const FullMatrix<number2> & full_matrix,
906 const bool elide_zero_values = true);
907
917 template <typename number2>
918 void
919 add(const size_type row,
920 const std::vector<size_type> &col_indices,
921 const std::vector<number2> & values,
922 const bool elide_zero_values = true);
923
933 template <typename number2>
934 void
935 add(const size_type row,
936 const size_type n_cols,
937 const size_type *col_indices,
938 const number2 * values,
939 const bool elide_zero_values = true,
940 const bool col_indices_are_sorted = false);
941
946 operator*=(const number factor);
947
952 operator/=(const number factor);
953
966 void
968
985 template <typename somenumber>
988
1005 template <typename ForwardIterator>
1006 void
1007 copy_from(const ForwardIterator begin, const ForwardIterator end);
1008
1018 template <typename somenumber>
1019 void
1021
1022#ifdef DEAL_II_WITH_TRILINOS
1034#endif
1035
1047 template <typename somenumber>
1048 void
1049 add(const number factor, const SparseMatrix<somenumber> &matrix);
1050
1070 const number &
1071 operator()(const size_type i, const size_type j) const;
1072
1076 number &
1077 operator()(const size_type i, const size_type j);
1078
1091 number
1092 el(const size_type i, const size_type j) const;
1093
1103 number
1104 diag_element(const size_type i) const;
1105
1110 number &
1112
1134 template <class OutVector, class InVector>
1135 void
1136 vmult(OutVector &dst, const InVector &src) const;
1137
1153 template <class OutVector, class InVector>
1154 void
1155 Tvmult(OutVector &dst, const InVector &src) const;
1156
1173 template <class OutVector, class InVector>
1174 void
1175 vmult_add(OutVector &dst, const InVector &src) const;
1176
1192 template <class OutVector, class InVector>
1193 void
1194 Tvmult_add(OutVector &dst, const InVector &src) const;
1195
1213 template <typename somenumber>
1214 somenumber
1216
1222 template <typename somenumber>
1223 somenumber
1225 const Vector<somenumber> &v) const;
1226
1236 template <typename somenumber>
1237 somenumber
1239 const Vector<somenumber> &x,
1240 const Vector<somenumber> &b) const;
1241
1277 template <typename numberB, typename numberC>
1278 void
1280 const SparseMatrix<numberB> &B,
1281 const Vector<number> & V = Vector<number>(),
1282 const bool rebuild_sparsity_pattern = true) const;
1283
1308 template <typename numberB, typename numberC>
1309 void
1311 const SparseMatrix<numberB> &B,
1312 const Vector<number> & V = Vector<number>(),
1313 const bool rebuild_sparsity_pattern = true) const;
1314
1328 real_type
1329 l1_norm() const;
1330
1338 real_type
1340
1345 real_type
1358 template <typename somenumber>
1359 void
1361 const Vector<somenumber> &src,
1362 const number omega = 1.) const;
1363
1370 template <typename somenumber>
1371 void
1373 const Vector<somenumber> & src,
1374 const number omega = 1.,
1375 const std::vector<std::size_t> &pos_right_of_diagonal =
1376 std::vector<std::size_t>()) const;
1377
1381 template <typename somenumber>
1382 void
1384 const Vector<somenumber> &src,
1385 const number omega = 1.) const;
1386
1390 template <typename somenumber>
1391 void
1393 const Vector<somenumber> &src,
1394 const number omega = 1.) const;
1395
1401 template <typename somenumber>
1402 void
1403 SSOR(Vector<somenumber> &v, const number omega = 1.) const;
1404
1409 template <typename somenumber>
1410 void
1411 SOR(Vector<somenumber> &v, const number omega = 1.) const;
1412
1417 template <typename somenumber>
1418 void
1419 TSOR(Vector<somenumber> &v, const number omega = 1.) const;
1420
1431 template <typename somenumber>
1432 void
1434 const std::vector<size_type> &permutation,
1435 const std::vector<size_type> &inverse_permutation,
1436 const number omega = 1.) const;
1437
1448 template <typename somenumber>
1449 void
1451 const std::vector<size_type> &permutation,
1452 const std::vector<size_type> &inverse_permutation,
1453 const number omega = 1.) const;
1454
1460 template <typename somenumber>
1461 void
1463 const Vector<somenumber> &b,
1464 const number omega = 1.) const;
1465
1470 template <typename somenumber>
1471 void
1473 const Vector<somenumber> &b,
1474 const number omega = 1.) const;
1475
1480 template <typename somenumber>
1481 void
1483 const Vector<somenumber> &b,
1484 const number omega = 1.) const;
1485
1490 template <typename somenumber>
1491 void
1493 const Vector<somenumber> &b,
1494 const number omega = 1.) const;
1508 begin() const;
1509
1513 iterator
1515
1520 end() const;
1521
1525 iterator
1527
1538 begin(const size_type r) const;
1539
1543 iterator
1545
1556 end(const size_type r) const;
1557
1561 iterator
1562 end(const size_type r);
1580 template <class StreamType>
1581 void
1582 print(StreamType &out,
1583 const bool across = false,
1584 const bool diagonal_first = true) const;
1585
1606 void
1607 print_formatted(std::ostream & out,
1608 const unsigned int precision = 3,
1609 const bool scientific = true,
1610 const unsigned int width = 0,
1611 const char * zero_string = " ",
1612 const double denominator = 1.) const;
1613
1619 void
1620 print_pattern(std::ostream &out, const double threshold = 0.) const;
1621
1630 void
1631 print_as_numpy_arrays(std::ostream & out,
1632 const unsigned int precision = 9) const;
1633
1644 void
1645 block_write(std::ostream &out) const;
1646
1663 void
1664 block_read(std::istream &in);
1675 int,
1676 int,
1677 << "You are trying to access the matrix entry with index <"
1678 << arg1 << ',' << arg2
1679 << ">, but this entry does not exist in the sparsity pattern "
1680 "of this matrix."
1681 "\n\n"
1682 "The most common cause for this problem is that you used "
1683 "a method to build the sparsity pattern that did not "
1684 "(completely) take into account all of the entries you "
1685 "will later try to write into. An example would be "
1686 "building a sparsity pattern that does not include "
1687 "the entries you will write into due to constraints "
1688 "on degrees of freedom such as hanging nodes or periodic "
1689 "boundary conditions. In such cases, building the "
1690 "sparsity pattern will succeed, but you will get errors "
1691 "such as the current one at one point or other when "
1692 "trying to write into the entries of the matrix.");
1697 "When copying one sparse matrix into another, "
1698 "or when adding one sparse matrix to another, "
1699 "both matrices need to refer to the same "
1700 "sparsity pattern.");
1705 int,
1706 int,
1707 << "The iterators denote a range of " << arg1
1708 << " elements, but the given number of rows was " << arg2);
1713 "You are attempting an operation on two vectors that "
1714 "are the same object, but the operation requires that the "
1715 "two objects are in fact different.");
1718protected:
1729 void
1731
1736 void
1738
1739private:
1746
1754 std::unique_ptr<number[]> val;
1755
1762 std::size_t max_len;
1763
1764 // make all other sparse matrices friends
1765 template <typename somenumber>
1766 friend class SparseMatrix;
1767 template <typename somenumber>
1769 template <typename>
1770 friend class SparseILU;
1771
1772 // To allow it calling private prepare_add() and prepare_set().
1773 template <typename>
1774 friend class BlockMatrixBase;
1775
1776 // Also give access to internal details to the iterator/accessor classes.
1777 template <typename, bool>
1779 template <typename, bool>
1781
1782#ifdef DEAL_II_WITH_MPI
1783 // Give access to internal datastructures to perform MPI operations.
1784 template <typename Number>
1785 friend void
1787 const MPI_Comm,
1789#endif
1790};
1791
1792#ifndef DOXYGEN
1793/*---------------------- Inline functions -----------------------------------*/
1794
1795
1796
1797template <typename number>
1798template <typename number2>
1799void
1801{
1802 this->reinit(sparse_matrix.get_sparsity_pattern());
1803}
1804
1805
1806
1807template <typename number>
1808inline typename SparseMatrix<number>::size_type
1810{
1811 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1812 return cols->rows;
1813}
1814
1815
1816template <typename number>
1817inline typename SparseMatrix<number>::size_type
1819{
1820 Assert(cols != nullptr, ExcNeedsSparsityPattern());
1821 return cols->cols;
1822}
1823
1824
1825// Inline the set() and add() functions, since they will be called frequently.
1826template <typename number>
1827inline void
1829 const size_type j,
1830 const number value)
1831{
1833
1834 const size_type index = cols->operator()(i, j);
1835
1836 // it is allowed to set elements of the matrix that are not part of the
1837 // sparsity pattern, if the value to which we set it is zero
1839 {
1840 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1841 ExcInvalidIndex(i, j));
1842 return;
1843 }
1844
1845 val[index] = value;
1846}
1847
1848
1849
1850template <typename number>
1851template <typename number2>
1852inline void
1853SparseMatrix<number>::set(const std::vector<size_type> &indices,
1854 const FullMatrix<number2> & values,
1855 const bool elide_zero_values)
1856{
1857 Assert(indices.size() == values.m(),
1858 ExcDimensionMismatch(indices.size(), values.m()));
1859 Assert(values.m() == values.n(), ExcNotQuadratic());
1860
1861 for (size_type i = 0; i < indices.size(); ++i)
1862 set(indices[i],
1863 indices.size(),
1864 indices.data(),
1865 &values(i, 0),
1866 elide_zero_values);
1867}
1868
1869
1870
1871template <typename number>
1872template <typename number2>
1873inline void
1874SparseMatrix<number>::set(const std::vector<size_type> &row_indices,
1875 const std::vector<size_type> &col_indices,
1876 const FullMatrix<number2> & values,
1877 const bool elide_zero_values)
1878{
1879 Assert(row_indices.size() == values.m(),
1880 ExcDimensionMismatch(row_indices.size(), values.m()));
1881 Assert(col_indices.size() == values.n(),
1882 ExcDimensionMismatch(col_indices.size(), values.n()));
1883
1884 for (size_type i = 0; i < row_indices.size(); ++i)
1885 set(row_indices[i],
1886 col_indices.size(),
1887 col_indices.data(),
1888 &values(i, 0),
1889 elide_zero_values);
1890}
1891
1892
1893
1894template <typename number>
1895template <typename number2>
1896inline void
1898 const std::vector<size_type> &col_indices,
1899 const std::vector<number2> & values,
1900 const bool elide_zero_values)
1901{
1902 Assert(col_indices.size() == values.size(),
1903 ExcDimensionMismatch(col_indices.size(), values.size()));
1904
1905 set(row,
1906 col_indices.size(),
1907 col_indices.data(),
1908 values.data(),
1909 elide_zero_values);
1910}
1911
1912
1913
1914template <typename number>
1915inline void
1917 const size_type j,
1918 const number value)
1919{
1921
1922 if (value == number())
1923 return;
1924
1925 const size_type index = cols->operator()(i, j);
1926
1927 // it is allowed to add elements to the matrix that are not part of the
1928 // sparsity pattern, if the value to which we set it is zero
1930 {
1931 Assert((index != SparsityPattern::invalid_entry) || (value == number()),
1932 ExcInvalidIndex(i, j));
1933 return;
1934 }
1935
1936 val[index] += value;
1937}
1938
1939
1940
1941template <typename number>
1942template <typename number2>
1943inline void
1944SparseMatrix<number>::add(const std::vector<size_type> &indices,
1945 const FullMatrix<number2> & values,
1946 const bool elide_zero_values)
1947{
1948 Assert(indices.size() == values.m(),
1949 ExcDimensionMismatch(indices.size(), values.m()));
1950 Assert(values.m() == values.n(), ExcNotQuadratic());
1951
1952 for (size_type i = 0; i < indices.size(); ++i)
1953 add(indices[i],
1954 indices.size(),
1955 indices.data(),
1956 &values(i, 0),
1957 elide_zero_values);
1958}
1959
1960
1961
1962template <typename number>
1963template <typename number2>
1964inline void
1965SparseMatrix<number>::add(const std::vector<size_type> &row_indices,
1966 const std::vector<size_type> &col_indices,
1967 const FullMatrix<number2> & values,
1968 const bool elide_zero_values)
1969{
1970 Assert(row_indices.size() == values.m(),
1971 ExcDimensionMismatch(row_indices.size(), values.m()));
1972 Assert(col_indices.size() == values.n(),
1973 ExcDimensionMismatch(col_indices.size(), values.n()));
1974
1975 for (size_type i = 0; i < row_indices.size(); ++i)
1976 add(row_indices[i],
1977 col_indices.size(),
1978 col_indices.data(),
1979 &values(i, 0),
1980 elide_zero_values);
1981}
1982
1983
1984
1985template <typename number>
1986template <typename number2>
1987inline void
1989 const std::vector<size_type> &col_indices,
1990 const std::vector<number2> & values,
1991 const bool elide_zero_values)
1992{
1993 Assert(col_indices.size() == values.size(),
1994 ExcDimensionMismatch(col_indices.size(), values.size()));
1995
1996 add(row,
1997 col_indices.size(),
1998 col_indices.data(),
1999 values.data(),
2000 elide_zero_values,
2001 std::is_sorted(col_indices.begin(), col_indices.end()));
2002}
2003
2004
2005
2006template <typename number>
2007inline SparseMatrix<number> &
2008SparseMatrix<number>::operator*=(const number factor)
2009{
2010 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2011 Assert(val != nullptr, ExcNotInitialized());
2012
2013 number * val_ptr = val.get();
2014 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2015
2016 while (val_ptr != end_ptr)
2017 *val_ptr++ *= factor;
2018
2019 return *this;
2020}
2021
2022
2023
2024template <typename number>
2025inline SparseMatrix<number> &
2026SparseMatrix<number>::operator/=(const number factor)
2027{
2028 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2029 Assert(val != nullptr, ExcNotInitialized());
2030 Assert(factor != number(), ExcDivideByZero());
2031
2032 const number factor_inv = number(1.) / factor;
2033
2034 number * val_ptr = val.get();
2035 const number *const end_ptr = val.get() + cols->n_nonzero_elements();
2036
2037 while (val_ptr != end_ptr)
2038 *val_ptr++ *= factor_inv;
2039
2040 return *this;
2041}
2042
2043
2044
2045template <typename number>
2046inline const number &
2048{
2049 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2050 Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2051 ExcInvalidIndex(i, j));
2052 return val[cols->operator()(i, j)];
2053}
2054
2055
2056
2057template <typename number>
2058inline number &
2060{
2061 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2062 Assert(cols->operator()(i, j) != SparsityPattern::invalid_entry,
2063 ExcInvalidIndex(i, j));
2064 return val[cols->operator()(i, j)];
2065}
2066
2067
2068
2069template <typename number>
2070inline number
2071SparseMatrix<number>::el(const size_type i, const size_type j) const
2072{
2073 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2074 const size_type index = cols->operator()(i, j);
2075
2077 return val[index];
2078 else
2079 return 0;
2080}
2081
2082
2083
2084template <typename number>
2085inline number
2087{
2088 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2089 Assert(m() == n(), ExcNotQuadratic());
2090 AssertIndexRange(i, m());
2091
2092 // Use that the first element in each row of a quadratic matrix is the main
2093 // diagonal
2094 return val[cols->rowstart[i]];
2095}
2096
2097
2098
2099template <typename number>
2100inline number &
2102{
2103 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2104 Assert(m() == n(), ExcNotQuadratic());
2105 AssertIndexRange(i, m());
2106
2107 // Use that the first element in each row of a quadratic matrix is the main
2108 // diagonal
2109 return val[cols->rowstart[i]];
2110}
2111
2112
2113
2114template <typename number>
2115template <typename ForwardIterator>
2116void
2117SparseMatrix<number>::copy_from(const ForwardIterator begin,
2118 const ForwardIterator end)
2119{
2120 Assert(static_cast<size_type>(std::distance(begin, end)) == m(),
2121 ExcIteratorRange(std::distance(begin, end), m()));
2122
2123 // for use in the inner loop, we define an alias to the type of the inner
2124 // iterators
2125 using inner_iterator =
2126 typename std::iterator_traits<ForwardIterator>::value_type::const_iterator;
2127 size_type row = 0;
2128 for (ForwardIterator i = begin; i != end; ++i, ++row)
2129 {
2130 const inner_iterator end_of_row = i->end();
2131 for (inner_iterator j = i->begin(); j != end_of_row; ++j)
2132 // write entries
2133 set(row, j->first, j->second);
2134 };
2135}
2136
2137
2138//---------------------------------------------------------------------------
2139
2140
2141namespace SparseMatrixIterators
2142{
2143 template <typename number>
2144 inline Accessor<number, true>::Accessor(const MatrixType *matrix,
2145 const std::size_t index_within_matrix)
2146 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(),
2147 index_within_matrix)
2148 , matrix(matrix)
2149 {}
2150
2151
2152
2153 template <typename number>
2154 inline Accessor<number, true>::Accessor(const MatrixType *matrix)
2155 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2156 , matrix(matrix)
2157 {}
2158
2159
2160
2161 template <typename number>
2162 inline Accessor<number, true>::Accessor(
2164 : SparsityPatternIterators::Accessor(a)
2165 , matrix(&a.get_matrix())
2166 {}
2167
2168
2169
2170 template <typename number>
2171 inline number
2172 Accessor<number, true>::value() const
2173 {
2174 AssertIndexRange(linear_index, matrix->n_nonzero_elements());
2175 return matrix->val[linear_index];
2176 }
2177
2178
2179
2180 template <typename number>
2181 inline const typename Accessor<number, true>::MatrixType &
2182 Accessor<number, true>::get_matrix() const
2183 {
2184 return *matrix;
2185 }
2186
2187
2188
2189 template <typename number>
2190 inline Accessor<number, false>::Reference::Reference(const Accessor *accessor,
2191 const bool)
2192 : accessor(accessor)
2193 {}
2194
2195
2196 template <typename number>
2197 inline Accessor<number, false>::Reference::operator number() const
2198 {
2199 AssertIndexRange(accessor->linear_index,
2200 accessor->matrix->n_nonzero_elements());
2201 return accessor->matrix->val[accessor->linear_index];
2202 }
2203
2204
2205
2206 template <typename number>
2207 inline const typename Accessor<number, false>::Reference &
2208 Accessor<number, false>::Reference::operator=(const number n) const
2209 {
2210 AssertIndexRange(accessor->linear_index,
2211 accessor->matrix->n_nonzero_elements());
2212 accessor->matrix->val[accessor->linear_index] = n;
2213 return *this;
2214 }
2215
2216
2217
2218 template <typename number>
2219 inline const typename Accessor<number, false>::Reference &
2220 Accessor<number, false>::Reference::operator+=(const number n) const
2221 {
2222 AssertIndexRange(accessor->linear_index,
2223 accessor->matrix->n_nonzero_elements());
2224 accessor->matrix->val[accessor->linear_index] += n;
2225 return *this;
2226 }
2227
2228
2229
2230 template <typename number>
2231 inline const typename Accessor<number, false>::Reference &
2232 Accessor<number, false>::Reference::operator-=(const number n) const
2233 {
2234 AssertIndexRange(accessor->linear_index,
2235 accessor->matrix->n_nonzero_elements());
2236 accessor->matrix->val[accessor->linear_index] -= n;
2237 return *this;
2238 }
2239
2240
2241
2242 template <typename number>
2243 inline const typename Accessor<number, false>::Reference &
2244 Accessor<number, false>::Reference::operator*=(const number n) const
2245 {
2246 AssertIndexRange(accessor->linear_index,
2247 accessor->matrix->n_nonzero_elements());
2248 accessor->matrix->val[accessor->linear_index] *= n;
2249 return *this;
2250 }
2251
2252
2253
2254 template <typename number>
2255 inline const typename Accessor<number, false>::Reference &
2256 Accessor<number, false>::Reference::operator/=(const number n) const
2257 {
2258 AssertIndexRange(accessor->linear_index,
2259 accessor->matrix->n_nonzero_elements());
2260 accessor->matrix->val[accessor->linear_index] /= n;
2261 return *this;
2262 }
2263
2264
2265
2266 template <typename number>
2267 inline Accessor<number, false>::Accessor(MatrixType * matrix,
2268 const std::size_t index)
2269 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern(), index)
2270 , matrix(matrix)
2271 {}
2272
2273
2274
2275 template <typename number>
2276 inline Accessor<number, false>::Accessor(MatrixType *matrix)
2277 : SparsityPatternIterators::Accessor(&matrix->get_sparsity_pattern())
2278 , matrix(matrix)
2279 {}
2280
2281
2282
2283 template <typename number>
2284 inline typename Accessor<number, false>::Reference
2285 Accessor<number, false>::value() const
2286 {
2287 return Reference(this, true);
2288 }
2289
2290
2291
2292 template <typename number>
2293 inline typename Accessor<number, false>::MatrixType &
2294 Accessor<number, false>::get_matrix() const
2295 {
2296 return *matrix;
2297 }
2298
2299
2300
2301 template <typename number, bool Constness>
2302 inline Iterator<number, Constness>::Iterator(MatrixType * matrix,
2303 const std::size_t index)
2304 : accessor(matrix, index)
2305 {}
2306
2307
2308
2309 template <typename number, bool Constness>
2310 inline Iterator<number, Constness>::Iterator(MatrixType *matrix)
2311 : accessor(matrix)
2312 {}
2313
2314
2315
2316 template <typename number, bool Constness>
2317 inline Iterator<number, Constness>::Iterator(
2319 : accessor(*i)
2320 {}
2321
2322
2323
2324 template <typename number, bool Constness>
2325 inline const Iterator<number, Constness> &
2326 Iterator<number, Constness>::operator=(
2328 {
2329 accessor = *i;
2330 return *this;
2331 }
2332
2333
2334
2335 template <typename number, bool Constness>
2336 inline Iterator<number, Constness> &
2337 Iterator<number, Constness>::operator++()
2338 {
2339 accessor.advance();
2340 return *this;
2341 }
2342
2343
2344 template <typename number, bool Constness>
2345 inline Iterator<number, Constness>
2346 Iterator<number, Constness>::operator++(int)
2347 {
2348 const Iterator iter = *this;
2349 accessor.advance();
2350 return iter;
2351 }
2352
2353
2354 template <typename number, bool Constness>
2355 inline const Accessor<number, Constness> &
2356 Iterator<number, Constness>::operator*() const
2357 {
2358 return accessor;
2359 }
2360
2361
2362 template <typename number, bool Constness>
2363 inline const Accessor<number, Constness> *
2364 Iterator<number, Constness>::operator->() const
2365 {
2366 return &accessor;
2367 }
2368
2369
2370 template <typename number, bool Constness>
2371 inline bool
2372 Iterator<number, Constness>::operator==(const Iterator &other) const
2373 {
2374 return (accessor == other.accessor);
2375 }
2376
2377
2378 template <typename number, bool Constness>
2379 inline bool
2380 Iterator<number, Constness>::operator!=(const Iterator &other) const
2381 {
2382 return !(*this == other);
2383 }
2384
2385
2386 template <typename number, bool Constness>
2387 inline bool
2388 Iterator<number, Constness>::operator<(const Iterator &other) const
2389 {
2390 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2392
2393 return (accessor < other.accessor);
2394 }
2395
2396
2397 template <typename number, bool Constness>
2398 inline bool
2399 Iterator<number, Constness>::operator>(const Iterator &other) const
2400 {
2401 return (other < *this);
2402 }
2403
2404
2405 template <typename number, bool Constness>
2406 inline int
2407 Iterator<number, Constness>::operator-(const Iterator &other) const
2408 {
2409 Assert(&accessor.get_matrix() == &other.accessor.get_matrix(),
2411
2412 return (*this)->linear_index - other->linear_index;
2413 }
2414
2415
2416
2417 template <typename number, bool Constness>
2418 inline Iterator<number, Constness>
2419 Iterator<number, Constness>::operator+(const size_type n) const
2420 {
2421 Iterator x = *this;
2422 for (size_type i = 0; i < n; ++i)
2423 ++x;
2424
2425 return x;
2426 }
2427
2428} // namespace SparseMatrixIterators
2429
2430
2431
2432template <typename number>
2435{
2436 return const_iterator(this, 0);
2437}
2438
2439
2440template <typename number>
2443{
2444 return const_iterator(this);
2445}
2446
2447
2448template <typename number>
2449inline typename SparseMatrix<number>::iterator
2451{
2452 return iterator(this, 0);
2453}
2454
2455
2456template <typename number>
2457inline typename SparseMatrix<number>::iterator
2459{
2460 return iterator(this, cols->rowstart[cols->rows]);
2461}
2462
2463
2464template <typename number>
2466SparseMatrix<number>::begin(const size_type r) const
2467{
2468 AssertIndexRange(r, m());
2469
2470 return const_iterator(this, cols->rowstart[r]);
2471}
2472
2473
2474
2475template <typename number>
2477SparseMatrix<number>::end(const size_type r) const
2478{
2479 AssertIndexRange(r, m());
2480
2481 return const_iterator(this, cols->rowstart[r + 1]);
2482}
2483
2484
2485
2486template <typename number>
2487inline typename SparseMatrix<number>::iterator
2488SparseMatrix<number>::begin(const size_type r)
2489{
2490 AssertIndexRange(r, m());
2491
2492 return iterator(this, cols->rowstart[r]);
2493}
2494
2495
2496
2497template <typename number>
2498inline typename SparseMatrix<number>::iterator
2499SparseMatrix<number>::end(const size_type r)
2500{
2501 AssertIndexRange(r, m());
2502
2503 return iterator(this, cols->rowstart[r + 1]);
2504}
2505
2506
2507
2508template <typename number>
2509template <class StreamType>
2510inline void
2511SparseMatrix<number>::print(StreamType &out,
2512 const bool across,
2513 const bool diagonal_first) const
2514{
2515 Assert(cols != nullptr, ExcNeedsSparsityPattern());
2516 Assert(val != nullptr, ExcNotInitialized());
2517
2518 bool hanging_diagonal = false;
2519 number diagonal = number();
2520
2521 for (size_type i = 0; i < cols->rows; ++i)
2522 {
2523 for (size_type j = cols->rowstart[i]; j < cols->rowstart[i + 1]; ++j)
2524 {
2525 if (!diagonal_first && i == cols->colnums[j])
2526 {
2527 diagonal = val[j];
2528 hanging_diagonal = true;
2529 }
2530 else
2531 {
2532 if (hanging_diagonal && cols->colnums[j] > i)
2533 {
2534 if (across)
2535 out << ' ' << i << ',' << i << ':' << diagonal;
2536 else
2537 out << '(' << i << ',' << i << ") " << diagonal
2538 << std::endl;
2539 hanging_diagonal = false;
2540 }
2541 if (across)
2542 out << ' ' << i << ',' << cols->colnums[j] << ':' << val[j];
2543 else
2544 out << '(' << i << ',' << cols->colnums[j] << ") " << val[j]
2545 << std::endl;
2546 }
2547 }
2548 if (hanging_diagonal)
2549 {
2550 if (across)
2551 out << ' ' << i << ',' << i << ':' << diagonal;
2552 else
2553 out << '(' << i << ',' << i << ") " << diagonal << std::endl;
2554 hanging_diagonal = false;
2555 }
2556 }
2557 if (across)
2558 out << std::endl;
2559}
2560
2561
2562template <typename number>
2563inline void
2565{
2566 // nothing to do here
2567}
2568
2569
2570
2571template <typename number>
2572inline void
2574{
2575 // nothing to do here
2576}
2577
2578#endif // DOXYGEN
2579
2580
2582
2583#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
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
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)
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
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
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
__global__ void set(Number *val, const Number s, const size_type N)
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:533
static ::ExceptionBase & ExcNeedsSparsityPattern()
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
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:82
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