Reference documentation for deal.II version GIT relicensing-136-gb80d0be4af 2024-03-18 08: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_sparse_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 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_sparse_matrix_h
16# define dealii_trilinos_sparse_matrix_h
17
18
19# include <deal.II/base/config.h>
20
21# ifdef DEAL_II_WITH_TRILINOS
22
26
35
36# include <Epetra_Comm.h>
37# include <Epetra_CrsGraph.h>
38# include <Epetra_Export.h>
39# include <Epetra_FECrsMatrix.h>
40# include <Epetra_Map.h>
41# include <Epetra_MpiComm.h>
42# include <Epetra_MultiVector.h>
43# include <Epetra_Operator.h>
44
45# include <cmath>
46# include <iterator>
47# include <memory>
48# include <type_traits>
49# include <vector>
50
52
53// forward declarations
54# ifndef DOXYGEN
55template <typename MatrixType>
56class BlockMatrixBase;
57
58template <typename number>
59class SparseMatrix;
60class SparsityPattern;
62
63namespace TrilinosWrappers
64{
65 class SparseMatrix;
66 class SparsityPattern;
67
68 namespace SparseMatrixIterators
69 {
70 template <bool Constness>
71 class Iterator;
72 }
73} // namespace TrilinosWrappers
74# endif
75
76namespace TrilinosWrappers
77{
82 {
87
92 std::size_t,
93 std::size_t,
94 std::size_t,
95 << "You tried to access row " << arg1
96 << " of a distributed sparsity pattern, "
97 << " but only rows " << arg2 << " through " << arg3
98 << " are stored locally and can be accessed.");
99
108 {
109 public:
114
119 const size_type row,
120 const size_type index);
121
126 row() const;
127
132 index() const;
133
138 column() const;
139
140 protected:
152
157
163 void
165
178 std::shared_ptr<std::vector<size_type>> colnum_cache;
179
183 std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
184 };
185
196 template <bool Constess>
197 class Accessor : public AccessorBase
198 {
203 value() const;
204
210 };
211
212
213
217 template <>
218 class Accessor<true> : public AccessorBase
219 {
220 public:
226
232
237 template <bool Other>
239
244 value() const;
245
246 private:
247 // Make iterator class a friend.
248 template <bool>
249 friend class Iterator;
250 };
251
255 template <>
256 class Accessor<false> : public AccessorBase
257 {
258 class Reference
259 {
260 public:
264 Reference(const Accessor<false> &accessor);
265
269 operator TrilinosScalar() const;
270
274 const Reference &
275 operator=(const TrilinosScalar n) const;
276
280 const Reference &
282
286 const Reference &
288
292 const Reference &
294
298 const Reference &
300
301 private:
307 };
308
309 public:
315
321
325 Reference
326 value() const;
327
328 private:
329 // Make iterator class a friend.
330 template <bool>
331 friend class Iterator;
332
333 // Make Reference object a friend.
334 friend class Reference;
335 };
336
350 template <bool Constness>
352 {
353 public:
358
364
370
376
381 Iterator(MatrixType *matrix, const size_type row, const size_type index);
382
386 template <bool Other>
388
394
400
404 const Accessor<Constness> &
405 operator*() const;
406
410 const Accessor<Constness> *
411 operator->() const;
412
417 template <bool OtherConstness>
418 bool
420
424 template <bool OtherConstness>
425 bool
427
433 template <bool OtherConstness>
434 bool
435 operator<(const Iterator<OtherConstness> &) const;
436
440 template <bool OtherConstness>
441 bool
443
448 size_type,
449 size_type,
450 << "Attempt to access element " << arg2 << " of row "
451 << arg1 << " which doesn't have that many elements.");
452
453 private:
458
459 template <bool Other>
460 friend class Iterator;
461 };
462
463 } // namespace SparseMatrixIterators
464} // namespace TrilinosWrappers
465
467
468namespace std
469{
470 template <bool Constness>
473 {
474 using iterator_category = forward_iterator_tag;
476 typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
477 Constness>::value_type;
479 typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
480 Constness>::difference_type;
481 };
482} // namespace std
483
485
486
487namespace TrilinosWrappers
488{
550 {
551 public:
556
561 std::size_t,
562 << "You tried to access row " << arg1
563 << " of a non-contiguous locally owned row set."
564 << " The row " << arg1
565 << " is not stored locally and can't be accessed.");
566
574 struct Traits
575 {
580 static const bool zero_addition_can_be_elided = true;
581 };
582
587
592
597
605 SparseMatrix();
606
615 const size_type n,
616 const unsigned int n_max_entries_per_row);
617
626 const size_type n,
627 const std::vector<unsigned int> &n_entries_per_row);
628
632 SparseMatrix(const SparsityPattern &InputSparsityPattern);
633
638 SparseMatrix(SparseMatrix &&other) noexcept;
639
643 SparseMatrix(const SparseMatrix &) = delete;
644
649 operator=(const SparseMatrix &) = delete;
650
654 virtual ~SparseMatrix() override = default;
655
671 template <typename SparsityPatternType>
672 void
673 reinit(const SparsityPatternType &sparsity_pattern);
674
687 void
688 reinit(const SparsityPattern &sparsity_pattern);
689
698 void
699 reinit(const SparseMatrix &sparse_matrix);
700
721 template <typename number>
722 void
723 reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
724 const double drop_tolerance = 1e-13,
725 const bool copy_values = true,
726 const ::SparsityPattern *use_this_sparsity = nullptr);
727
733 void
734 reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
753 SparseMatrix(const IndexSet &parallel_partitioning,
754 const MPI_Comm communicator = MPI_COMM_WORLD,
755 const unsigned int n_max_entries_per_row = 0);
756
764 SparseMatrix(const IndexSet &parallel_partitioning,
765 const MPI_Comm communicator,
766 const std::vector<unsigned int> &n_entries_per_row);
767
782 SparseMatrix(const IndexSet &row_parallel_partitioning,
783 const IndexSet &col_parallel_partitioning,
784 const MPI_Comm communicator = MPI_COMM_WORLD,
785 const size_type n_max_entries_per_row = 0);
786
801 SparseMatrix(const IndexSet &row_parallel_partitioning,
802 const IndexSet &col_parallel_partitioning,
803 const MPI_Comm communicator,
804 const std::vector<unsigned int> &n_entries_per_row);
805
826 template <typename SparsityPatternType>
827 void
828 reinit(const IndexSet &parallel_partitioning,
829 const SparsityPatternType &sparsity_pattern,
830 const MPI_Comm communicator = MPI_COMM_WORLD,
831 const bool exchange_data = false);
832
845 template <typename SparsityPatternType>
846 std::enable_if_t<
847 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
848 reinit(const IndexSet &row_parallel_partitioning,
849 const IndexSet &col_parallel_partitioning,
850 const SparsityPatternType &sparsity_pattern,
851 const MPI_Comm communicator = MPI_COMM_WORLD,
852 const bool exchange_data = false);
853
870 template <typename number>
871 void
872 reinit(const IndexSet &parallel_partitioning,
873 const ::SparseMatrix<number> &dealii_sparse_matrix,
874 const MPI_Comm communicator = MPI_COMM_WORLD,
875 const double drop_tolerance = 1e-13,
876 const bool copy_values = true,
877 const ::SparsityPattern *use_this_sparsity = nullptr);
878
892 template <typename number>
893 void
894 reinit(const IndexSet &row_parallel_partitioning,
895 const IndexSet &col_parallel_partitioning,
896 const ::SparseMatrix<number> &dealii_sparse_matrix,
897 const MPI_Comm communicator = MPI_COMM_WORLD,
898 const double drop_tolerance = 1e-13,
899 const bool copy_values = true,
900 const ::SparsityPattern *use_this_sparsity = nullptr);
911 m() const;
912
917 n() const;
918
927 unsigned int
928 local_size() const;
929
938 std::pair<size_type, size_type>
939 local_range() const;
940
945 bool
946 in_local_range(const size_type index) const;
947
952 std::uint64_t
954
958 unsigned int
959 row_length(const size_type row) const;
960
967 bool
969
976 memory_consumption() const;
977
982 get_mpi_communicator() const;
983
1000 operator=(const double d);
1001
1009 void
1010 clear();
1011
1039 void
1041
1063 void
1064 set(const size_type i, const size_type j, const TrilinosScalar value);
1065
1098 void
1099 set(const std::vector<size_type> &indices,
1100 const FullMatrix<TrilinosScalar> &full_matrix,
1101 const bool elide_zero_values = false);
1102
1108 void
1109 set(const std::vector<size_type> &row_indices,
1110 const std::vector<size_type> &col_indices,
1111 const FullMatrix<TrilinosScalar> &full_matrix,
1112 const bool elide_zero_values = false);
1113
1141 void
1142 set(const size_type row,
1143 const std::vector<size_type> &col_indices,
1144 const std::vector<TrilinosScalar> &values,
1145 const bool elide_zero_values = false);
1146
1174 template <typename Number>
1175 void
1176 set(const size_type row,
1177 const size_type n_cols,
1178 const size_type *col_indices,
1179 const Number *values,
1180 const bool elide_zero_values = false);
1181
1191 void
1192 add(const size_type i, const size_type j, const TrilinosScalar value);
1193
1212 void
1213 add(const std::vector<size_type> &indices,
1214 const FullMatrix<TrilinosScalar> &full_matrix,
1215 const bool elide_zero_values = true);
1216
1222 void
1223 add(const std::vector<size_type> &row_indices,
1224 const std::vector<size_type> &col_indices,
1225 const FullMatrix<TrilinosScalar> &full_matrix,
1226 const bool elide_zero_values = true);
1227
1241 void
1242 add(const size_type row,
1243 const std::vector<size_type> &col_indices,
1244 const std::vector<TrilinosScalar> &values,
1245 const bool elide_zero_values = true);
1246
1260 void
1261 add(const size_type row,
1262 const size_type n_cols,
1263 const size_type *col_indices,
1264 const TrilinosScalar *values,
1265 const bool elide_zero_values = true,
1266 const bool col_indices_are_sorted = false);
1267
1271 SparseMatrix &
1272 operator*=(const TrilinosScalar factor);
1273
1277 SparseMatrix &
1278 operator/=(const TrilinosScalar factor);
1279
1283 void
1284 copy_from(const SparseMatrix &source);
1285
1293 void
1294 add(const TrilinosScalar factor, const SparseMatrix &matrix);
1295
1322 void
1323 clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1324
1345 void
1347 const TrilinosScalar new_diag_value = 0);
1348
1358 void
1359 transpose();
1360
1376 operator()(const size_type i, const size_type j) const;
1377
1395 el(const size_type i, const size_type j) const;
1396
1404 diag_element(const size_type i) const;
1405
1442 template <typename VectorType>
1443 std::enable_if_t<
1444 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1445 vmult(VectorType &dst, const VectorType &src) const;
1446
1453 template <typename VectorType>
1454 std::enable_if_t<
1455 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1456 vmult(VectorType &dst, const VectorType &src) const;
1457
1472 template <typename VectorType>
1473 std::enable_if_t<
1474 std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1475 Tvmult(VectorType &dst, const VectorType &src) const;
1476
1483 template <typename VectorType>
1484 std::enable_if_t<
1485 !std::is_same_v<typename VectorType::value_type, TrilinosScalar>>
1486 Tvmult(VectorType &dst, const VectorType &src) const;
1487
1497 template <typename VectorType>
1498 void
1499 vmult_add(VectorType &dst, const VectorType &src) const;
1500
1511 template <typename VectorType>
1512 void
1513 Tvmult_add(VectorType &dst, const VectorType &src) const;
1514
1537 matrix_norm_square(const MPI::Vector &v) const;
1538
1559 matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1560
1577 template <typename VectorType>
1579 residual(VectorType &dst, const VectorType &x, const VectorType &b) const;
1580
1595 void
1597 const SparseMatrix &B,
1598 const MPI::Vector &V = MPI::Vector()) const;
1599
1600
1617 void
1619 const SparseMatrix &B,
1620 const MPI::Vector &V = MPI::Vector()) const;
1621
1636 l1_norm() const;
1637
1647 linfty_norm() const;
1648
1654 frobenius_norm() const;
1655
1666 const Epetra_CrsMatrix &
1668
1673 const Epetra_CrsGraph &
1675
1687 IndexSet
1689
1695 IndexSet
1697
1724 begin() const;
1725
1729 iterator
1731
1737 end() const;
1738
1742 iterator
1744
1774 begin(const size_type r) const;
1775
1779 iterator
1781
1792 end(const size_type r) const;
1793
1797 iterator
1798 end(const size_type r);
1799
1811 void
1812 write_ascii();
1813
1821 void
1822 print(std::ostream &out,
1823 const bool write_extended_trilinos_info = false) const;
1824
1835 int,
1836 << "An error with error number " << arg1
1837 << " occurred while calling a Trilinos function. "
1838 "\n\n"
1839 "For historical reasons, many Trilinos functions express "
1840 "errors by returning specific integer values to indicate "
1841 "certain errors. Unfortunately, different Trilinos functions "
1842 "often use the same integer values for different kinds of "
1843 "errors, and in most cases it is also not documented what "
1844 "each error code actually means. As a consequence, it is often "
1845 "difficult to say what a particular error (in this case, "
1846 "the error with integer code '"
1847 << arg1
1848 << "') represents and how one should fix a code to avoid it. "
1849 "The best one can often do is to look up the call stack to "
1850 "see which deal.II function generated the error, and which "
1851 "Trilinos function the error code had originated from; "
1852 "then look up the Trilinos source code of that function (for "
1853 "example on github) to see what code path set that error "
1854 "code. Short of going through all of that, the only other "
1855 "option is to guess the cause of the error from "
1856 "the context in which the error appeared.");
1857
1858
1863 size_type,
1864 size_type,
1865 << "The entry with index <" << arg1 << ',' << arg2
1866 << "> does not exist.");
1867
1872 "You are attempting an operation on two vectors that "
1873 "are the same object, but the operation requires that the "
1874 "two objects are in fact different.");
1875
1880
1885 size_type,
1886 size_type,
1887 size_type,
1888 size_type,
1889 << "You tried to access element (" << arg1 << '/' << arg2
1890 << ')'
1891 << " of a distributed matrix, but only rows in range ["
1892 << arg3 << ',' << arg4
1893 << "] are stored locally and can be accessed.");
1894
1899 size_type,
1900 size_type,
1901 << "You tried to access element (" << arg1 << '/' << arg2
1902 << ')' << " of a sparse matrix, but it appears to not"
1903 << " exist in the Trilinos sparsity pattern.");
1908 protected:
1909 private:
1914 std::unique_ptr<Epetra_Map> column_space_map;
1915
1921 std::unique_ptr<Epetra_FECrsMatrix> matrix;
1922
1928 std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1929
1933 std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1934
1946 Epetra_CombineMode last_action;
1947
1953
1966 void
1968
1976 void
1978
1979 // To allow calling protected prepare_add() and prepare_set().
1980 friend class BlockMatrixBase<SparseMatrix>;
1981 };
1982
1983
1984
1985 // forwards declarations
1986 class SolverBase;
1987 class PreconditionBase;
1988
1989 namespace internal
1990 {
1991 inline void
1992 check_vector_map_equality(const Epetra_CrsMatrix &mtrx,
1993 const Epetra_MultiVector &src,
1994 const Epetra_MultiVector &dst,
1995 const bool transpose)
1996 {
1997 if (transpose == false)
1998 {
1999 Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
2000 ExcMessage(
2001 "Column map of matrix does not fit with vector map!"));
2002 Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
2003 ExcMessage("Row map of matrix does not fit with vector map!"));
2004 }
2005 else
2006 {
2007 Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2008 ExcMessage(
2009 "Column map of matrix does not fit with vector map!"));
2010 Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2011 ExcMessage("Row map of matrix does not fit with vector map!"));
2012 }
2013 (void)mtrx; // removes -Wunused-variable in optimized mode
2014 (void)src;
2015 (void)dst;
2016 }
2017
2018 inline void
2020 const Epetra_MultiVector &src,
2021 const Epetra_MultiVector &dst,
2022 const bool transpose)
2023 {
2024 if (transpose == false)
2025 {
2026 Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2027 ExcMessage(
2028 "Column map of operator does not fit with vector map!"));
2029 Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2030 ExcMessage(
2031 "Row map of operator does not fit with vector map!"));
2032 }
2033 else
2034 {
2035 Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2036 ExcMessage(
2037 "Column map of operator does not fit with vector map!"));
2038 Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2039 ExcMessage(
2040 "Row map of operator does not fit with vector map!"));
2041 }
2042 (void)op; // removes -Wunused-variable in optimized mode
2043 (void)src;
2044 (void)dst;
2045 }
2046
2047 namespace LinearOperatorImplementation
2048 {
2069 {
2070 public:
2074 using VectorType = Epetra_MultiVector;
2075
2080
2085
2099
2103 TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2104 const TrilinosWrappers::SparseMatrix &matrix);
2105
2109 TrilinosPayload(const TrilinosPayload &payload_exemplar,
2110 const TrilinosWrappers::SparseMatrix &matrix);
2111
2116 const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2117 const TrilinosWrappers::PreconditionBase &preconditioner);
2118
2123 const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2124 const TrilinosWrappers::PreconditionBase &preconditioner);
2125
2130 const TrilinosPayload &payload_exemplar,
2131 const TrilinosWrappers::PreconditionBase &preconditioner);
2132
2136 TrilinosPayload(const TrilinosPayload &payload);
2137
2145 TrilinosPayload(const TrilinosPayload &first_op,
2146 const TrilinosPayload &second_op);
2147
2151 virtual ~TrilinosPayload() override = default;
2152
2157 identity_payload() const;
2158
2163 null_payload() const;
2164
2169 transpose_payload() const;
2170
2187 template <typename Solver, typename Preconditioner>
2188 std::enable_if_t<
2189 std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
2190 std::is_base_of_v<TrilinosWrappers::PreconditionBase,
2191 Preconditioner>,
2193 inverse_payload(Solver &, const Preconditioner &) const;
2194
2212 template <typename Solver, typename Preconditioner>
2213 std::enable_if_t<
2214 !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
2215 std::is_base_of_v<TrilinosWrappers::PreconditionBase,
2216 Preconditioner>),
2218 inverse_payload(Solver &, const Preconditioner &) const;
2219
2232 IndexSet
2234
2240 IndexSet
2242
2246 MPI_Comm
2247 get_mpi_communicator() const;
2248
2255 void
2256 transpose();
2257
2265 std::function<void(VectorType &, const VectorType &)> vmult;
2266
2274 std::function<void(VectorType &, const VectorType &)> Tvmult;
2275
2284 std::function<void(VectorType &, const VectorType &)> inv_vmult;
2285
2294 std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2295
2309 virtual bool
2310 UseTranspose() const override;
2311
2327 virtual int
2328 SetUseTranspose(bool UseTranspose) override;
2329
2341 virtual int
2342 Apply(const VectorType &X, VectorType &Y) const override;
2343
2362 virtual int
2363 ApplyInverse(const VectorType &Y, VectorType &X) const override;
2377 virtual const char *
2378 Label() const override;
2379
2387 virtual const Epetra_Comm &
2388 Comm() const override;
2389
2397 virtual const Epetra_Map &
2398 OperatorDomainMap() const override;
2399
2408 virtual const Epetra_Map &
2409 OperatorRangeMap() const override;
2412 private:
2422 template <typename EpetraOpType>
2423 TrilinosPayload(EpetraOpType &op,
2424 const bool supports_inverse_operations,
2425 const bool use_transpose,
2426 const MPI_Comm mpi_communicator,
2429
2435
2440 Epetra_MpiComm communicator;
2441
2446 Epetra_Map domain_map;
2447
2452 Epetra_Map range_map;
2453
2462 virtual bool
2463 HasNormInf() const override;
2464
2472 virtual double
2473 NormInf() const override;
2474 };
2475
2481 operator+(const TrilinosPayload &first_op,
2482 const TrilinosPayload &second_op);
2483
2489 operator*(const TrilinosPayload &first_op,
2490 const TrilinosPayload &second_op);
2491
2492 } // namespace LinearOperatorImplementation
2493 } /* namespace internal */
2494
2495
2496
2497 // ----------------------- inline and template functions --------------------
2498
2499# ifndef DOXYGEN
2500
2501 namespace SparseMatrixIterators
2502 {
2504 size_type row,
2505 size_type index)
2506 : matrix(matrix)
2507 , a_row(row)
2508 , a_index(index)
2509 {
2510 visit_present_row();
2511 }
2512
2513
2514 inline AccessorBase::size_type
2515 AccessorBase::row() const
2516 {
2517 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2518 return a_row;
2519 }
2520
2521
2522 inline AccessorBase::size_type
2523 AccessorBase::column() const
2524 {
2525 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2526 return (*colnum_cache)[a_index];
2527 }
2528
2529
2530 inline AccessorBase::size_type
2531 AccessorBase::index() const
2532 {
2533 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2534 return a_index;
2535 }
2536
2537
2538 inline Accessor<true>::Accessor(MatrixType *matrix,
2539 const size_type row,
2540 const size_type index)
2541 : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2542 {}
2543
2544
2545 template <bool Other>
2546 inline Accessor<true>::Accessor(const Accessor<Other> &other)
2547 : AccessorBase(other)
2548 {}
2549
2550
2551 inline TrilinosScalar
2552 Accessor<true>::value() const
2553 {
2554 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2555 return (*value_cache)[a_index];
2556 }
2557
2558
2559 inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2560 : accessor(const_cast<Accessor<false> &>(acc))
2561 {}
2562
2563
2564 inline Accessor<false>::Reference::operator TrilinosScalar() const
2565 {
2566 return (*accessor.value_cache)[accessor.a_index];
2567 }
2568
2569 inline const Accessor<false>::Reference &
2570 Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2571 {
2572 (*accessor.value_cache)[accessor.a_index] = n;
2573 accessor.matrix->set(accessor.row(),
2574 accessor.column(),
2575 static_cast<TrilinosScalar>(*this));
2576 return *this;
2577 }
2578
2579
2580 inline const Accessor<false>::Reference &
2581 Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2582 {
2583 (*accessor.value_cache)[accessor.a_index] += n;
2584 accessor.matrix->set(accessor.row(),
2585 accessor.column(),
2586 static_cast<TrilinosScalar>(*this));
2587 return *this;
2588 }
2589
2590
2591 inline const Accessor<false>::Reference &
2592 Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2593 {
2594 (*accessor.value_cache)[accessor.a_index] -= n;
2595 accessor.matrix->set(accessor.row(),
2596 accessor.column(),
2597 static_cast<TrilinosScalar>(*this));
2598 return *this;
2599 }
2600
2601
2602 inline const Accessor<false>::Reference &
2603 Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2604 {
2605 (*accessor.value_cache)[accessor.a_index] *= n;
2606 accessor.matrix->set(accessor.row(),
2607 accessor.column(),
2608 static_cast<TrilinosScalar>(*this));
2609 return *this;
2610 }
2611
2612
2613 inline const Accessor<false>::Reference &
2614 Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2615 {
2616 (*accessor.value_cache)[accessor.a_index] /= n;
2617 accessor.matrix->set(accessor.row(),
2618 accessor.column(),
2619 static_cast<TrilinosScalar>(*this));
2620 return *this;
2621 }
2622
2623
2624 inline Accessor<false>::Accessor(MatrixType *matrix,
2625 const size_type row,
2626 const size_type index)
2627 : AccessorBase(matrix, row, index)
2628 {}
2629
2630
2631 inline Accessor<false>::Reference
2632 Accessor<false>::value() const
2633 {
2634 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2635 return {*this};
2636 }
2637
2638
2639
2640 template <bool Constness>
2641 inline Iterator<Constness>::Iterator(MatrixType *matrix,
2642 const size_type row,
2643 const size_type index)
2644 : accessor(matrix, row, index)
2645 {}
2646
2647
2648 template <bool Constness>
2649 template <bool Other>
2650 inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2651 : accessor(other.accessor)
2652 {}
2653
2654
2655 template <bool Constness>
2656 inline Iterator<Constness> &
2657 Iterator<Constness>::operator++()
2658 {
2659 Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2660
2661 ++accessor.a_index;
2662
2663 // If at end of line: do one
2664 // step, then cycle until we
2665 // find a row with a nonzero
2666 // number of entries.
2667 if (accessor.a_index >= accessor.colnum_cache->size())
2668 {
2669 accessor.a_index = 0;
2670 ++accessor.a_row;
2671
2672 while ((accessor.a_row < accessor.matrix->m()) &&
2673 ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2674 (accessor.matrix->row_length(accessor.a_row) == 0)))
2675 ++accessor.a_row;
2676
2677 accessor.visit_present_row();
2678 }
2679 return *this;
2680 }
2681
2682
2683 template <bool Constness>
2684 inline Iterator<Constness>
2685 Iterator<Constness>::operator++(int)
2686 {
2687 const Iterator<Constness> old_state = *this;
2688 ++(*this);
2689 return old_state;
2690 }
2691
2692
2693
2694 template <bool Constness>
2695 inline const Accessor<Constness> &
2696 Iterator<Constness>::operator*() const
2697 {
2698 return accessor;
2699 }
2700
2701
2702
2703 template <bool Constness>
2704 inline const Accessor<Constness> *
2705 Iterator<Constness>::operator->() const
2706 {
2707 return &accessor;
2708 }
2709
2710
2711
2712 template <bool Constness>
2713 template <bool OtherConstness>
2714 inline bool
2715 Iterator<Constness>::operator==(const Iterator<OtherConstness> &other) const
2716 {
2717 return (accessor.a_row == other.accessor.a_row &&
2718 accessor.a_index == other.accessor.a_index);
2719 }
2720
2721
2722
2723 template <bool Constness>
2724 template <bool OtherConstness>
2725 inline bool
2726 Iterator<Constness>::operator!=(const Iterator<OtherConstness> &other) const
2727 {
2728 return !(*this == other);
2729 }
2730
2731
2732
2733 template <bool Constness>
2734 template <bool OtherConstness>
2735 inline bool
2736 Iterator<Constness>::operator<(const Iterator<OtherConstness> &other) const
2737 {
2738 return (accessor.row() < other.accessor.row() ||
2739 (accessor.row() == other.accessor.row() &&
2740 accessor.index() < other.accessor.index()));
2741 }
2742
2743
2744 template <bool Constness>
2745 template <bool OtherConstness>
2746 inline bool
2747 Iterator<Constness>::operator>(const Iterator<OtherConstness> &other) const
2748 {
2749 return (other < *this);
2750 }
2751
2752 } // namespace SparseMatrixIterators
2753
2754
2755
2757 SparseMatrix::begin() const
2758 {
2759 return begin(0);
2760 }
2761
2762
2763
2765 SparseMatrix::end() const
2766 {
2767 return const_iterator(this, m(), 0);
2768 }
2769
2770
2771
2773 SparseMatrix::begin(const size_type r) const
2774 {
2775 AssertIndexRange(r, m());
2776 if (in_local_range(r) && (row_length(r) > 0))
2777 return const_iterator(this, r, 0);
2778 else
2779 return end(r);
2780 }
2781
2782
2783
2785 SparseMatrix::end(const size_type r) const
2786 {
2787 AssertIndexRange(r, m());
2788
2789 // place the iterator on the first entry
2790 // past this line, or at the end of the
2791 // matrix
2792 for (size_type i = r + 1; i < m(); ++i)
2793 if (in_local_range(i) && (row_length(i) > 0))
2794 return const_iterator(this, i, 0);
2795
2796 // if there is no such line, then take the
2797 // end iterator of the matrix
2798 return end();
2799 }
2800
2801
2802
2805 {
2806 return begin(0);
2807 }
2808
2809
2810
2813 {
2814 return iterator(this, m(), 0);
2815 }
2816
2817
2818
2820 SparseMatrix::begin(const size_type r)
2821 {
2822 AssertIndexRange(r, m());
2823 if (in_local_range(r) && (row_length(r) > 0))
2824 return iterator(this, r, 0);
2825 else
2826 return end(r);
2827 }
2828
2829
2830
2832 SparseMatrix::end(const size_type r)
2833 {
2834 AssertIndexRange(r, m());
2835
2836 // place the iterator on the first entry
2837 // past this line, or at the end of the
2838 // matrix
2839 for (size_type i = r + 1; i < m(); ++i)
2840 if (in_local_range(i) && (row_length(i) > 0))
2841 return iterator(this, i, 0);
2842
2843 // if there is no such line, then take the
2844 // end iterator of the matrix
2845 return end();
2846 }
2847
2848
2849
2850 inline bool
2851 SparseMatrix::in_local_range(const size_type index) const
2852 {
2854# ifndef DEAL_II_WITH_64BIT_INDICES
2855 begin = matrix->RowMap().MinMyGID();
2856 end = matrix->RowMap().MaxMyGID() + 1;
2857# else
2858 begin = matrix->RowMap().MinMyGID64();
2859 end = matrix->RowMap().MaxMyGID64() + 1;
2860# endif
2861
2862 return ((index >= static_cast<size_type>(begin)) &&
2863 (index < static_cast<size_type>(end)));
2864 }
2865
2866
2867
2868 inline bool
2870 {
2871 return compressed;
2872 }
2873
2874
2875
2876 // Inline the set() and add() functions, since they will be called
2877 // frequently, and the compiler can optimize away some unnecessary loops
2878 // when the sizes are given at compile time.
2879 template <>
2880 void
2881 SparseMatrix::set<TrilinosScalar>(const size_type row,
2882 const size_type n_cols,
2883 const size_type *col_indices,
2884 const TrilinosScalar *values,
2885 const bool elide_zero_values);
2886
2887
2888
2889 template <typename Number>
2890 void
2891 SparseMatrix::set(const size_type row,
2892 const size_type n_cols,
2893 const size_type *col_indices,
2894 const Number *values,
2895 const bool elide_zero_values)
2896 {
2897 std::vector<TrilinosScalar> trilinos_values(n_cols);
2898 std::copy(values, values + n_cols, trilinos_values.begin());
2899 this->set(
2900 row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2901 }
2902
2903
2904
2905 inline void
2906 SparseMatrix::set(const size_type i,
2907 const size_type j,
2908 const TrilinosScalar value)
2909 {
2910 AssertIsFinite(value);
2911
2912 set(i, 1, &j, &value, false);
2913 }
2914
2915
2916
2917 inline void
2918 SparseMatrix::set(const std::vector<size_type> &indices,
2919 const FullMatrix<TrilinosScalar> &values,
2920 const bool elide_zero_values)
2921 {
2922 Assert(indices.size() == values.m(),
2923 ExcDimensionMismatch(indices.size(), values.m()));
2924 Assert(values.m() == values.n(), ExcNotQuadratic());
2925
2926 for (size_type i = 0; i < indices.size(); ++i)
2927 set(indices[i],
2928 indices.size(),
2929 indices.data(),
2930 &values(i, 0),
2931 elide_zero_values);
2932 }
2933
2934
2935
2936 inline void
2937 SparseMatrix::add(const size_type i,
2938 const size_type j,
2939 const TrilinosScalar value)
2940 {
2941 AssertIsFinite(value);
2942
2943 if (value == 0)
2944 {
2945 // we have to check after Insert/Add in any case to be consistent
2946 // with the MPI communication model, but we can save some
2947 // work if the addend is zero. However, these actions are done in case
2948 // we pass on to the other function.
2949
2950 // TODO: fix this (do not run compress here, but fail)
2951 if (last_action == Insert)
2952 {
2953 const int ierr = matrix->GlobalAssemble(*column_space_map,
2954 matrix->RowMap(),
2955 false);
2956
2957 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2958 }
2959
2960 last_action = Add;
2961
2962 return;
2963 }
2964 else
2965 add(i, 1, &j, &value, false);
2966 }
2967
2968
2969
2970 // inline "simple" functions that are called frequently and do only involve
2971 // a call to some Trilinos function.
2973 SparseMatrix::m() const
2974 {
2975# ifndef DEAL_II_WITH_64BIT_INDICES
2976 return matrix->NumGlobalRows();
2977# else
2978 return matrix->NumGlobalRows64();
2979# endif
2980 }
2981
2982
2983
2985 SparseMatrix::n() const
2986 {
2987 // If the matrix structure has not been fixed (i.e., we did not have a
2988 // sparsity pattern), it does not know about the number of columns so we
2989 // must always take this from the additional column space map
2990 Assert(column_space_map.get() != nullptr, ExcInternalError());
2991 return n_global_elements(*column_space_map);
2992 }
2993
2994
2995
2996 inline unsigned int
2998 {
2999 return matrix->NumMyRows();
3000 }
3001
3002
3003
3004 inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
3006 {
3008# ifndef DEAL_II_WITH_64BIT_INDICES
3009 begin = matrix->RowMap().MinMyGID();
3010 end = matrix->RowMap().MaxMyGID() + 1;
3011# else
3012 begin = matrix->RowMap().MinMyGID64();
3013 end = matrix->RowMap().MaxMyGID64() + 1;
3014# endif
3015
3016 return std::make_pair(begin, end);
3017 }
3018
3019
3020
3021 inline std::uint64_t
3023 {
3024 // Trilinos uses 64bit functions internally for attribute access, which
3025 // return `long long`. They also offer 32bit variants that return `int`,
3026 // however those call the 64bit version and convert the values to 32bit.
3027 // There is no necessity in using the 32bit versions at all.
3028 return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
3029 }
3030
3031
3032
3033 template <typename SparsityPatternType>
3034 inline void
3035 SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3036 const SparsityPatternType &sparsity_pattern,
3037 const MPI_Comm communicator,
3038 const bool exchange_data)
3039 {
3040 reinit(parallel_partitioning,
3041 parallel_partitioning,
3042 sparsity_pattern,
3043 communicator,
3044 exchange_data);
3045 }
3046
3047
3048
3049 template <typename number>
3050 inline void
3051 SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3052 const ::SparseMatrix<number> &sparse_matrix,
3053 const MPI_Comm communicator,
3054 const double drop_tolerance,
3055 const bool copy_values,
3056 const ::SparsityPattern *use_this_sparsity)
3057 {
3058 Epetra_Map map =
3059 parallel_partitioning.make_trilinos_map(communicator, false);
3060 reinit(parallel_partitioning,
3061 parallel_partitioning,
3062 sparse_matrix,
3063 drop_tolerance,
3064 copy_values,
3065 use_this_sparsity);
3066 }
3067
3068
3069
3070 inline const Epetra_CrsMatrix &
3072 {
3073 return static_cast<const Epetra_CrsMatrix &>(*matrix);
3074 }
3075
3076
3077
3078 inline const Epetra_CrsGraph &
3080 {
3081 return matrix->Graph();
3082 }
3083
3084
3085
3086 inline IndexSet
3088 {
3089 return IndexSet(matrix->DomainMap());
3090 }
3091
3092
3093
3094 inline IndexSet
3096 {
3097 return IndexSet(matrix->RangeMap());
3098 }
3099
3100
3101
3102 inline void
3104 {
3105 // nothing to do here
3106 }
3107
3108
3109
3110 inline void
3112 {
3113 // nothing to do here
3114 }
3115
3116
3117
3118 template <typename VectorType>
3119 inline TrilinosScalar
3120 SparseMatrix::residual(VectorType &dst,
3121 const VectorType &x,
3122 const VectorType &b) const
3123 {
3124 vmult(dst, x);
3125 dst -= b;
3126 dst *= -1.;
3127
3128 return dst.l2_norm();
3129 }
3130
3131
3132 namespace internal
3133 {
3134 namespace LinearOperatorImplementation
3135 {
3136 template <typename EpetraOpType>
3137 TrilinosPayload::TrilinosPayload(
3138 EpetraOpType &op,
3139 const bool supports_inverse_operations,
3140 const bool use_transpose,
3141 const MPI_Comm mpi_communicator,
3142 const IndexSet &locally_owned_domain_indices,
3143 const IndexSet &locally_owned_range_indices)
3144 : use_transpose(use_transpose)
3145 , communicator(mpi_communicator)
3146 , domain_map(
3147 locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3148 , range_map(
3149 locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3150 {
3151 vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3152 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3153 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3154 Assert(&tril_src != &tril_dst,
3157 tril_src,
3158 tril_dst,
3159 op.UseTranspose());
3160
3161 const int ierr = op.Apply(tril_src, tril_dst);
3162 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3163 };
3164
3165 Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3166 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3167 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3168 Assert(&tril_src != &tril_dst,
3171 tril_src,
3172 tril_dst,
3173 !op.UseTranspose());
3174
3175 op.SetUseTranspose(!op.UseTranspose());
3176 const int ierr = op.Apply(tril_src, tril_dst);
3177 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3178 op.SetUseTranspose(!op.UseTranspose());
3179 };
3180
3181 if (supports_inverse_operations)
3182 {
3183 inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3184 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3185 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3186 Assert(
3187 &tril_src != &tril_dst,
3190 tril_src,
3191 tril_dst,
3192 !op.UseTranspose());
3193
3194 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3195 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3196 };
3197
3198 inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
3199 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3200 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3201 Assert(
3202 &tril_src != &tril_dst,
3205 tril_src,
3206 tril_dst,
3207 op.UseTranspose());
3208
3209 op.SetUseTranspose(!op.UseTranspose());
3210 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3211 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3212 op.SetUseTranspose(!op.UseTranspose());
3213 };
3214 }
3215 else
3216 {
3217 inv_vmult = [](Domain &, const Range &) {
3218 Assert(false,
3219 ExcMessage(
3220 "Uninitialized TrilinosPayload::inv_vmult called. "
3221 "The operator does not support inverse operations."));
3222 };
3223
3224 inv_Tvmult = [](Range &, const Domain &) {
3225 Assert(false,
3226 ExcMessage(
3227 "Uninitialized TrilinosPayload::inv_Tvmult called. "
3228 "The operator does not support inverse operations."));
3229 };
3230 }
3231 }
3232
3233
3234 template <typename Solver, typename Preconditioner>
3235 std::enable_if_t<
3236 std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
3237 std::is_base_of_v<TrilinosWrappers::PreconditionBase, Preconditioner>,
3238 TrilinosPayload>
3240 Solver &solver,
3241 const Preconditioner &preconditioner) const
3242 {
3243 const auto &payload = *this;
3244
3245 TrilinosPayload return_op(payload);
3246
3247 // Capture by copy so the payloads are always valid
3248
3249 return_op.inv_vmult = [payload, &solver, &preconditioner](
3250 TrilinosPayload::Domain &tril_dst,
3251 const TrilinosPayload::Range &tril_src) {
3252 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3253 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3254 Assert(&tril_src != &tril_dst,
3257 tril_src,
3258 tril_dst,
3259 !payload.UseTranspose());
3260 solver.solve(payload, tril_dst, tril_src, preconditioner);
3261 };
3262
3263 return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3264 TrilinosPayload::Range &tril_dst,
3265 const TrilinosPayload::Domain &tril_src) {
3266 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3267 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3268 Assert(&tril_src != &tril_dst,
3271 tril_src,
3272 tril_dst,
3273 payload.UseTranspose());
3274
3275 const_cast<TrilinosPayload &>(payload).transpose();
3276 solver.solve(payload, tril_dst, tril_src, preconditioner);
3277 const_cast<TrilinosPayload &>(payload).transpose();
3278 };
3279
3280 // If the input operator is already setup for transpose operations, then
3281 // we must do similar with its inverse.
3282 if (return_op.UseTranspose() == true)
3283 std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3284
3285 return return_op;
3286 }
3287
3288 template <typename Solver, typename Preconditioner>
3289 std::enable_if_t<
3290 !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
3291 std::is_base_of_v<TrilinosWrappers::PreconditionBase,
3292 Preconditioner>),
3293 TrilinosPayload>
3294 TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3295 {
3296 TrilinosPayload return_op(*this);
3297
3298 return_op.inv_vmult = [](TrilinosPayload::Domain &,
3299 const TrilinosPayload::Range &) {
3300 AssertThrow(false,
3301 ExcMessage("Payload inv_vmult disabled because of "
3302 "incompatible solver/preconditioner choice."));
3303 };
3304
3305 return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3306 const TrilinosPayload::Domain &) {
3307 AssertThrow(false,
3308 ExcMessage("Payload inv_vmult disabled because of "
3309 "incompatible solver/preconditioner choice."));
3310 };
3311
3312 return return_op;
3313 }
3314 } // namespace LinearOperatorImplementation
3315 } // namespace internal
3316
3317 template <>
3318 void
3319 SparseMatrix::set<TrilinosScalar>(const size_type row,
3320 const size_type n_cols,
3321 const size_type *col_indices,
3322 const TrilinosScalar *values,
3323 const bool elide_zero_values);
3324# endif // DOXYGEN
3325
3326} /* namespace TrilinosWrappers */
3327
3328
3330
3331
3332# endif // DEAL_II_WITH_TRILINOS
3333
3334
3335/*----------------------- trilinos_sparse_matrix.h --------------------*/
3336
3337#endif
3338/*----------------------- trilinos_sparse_matrix.h --------------------*/
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
AccessorBase(SparseMatrix *matrix, const size_type row, const size_type index)
std::shared_ptr< std::vector< TrilinosScalar > > value_cache
std::shared_ptr< std::vector< size_type > > colnum_cache
const Reference & operator-=(const TrilinosScalar n) const
const Reference & operator*=(const TrilinosScalar n) const
const Reference & operator+=(const TrilinosScalar n) const
const Reference & operator/=(const TrilinosScalar n) const
const Reference & operator=(const TrilinosScalar n) const
Accessor(MatrixType *matrix, const size_type row, const size_type index)
Accessor(MatrixType *matrix, const size_type row, const size_type index)
const Accessor< Constness > & operator*() const
const Accessor< Constness > * operator->() const
typename Accessor< Constness >::MatrixType MatrixType
bool operator==(const Iterator< OtherConstness > &) const
bool operator!=(const Iterator< OtherConstness > &) const
bool operator<(const Iterator< OtherConstness > &) const
Iterator(const Iterator< Other > &other)
bool operator>(const Iterator< OtherConstness > &) const
Iterator(MatrixType *matrix, const size_type row, const size_type index)
void set(const size_type i, const size_type j, const TrilinosScalar value)
std::unique_ptr< Epetra_Map > column_space_map
TrilinosScalar residual(VectorType &dst, const VectorType &x, const VectorType &b) const
SparseMatrix & operator*=(const TrilinosScalar factor)
void mmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::unique_ptr< Epetra_Export > nonlocal_matrix_exporter
void compress(VectorOperation::values operation)
std::unique_ptr< Epetra_FECrsMatrix > matrix
const Epetra_CrsMatrix & trilinos_matrix() const
TrilinosScalar matrix_norm_square(const MPI::Vector &v) const
IndexSet locally_owned_range_indices() const
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const MPI::Vector &V=MPI::Vector()) const
std::enable_if_t< std::is_same_v< typename VectorType::value_type, TrilinosScalar > > Tvmult(VectorType &dst, const VectorType &src) const
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
void clear_rows(const ArrayView< const size_type > &rows, const TrilinosScalar new_diag_value=0)
const_iterator begin() const
void reinit(const SparsityPatternType &sparsity_pattern)
void vmult_add(VectorType &dst, const VectorType &src) const
SparseMatrix & operator=(const SparseMatrix &)=delete
const_iterator end(const size_type r) const
SparseMatrix & operator/=(const TrilinosScalar factor)
void Tvmult_add(VectorType &dst, const VectorType &src) const
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)
TrilinosScalar el(const size_type i, const size_type j) const
IndexSet locally_owned_domain_indices() const
bool in_local_range(const size_type index) const
void copy_from(const SparseMatrix &source)
unsigned int row_length(const size_type row) const
std::uint64_t n_nonzero_elements() const
std::enable_if_t< std::is_same_v< typename VectorType::value_type, TrilinosScalar > > vmult(VectorType &dst, const VectorType &src) const
const_iterator begin(const size_type r) const
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
SparseMatrix(const SparseMatrix &)=delete
iterator end(const size_type r)
TrilinosScalar diag_element(const size_type i) const
void add(const size_type i, const size_type j, const TrilinosScalar value)
unsigned int local_size() const
iterator begin(const size_type r)
TrilinosScalar operator()(const size_type i, const size_type j) const
void reinit(const IndexSet &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)
virtual ~SparseMatrix() override=default
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)
TrilinosScalar matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const
const_iterator end() const
void set(const std::vector< size_type > &indices, const FullMatrix< TrilinosScalar > &full_matrix, const bool elide_zero_values=false)
std::pair< size_type, size_type > local_range() const
void reinit(const IndexSet &parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator=MPI_COMM_WORLD, const bool exchange_data=false)
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
std::enable_if_t< std::is_base_of_v< TrilinosWrappers::SolverBase, Solver > &&std::is_base_of_v< TrilinosWrappers::PreconditionBase, Preconditioner >, TrilinosPayload > inverse_payload(Solver &, const Preconditioner &) const
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
virtual int Apply(const VectorType &X, VectorType &Y) const override
TrilinosPayload(EpetraOpType &op, const bool supports_inverse_operations, const bool use_transpose, const MPI_Comm mpi_communicator, const IndexSet &locally_owned_domain_indices, const IndexSet &locally_owned_range_indices)
std::enable_if_t< !(std::is_base_of_v< TrilinosWrappers::SolverBase, Solver > &&std::is_base_of_v< TrilinosWrappers::PreconditionBase, Preconditioner >), TrilinosPayload > inverse_payload(Solver &, const Preconditioner &) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
#define DeclException0(Exception0)
Definition exceptions.h:471
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:585
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define AssertIsFinite(number)
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
static ::ExceptionBase & ExcBeyondEndOfMatrix()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition exceptions.h:562
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
static ::ExceptionBase & ExcMatrixNotCompressed()
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertThrow(cond, exc)
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
TrilinosPayload operator+(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
TrilinosPayload operator*(const TrilinosPayload &first_op, const TrilinosPayload &second_op)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
void check_vector_map_equality(const Epetra_CrsMatrix &mtrx, const Epetra_MultiVector &src, const Epetra_MultiVector &dst, const bool transpose)
TrilinosWrappers::types::int64_type n_global_elements(const Epetra_BlockMap &map)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
STL namespace.
unsigned int global_dof_index
Definition types.h:81
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::value_type value_type
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::difference_type difference_type
double TrilinosScalar
Definition types.h:178