deal.II version GIT relicensing-1982-gbb401394a0 2024-10-13 08:40:01+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
170
175
184 const size_type n,
185 const unsigned int n_max_entries_per_row);
186
195 const size_type n,
196 const std::vector<unsigned int> &n_entries_per_row);
197
203
208
214
220
225 virtual ~SparseMatrix() override = default;
226
242 template <typename SparsityPatternType>
243 void
244 reinit(const SparsityPatternType &sparsity_pattern);
245
255 void
256 reinit(const SparsityPattern<MemorySpace> &sparsity_pattern);
275 SparseMatrix(const IndexSet &parallel_partitioning,
276 const MPI_Comm communicator = MPI_COMM_WORLD,
277 const unsigned int n_max_entries_per_row = 0);
278
287 SparseMatrix(const IndexSet &parallel_partitioning,
288 const MPI_Comm communicator,
289 const std::vector<unsigned int> &n_entries_per_row);
290
305 SparseMatrix(const IndexSet &row_parallel_partitioning,
306 const IndexSet &col_parallel_partitioning,
307 const MPI_Comm communicator = MPI_COMM_WORLD,
308 const size_type n_max_entries_per_row = 0);
309
318 SparseMatrix(const IndexSet &row_parallel_partitioning,
319 const IndexSet &col_parallel_partitioning,
320 const MPI_Comm communicator,
321 const std::vector<unsigned int> &n_entries_per_row);
322
342 template <typename SparsityPatternType>
343 std::enable_if_t<
344 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
345 reinit(const IndexSet &parallel_partitioning,
346 const SparsityPatternType &sparsity_pattern,
347 const MPI_Comm communicator = MPI_COMM_WORLD,
348 const bool exchange_data = false);
349
362 template <typename SparsityPatternType>
363 std::enable_if_t<
364 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
365 reinit(const IndexSet &row_parallel_partitioning,
366 const IndexSet &col_parallel_partitioning,
367 const SparsityPatternType &sparsity_pattern,
368 const MPI_Comm communicator = MPI_COMM_WORLD,
369 const bool exchange_data = false);
370
387 void
388 reinit(const IndexSet &row_parallel_partitioning,
389 const IndexSet &col_parallel_partitioning,
390 const ::SparseMatrix<Number> &dealii_sparse_matrix,
391 const MPI_Comm communicator = MPI_COMM_WORLD,
392 const double drop_tolerance = 1e-13,
393 const bool copy_values = true,
394 const ::SparsityPattern *use_this_sparsity = nullptr);
395
406 m() const;
407
412 n() const;
413
414
423 unsigned int
424 local_size() const;
425
434 std::pair<size_type, size_type>
435 local_range() const;
436
441 bool
442 in_local_range(const size_type index) const;
443
448 size_t
450
454 unsigned int
455 row_length(const size_type row) const;
456
463 bool
465
487 operator=(const double d);
488
493 operator*=(const Number factor);
494
499 operator/=(const Number factor);
500
504 void
506
516 void
517 add(const size_type i, const size_type j, const Number value);
518
532 void
533 add(const size_type row,
534 const size_type n_cols,
535 const size_type *col_indices,
536 const Number *values,
537 const bool elide_zero_values = true,
538 const bool col_indices_are_sorted = false);
539
547 void
548 add(const Number factor, const SparseMatrix<Number, MemorySpace> &matrix);
549
571 void
572 set(const size_type i, const size_type j, const Number value);
573
606 void
607 set(const std::vector<size_type> &indices,
608 const FullMatrix<Number> &full_matrix,
609 const bool elide_zero_values = false);
610
616 void
617 set(const std::vector<size_type> &row_indices,
618 const std::vector<size_type> &col_indices,
619 const FullMatrix<Number> &full_matrix,
620 const bool elide_zero_values = false);
621
649 void
650 set(const size_type row,
651 const std::vector<size_type> &col_indices,
652 const std::vector<Number> &values,
653 const bool elide_zero_values = false);
654
682 void
683 set(const size_type row,
684 const size_type n_cols,
685 const size_type *col_indices,
686 const Number *values,
687 const bool elide_zero_values = false);
688
715 void
716 clear_row(const size_type row, const Number new_diag_value = 0);
717
738 void
740 const Number new_diag_value = 0);
741
749 void
751
766 Number
767 operator()(const size_type i, const size_type j) const;
768
785 Number
786 el(const size_type i, const size_type j) const;
787
793 Number
794 diag_element(const size_type i) const;
795
801 /*
802 * Matrix-vector multiplication: let <i>dst = M*src</i> with <i>M</i>
803 * being this matrix.
804 *
805 * Source and destination must not be the same vector.
806 *
807 * The vector @p dst has to be initialized with the same IndexSet that was
808 * used for the row indices of the matrix and the vector @p src has to be
809 * initialized with the same IndexSet that was used for the column indices
810 * of the matrix.
811 */
812 template <typename InputVectorType>
813 void
814 vmult(InputVectorType &dst, const InputVectorType &src) const;
815
816 /*
817 * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
818 * <i>M</i> being this matrix. This function does the same as vmult() but
819 * takes the transposed matrix.
820 *
821 * Source and destination must not be the same vector.
822 */
823 template <typename InputVectorType>
824 void
825 Tvmult(InputVectorType &dst, const InputVectorType &src) const;
826
833 template <typename InputVectorType>
834 void
835 vmult_add(InputVectorType &dst, const InputVectorType &src) const;
836
844 template <typename InputVectorType>
845 void
846 Tvmult_add(InputVectorType &dst, const InputVectorType &src) const;
847
869 Number
871
891 Number
893 const Vector<Number, MemorySpace> &v) const;
894
911 Number
914 const Vector<Number, MemorySpace> &b) const;
915
927 Number
929
943 void
944 print(std::ostream &out,
945 const bool print_detailed_trilinos_information = false) const;
946
978 void
980
987 void
989
998
1007
1017 Teuchos::RCP<const TpetraTypes::MatrixType<Number, MemorySpace>>
1019
1029 Teuchos::RCP<TpetraTypes::MatrixType<Number, MemorySpace>>
1042 IndexSet
1044
1050 IndexSet
1052
1079 begin() const;
1080
1084 iterator
1086
1092 end() const;
1093
1097 iterator
1099
1129 begin(const size_type r) const;
1130
1134 iterator
1136
1147 end(const size_type r) const;
1148
1152 iterator
1153 end(const size_type r);
1154
1165
1171 "You are attempting an operation on two vectors that "
1172 "are the same object, but the operation requires that the "
1173 "two objects are in fact different.");
1174
1175 /*
1176 * Exception
1177 */
1179 "The column partitioning of a matrix does not match "
1180 "the partitioning of a vector you are trying to "
1181 "multiply it with. Are you multiplying the "
1182 "matrix with a vector that has ghost elements?");
1183
1184 /*
1185 * Exception
1186 */
1188 "The row partitioning of a matrix does not match "
1189 "the partitioning of a vector you are trying to "
1190 "put the result of a matrix-vector product in. "
1191 "Are you trying to put the product of the "
1192 "matrix with a vector into a vector that has "
1193 "ghost elements?");
1194
1199 size_type,
1200 size_type,
1201 << "The entry with index <" << arg1 << ',' << arg2
1202 << "> does not exist.");
1203
1208 size_type,
1209 size_type,
1210 size_type,
1211 size_type,
1212 << "You tried to access element (" << arg1 << '/' << arg2
1213 << ')'
1214 << " of a distributed matrix, but only rows in range ["
1215 << arg3 << ',' << arg4
1216 << "] are stored locally and can be accessed.");
1217
1220 private:
1224 Number
1225 element(const size_type i, const size_type j, const bool no_error) const;
1226
1236 Teuchos::RCP<TpetraTypes::MapType<MemorySpace>> column_space_map;
1237
1243 Teuchos::RCP<TpetraTypes::MatrixType<Number, MemorySpace>> matrix;
1244
1250
1263 void
1265
1273 void
1275
1276 // To allow calling protected prepare_add() and prepare_set().
1277 friend class BlockMatrixBase<SparseMatrix<Number, MemorySpace>>;
1278 }; // class SparseMatrix
1279
1284 {
1289
1297 template <typename Number, typename MemorySpace>
1299 {
1300 public:
1305
1310 const size_type row,
1311 const size_type index);
1312
1316 size_type
1317 row() const;
1318
1322 size_type
1323 index() const;
1324
1328 size_type
1329 column() const;
1330
1331 protected:
1343
1348
1354 void
1356
1369 std::shared_ptr<std::vector<::types::signed_global_dof_index>>
1371
1375 std::shared_ptr<std::vector<Number>> value_cache;
1376
1377 private:
1378 friend class Iterator<Number, MemorySpace, false>;
1379 friend class Iterator<Number, MemorySpace, true>;
1380 };
1381
1392 template <typename Number, typename MemorySpace, bool Constness>
1393 class Accessor : public AccessorBase<Number, MemorySpace>
1394 {
1398 Number
1399 value() const;
1400
1404 Number &
1406 };
1407
1411 template <typename Number, typename MemorySpace>
1412 class Accessor<Number, MemorySpace, true>
1413 : public AccessorBase<Number, MemorySpace>
1414 {
1415 public:
1421
1426
1432 const size_type row,
1433 const size_type index);
1434
1439 template <bool Other>
1441
1445 Number
1446 value() const;
1447 };
1448
1452 template <typename Number, typename MemorySpace>
1453 class Accessor<Number, MemorySpace, false>
1454 : public AccessorBase<Number, MemorySpace>
1455 {
1456 class Reference
1457 {
1458 public:
1462 Reference(const Accessor<Number, MemorySpace, false> &accessor);
1463
1467 operator Number() const;
1468
1472 const Reference &
1473 operator=(const Number n) const;
1474
1478 const Reference &
1479 operator+=(const Number n) const;
1480
1484 const Reference &
1485 operator-=(const Number n) const;
1486
1490 const Reference &
1491 operator*=(const Number n) const;
1492
1496 const Reference &
1497 operator/=(const Number n) const;
1498
1499 private:
1505 };
1506
1507 public:
1513
1518
1524 const size_type row,
1525 const size_type index);
1526
1530 Reference
1531 value() const;
1532
1533 private:
1534 // Make Reference object a friend.
1535 friend class Reference;
1536 };
1537
1551 template <typename Number, typename MemorySpace, bool Constness>
1553 {
1554 public:
1559
1565
1570 using value_type = Number;
1571
1578
1583 Iterator(MatrixType *matrix,
1584 const size_type row,
1585 const size_type index);
1586
1590 template <bool Other>
1592
1597 operator++();
1598
1603 operator++(int);
1604
1609 operator*() const;
1610
1615 operator->() const;
1616
1621 template <bool OtherConstness>
1622 bool
1624
1628 template <bool OtherConstness>
1629 bool
1631
1637 template <bool OtherConstness>
1638 bool
1640
1644 template <bool OtherConstness>
1645 bool
1647
1652 size_type,
1653 size_type,
1654 << "Attempt to access element " << arg2 << " of row "
1655 << arg1 << " which doesn't have that many elements.");
1656
1657 private:
1662
1663 friend class Iterator<Number, MemorySpace, true>;
1664 friend class Iterator<Number, MemorySpace, false>;
1665 };
1666 } // namespace SparseMatrixIterators
1667
1668 } // namespace TpetraWrappers
1669
1670} // namespace LinearAlgebra
1671
1672
1674
1675namespace std
1676{
1677 template <typename Number, typename MemorySpace, bool Constness>
1680 Iterator<Number, MemorySpace, Constness>>
1681 {
1682 using iterator_category = forward_iterator_tag;
1684 typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
1685 Iterator<Number, MemorySpace, Constness>::value_type;
1687 typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
1688 Iterator<Number, MemorySpace, Constness>::difference_type;
1689 };
1690} // namespace std
1691
1692/* ------------------------- Inline functions ---------------------- */
1693
1694
1696
1697namespace LinearAlgebra
1698{
1699
1700 namespace TpetraWrappers
1701 {
1702 template <typename Number, typename MemorySpace>
1703 inline void
1705 const size_type j,
1706 const Number value)
1707 {
1708 set(i, 1, &j, &value, false);
1709 }
1710
1711
1712
1713 template <typename Number, typename MemorySpace>
1714 inline void
1716 const size_type j,
1717 const Number value)
1718 {
1719 add(i, 1, &j, &value, false);
1720 }
1721
1722
1723
1724 template <typename Number, typename MemorySpace>
1725 inline Number
1729 const Vector<Number, MemorySpace> &b) const
1730 {
1731 vmult(dst, x);
1732 dst -= b;
1733 dst *= -1.;
1734
1735 return dst.l2_norm();
1736 }
1737
1738
1739
1740 template <typename Number, typename MemorySpace>
1743 {
1744 return matrix->getRowMap()->getGlobalNumElements();
1745 }
1746
1747
1748
1749 template <typename Number, typename MemorySpace>
1752 {
1753 // If the matrix structure has not been fixed (i.e., we did not have a
1754 // sparsity pattern), it does not know about the number of columns, so we
1755 // must always take this from the additional column space map
1756 Assert(column_space_map.get() != nullptr, ExcInternalError());
1757 return column_space_map->getGlobalNumElements();
1758 }
1759
1760
1761
1762 template <typename Number, typename MemorySpace>
1763 inline bool
1765 {
1766 return compressed;
1767 }
1768
1769
1770
1771 template <typename Number, typename MemorySpace>
1774 {
1775 return begin(0);
1776 }
1777
1778
1779
1780 template <typename Number, typename MemorySpace>
1783 {
1784 return const_iterator(this, m(), 0);
1785 }
1786
1787
1788
1789 template <typename Number, typename MemorySpace>
1792 {
1793 AssertIndexRange(r, m());
1794 if (in_local_range(r) && (row_length(r) > 0))
1795 return const_iterator(this, r, 0);
1796 else
1797 return end(r);
1798 }
1799
1800
1801 template <typename Number, typename MemorySpace>
1804 {
1805 AssertIndexRange(r, m());
1806
1807 // place the iterator on the first entry
1808 // past this line, or at the end of the
1809 // matrix
1810 for (size_type i = r + 1; i < m(); ++i)
1811 if (in_local_range(i) && (row_length(i) > 0))
1812 return const_iterator(this, i, 0);
1813
1814 // if there is no such line, then take the
1815 // end iterator of the matrix
1816 return end();
1817 }
1818
1819
1820
1821 template <typename Number, typename MemorySpace>
1824 {
1825 return begin(0);
1826 }
1827
1828
1829
1830 template <typename Number, typename MemorySpace>
1833 {
1834 return iterator(this, m(), 0);
1835 }
1836
1837
1838
1839 template <typename Number, typename MemorySpace>
1842 {
1843 AssertIndexRange(r, m());
1844 if (in_local_range(r) && (row_length(r) > 0))
1845 return iterator(this, r, 0);
1846 else
1847 return end(r);
1848 }
1849
1850
1851
1852 template <typename Number, typename MemorySpace>
1855 {
1856 AssertIndexRange(r, m());
1857
1858 // place the iterator on the first entry
1859 // past this line, or at the end of the
1860 // matrix
1861 for (size_type i = r + 1; i < m(); ++i)
1862 if (in_local_range(i) && (row_length(i) > 0))
1863 return iterator(this, i, 0);
1864
1865 // if there is no such line, then take the
1866 // end iterator of the matrix
1867 return end();
1868 }
1869
1870
1871
1872 template <typename Number, typename MemorySpace>
1873 inline bool
1875 const size_type index) const
1876 {
1877 const size_type begin = matrix->getRowMap()->getMinGlobalIndex();
1878 const size_type end = matrix->getRowMap()->getMaxGlobalIndex() + 1;
1879
1880 return ((index >= begin) && (index < end));
1881 }
1882
1883
1884
1885 template <typename Number, typename MemorySpace>
1886 unsigned int
1888 {
1889 auto n_entries = matrix->getNumEntriesInGlobalRow(row);
1890 Assert(n_entries !=
1891 Teuchos::OrdinalTraits<decltype(n_entries)>::invalid(),
1892 ExcAccessToNonlocalRow(row));
1893
1894 return n_entries;
1895 }
1896
1897
1898
1899 template <typename Number, typename MemorySpace>
1900 inline void
1902 {
1903 // nothing to do here
1904 }
1905
1906
1907
1908 template <typename Number, typename MemorySpace>
1909 inline void
1911 {
1912 // nothing to do here
1913 }
1914
1915
1916
1917 template <typename Number, typename MemorySpace>
1920 {
1921 return *matrix;
1922 }
1923
1924
1925
1926 template <typename Number, typename MemorySpace>
1929 {
1930 return *matrix;
1931 }
1932
1933
1934
1935 template <typename Number, typename MemorySpace>
1936 inline Teuchos::RCP<const TpetraTypes::MatrixType<Number, MemorySpace>>
1938 {
1939 return matrix.getConst();
1940 }
1941
1942
1943
1944 template <typename Number, typename MemorySpace>
1945 inline Teuchos::RCP<TpetraTypes::MatrixType<Number, MemorySpace>>
1947 {
1948 return matrix;
1949 }
1950
1951
1952
1953 template <typename Number, typename MemorySpace>
1954 inline IndexSet
1956 {
1957 return IndexSet(matrix->getDomainMap());
1958 }
1959
1960
1961
1962 template <typename Number, typename MemorySpace>
1963 inline IndexSet
1965 {
1966 return IndexSet(matrix->getRangeMap());
1967 }
1968
1969
1970 namespace SparseMatrixIterators
1971 {
1972 template <typename Number, typename MemorySpace>
1975 size_type row,
1976 size_type index)
1977 : matrix(matrix)
1978 , a_row(row)
1979 , a_index(index)
1980 {
1982 }
1983
1984
1985 template <typename Number, typename MemorySpace>
1988 {
1989 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1990 return a_row;
1991 }
1992
1993
1994 template <typename Number, typename MemorySpace>
1997 {
1998 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
1999 return (*colnum_cache)[a_index];
2000 }
2001
2002
2003 template <typename Number, typename MemorySpace>
2006 {
2007 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2008 return a_index;
2009 }
2010
2011
2012
2013 template <typename Number, typename MemorySpace>
2014 void
2016 {
2017 // if we are asked to visit the past-the-end line, then simply
2018 // release all our caches and go on with life.
2019 //
2020 // do the same if the row we're supposed to visit is not locally
2021 // owned. this is simply going to make non-locally owned rows
2022 // look like they're empty
2023 if ((this->a_row == matrix->m()) ||
2024 (matrix->in_local_range(this->a_row) == false))
2025 {
2026 colnum_cache.reset();
2027 value_cache.reset();
2028
2029 return;
2030 }
2031
2032 // get a representation of the present row
2033 size_t ncols;
2035 matrix->row_length(this->a_row);
2036 if (value_cache.get() == nullptr)
2037 {
2038 value_cache = std::make_shared<std::vector<Number>>(colnums);
2039 colnum_cache = std::make_shared<
2040 std::vector<::types::signed_global_dof_index>>(colnums);
2041 }
2042 else
2043 {
2044 value_cache->resize(colnums);
2045 colnum_cache->resize(colnums);
2046 }
2047
2048# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
2050 nonconst_global_inds_host_view_type col_indices(colnum_cache->data(),
2051 colnums);
2053 nonconst_values_host_view_type values(value_cache->data(), colnums);
2054# else
2055 Teuchos::ArrayView<::types::signed_global_dof_index> col_indices(
2056 *colnum_cache);
2057 Teuchos::ArrayView<Number> values(*value_cache);
2058# endif
2059 matrix->trilinos_matrix().getGlobalRowCopy(this->a_row,
2060 col_indices,
2061 values,
2062 ncols);
2063
2064 AssertDimension(ncols, colnums);
2065
2066 // copy it into our caches if the
2067 // line isn't empty. if it is, then
2068 // we've done something wrong, since
2069 // we shouldn't have initialized an
2070 // iterator for an empty line (what
2071 // would it point to?)
2072 }
2073
2074
2075
2076 template <typename Number, typename MemorySpace>
2078 MatrixType *matrix,
2079 const size_type row,
2080 const size_type index)
2081 : AccessorBase<Number, MemorySpace>(
2082 const_cast<SparseMatrix<Number, MemorySpace> *>(matrix),
2083 row,
2084 index)
2085 {}
2086
2087
2088
2089 template <typename Number, typename MemorySpace>
2090 template <bool Other>
2095
2096
2097
2098 template <typename Number, typename MemorySpace>
2099 inline Number
2108
2109
2110
2111 template <typename Number, typename MemorySpace>
2114 : accessor(const_cast<Accessor<Number, MemorySpace, false> &>(acc))
2115 {}
2116
2117
2118
2119 template <typename Number, typename MemorySpace>
2121 const
2122 {
2123 return (*accessor.value_cache)[accessor.a_index];
2124 }
2125
2126
2127
2128 template <typename Number, typename MemorySpace>
2131 const Number n) const
2132 {
2133 (*accessor.value_cache)[accessor.a_index] = n;
2134 accessor.matrix->set(accessor.row(),
2135 accessor.column(),
2136 static_cast<Number>(*this));
2137 return *this;
2138 }
2139
2140
2141
2142 template <typename Number, typename MemorySpace>
2145 const Number n) const
2146 {
2147 (*accessor.value_cache)[accessor.a_index] += n;
2148 accessor.matrix->set(accessor.row(),
2149 accessor.column(),
2150 static_cast<Number>(*this));
2151 return *this;
2152 }
2153
2154
2155
2156 template <typename Number, typename MemorySpace>
2159 const Number n) const
2160 {
2161 (*accessor.value_cache)[accessor.a_index] -= n;
2162 accessor.matrix->set(accessor.row(),
2163 accessor.column(),
2164 static_cast<Number>(*this));
2165 return *this;
2166 }
2167
2168
2169 template <typename Number, typename MemorySpace>
2172 const Number n) const
2173 {
2174 (*accessor.value_cache)[accessor.a_index] *= n;
2175 accessor.matrix->set(accessor.row(),
2176 accessor.column(),
2177 static_cast<Number>(*this));
2178 return *this;
2179 }
2180
2181
2182 template <typename Number, typename MemorySpace>
2185 const Number n) const
2186 {
2187 (*accessor.value_cache)[accessor.a_index] /= n;
2188 accessor.matrix->set(accessor.row(),
2189 accessor.column(),
2190 static_cast<Number>(*this));
2191 return *this;
2192 }
2193
2194
2195 template <typename Number, typename MemorySpace>
2202
2203
2204 template <typename Number, typename MemorySpace>
2213
2214
2215
2216 template <typename Number, typename MemorySpace, bool Constness>
2218 MatrixType *matrix,
2219 const size_type row,
2220 const size_type index)
2221 : accessor(matrix, row, index)
2222 {}
2223
2224
2225 template <typename Number, typename MemorySpace, bool Constness>
2226 template <bool Other>
2229 : accessor(other.accessor)
2230 {}
2231
2232
2233
2234 template <typename Number, typename MemorySpace, bool Constness>
2237 {
2238 Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2239
2240 ++accessor.a_index;
2241
2242 // If at end of line: do one
2243 // step, then cycle until we
2244 // find a row with a nonzero
2245 // number of entries.
2246 if (accessor.a_index >= accessor.colnum_cache->size())
2247 {
2248 accessor.a_index = 0;
2249 ++accessor.a_row;
2250
2251 while (
2252 (accessor.a_row < accessor.matrix->m()) &&
2253 ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2254 (accessor.matrix->row_length(accessor.a_row) == 0)))
2255 ++accessor.a_row;
2256
2257 accessor.visit_present_row();
2258 }
2259 return *this;
2260 }
2261
2262
2263 template <typename Number, typename MemorySpace, bool Constness>
2266 {
2267 const Iterator<Number, MemorySpace, Constness> old_state = *this;
2268 ++(*this);
2269 return old_state;
2270 }
2271
2272
2273
2274 template <typename Number, typename MemorySpace, bool Constness>
2277 {
2278 return accessor;
2279 }
2280
2281
2282
2283 template <typename Number, typename MemorySpace, bool Constness>
2286 {
2287 return &accessor;
2288 }
2289
2290
2291
2292 template <typename Number, typename MemorySpace, bool Constness>
2293 template <bool OtherConstness>
2294 inline bool
2297 {
2298 return (accessor.a_row == other.accessor.a_row &&
2299 accessor.a_index == other.accessor.a_index);
2300 }
2301
2302
2303
2304 template <typename Number, typename MemorySpace, bool Constness>
2305 template <bool OtherConstness>
2306 inline bool
2309 {
2310 return !(*this == other);
2311 }
2312
2313
2314
2315 template <typename Number, typename MemorySpace, bool Constness>
2316 template <bool OtherConstness>
2317 inline bool
2320 {
2321 return (accessor.row() < other.accessor.row() ||
2322 (accessor.row() == other.accessor.row() &&
2323 accessor.index() < other.accessor.index()));
2324 }
2325
2326
2327 template <typename Number, typename MemorySpace, bool Constness>
2328 template <bool OtherConstness>
2329 inline bool
2332 {
2333 return (other < *this);
2334 }
2335
2336 } // namespace SparseMatrixIterators
2337 } // namespace TpetraWrappers
2338
2339} // namespace LinearAlgebra
2340
2342
2343#endif // DEAL_II_TRILINOS_WITH_TPETRA
2344
2345#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
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
SparseMatrix< Number, MemorySpace > & operator=(SparseMatrix< Number, MemorySpace > &&other) noexcept
TpetraTypes::MatrixType< Number, MemorySpace > & trilinos_matrix()
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
Number matrix_scalar_product(const Vector< Number, MemorySpace > &u, const Vector< Number, MemorySpace > &v) const
Teuchos::RCP< const TpetraTypes::MatrixType< Number, MemorySpace > > trilinos_rcp() const
Teuchos::RCP< TpetraTypes::MatrixType< Number, MemorySpace > > matrix
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< TpetraTypes::MapType< MemorySpace > > column_space_map
SparseMatrix & operator*=(const Number factor)
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
Teuchos::RCP< TpetraTypes::MatrixType< Number, MemorySpace > > trilinos_rcp()
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
const TpetraTypes::MatrixType< Number, MemorySpace > & trilinos_matrix() 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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
static ::ExceptionBase & ExcDomainMapMismatch()
static ::ExceptionBase & ExcColMapMismatch()
#define DeclException0(Exception0)
Definition exceptions.h:466
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:580
#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:534
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:489
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:511
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcSourceEqualsDestination()
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
STL namespace.
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