Reference documentation for deal.II version GIT relicensing-487-ge9eb5ab491 2024-04-25 07:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
trilinos_tpetra_sparse_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2024 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_trilinos_tpetra_sparse_matrix_h
16#define dealii_trilinos_tpetra_sparse_matrix_h
17
18#include <deal.II/base/config.h>
19
20#ifdef DEAL_II_TRILINOS_WITH_TPETRA
21
25
30# include <deal.II/lac/vector.h>
31
32// Tpetra includes
33# include <Tpetra_Core.hpp>
34# include <Tpetra_CrsMatrix.hpp>
35
36# include <type_traits>
37
38
40
41// forward declarations
42
43template <typename Number>
44class SparseMatrix;
45
46# ifndef DOXYGEN
47template <typename MatrixType>
48class BlockMatrixBase;
49
50namespace LinearAlgebra
51{
52 namespace TpetraWrappers
53 {
54 template <typename MemorySpace>
55 class SparsityPattern;
56
57 namespace SparseMatrixIterators
58 {
59 template <typename Number, typename MemorySpace, bool Constness>
60 class Iterator;
61 }
62 } // namespace TpetraWrappers
63} // namespace LinearAlgebra
64# endif
65
66namespace LinearAlgebra
67{
68
69 namespace TpetraWrappers
70 {
109 template <typename Number, typename MemorySpace = ::MemorySpace::Host>
111 {
112 public:
117
122 std::size_t,
123 << "You tried to access row " << arg1
124 << " of a non-contiguous locally owned row set."
125 << " The row " << arg1
126 << " is not stored locally and can't be accessed.");
127
135 struct Traits
136 {
141 static const bool zero_addition_can_be_elided = true;
142 };
143
147 using iterator =
149
155
160 using value_type = Number;
161
165# if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0)
166 using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode<
167 typename MemorySpace::kokkos_space::execution_space,
168 typename MemorySpace::kokkos_space>;
169# else
170 using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode<
171 typename MemorySpace::kokkos_space::execution_space,
172 typename MemorySpace::kokkos_space>;
173# endif
174
179 Tpetra::CrsMatrix<Number,
180 int,
182 NodeType>;
183
187 using MapType =
188 Tpetra::Map<int, ::types::signed_global_dof_index, NodeType>;
189
193 using GraphType =
194 Tpetra::CrsGraph<int, ::types::signed_global_dof_index, NodeType>;
195
199 using VectorType = Tpetra::
200 Vector<Number, int, ::types::signed_global_dof_index, NodeType>;
201
210
215
224 const size_type n,
225 const unsigned int n_max_entries_per_row);
226
235 const size_type n,
236 const std::vector<unsigned int> &n_entries_per_row);
237
243
248
254
260
265 virtual ~SparseMatrix() override = default;
266
282 template <typename SparsityPatternType>
283 void
284 reinit(const SparsityPatternType &sparsity_pattern);
285
295 void
296 reinit(const SparsityPattern<MemorySpace> &sparsity_pattern);
315 SparseMatrix(const IndexSet &parallel_partitioning,
316 const MPI_Comm communicator = MPI_COMM_WORLD,
317 const unsigned int n_max_entries_per_row = 0);
318
327 SparseMatrix(const IndexSet &parallel_partitioning,
328 const MPI_Comm communicator,
329 const std::vector<unsigned int> &n_entries_per_row);
330
345 SparseMatrix(const IndexSet &row_parallel_partitioning,
346 const IndexSet &col_parallel_partitioning,
347 const MPI_Comm communicator = MPI_COMM_WORLD,
348 const size_type n_max_entries_per_row = 0);
349
358 SparseMatrix(const IndexSet &row_parallel_partitioning,
359 const IndexSet &col_parallel_partitioning,
360 const MPI_Comm communicator,
361 const std::vector<unsigned int> &n_entries_per_row);
362
382 template <typename SparsityPatternType>
383 std::enable_if_t<
384 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
385 reinit(const IndexSet &parallel_partitioning,
386 const SparsityPatternType &sparsity_pattern,
387 const MPI_Comm communicator = MPI_COMM_WORLD,
388 const bool exchange_data = false);
389
402 template <typename SparsityPatternType>
403 std::enable_if_t<
404 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
405 reinit(const IndexSet &row_parallel_partitioning,
406 const IndexSet &col_parallel_partitioning,
407 const SparsityPatternType &sparsity_pattern,
408 const MPI_Comm communicator = MPI_COMM_WORLD,
409 const bool exchange_data = false);
410
427 void
428 reinit(const IndexSet &row_parallel_partitioning,
429 const IndexSet &col_parallel_partitioning,
430 const ::SparseMatrix<Number> &dealii_sparse_matrix,
431 const MPI_Comm communicator = MPI_COMM_WORLD,
432 const double drop_tolerance = 1e-13,
433 const bool copy_values = true,
434 const ::SparsityPattern *use_this_sparsity = nullptr);
435
446 m() const;
447
452 n() const;
453
454
463 unsigned int
464 local_size() const;
465
474 std::pair<size_type, size_type>
475 local_range() const;
476
481 bool
482 in_local_range(const size_type index) const;
483
488 size_t
490
494 unsigned int
495 row_length(const size_type row) const;
496
503 bool
505
527 operator=(const double d);
528
533 operator*=(const Number factor);
534
539 operator/=(const Number factor);
540
544 void
546
556 void
557 add(const size_type i, const size_type j, const Number value);
558
572 void
573 add(const size_type row,
574 const size_type n_cols,
575 const size_type *col_indices,
576 const Number *values,
577 const bool elide_zero_values = true,
578 const bool col_indices_are_sorted = false);
579
587 void
588 add(const Number factor, const SparseMatrix<Number, MemorySpace> &matrix);
589
611 void
612 set(const size_type i, const size_type j, const Number value);
613
646 void
647 set(const std::vector<size_type> &indices,
648 const FullMatrix<Number> &full_matrix,
649 const bool elide_zero_values = false);
650
656 void
657 set(const std::vector<size_type> &row_indices,
658 const std::vector<size_type> &col_indices,
659 const FullMatrix<Number> &full_matrix,
660 const bool elide_zero_values = false);
661
689 void
690 set(const size_type row,
691 const std::vector<size_type> &col_indices,
692 const std::vector<Number> &values,
693 const bool elide_zero_values = false);
694
722 void
723 set(const size_type row,
724 const size_type n_cols,
725 const size_type *col_indices,
726 const Number *values,
727 const bool elide_zero_values = false);
728
755 void
756 clear_row(const size_type row, const Number new_diag_value = 0);
757
778 void
780 const Number new_diag_value = 0);
781
789 void
791
806 Number
807 operator()(const size_type i, const size_type j) const;
808
825 Number
826 el(const size_type i, const size_type j) const;
827
833 Number
834 diag_element(const size_type i) const;
835
841 /*
842 * Matrix-vector multiplication: let <i>dst = M*src</i> with <i>M</i>
843 * being this matrix.
844 *
845 * Source and destination must not be the same vector.
846 *
847 * The vector @p dst has to be initialized with the same IndexSet that was
848 * used for the row indices of the matrix and the vector @p src has to be
849 * initialized with the same IndexSet that was used for the column indices
850 * of the matrix.
851 */
852 template <typename InputVectorType>
853 void
854 vmult(InputVectorType &dst, const InputVectorType &src) const;
855
856 /*
857 * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
858 * <i>M</i> being this matrix. This function does the same as vmult() but
859 * takes the transposed matrix.
860 *
861 * Source and destination must not be the same vector.
862 */
863 template <typename InputVectorType>
864 void
865 Tvmult(InputVectorType &dst, const InputVectorType &src) const;
866
873 template <typename InputVectorType>
874 void
875 vmult_add(InputVectorType &dst, const InputVectorType &src) const;
876
884 template <typename InputVectorType>
885 void
886 Tvmult_add(InputVectorType &dst, const InputVectorType &src) const;
887
909 Number
911
931 Number
933 const Vector<Number, MemorySpace> &v) const;
934
951 Number
954 const Vector<Number, MemorySpace> &b) const;
955
967 Number
969
983 void
984 print(std::ostream &out,
985 const bool print_detailed_trilinos_information = false) const;
986
1018 void
1020
1027 void
1029
1036 const MatrixType &
1038
1045 MatrixType &
1047
1057 Teuchos::RCP<const MatrixType>
1059
1069 Teuchos::RCP<MatrixType>
1082 IndexSet
1084
1090 IndexSet
1092
1119 begin() const;
1120
1124 iterator
1126
1132 end() const;
1133
1137 iterator
1139
1169 begin(const size_type r) const;
1170
1174 iterator
1176
1187 end(const size_type r) const;
1188
1192 iterator
1193 end(const size_type r);
1194
1205
1211 "You are attempting an operation on two vectors that "
1212 "are the same object, but the operation requires that the "
1213 "two objects are in fact different.");
1214
1215 /*
1216 * Exception
1217 */
1219 "The column partitioning of a matrix does not match "
1220 "the partitioning of a vector you are trying to "
1221 "multiply it with. Are you multiplying the "
1222 "matrix with a vector that has ghost elements?");
1223
1224 /*
1225 * Exception
1226 */
1228 "The row partitioning of a matrix does not match "
1229 "the partitioning of a vector you are trying to "
1230 "put the result of a matrix-vector product in. "
1231 "Are you trying to put the product of the "
1232 "matrix with a vector into a vector that has "
1233 "ghost elements?");
1234
1239 size_type,
1240 size_type,
1241 << "The entry with index <" << arg1 << ',' << arg2
1242 << "> does not exist.");
1243
1248 size_type,
1249 size_type,
1250 size_type,
1251 size_type,
1252 << "You tried to access element (" << arg1 << '/' << arg2
1253 << ')'
1254 << " of a distributed matrix, but only rows in range ["
1255 << arg3 << ',' << arg4
1256 << "] are stored locally and can be accessed.");
1257
1260 private:
1264 Number
1265 element(const size_type i, const size_type j, const bool no_error) const;
1266
1276 Teuchos::RCP<MapType> column_space_map;
1277
1283 Teuchos::RCP<MatrixType> matrix;
1284
1290
1303 void
1305
1313 void
1315
1316 // To allow calling protected prepare_add() and prepare_set().
1317 friend class BlockMatrixBase<SparseMatrix<Number, MemorySpace>>;
1318 }; // class SparseMatrix
1319
1324 {
1329
1337 template <typename Number, typename MemorySpace>
1339 {
1340 public:
1345
1350 const size_type row,
1351 const size_type index);
1352
1356 size_type
1357 row() const;
1358
1362 size_type
1363 index() const;
1364
1368 size_type
1369 column() const;
1370
1371 protected:
1383
1388
1394 void
1396
1409 std::shared_ptr<std::vector<::types::signed_global_dof_index>>
1411
1415 std::shared_ptr<std::vector<Number>> value_cache;
1416
1417 private:
1418 friend class Iterator<Number, MemorySpace, false>;
1419 friend class Iterator<Number, MemorySpace, true>;
1420 };
1421
1432 template <typename Number, typename MemorySpace, bool Constness>
1433 class Accessor : public AccessorBase<Number, MemorySpace>
1434 {
1438 Number
1439 value() const;
1440
1444 Number &
1446 };
1447
1451 template <typename Number, typename MemorySpace>
1452 class Accessor<Number, MemorySpace, true>
1453 : public AccessorBase<Number, MemorySpace>
1454 {
1455 public:
1461
1466
1472 const size_type row,
1473 const size_type index);
1474
1479 template <bool Other>
1481
1485 Number
1486 value() const;
1487 };
1488
1492 template <typename Number, typename MemorySpace>
1493 class Accessor<Number, MemorySpace, false>
1494 : public AccessorBase<Number, MemorySpace>
1495 {
1496 class Reference
1497 {
1498 public:
1502 Reference(const Accessor<Number, MemorySpace, false> &accessor);
1503
1507 operator Number() const;
1508
1512 const Reference &
1513 operator=(const Number n) const;
1514
1518 const Reference &
1519 operator+=(const Number n) const;
1520
1524 const Reference &
1525 operator-=(const Number n) const;
1526
1530 const Reference &
1531 operator*=(const Number n) const;
1532
1536 const Reference &
1537 operator/=(const Number n) const;
1538
1539 private:
1545 };
1546
1547 public:
1553
1558
1564 const size_type row,
1565 const size_type index);
1566
1570 Reference
1571 value() const;
1572
1573 private:
1574 // Make Reference object a friend.
1575 friend class Reference;
1576 };
1577
1591 template <typename Number, typename MemorySpace, bool Constness>
1593 {
1594 public:
1599
1605
1610 using value_type = Number;
1611
1618
1623 Iterator(MatrixType *matrix,
1624 const size_type row,
1625 const size_type index);
1626
1630 template <bool Other>
1632
1637 operator++();
1638
1643 operator++(int);
1644
1649 operator*() const;
1650
1655 operator->() const;
1656
1661 template <bool OtherConstness>
1662 bool
1664
1668 template <bool OtherConstness>
1669 bool
1671
1677 template <bool OtherConstness>
1678 bool
1680
1684 template <bool OtherConstness>
1685 bool
1687
1692 size_type,
1693 size_type,
1694 << "Attempt to access element " << arg2 << " of row "
1695 << arg1 << " which doesn't have that many elements.");
1696
1697 private:
1702
1703 friend class Iterator<Number, MemorySpace, true>;
1704 friend class Iterator<Number, MemorySpace, false>;
1705 };
1706 } // namespace SparseMatrixIterators
1707
1708 } // namespace TpetraWrappers
1709
1710} // namespace LinearAlgebra
1711
1712
1714
1715namespace std
1716{
1717 template <typename Number, typename MemorySpace, bool Constness>
1720 Iterator<Number, MemorySpace, Constness>>
1721 {
1722 using iterator_category = forward_iterator_tag;
1724 typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
1725 Iterator<Number, MemorySpace, Constness>::value_type;
1727 typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
1728 Iterator<Number, MemorySpace, Constness>::difference_type;
1729 };
1730} // namespace std
1731
1732/* ------------------------- Inline functions ---------------------- */
1733
1734
1736
1737namespace LinearAlgebra
1738{
1739
1740 namespace TpetraWrappers
1741 {
1742 template <typename Number, typename MemorySpace>
1743 inline void
1745 const size_type j,
1746 const Number value)
1747 {
1748 set(i, 1, &j, &value, false);
1749 }
1750
1751
1752
1753 template <typename Number, typename MemorySpace>
1754 inline void
1756 const size_type j,
1757 const Number value)
1758 {
1759 add(i, 1, &j, &value, false);
1760 }
1761
1762
1763
1764 template <typename Number, typename MemorySpace>
1765 inline Number
1769 const Vector<Number, MemorySpace> &b) const
1770 {
1771 vmult(dst, x);
1772 dst -= b;
1773 dst *= -1.;
1774
1775 return dst.l2_norm();
1776 }
1777
1778
1779
1780 template <typename Number, typename MemorySpace>
1783 {
1784 return matrix->getRowMap()->getGlobalNumElements();
1785 }
1786
1787
1788
1789 template <typename Number, typename MemorySpace>
1792 {
1793 // If the matrix structure has not been fixed (i.e., we did not have a
1794 // sparsity pattern), it does not know about the number of columns, so we
1795 // must always take this from the additional column space map
1796 Assert(column_space_map.get() != nullptr, ExcInternalError());
1797 return column_space_map->getGlobalNumElements();
1798 }
1799
1800
1801
1802 template <typename Number, typename MemorySpace>
1803 inline bool
1805 {
1806 return compressed;
1807 }
1808
1809
1810
1811 template <typename Number, typename MemorySpace>
1814 {
1815 return begin(0);
1816 }
1817
1818
1819
1820 template <typename Number, typename MemorySpace>
1823 {
1824 return const_iterator(this, m(), 0);
1825 }
1826
1827
1828
1829 template <typename Number, typename MemorySpace>
1832 {
1833 AssertIndexRange(r, m());
1834 if (in_local_range(r) && (row_length(r) > 0))
1835 return const_iterator(this, r, 0);
1836 else
1837 return end(r);
1838 }
1839
1840
1841 template <typename Number, typename MemorySpace>
1844 {
1845 AssertIndexRange(r, m());
1846
1847 // place the iterator on the first entry
1848 // past this line, or at the end of the
1849 // matrix
1850 for (size_type i = r + 1; i < m(); ++i)
1851 if (in_local_range(i) && (row_length(i) > 0))
1852 return const_iterator(this, i, 0);
1853
1854 // if there is no such line, then take the
1855 // end iterator of the matrix
1856 return end();
1857 }
1858
1859
1860
1861 template <typename Number, typename MemorySpace>
1864 {
1865 return begin(0);
1866 }
1867
1868
1869
1870 template <typename Number, typename MemorySpace>
1873 {
1874 return iterator(this, m(), 0);
1875 }
1876
1877
1878
1879 template <typename Number, typename MemorySpace>
1882 {
1883 AssertIndexRange(r, m());
1884 if (in_local_range(r) && (row_length(r) > 0))
1885 return iterator(this, r, 0);
1886 else
1887 return end(r);
1888 }
1889
1890
1891
1892 template <typename Number, typename MemorySpace>
1895 {
1896 AssertIndexRange(r, m());
1897
1898 // place the iterator on the first entry
1899 // past this line, or at the end of the
1900 // matrix
1901 for (size_type i = r + 1; i < m(); ++i)
1902 if (in_local_range(i) && (row_length(i) > 0))
1903 return iterator(this, i, 0);
1904
1905 // if there is no such line, then take the
1906 // end iterator of the matrix
1907 return end();
1908 }
1909
1910
1911
1912 template <typename Number, typename MemorySpace>
1913 inline bool
1915 const size_type index) const
1916 {
1917 const size_type begin = matrix->getRowMap()->getMinGlobalIndex();
1918 const size_type end = matrix->getRowMap()->getMaxGlobalIndex() + 1;
1919
1920 return ((index >= begin) && (index < end));
1921 }
1922
1923
1924
1925 template <typename Number, typename MemorySpace>
1926 unsigned int
1928 {
1929 auto n_entries = matrix->getNumEntriesInGlobalRow(row);
1930 Assert(n_entries !=
1931 Teuchos::OrdinalTraits<decltype(n_entries)>::invalid(),
1932 ExcAccessToNonlocalRow(row));
1933
1934 return n_entries;
1935 }
1936
1937
1938
1939 template <typename Number, typename MemorySpace>
1940 inline void
1942 {
1943 // nothing to do here
1944 }
1945
1946
1947
1948 template <typename Number, typename MemorySpace>
1949 inline void
1951 {
1952 // nothing to do here
1953 }
1954
1955
1956
1957 template <typename Number, typename MemorySpace>
1958 inline const Tpetra::CrsMatrix<
1959 Number,
1960 int,
1964 {
1965 return *matrix;
1966 }
1967
1968
1969
1970 template <typename Number, typename MemorySpace>
1971 inline Tpetra::CrsMatrix<
1972 Number,
1973 int,
1977 {
1978 return *matrix;
1979 }
1980
1981
1982
1983 template <typename Number, typename MemorySpace>
1984 inline Teuchos::RCP<const Tpetra::CrsMatrix<
1985 Number,
1986 int,
1990 {
1991 return matrix.getConst();
1992 }
1993
1994
1995
1996 template <typename Number, typename MemorySpace>
1997 inline Teuchos::RCP<
1998 Tpetra::CrsMatrix<Number,
1999 int,
2003 {
2004 return matrix;
2005 }
2006
2007
2008
2009 template <typename Number, typename MemorySpace>
2010 inline IndexSet
2012 {
2013 return IndexSet(matrix->getDomainMap());
2014 }
2015
2016
2017
2018 template <typename Number, typename MemorySpace>
2019 inline IndexSet
2021 {
2022 return IndexSet(matrix->getRangeMap());
2023 }
2024
2025
2026 namespace SparseMatrixIterators
2027 {
2028 template <typename Number, typename MemorySpace>
2031 size_type row,
2032 size_type index)
2033 : matrix(matrix)
2034 , a_row(row)
2035 , a_index(index)
2036 {
2038 }
2039
2040
2041 template <typename Number, typename MemorySpace>
2044 {
2045 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2046 return a_row;
2047 }
2048
2049
2050 template <typename Number, typename MemorySpace>
2053 {
2054 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2055 return (*colnum_cache)[a_index];
2056 }
2057
2058
2059 template <typename Number, typename MemorySpace>
2062 {
2063 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2064 return a_index;
2065 }
2066
2067
2068
2069 template <typename Number, typename MemorySpace>
2070 void
2072 {
2073 // if we are asked to visit the past-the-end line, then simply
2074 // release all our caches and go on with life.
2075 //
2076 // do the same if the row we're supposed to visit is not locally
2077 // owned. this is simply going to make non-locally owned rows
2078 // look like they're empty
2079 if ((this->a_row == matrix->m()) ||
2080 (matrix->in_local_range(this->a_row) == false))
2081 {
2082 colnum_cache.reset();
2083 value_cache.reset();
2084
2085 return;
2086 }
2087
2088 // get a representation of the present row
2089 size_t ncols;
2091 matrix->row_length(this->a_row);
2092 if (value_cache.get() == nullptr)
2093 {
2094 value_cache = std::make_shared<std::vector<Number>>(colnums);
2095 colnum_cache = std::make_shared<
2096 std::vector<::types::signed_global_dof_index>>(colnums);
2097 }
2098 else
2099 {
2100 value_cache->resize(colnums);
2101 colnum_cache->resize(colnums);
2102 }
2103
2104# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
2106 nonconst_global_inds_host_view_type col_indices(colnum_cache->data(),
2107 colnums);
2109 nonconst_values_host_view_type values(value_cache->data(), colnums);
2110# else
2111 Teuchos::ArrayView<::types::signed_global_dof_index> col_indices(
2112 *colnum_cache);
2113 Teuchos::ArrayView<Number> values(*value_cache);
2114# endif
2115 matrix->trilinos_matrix().getGlobalRowCopy(this->a_row,
2116 col_indices,
2117 values,
2118 ncols);
2119
2120 AssertDimension(ncols, colnums);
2121
2122 // copy it into our caches if the
2123 // line isn't empty. if it is, then
2124 // we've done something wrong, since
2125 // we shouldn't have initialized an
2126 // iterator for an empty line (what
2127 // would it point to?)
2128 }
2129
2130
2131
2132 template <typename Number, typename MemorySpace>
2134 MatrixType *matrix,
2135 const size_type row,
2136 const size_type index)
2137 : AccessorBase<Number, MemorySpace>(
2138 const_cast<SparseMatrix<Number, MemorySpace> *>(matrix),
2139 row,
2140 index)
2141 {}
2142
2143
2144
2145 template <typename Number, typename MemorySpace>
2146 template <bool Other>
2151
2152
2153
2154 template <typename Number, typename MemorySpace>
2155 inline Number
2164
2165
2166
2167 template <typename Number, typename MemorySpace>
2170 : accessor(const_cast<Accessor<Number, MemorySpace, false> &>(acc))
2171 {}
2172
2173
2174
2175 template <typename Number, typename MemorySpace>
2177 const
2178 {
2179 return (*accessor.value_cache)[accessor.a_index];
2180 }
2181
2182
2183
2184 template <typename Number, typename MemorySpace>
2187 const Number n) const
2188 {
2189 (*accessor.value_cache)[accessor.a_index] = n;
2190 accessor.matrix->set(accessor.row(),
2191 accessor.column(),
2192 static_cast<Number>(*this));
2193 return *this;
2194 }
2195
2196
2197
2198 template <typename Number, typename MemorySpace>
2201 const Number n) const
2202 {
2203 (*accessor.value_cache)[accessor.a_index] += n;
2204 accessor.matrix->set(accessor.row(),
2205 accessor.column(),
2206 static_cast<Number>(*this));
2207 return *this;
2208 }
2209
2210
2211
2212 template <typename Number, typename MemorySpace>
2215 const Number n) const
2216 {
2217 (*accessor.value_cache)[accessor.a_index] -= n;
2218 accessor.matrix->set(accessor.row(),
2219 accessor.column(),
2220 static_cast<Number>(*this));
2221 return *this;
2222 }
2223
2224
2225 template <typename Number, typename MemorySpace>
2228 const Number n) const
2229 {
2230 (*accessor.value_cache)[accessor.a_index] *= n;
2231 accessor.matrix->set(accessor.row(),
2232 accessor.column(),
2233 static_cast<Number>(*this));
2234 return *this;
2235 }
2236
2237
2238 template <typename Number, typename MemorySpace>
2241 const Number n) const
2242 {
2243 (*accessor.value_cache)[accessor.a_index] /= n;
2244 accessor.matrix->set(accessor.row(),
2245 accessor.column(),
2246 static_cast<Number>(*this));
2247 return *this;
2248 }
2249
2250
2251 template <typename Number, typename MemorySpace>
2258
2259
2260 template <typename Number, typename MemorySpace>
2269
2270
2271
2272 template <typename Number, typename MemorySpace, bool Constness>
2274 MatrixType *matrix,
2275 const size_type row,
2276 const size_type index)
2277 : accessor(matrix, row, index)
2278 {}
2279
2280
2281 template <typename Number, typename MemorySpace, bool Constness>
2282 template <bool Other>
2285 : accessor(other.accessor)
2286 {}
2287
2288
2289
2290 template <typename Number, typename MemorySpace, bool Constness>
2293 {
2294 Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2295
2296 ++accessor.a_index;
2297
2298 // If at end of line: do one
2299 // step, then cycle until we
2300 // find a row with a nonzero
2301 // number of entries.
2302 if (accessor.a_index >= accessor.colnum_cache->size())
2303 {
2304 accessor.a_index = 0;
2305 ++accessor.a_row;
2306
2307 while (
2308 (accessor.a_row < accessor.matrix->m()) &&
2309 ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2310 (accessor.matrix->row_length(accessor.a_row) == 0)))
2311 ++accessor.a_row;
2312
2313 accessor.visit_present_row();
2314 }
2315 return *this;
2316 }
2317
2318
2319 template <typename Number, typename MemorySpace, bool Constness>
2322 {
2323 const Iterator<Number, MemorySpace, Constness> old_state = *this;
2324 ++(*this);
2325 return old_state;
2326 }
2327
2328
2329
2330 template <typename Number, typename MemorySpace, bool Constness>
2333 {
2334 return accessor;
2335 }
2336
2337
2338
2339 template <typename Number, typename MemorySpace, bool Constness>
2342 {
2343 return &accessor;
2344 }
2345
2346
2347
2348 template <typename Number, typename MemorySpace, bool Constness>
2349 template <bool OtherConstness>
2350 inline bool
2353 {
2354 return (accessor.a_row == other.accessor.a_row &&
2355 accessor.a_index == other.accessor.a_index);
2356 }
2357
2358
2359
2360 template <typename Number, typename MemorySpace, bool Constness>
2361 template <bool OtherConstness>
2362 inline bool
2365 {
2366 return !(*this == other);
2367 }
2368
2369
2370
2371 template <typename Number, typename MemorySpace, bool Constness>
2372 template <bool OtherConstness>
2373 inline bool
2376 {
2377 return (accessor.row() < other.accessor.row() ||
2378 (accessor.row() == other.accessor.row() &&
2379 accessor.index() < other.accessor.index()));
2380 }
2381
2382
2383 template <typename Number, typename MemorySpace, bool Constness>
2384 template <bool OtherConstness>
2385 inline bool
2388 {
2389 return (other < *this);
2390 }
2391
2392 } // namespace SparseMatrixIterators
2393 } // namespace TpetraWrappers
2394
2395} // namespace LinearAlgebra
2396
2398
2399#endif // DEAL_II_TRILINOS_WITH_TPETRA
2400
2401#endif // dealii_trilinos_tpetra_sparse_matrix_h
std::shared_ptr< std::vector<::types::signed_global_dof_index > > colnum_cache
AccessorBase(SparseMatrix< Number, MemorySpace > *matrix, const size_type row, const size_type index)
bool operator!=(const Iterator< Number, MemorySpace, OtherConstness > &) const
bool operator<(const Iterator< Number, MemorySpace, OtherConstness > &) const
bool operator>(const Iterator< Number, MemorySpace, OtherConstness > &) const
Iterator(MatrixType *matrix, const size_type row, const size_type index)
const Accessor< Number, MemorySpace, Constness > * operator->() const
typename Accessor< Number, MemorySpace, Constness >::MatrixType MatrixType
const Accessor< Number, MemorySpace, Constness > & operator*() const
bool operator==(const Iterator< Number, MemorySpace, OtherConstness > &) const
virtual ~SparseMatrix() override=default
Tpetra::CrsGraph< int, ::types::signed_global_dof_index, NodeType > GraphType
unsigned int row_length(const size_type row) const
SparseMatrix(const SparseMatrix< Number, MemorySpace > &)=delete
void clear_row(const size_type row, const Number new_diag_value=0)
SparseMatrix(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const MPI_Comm communicator=MPI_COMM_WORLD, const size_type n_max_entries_per_row=0)
void print(std::ostream &out, const bool print_detailed_trilinos_information=false) const
Tpetra::Vector< Number, int, ::types::signed_global_dof_index, NodeType > VectorType
SparseMatrix< Number, MemorySpace > & operator=(SparseMatrix< Number, MemorySpace > &&other) noexcept
void add(const size_type i, const size_type j, const Number value)
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< Number > &full_matrix, const bool elide_zero_values=false)
SparseMatrix(SparseMatrix< Number, MemorySpace > &&other) noexcept
void clear_rows(const ArrayView< const size_type > &rows, const Number new_diag_value=0)
void set(const std::vector< size_type > &indices, const FullMatrix< Number > &full_matrix, const bool elide_zero_values=false)
Number residual(Vector< Number, MemorySpace > &dst, const Vector< Number, MemorySpace > &x, const Vector< Number, MemorySpace > &b) const
std::pair< size_type, size_type > local_range() const
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< Number > &values, const bool elide_zero_values=false)
void Tvmult(InputVectorType &dst, const InputVectorType &src) const
Tpetra::CrsMatrix< Number, int, ::types::signed_global_dof_index, NodeType > MatrixType
Number matrix_scalar_product(const Vector< Number, MemorySpace > &u, const Vector< Number, MemorySpace > &v) const
void reinit(const SparsityPatternType &sparsity_pattern)
SparseMatrix & operator=(const double d)
Number element(const size_type i, const size_type j, const bool no_error) const
SparseMatrix(const size_type m, const size_type n, const unsigned int n_max_entries_per_row)
SparseMatrix & operator/=(const Number factor)
Teuchos::RCP< const MatrixType > trilinos_rcp() const
SparseMatrix & operator*=(const Number factor)
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const Number *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
SparseMatrix(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const MPI_Comm communicator, const std::vector< unsigned int > &n_entries_per_row)
SparseMatrix(const SparsityPattern< MemorySpace > &sparsity_pattern)
SparseMatrix(const IndexSet &parallel_partitioning, const MPI_Comm communicator=MPI_COMM_WORLD, const unsigned int n_max_entries_per_row=0)
Number diag_element(const size_type i) const
Number operator()(const size_type i, const size_type j) const
void vmult(InputVectorType &dst, const InputVectorType &src) const
SparseMatrix(const size_type m, const size_type n, const std::vector< unsigned int > &n_entries_per_row)
void set(const size_type i, const size_type j, const Number value)
void copy_from(const SparseMatrix< Number, MemorySpace > &source)
SparseMatrix(const IndexSet &parallel_partitioning, const MPI_Comm communicator, const std::vector< unsigned int > &n_entries_per_row)
std::enable_if_t< !std::is_same_v< SparsityPatternType, ::SparseMatrix< double > > > reinit(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator=MPI_COMM_WORLD, const bool exchange_data=false)
SparseMatrix< Number, MemorySpace > & operator=(const SparseMatrix< Number, MemorySpace > &)=delete
Tpetra::Map< int, ::types::signed_global_dof_index, NodeType > MapType
void Tvmult_add(InputVectorType &dst, const InputVectorType &src) const
void reinit(const SparsityPattern< MemorySpace > &sparsity_pattern)
Number el(const size_type i, const size_type j) const
void vmult_add(InputVectorType &dst, const InputVectorType &src) const
void compress(VectorOperation::values operation)
std::enable_if_t< !std::is_same_v< SparsityPatternType, ::SparseMatrix< double > > > reinit(const IndexSet &parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator=MPI_COMM_WORLD, const bool exchange_data=false)
void add(const Number factor, const SparseMatrix< Number, MemorySpace > &matrix)
void reinit(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const ::SparseMatrix< Number > &dealii_sparse_matrix, const MPI_Comm communicator=MPI_COMM_WORLD, const double drop_tolerance=1e-13, const bool copy_values=true, const ::SparsityPattern *use_this_sparsity=nullptr)
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const Number *values, const bool elide_zero_values=false)
Number matrix_norm_square(const Vector< Number, MemorySpace > &v) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcDomainMapMismatch()
static ::ExceptionBase & ExcColMapMismatch()
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:585
#define Assert(cond, exc)
static ::ExceptionBase & ExcMatrixNotCompressed()
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcSourceEqualsDestination()
STL namespace.
int signed_global_dof_index
Definition types.h:92
unsigned int global_dof_index
Definition types.h:81
typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::Iterator< Number, MemorySpace, Constness >::difference_type difference_type
typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::Iterator< Number, MemorySpace, Constness >::value_type value_type