Reference documentation for deal.II version GIT relicensing-891-g13379fae33 2024-06-22 22:40: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
2048 namespace LinearOperatorImplementation
2049 {
2070 {
2071 public:
2075 using VectorType = Epetra_MultiVector;
2076
2081
2086
2091
2105
2109 TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2110 const TrilinosWrappers::SparseMatrix &matrix);
2111
2115 TrilinosPayload(const TrilinosPayload &payload_exemplar,
2116 const TrilinosWrappers::SparseMatrix &matrix);
2117
2122 const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2123 const TrilinosWrappers::PreconditionBase &preconditioner);
2124
2129 const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2130 const TrilinosWrappers::PreconditionBase &preconditioner);
2131
2136 const TrilinosPayload &payload_exemplar,
2137 const TrilinosWrappers::PreconditionBase &preconditioner);
2138
2142 TrilinosPayload(const TrilinosPayload &payload);
2143
2151 TrilinosPayload(const TrilinosPayload &first_op,
2152 const TrilinosPayload &second_op);
2153
2157 virtual ~TrilinosPayload() override = default;
2158
2163 identity_payload() const;
2164
2169 null_payload() const;
2170
2175 transpose_payload() const;
2176
2193 template <typename Solver, typename Preconditioner>
2194 std::enable_if_t<
2195 std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
2196 std::is_base_of_v<TrilinosWrappers::PreconditionBase,
2197 Preconditioner>,
2199 inverse_payload(Solver &, const Preconditioner &) const;
2200
2218 template <typename Solver, typename Preconditioner>
2219 std::enable_if_t<
2220 !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
2221 std::is_base_of_v<TrilinosWrappers::PreconditionBase,
2222 Preconditioner>),
2224 inverse_payload(Solver &, const Preconditioner &) const;
2225
2238 IndexSet
2240
2246 IndexSet
2248
2252 MPI_Comm
2253 get_mpi_communicator() const;
2254
2261 void
2262 transpose();
2263
2271 std::function<void(VectorType &, const VectorType &)> vmult;
2272
2280 std::function<void(VectorType &, const VectorType &)> Tvmult;
2281
2290 std::function<void(VectorType &, const VectorType &)> inv_vmult;
2291
2300 std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2301
2315 virtual bool
2316 UseTranspose() const override;
2317
2333 virtual int
2334 SetUseTranspose(bool UseTranspose) override;
2335
2347 virtual int
2348 Apply(const VectorType &X, VectorType &Y) const override;
2349
2368 virtual int
2369 ApplyInverse(const VectorType &Y, VectorType &X) const override;
2383 virtual const char *
2384 Label() const override;
2385
2393 virtual const Epetra_Comm &
2394 Comm() const override;
2395
2403 virtual const Epetra_Map &
2404 OperatorDomainMap() const override;
2405
2414 virtual const Epetra_Map &
2415 OperatorRangeMap() const override;
2418 private:
2428 template <typename EpetraOpType>
2429 TrilinosPayload(EpetraOpType &op,
2430 const bool supports_inverse_operations,
2431 const bool use_transpose,
2432 const MPI_Comm mpi_communicator,
2435
2441
2446 Epetra_MpiComm communicator;
2447
2452 Epetra_Map domain_map;
2453
2458 Epetra_Map range_map;
2459
2468 virtual bool
2469 HasNormInf() const override;
2470
2478 virtual double
2479 NormInf() const override;
2480 };
2481
2487 operator+(const TrilinosPayload &first_op,
2488 const TrilinosPayload &second_op);
2489
2495 operator*(const TrilinosPayload &first_op,
2496 const TrilinosPayload &second_op);
2497
2498 } // namespace LinearOperatorImplementation
2499 } /* namespace internal */
2500
2501
2502
2503 // ----------------------- inline and template functions --------------------
2504
2505# ifndef DOXYGEN
2506
2507 namespace SparseMatrixIterators
2508 {
2510 size_type row,
2511 size_type index)
2512 : matrix(matrix)
2513 , a_row(row)
2514 , a_index(index)
2515 {
2516 visit_present_row();
2517 }
2518
2519
2520 inline AccessorBase::size_type
2521 AccessorBase::row() const
2522 {
2523 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2524 return a_row;
2525 }
2526
2527
2528 inline AccessorBase::size_type
2529 AccessorBase::column() const
2530 {
2531 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2532 return (*colnum_cache)[a_index];
2533 }
2534
2535
2536 inline AccessorBase::size_type
2537 AccessorBase::index() const
2538 {
2539 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2540 return a_index;
2541 }
2542
2543
2544 inline Accessor<true>::Accessor(MatrixType *matrix,
2545 const size_type row,
2546 const size_type index)
2547 : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2548 {}
2549
2550
2551 template <bool Other>
2552 inline Accessor<true>::Accessor(const Accessor<Other> &other)
2553 : AccessorBase(other)
2554 {}
2555
2556
2557 inline TrilinosScalar
2558 Accessor<true>::value() const
2559 {
2560 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2561 return (*value_cache)[a_index];
2562 }
2563
2564
2565 inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2566 : accessor(const_cast<Accessor<false> &>(acc))
2567 {}
2568
2569
2570 inline Accessor<false>::Reference::operator TrilinosScalar() const
2571 {
2572 return (*accessor.value_cache)[accessor.a_index];
2573 }
2574
2575 inline const Accessor<false>::Reference &
2576 Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2577 {
2578 (*accessor.value_cache)[accessor.a_index] = n;
2579 accessor.matrix->set(accessor.row(),
2580 accessor.column(),
2581 static_cast<TrilinosScalar>(*this));
2582 return *this;
2583 }
2584
2585
2586 inline const Accessor<false>::Reference &
2587 Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2588 {
2589 (*accessor.value_cache)[accessor.a_index] += n;
2590 accessor.matrix->set(accessor.row(),
2591 accessor.column(),
2592 static_cast<TrilinosScalar>(*this));
2593 return *this;
2594 }
2595
2596
2597 inline const Accessor<false>::Reference &
2598 Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2599 {
2600 (*accessor.value_cache)[accessor.a_index] -= n;
2601 accessor.matrix->set(accessor.row(),
2602 accessor.column(),
2603 static_cast<TrilinosScalar>(*this));
2604 return *this;
2605 }
2606
2607
2608 inline const Accessor<false>::Reference &
2609 Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2610 {
2611 (*accessor.value_cache)[accessor.a_index] *= n;
2612 accessor.matrix->set(accessor.row(),
2613 accessor.column(),
2614 static_cast<TrilinosScalar>(*this));
2615 return *this;
2616 }
2617
2618
2619 inline const Accessor<false>::Reference &
2620 Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2621 {
2622 (*accessor.value_cache)[accessor.a_index] /= n;
2623 accessor.matrix->set(accessor.row(),
2624 accessor.column(),
2625 static_cast<TrilinosScalar>(*this));
2626 return *this;
2627 }
2628
2629
2630 inline Accessor<false>::Accessor(MatrixType *matrix,
2631 const size_type row,
2632 const size_type index)
2633 : AccessorBase(matrix, row, index)
2634 {}
2635
2636
2637 inline Accessor<false>::Reference
2638 Accessor<false>::value() const
2639 {
2640 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2641 return {*this};
2642 }
2643
2644
2645
2646 template <bool Constness>
2647 inline Iterator<Constness>::Iterator(MatrixType *matrix,
2648 const size_type row,
2649 const size_type index)
2650 : accessor(matrix, row, index)
2651 {}
2652
2653
2654 template <bool Constness>
2655 template <bool Other>
2656 inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2657 : accessor(other.accessor)
2658 {}
2659
2660
2661 template <bool Constness>
2662 inline Iterator<Constness> &
2663 Iterator<Constness>::operator++()
2664 {
2665 Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2666
2667 ++accessor.a_index;
2668
2669 // If at end of line: do one
2670 // step, then cycle until we
2671 // find a row with a nonzero
2672 // number of entries.
2673 if (accessor.a_index >= accessor.colnum_cache->size())
2674 {
2675 accessor.a_index = 0;
2676 ++accessor.a_row;
2677
2678 while ((accessor.a_row < accessor.matrix->m()) &&
2679 ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2680 (accessor.matrix->row_length(accessor.a_row) == 0)))
2681 ++accessor.a_row;
2682
2683 accessor.visit_present_row();
2684 }
2685 return *this;
2686 }
2687
2688
2689 template <bool Constness>
2690 inline Iterator<Constness>
2691 Iterator<Constness>::operator++(int)
2692 {
2693 const Iterator<Constness> old_state = *this;
2694 ++(*this);
2695 return old_state;
2696 }
2697
2698
2699
2700 template <bool Constness>
2701 inline const Accessor<Constness> &
2702 Iterator<Constness>::operator*() const
2703 {
2704 return accessor;
2705 }
2706
2707
2708
2709 template <bool Constness>
2710 inline const Accessor<Constness> *
2711 Iterator<Constness>::operator->() const
2712 {
2713 return &accessor;
2714 }
2715
2716
2717
2718 template <bool Constness>
2719 template <bool OtherConstness>
2720 inline bool
2721 Iterator<Constness>::operator==(const Iterator<OtherConstness> &other) const
2722 {
2723 return (accessor.a_row == other.accessor.a_row &&
2724 accessor.a_index == other.accessor.a_index);
2725 }
2726
2727
2728
2729 template <bool Constness>
2730 template <bool OtherConstness>
2731 inline bool
2732 Iterator<Constness>::operator!=(const Iterator<OtherConstness> &other) const
2733 {
2734 return !(*this == other);
2735 }
2736
2737
2738
2739 template <bool Constness>
2740 template <bool OtherConstness>
2741 inline bool
2742 Iterator<Constness>::operator<(const Iterator<OtherConstness> &other) const
2743 {
2744 return (accessor.row() < other.accessor.row() ||
2745 (accessor.row() == other.accessor.row() &&
2746 accessor.index() < other.accessor.index()));
2747 }
2748
2749
2750 template <bool Constness>
2751 template <bool OtherConstness>
2752 inline bool
2753 Iterator<Constness>::operator>(const Iterator<OtherConstness> &other) const
2754 {
2755 return (other < *this);
2756 }
2757
2758 } // namespace SparseMatrixIterators
2759
2760
2761
2763 SparseMatrix::begin() const
2764 {
2765 return begin(0);
2766 }
2767
2768
2769
2771 SparseMatrix::end() const
2772 {
2773 return const_iterator(this, m(), 0);
2774 }
2775
2776
2777
2779 SparseMatrix::begin(const size_type r) const
2780 {
2781 AssertIndexRange(r, m());
2782 if (in_local_range(r) && (row_length(r) > 0))
2783 return const_iterator(this, r, 0);
2784 else
2785 return end(r);
2786 }
2787
2788
2789
2791 SparseMatrix::end(const size_type r) const
2792 {
2793 AssertIndexRange(r, m());
2794
2795 // place the iterator on the first entry
2796 // past this line, or at the end of the
2797 // matrix
2798 for (size_type i = r + 1; i < m(); ++i)
2799 if (in_local_range(i) && (row_length(i) > 0))
2800 return const_iterator(this, i, 0);
2801
2802 // if there is no such line, then take the
2803 // end iterator of the matrix
2804 return end();
2805 }
2806
2807
2808
2811 {
2812 return begin(0);
2813 }
2814
2815
2816
2819 {
2820 return iterator(this, m(), 0);
2821 }
2822
2823
2824
2826 SparseMatrix::begin(const size_type r)
2827 {
2828 AssertIndexRange(r, m());
2829 if (in_local_range(r) && (row_length(r) > 0))
2830 return iterator(this, r, 0);
2831 else
2832 return end(r);
2833 }
2834
2835
2836
2838 SparseMatrix::end(const size_type r)
2839 {
2840 AssertIndexRange(r, m());
2841
2842 // place the iterator on the first entry
2843 // past this line, or at the end of the
2844 // matrix
2845 for (size_type i = r + 1; i < m(); ++i)
2846 if (in_local_range(i) && (row_length(i) > 0))
2847 return iterator(this, i, 0);
2848
2849 // if there is no such line, then take the
2850 // end iterator of the matrix
2851 return end();
2852 }
2853
2854
2855
2856 inline bool
2857 SparseMatrix::in_local_range(const size_type index) const
2858 {
2860# ifndef DEAL_II_WITH_64BIT_INDICES
2861 begin = matrix->RowMap().MinMyGID();
2862 end = matrix->RowMap().MaxMyGID() + 1;
2863# else
2864 begin = matrix->RowMap().MinMyGID64();
2865 end = matrix->RowMap().MaxMyGID64() + 1;
2866# endif
2867
2868 return ((index >= static_cast<size_type>(begin)) &&
2869 (index < static_cast<size_type>(end)));
2870 }
2871
2872
2873
2874 inline bool
2876 {
2877 return compressed;
2878 }
2879
2880
2881
2882 // Inline the set() and add() functions, since they will be called
2883 // frequently, and the compiler can optimize away some unnecessary loops
2884 // when the sizes are given at compile time.
2885 template <>
2886 void
2887 SparseMatrix::set<TrilinosScalar>(const size_type row,
2888 const size_type n_cols,
2889 const size_type *col_indices,
2890 const TrilinosScalar *values,
2891 const bool elide_zero_values);
2892
2893
2894
2895 template <typename Number>
2896 void
2897 SparseMatrix::set(const size_type row,
2898 const size_type n_cols,
2899 const size_type *col_indices,
2900 const Number *values,
2901 const bool elide_zero_values)
2902 {
2903 std::vector<TrilinosScalar> trilinos_values(n_cols);
2904 std::copy(values, values + n_cols, trilinos_values.begin());
2905 this->set(
2906 row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2907 }
2908
2909
2910
2911 inline void
2912 SparseMatrix::set(const size_type i,
2913 const size_type j,
2914 const TrilinosScalar value)
2915 {
2916 AssertIsFinite(value);
2917
2918 set(i, 1, &j, &value, false);
2919 }
2920
2921
2922
2923 inline void
2924 SparseMatrix::set(const std::vector<size_type> &indices,
2925 const FullMatrix<TrilinosScalar> &values,
2926 const bool elide_zero_values)
2927 {
2928 Assert(indices.size() == values.m(),
2929 ExcDimensionMismatch(indices.size(), values.m()));
2930 Assert(values.m() == values.n(), ExcNotQuadratic());
2931
2932 for (size_type i = 0; i < indices.size(); ++i)
2933 set(indices[i],
2934 indices.size(),
2935 indices.data(),
2936 &values(i, 0),
2937 elide_zero_values);
2938 }
2939
2940
2941
2942 inline void
2943 SparseMatrix::add(const size_type i,
2944 const size_type j,
2945 const TrilinosScalar value)
2946 {
2947 AssertIsFinite(value);
2948
2949 if (value == 0)
2950 {
2951 // we have to check after Insert/Add in any case to be consistent
2952 // with the MPI communication model, but we can save some
2953 // work if the addend is zero. However, these actions are done in case
2954 // we pass on to the other function.
2955
2956 // TODO: fix this (do not run compress here, but fail)
2957 if (last_action == Insert)
2958 {
2959 const int ierr = matrix->GlobalAssemble(*column_space_map,
2960 matrix->RowMap(),
2961 false);
2962
2963 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2964 }
2965
2966 last_action = Add;
2967
2968 return;
2969 }
2970 else
2971 add(i, 1, &j, &value, false);
2972 }
2973
2974
2975
2976 // inline "simple" functions that are called frequently and do only involve
2977 // a call to some Trilinos function.
2979 SparseMatrix::m() const
2980 {
2981# ifndef DEAL_II_WITH_64BIT_INDICES
2982 return matrix->NumGlobalRows();
2983# else
2984 return matrix->NumGlobalRows64();
2985# endif
2986 }
2987
2988
2989
2991 SparseMatrix::n() const
2992 {
2993 // If the matrix structure has not been fixed (i.e., we did not have a
2994 // sparsity pattern), it does not know about the number of columns so we
2995 // must always take this from the additional column space map
2996 Assert(column_space_map.get() != nullptr, ExcInternalError());
2997 return n_global_elements(*column_space_map);
2998 }
2999
3000
3001
3002 inline unsigned int
3004 {
3005 return matrix->NumMyRows();
3006 }
3007
3008
3009
3010 inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
3012 {
3014# ifndef DEAL_II_WITH_64BIT_INDICES
3015 begin = matrix->RowMap().MinMyGID();
3016 end = matrix->RowMap().MaxMyGID() + 1;
3017# else
3018 begin = matrix->RowMap().MinMyGID64();
3019 end = matrix->RowMap().MaxMyGID64() + 1;
3020# endif
3021
3022 return std::make_pair(begin, end);
3023 }
3024
3025
3026
3027 inline std::uint64_t
3029 {
3030 // Trilinos uses 64bit functions internally for attribute access, which
3031 // return `long long`. They also offer 32bit variants that return `int`,
3032 // however those call the 64bit version and convert the values to 32bit.
3033 // There is no necessity in using the 32bit versions at all.
3034 return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
3035 }
3036
3037
3038
3039 template <typename SparsityPatternType>
3040 inline void
3041 SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3042 const SparsityPatternType &sparsity_pattern,
3043 const MPI_Comm communicator,
3044 const bool exchange_data)
3045 {
3046 reinit(parallel_partitioning,
3047 parallel_partitioning,
3048 sparsity_pattern,
3049 communicator,
3050 exchange_data);
3051 }
3052
3053
3054
3055 template <typename number>
3056 inline void
3057 SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3058 const ::SparseMatrix<number> &sparse_matrix,
3059 const MPI_Comm communicator,
3060 const double drop_tolerance,
3061 const bool copy_values,
3062 const ::SparsityPattern *use_this_sparsity)
3063 {
3064 Epetra_Map map =
3065 parallel_partitioning.make_trilinos_map(communicator, false);
3066 reinit(parallel_partitioning,
3067 parallel_partitioning,
3068 sparse_matrix,
3069 drop_tolerance,
3070 copy_values,
3071 use_this_sparsity);
3072 }
3073
3074
3075
3076 inline const Epetra_CrsMatrix &
3078 {
3079 return static_cast<const Epetra_CrsMatrix &>(*matrix);
3080 }
3081
3082
3083
3084 inline const Epetra_CrsGraph &
3086 {
3087 return matrix->Graph();
3088 }
3089
3090
3091
3092 inline IndexSet
3094 {
3095 return IndexSet(matrix->DomainMap());
3096 }
3097
3098
3099
3100 inline IndexSet
3102 {
3103 return IndexSet(matrix->RangeMap());
3104 }
3105
3106
3107
3108 inline void
3110 {
3111 // nothing to do here
3112 }
3113
3114
3115
3116 inline void
3118 {
3119 // nothing to do here
3120 }
3121
3122
3123
3124 template <typename VectorType>
3125 inline TrilinosScalar
3126 SparseMatrix::residual(VectorType &dst,
3127 const VectorType &x,
3128 const VectorType &b) const
3129 {
3130 vmult(dst, x);
3131 dst -= b;
3132 dst *= -1.;
3133
3134 return dst.l2_norm();
3135 }
3136
3137
3138 namespace internal
3139 {
3140 namespace LinearOperatorImplementation
3141 {
3142 template <typename EpetraOpType>
3143 TrilinosPayload::TrilinosPayload(
3144 EpetraOpType &op,
3145 const bool supports_inverse_operations,
3146 const bool use_transpose,
3147 const MPI_Comm mpi_communicator,
3148 const IndexSet &locally_owned_domain_indices,
3149 const IndexSet &locally_owned_range_indices)
3150 : use_transpose(use_transpose)
3151 , communicator(mpi_communicator)
3152 , domain_map(
3153 locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3154 , range_map(
3155 locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3156 {
3157 vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3158 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3159 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3160 Assert(&tril_src != &tril_dst,
3163 tril_src,
3164 tril_dst,
3165 op.UseTranspose());
3166
3167 const int ierr = op.Apply(tril_src, tril_dst);
3168 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3169 };
3170
3171 Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3172 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3173 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3174 Assert(&tril_src != &tril_dst,
3177 tril_src,
3178 tril_dst,
3179 !op.UseTranspose());
3180
3181 op.SetUseTranspose(!op.UseTranspose());
3182 const int ierr = op.Apply(tril_src, tril_dst);
3183 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3184 op.SetUseTranspose(!op.UseTranspose());
3185 };
3186
3187 if (supports_inverse_operations)
3188 {
3189 inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3190 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3191 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3192 Assert(
3193 &tril_src != &tril_dst,
3196 tril_src,
3197 tril_dst,
3198 !op.UseTranspose());
3199
3200 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3201 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3202 };
3203
3204 inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
3205 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3206 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3207 Assert(
3208 &tril_src != &tril_dst,
3211 tril_src,
3212 tril_dst,
3213 op.UseTranspose());
3214
3215 op.SetUseTranspose(!op.UseTranspose());
3216 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3217 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3218 op.SetUseTranspose(!op.UseTranspose());
3219 };
3220 }
3221 else
3222 {
3223 inv_vmult = [](Domain &, const Range &) {
3224 Assert(false,
3225 ExcMessage(
3226 "Uninitialized TrilinosPayload::inv_vmult called. "
3227 "The operator does not support inverse operations."));
3228 };
3229
3230 inv_Tvmult = [](Range &, const Domain &) {
3231 Assert(false,
3232 ExcMessage(
3233 "Uninitialized TrilinosPayload::inv_Tvmult called. "
3234 "The operator does not support inverse operations."));
3235 };
3236 }
3237 }
3238
3239
3240 template <typename Solver, typename Preconditioner>
3241 std::enable_if_t<
3242 std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
3243 std::is_base_of_v<TrilinosWrappers::PreconditionBase, Preconditioner>,
3244 TrilinosPayload>
3246 Solver &solver,
3247 const Preconditioner &preconditioner) const
3248 {
3249 const auto &payload = *this;
3250
3251 TrilinosPayload return_op(payload);
3252
3253 // Capture by copy so the payloads are always valid
3254
3255 return_op.inv_vmult = [payload, &solver, &preconditioner](
3256 TrilinosPayload::Domain &tril_dst,
3257 const TrilinosPayload::Range &tril_src) {
3258 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3259 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3260 Assert(&tril_src != &tril_dst,
3263 tril_src,
3264 tril_dst,
3265 !payload.UseTranspose());
3266 solver.solve(payload, tril_dst, tril_src, preconditioner);
3267 };
3268
3269 return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3270 TrilinosPayload::Range &tril_dst,
3271 const TrilinosPayload::Domain &tril_src) {
3272 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3273 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3274 Assert(&tril_src != &tril_dst,
3277 tril_src,
3278 tril_dst,
3279 payload.UseTranspose());
3280
3281 const_cast<TrilinosPayload &>(payload).transpose();
3282 solver.solve(payload, tril_dst, tril_src, preconditioner);
3283 const_cast<TrilinosPayload &>(payload).transpose();
3284 };
3285
3286 // If the input operator is already setup for transpose operations, then
3287 // we must do similar with its inverse.
3288 if (return_op.UseTranspose() == true)
3289 std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3290
3291 return return_op;
3292 }
3293
3294 template <typename Solver, typename Preconditioner>
3295 std::enable_if_t<
3296 !(std::is_base_of_v<TrilinosWrappers::SolverBase, Solver> &&
3297 std::is_base_of_v<TrilinosWrappers::PreconditionBase,
3298 Preconditioner>),
3299 TrilinosPayload>
3300 TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3301 {
3302 TrilinosPayload return_op(*this);
3303
3304 return_op.inv_vmult = [](TrilinosPayload::Domain &,
3305 const TrilinosPayload::Range &) {
3306 AssertThrow(false,
3307 ExcMessage("Payload inv_vmult disabled because of "
3308 "incompatible solver/preconditioner choice."));
3309 };
3310
3311 return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3312 const TrilinosPayload::Domain &) {
3313 AssertThrow(false,
3314 ExcMessage("Payload inv_vmult disabled because of "
3315 "incompatible solver/preconditioner choice."));
3316 };
3317
3318 return return_op;
3319 }
3320 } // namespace LinearOperatorImplementation
3321 } // namespace internal
3322
3323 template <>
3324 void
3325 SparseMatrix::set<TrilinosScalar>(const size_type row,
3326 const size_type n_cols,
3327 const size_type *col_indices,
3328 const TrilinosScalar *values,
3329 const bool elide_zero_values);
3330# endif // DOXYGEN
3331
3332} /* namespace TrilinosWrappers */
3333
3334
3336
3337
3338# endif // DEAL_II_WITH_TRILINOS
3339
3340
3341/*----------------------- trilinos_sparse_matrix.h --------------------*/
3342
3343#endif
3344/*----------------------- 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