Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_sparse_matrix.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_trilinos_sparse_matrix_h
17# define dealii_trilinos_sparse_matrix_h
18
19
20# include <deal.II/base/config.h>
21
22# ifdef DEAL_II_WITH_TRILINOS
23
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# include <mpi.h>
45
46# include <cmath>
47# include <iterator>
48# include <memory>
49# include <type_traits>
50# include <vector>
51
53
54// forward declarations
55# ifndef DOXYGEN
56template <typename MatrixType>
57class BlockMatrixBase;
58
59template <typename number>
60class SparseMatrix;
61class SparsityPattern;
63
64namespace TrilinosWrappers
65{
66 class SparseMatrix;
67 class SparsityPattern;
68
69 namespace SparseMatrixIterators
70 {
71 template <bool Constness>
72 class Iterator;
73 }
74} // namespace TrilinosWrappers
75# endif
76
77namespace TrilinosWrappers
78{
83 {
88
93 std::size_t,
94 std::size_t,
95 std::size_t,
96 << "You tried to access row " << arg1
97 << " of a distributed sparsity pattern, "
98 << " but only rows " << arg2 << " through " << arg3
99 << " are stored locally and can be accessed.");
100
109 {
110 public:
115
120 const size_type row,
121 const size_type index);
122
127 row() const;
128
133 index() const;
134
139 column() const;
140
141 protected:
153
158
164 void
166
179 std::shared_ptr<std::vector<size_type>> colnum_cache;
180
184 std::shared_ptr<std::vector<TrilinosScalar>> value_cache;
185 };
186
197 template <bool Constess>
198 class Accessor : public AccessorBase
199 {
204 value() const;
205
211 };
212
213
214
218 template <>
219 class Accessor<true> : public AccessorBase
220 {
221 public:
227
233
238 template <bool Other>
240
245 value() const;
246
247 private:
248 // Make iterator class a friend.
249 template <bool>
250 friend class Iterator;
251 };
252
256 template <>
257 class Accessor<false> : public AccessorBase
258 {
259 class Reference
260 {
261 public:
265 Reference(const Accessor<false> &accessor);
266
270 operator TrilinosScalar() const;
271
275 const Reference &
276 operator=(const TrilinosScalar n) const;
277
281 const Reference &
283
287 const Reference &
289
293 const Reference &
295
299 const Reference &
301
302 private:
308 };
309
310 public:
316
322
326 Reference
327 value() const;
328
329 private:
330 // Make iterator class a friend.
331 template <bool>
332 friend class Iterator;
333
334 // Make Reference object a friend.
335 friend class Reference;
336 };
337
351 template <bool Constness>
353 {
354 public:
359
365
371
377
382 Iterator(MatrixType *matrix, const size_type row, const size_type index);
383
387 template <bool Other>
389
395
401
405 const Accessor<Constness> &
406 operator*() const;
407
411 const Accessor<Constness> *
412 operator->() const;
413
418 bool
420
424 bool
426
432 bool
433 operator<(const Iterator<Constness> &) const;
434
438 bool
440
445 size_type,
446 size_type,
447 << "Attempt to access element " << arg2 << " of row "
448 << arg1 << " which doesn't have that many elements.");
449
450 private:
455
456 template <bool Other>
457 friend class Iterator;
458 };
459
460 } // namespace SparseMatrixIterators
461} // namespace TrilinosWrappers
462
464
465namespace std
466{
467 template <bool Constness>
468 struct iterator_traits<
470 {
471 using iterator_category = forward_iterator_tag;
473 typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
474 Constness>::value_type;
476 typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
477 Constness>::difference_type;
478 };
479} // namespace std
480
482
483
484namespace TrilinosWrappers
485{
547 {
548 public:
553
558 std::size_t,
559 << "You tried to access row " << arg1
560 << " of a non-contiguous locally owned row set."
561 << " The row " << arg1
562 << " is not stored locally and can't be accessed.");
563
571 struct Traits
572 {
577 static const bool zero_addition_can_be_elided = true;
578 };
579
584
589
594
602 SparseMatrix();
603
612 const size_type n,
613 const unsigned int n_max_entries_per_row);
614
623 const size_type n,
624 const std::vector<unsigned int> &n_entries_per_row);
625
629 SparseMatrix(const SparsityPattern &InputSparsityPattern);
630
635 SparseMatrix(SparseMatrix &&other) noexcept;
636
640 SparseMatrix(const SparseMatrix &) = delete;
641
646 operator=(const SparseMatrix &) = delete;
647
651 virtual ~SparseMatrix() override = default;
652
668 template <typename SparsityPatternType>
669 void
670 reinit(const SparsityPatternType &sparsity_pattern);
671
684 void
685 reinit(const SparsityPattern &sparsity_pattern);
686
695 void
696 reinit(const SparseMatrix &sparse_matrix);
697
718 template <typename number>
719 void
720 reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
721 const double drop_tolerance = 1e-13,
722 const bool copy_values = true,
723 const ::SparsityPattern * use_this_sparsity = nullptr);
724
730 void
731 reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
733
750 SparseMatrix(const IndexSet & parallel_partitioning,
751 const MPI_Comm & communicator = MPI_COMM_WORLD,
752 const unsigned int n_max_entries_per_row = 0);
753
761 SparseMatrix(const IndexSet & parallel_partitioning,
762 const MPI_Comm & communicator,
763 const std::vector<unsigned int> &n_entries_per_row);
764
779 SparseMatrix(const IndexSet &row_parallel_partitioning,
780 const IndexSet &col_parallel_partitioning,
781 const MPI_Comm &communicator = MPI_COMM_WORLD,
782 const size_type n_max_entries_per_row = 0);
783
798 SparseMatrix(const IndexSet & row_parallel_partitioning,
799 const IndexSet & col_parallel_partitioning,
800 const MPI_Comm & communicator,
801 const std::vector<unsigned int> &n_entries_per_row);
802
823 template <typename SparsityPatternType>
824 void
825 reinit(const IndexSet & parallel_partitioning,
826 const SparsityPatternType &sparsity_pattern,
827 const MPI_Comm & communicator = MPI_COMM_WORLD,
828 const bool exchange_data = false);
829
842 template <typename SparsityPatternType>
843 typename std::enable_if<
844 !std::is_same<SparsityPatternType,
845 ::SparseMatrix<double>>::value>::type
846 reinit(const IndexSet & row_parallel_partitioning,
847 const IndexSet & col_parallel_partitioning,
848 const SparsityPatternType &sparsity_pattern,
849 const MPI_Comm & communicator = MPI_COMM_WORLD,
850 const bool exchange_data = false);
851
868 template <typename number>
869 void
870 reinit(const IndexSet & parallel_partitioning,
871 const ::SparseMatrix<number> &dealii_sparse_matrix,
872 const MPI_Comm & communicator = MPI_COMM_WORLD,
873 const double drop_tolerance = 1e-13,
874 const bool copy_values = true,
875 const ::SparsityPattern * use_this_sparsity = nullptr);
876
890 template <typename number>
891 void
892 reinit(const IndexSet & row_parallel_partitioning,
893 const IndexSet & col_parallel_partitioning,
894 const ::SparseMatrix<number> &dealii_sparse_matrix,
895 const MPI_Comm & communicator = MPI_COMM_WORLD,
896 const double drop_tolerance = 1e-13,
897 const bool copy_values = true,
898 const ::SparsityPattern * use_this_sparsity = nullptr);
900
904
909 m() const;
910
915 n() const;
916
925 unsigned int
926 local_size() const;
927
936 std::pair<size_type, size_type>
937 local_range() const;
938
943 bool
944 in_local_range(const size_type index) const;
945
950 std::uint64_t
952
956 unsigned int
957 row_length(const size_type row) const;
958
965 bool
967
974 memory_consumption() const;
975
980 get_mpi_communicator() const;
981
983
987
998 operator=(const double d);
999
1007 void
1008 clear();
1009
1037 void
1039
1061 void
1062 set(const size_type i, const size_type j, const TrilinosScalar value);
1063
1096 void
1097 set(const std::vector<size_type> & indices,
1098 const FullMatrix<TrilinosScalar> &full_matrix,
1099 const bool elide_zero_values = false);
1100
1106 void
1107 set(const std::vector<size_type> & row_indices,
1108 const std::vector<size_type> & col_indices,
1109 const FullMatrix<TrilinosScalar> &full_matrix,
1110 const bool elide_zero_values = false);
1111
1139 void
1140 set(const size_type row,
1141 const std::vector<size_type> & col_indices,
1142 const std::vector<TrilinosScalar> &values,
1143 const bool elide_zero_values = false);
1144
1172 template <typename Number>
1173 void
1174 set(const size_type row,
1175 const size_type n_cols,
1176 const size_type *col_indices,
1177 const Number * values,
1178 const bool elide_zero_values = false);
1179
1189 void
1190 add(const size_type i, const size_type j, const TrilinosScalar value);
1191
1210 void
1211 add(const std::vector<size_type> & indices,
1212 const FullMatrix<TrilinosScalar> &full_matrix,
1213 const bool elide_zero_values = true);
1214
1220 void
1221 add(const std::vector<size_type> & row_indices,
1222 const std::vector<size_type> & col_indices,
1223 const FullMatrix<TrilinosScalar> &full_matrix,
1224 const bool elide_zero_values = true);
1225
1239 void
1240 add(const size_type row,
1241 const std::vector<size_type> & col_indices,
1242 const std::vector<TrilinosScalar> &values,
1243 const bool elide_zero_values = true);
1244
1258 void
1259 add(const size_type row,
1260 const size_type n_cols,
1261 const size_type * col_indices,
1262 const TrilinosScalar *values,
1263 const bool elide_zero_values = true,
1264 const bool col_indices_are_sorted = false);
1265
1269 SparseMatrix &
1270 operator*=(const TrilinosScalar factor);
1271
1275 SparseMatrix &
1276 operator/=(const TrilinosScalar factor);
1277
1281 void
1282 copy_from(const SparseMatrix &source);
1283
1291 void
1292 add(const TrilinosScalar factor, const SparseMatrix &matrix);
1293
1320 void
1321 clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1322
1343 void
1344 clear_rows(const std::vector<size_type> &rows,
1345 const TrilinosScalar new_diag_value = 0);
1346
1356 void
1357 transpose();
1358
1360
1364
1374 operator()(const size_type i, const size_type j) const;
1375
1393 el(const size_type i, const size_type j) const;
1394
1402 diag_element(const size_type i) const;
1403
1405
1409
1440 template <typename VectorType>
1441 typename std::enable_if<std::is_same<typename VectorType::value_type,
1442 TrilinosScalar>::value>::type
1443 vmult(VectorType &dst, const VectorType &src) const;
1444
1451 template <typename VectorType>
1452 typename std::enable_if<!std::is_same<typename VectorType::value_type,
1453 TrilinosScalar>::value>::type
1454 vmult(VectorType &dst, const VectorType &src) const;
1455
1470 template <typename VectorType>
1471 typename std::enable_if<std::is_same<typename VectorType::value_type,
1472 TrilinosScalar>::value>::type
1473 Tvmult(VectorType &dst, const VectorType &src) const;
1474
1481 template <typename VectorType>
1482 typename std::enable_if<!std::is_same<typename VectorType::value_type,
1483 TrilinosScalar>::value>::type
1484 Tvmult(VectorType &dst, const VectorType &src) const;
1485
1495 template <typename VectorType>
1496 void
1497 vmult_add(VectorType &dst, const VectorType &src) const;
1498
1509 template <typename VectorType>
1510 void
1511 Tvmult_add(VectorType &dst, const VectorType &src) const;
1512
1535 matrix_norm_square(const MPI::Vector &v) const;
1536
1557 matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1558
1576 residual(MPI::Vector & dst,
1577 const MPI::Vector &x,
1578 const MPI::Vector &b) const;
1579
1594 void
1596 const SparseMatrix &B,
1597 const MPI::Vector & V = MPI::Vector()) const;
1598
1599
1616 void
1618 const SparseMatrix &B,
1619 const MPI::Vector & V = MPI::Vector()) const;
1620
1622
1626
1635 l1_norm() const;
1636
1646 linfty_norm() const;
1647
1653 frobenius_norm() const;
1654
1656
1660
1665 const Epetra_CrsMatrix &
1667
1672 const Epetra_CrsGraph &
1674
1676
1681
1686 IndexSet
1688
1694 IndexSet
1696
1698
1703
1723 begin() const;
1724
1728 iterator
1730
1736 end() const;
1737
1741 iterator
1743
1773 begin(const size_type r) const;
1774
1778 iterator
1780
1791 end(const size_type r) const;
1792
1796 iterator
1797 end(const size_type r);
1798
1800
1804
1810 void
1811 write_ascii();
1812
1820 void
1821 print(std::ostream &out,
1822 const bool write_extended_trilinos_info = false) const;
1823
1825
1834 int,
1835 << "An error with error number " << arg1
1836 << " occurred while calling a Trilinos function. "
1837 "\n\n"
1838 "For historical reasons, many Trilinos functions express "
1839 "errors by returning specific integer values to indicate "
1840 "certain errors. Unfortunately, different Trilinos functions "
1841 "often use the same integer values for different kinds of "
1842 "errors, and in most cases it is also not documented what "
1843 "each error code actually means. As a consequence, it is often "
1844 "difficult to say what a particular error (in this case, "
1845 "the error with integer code '"
1846 << arg1
1847 << "') represents and how one should fix a code to avoid it. "
1848 "The best one can often do is to look up the call stack to "
1849 "see which deal.II function generated the error, and which "
1850 "Trilinos function the error code had originated from; "
1851 "then look up the Trilinos source code of that function (for "
1852 "example on github) to see what code path set that error "
1853 "code. Short of going through all of that, the only other "
1854 "option is to guess the cause of the error from "
1855 "the context in which the error appeared.");
1856
1857
1862 size_type,
1863 size_type,
1864 << "The entry with index <" << arg1 << ',' << arg2
1865 << "> does not exist.");
1866
1871 "You are attempting an operation on two matrices that "
1872 "are the same object, but the operation requires that the "
1873 "two objects are in fact different.");
1874
1879
1884 size_type,
1885 size_type,
1886 size_type,
1887 size_type,
1888 << "You tried to access element (" << arg1 << '/' << arg2
1889 << ')'
1890 << " of a distributed matrix, but only rows in range ["
1891 << arg3 << ',' << arg4
1892 << "] are stored locally and can be accessed.");
1893
1898 size_type,
1899 size_type,
1900 << "You tried to access element (" << arg1 << '/' << arg2
1901 << ')' << " of a sparse matrix, but it appears to not"
1902 << " exist in the Trilinos sparsity pattern.");
1904
1905
1906
1907 protected:
1918 void
1920
1926 void
1928
1929
1930
1931 private:
1936 std::unique_ptr<Epetra_Map> column_space_map;
1937
1943 std::unique_ptr<Epetra_FECrsMatrix> matrix;
1944
1950 std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1951
1955 std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1956
1968 Epetra_CombineMode last_action;
1969
1975
1976 // To allow calling protected prepare_add() and prepare_set().
1977 friend class BlockMatrixBase<SparseMatrix>;
1978 };
1979
1980
1981
1982 // forwards declarations
1983 class SolverBase;
1984 class PreconditionBase;
1985
1986 namespace internal
1987 {
1988 inline void
1989 check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1990 const Epetra_MultiVector &src,
1991 const Epetra_MultiVector &dst,
1992 const bool transpose)
1993 {
1994 if (transpose == false)
1995 {
1996 Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
1997 ExcMessage(
1998 "Column map of matrix does not fit with vector map!"));
1999 Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
2000 ExcMessage("Row map of matrix does not fit with vector map!"));
2001 }
2002 else
2003 {
2004 Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2005 ExcMessage(
2006 "Column map of matrix does not fit with vector map!"));
2007 Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2008 ExcMessage("Row map of matrix does not fit with vector map!"));
2009 }
2010 (void)mtrx; // removes -Wunused-variable in optimized mode
2011 (void)src;
2012 (void)dst;
2013 }
2014
2015 inline void
2017 const Epetra_MultiVector &src,
2018 const Epetra_MultiVector &dst,
2019 const bool transpose)
2020 {
2021 if (transpose == false)
2022 {
2023 Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2024 ExcMessage(
2025 "Column map of operator does not fit with vector map!"));
2026 Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2027 ExcMessage(
2028 "Row map of operator does not fit with vector map!"));
2029 }
2030 else
2031 {
2032 Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2033 ExcMessage(
2034 "Column map of operator does not fit with vector map!"));
2035 Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2036 ExcMessage(
2037 "Row map of operator does not fit with vector map!"));
2038 }
2039 (void)op; // removes -Wunused-variable in optimized mode
2040 (void)src;
2041 (void)dst;
2042 }
2043
2044 namespace LinearOperatorImplementation
2045 {
2066 {
2067 public:
2071 using VectorType = Epetra_MultiVector;
2072
2077
2082
2087
2096
2100 TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2102
2106 TrilinosPayload(const TrilinosPayload & payload_exemplar,
2108
2113 const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2114 const TrilinosWrappers::PreconditionBase &preconditioner);
2115
2120 const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2121 const TrilinosWrappers::PreconditionBase &preconditioner);
2122
2127 const TrilinosPayload & payload_exemplar,
2128 const TrilinosWrappers::PreconditionBase &preconditioner);
2129
2133 TrilinosPayload(const TrilinosPayload &payload);
2134
2142 TrilinosPayload(const TrilinosPayload &first_op,
2143 const TrilinosPayload &second_op);
2144
2148 virtual ~TrilinosPayload() override = default;
2149
2154 identity_payload() const;
2155
2160 null_payload() const;
2161
2166 transpose_payload() const;
2167
2184 template <typename Solver, typename Preconditioner>
2185 typename std::enable_if<
2186 std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2187 std::is_base_of<TrilinosWrappers::PreconditionBase,
2188 Preconditioner>::value,
2189 TrilinosPayload>::type
2190 inverse_payload(Solver &, const Preconditioner &) const;
2191
2209 template <typename Solver, typename Preconditioner>
2210 typename std::enable_if<
2211 !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2212 std::is_base_of<TrilinosWrappers::PreconditionBase,
2213 Preconditioner>::value),
2214 TrilinosPayload>::type
2215 inverse_payload(Solver &, const Preconditioner &) const;
2216
2218
2223
2229 IndexSet
2231
2237 IndexSet
2239
2243 MPI_Comm
2244 get_mpi_communicator() const;
2245
2252 void
2253 transpose();
2254
2262 std::function<void(VectorType &, const VectorType &)> vmult;
2263
2271 std::function<void(VectorType &, const VectorType &)> Tvmult;
2272
2281 std::function<void(VectorType &, const VectorType &)> inv_vmult;
2282
2291 std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2292
2294
2299
2306 virtual bool
2307 UseTranspose() const override;
2308
2324 virtual int
2325 SetUseTranspose(bool UseTranspose) override;
2326
2338 virtual int
2339 Apply(const VectorType &X, VectorType &Y) const override;
2340
2359 virtual int
2360 ApplyInverse(const VectorType &Y, VectorType &X) const override;
2362
2367
2374 virtual const char *
2375 Label() const override;
2376
2384 virtual const Epetra_Comm &
2385 Comm() const override;
2386
2394 virtual const Epetra_Map &
2395 OperatorDomainMap() const override;
2396
2405 virtual const Epetra_Map &
2406 OperatorRangeMap() const override;
2408
2409 private:
2419 template <typename EpetraOpType>
2420 TrilinosPayload(EpetraOpType & op,
2421 const bool supports_inverse_operations,
2422 const bool use_transpose,
2423 const MPI_Comm &mpi_communicator,
2426
2432
2437 Epetra_MpiComm communicator;
2438
2443 Epetra_Map domain_map;
2444
2449 Epetra_Map range_map;
2450
2459 virtual bool
2460 HasNormInf() const override;
2461
2469 virtual double
2470 NormInf() const override;
2471 };
2472
2478 operator+(const TrilinosPayload &first_op,
2479 const TrilinosPayload &second_op);
2480
2486 operator*(const TrilinosPayload &first_op,
2487 const TrilinosPayload &second_op);
2488
2489 } // namespace LinearOperatorImplementation
2490 } /* namespace internal */
2491
2492
2493
2494 // ----------------------- inline and template functions --------------------
2495
2496# ifndef DOXYGEN
2497
2498 namespace SparseMatrixIterators
2499 {
2501 size_type row,
2502 size_type index)
2503 : matrix(matrix)
2504 , a_row(row)
2505 , a_index(index)
2506 {
2507 visit_present_row();
2508 }
2509
2510
2512 AccessorBase::row() const
2513 {
2514 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2515 return a_row;
2516 }
2517
2518
2520 AccessorBase::column() const
2521 {
2522 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2523 return (*colnum_cache)[a_index];
2524 }
2525
2526
2528 AccessorBase::index() const
2529 {
2530 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2531 return a_index;
2532 }
2533
2534
2535 inline Accessor<true>::Accessor(MatrixType * matrix,
2536 const size_type row,
2537 const size_type index)
2538 : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2539 {}
2540
2541
2542 template <bool Other>
2543 inline Accessor<true>::Accessor(const Accessor<Other> &other)
2544 : AccessorBase(other)
2545 {}
2546
2547
2548 inline TrilinosScalar
2549 Accessor<true>::value() const
2550 {
2551 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2552 return (*value_cache)[a_index];
2553 }
2554
2555
2556 inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2557 : accessor(const_cast<Accessor<false> &>(acc))
2558 {}
2559
2560
2561 inline Accessor<false>::Reference::operator TrilinosScalar() const
2562 {
2563 return (*accessor.value_cache)[accessor.a_index];
2564 }
2565
2566 inline const Accessor<false>::Reference &
2567 Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2568 {
2569 (*accessor.value_cache)[accessor.a_index] = n;
2570 accessor.matrix->set(accessor.row(),
2571 accessor.column(),
2572 static_cast<TrilinosScalar>(*this));
2573 return *this;
2574 }
2575
2576
2577 inline const Accessor<false>::Reference &
2578 Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2579 {
2580 (*accessor.value_cache)[accessor.a_index] += n;
2581 accessor.matrix->set(accessor.row(),
2582 accessor.column(),
2583 static_cast<TrilinosScalar>(*this));
2584 return *this;
2585 }
2586
2587
2588 inline const Accessor<false>::Reference &
2589 Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2590 {
2591 (*accessor.value_cache)[accessor.a_index] -= n;
2592 accessor.matrix->set(accessor.row(),
2593 accessor.column(),
2594 static_cast<TrilinosScalar>(*this));
2595 return *this;
2596 }
2597
2598
2599 inline const Accessor<false>::Reference &
2600 Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2601 {
2602 (*accessor.value_cache)[accessor.a_index] *= n;
2603 accessor.matrix->set(accessor.row(),
2604 accessor.column(),
2605 static_cast<TrilinosScalar>(*this));
2606 return *this;
2607 }
2608
2609
2610 inline const Accessor<false>::Reference &
2611 Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2612 {
2613 (*accessor.value_cache)[accessor.a_index] /= n;
2614 accessor.matrix->set(accessor.row(),
2615 accessor.column(),
2616 static_cast<TrilinosScalar>(*this));
2617 return *this;
2618 }
2619
2620
2621 inline Accessor<false>::Accessor(MatrixType * matrix,
2622 const size_type row,
2623 const size_type index)
2624 : AccessorBase(matrix, row, index)
2625 {}
2626
2627
2628 inline Accessor<false>::Reference
2629 Accessor<false>::value() const
2630 {
2631 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2632 return {*this};
2633 }
2634
2635
2636
2637 template <bool Constness>
2638 inline Iterator<Constness>::Iterator(MatrixType * matrix,
2639 const size_type row,
2640 const size_type index)
2641 : accessor(matrix, row, index)
2642 {}
2643
2644
2645 template <bool Constness>
2646 template <bool Other>
2647 inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2648 : accessor(other.accessor)
2649 {}
2650
2651
2652 template <bool Constness>
2653 inline Iterator<Constness> &
2655 {
2656 Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2657
2658 ++accessor.a_index;
2659
2660 // If at end of line: do one
2661 // step, then cycle until we
2662 // find a row with a nonzero
2663 // number of entries.
2664 if (accessor.a_index >= accessor.colnum_cache->size())
2665 {
2666 accessor.a_index = 0;
2667 ++accessor.a_row;
2668
2669 while ((accessor.a_row < accessor.matrix->m()) &&
2670 ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2671 (accessor.matrix->row_length(accessor.a_row) == 0)))
2672 ++accessor.a_row;
2673
2674 accessor.visit_present_row();
2675 }
2676 return *this;
2677 }
2678
2679
2680 template <bool Constness>
2681 inline Iterator<Constness>
2683 {
2684 const Iterator<Constness> old_state = *this;
2685 ++(*this);
2686 return old_state;
2687 }
2688
2689
2690
2691 template <bool Constness>
2692 inline const Accessor<Constness> &
2694 {
2695 return accessor;
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 bool
2711 Iterator<Constness>::operator==(const Iterator<Constness> &other) const
2712 {
2713 return (accessor.a_row == other.accessor.a_row &&
2714 accessor.a_index == other.accessor.a_index);
2715 }
2716
2717
2718
2719 template <bool Constness>
2720 inline bool
2721 Iterator<Constness>::operator!=(const Iterator<Constness> &other) const
2722 {
2723 return !(*this == other);
2724 }
2725
2726
2727
2728 template <bool Constness>
2729 inline bool
2730 Iterator<Constness>::operator<(const Iterator<Constness> &other) const
2731 {
2732 return (accessor.row() < other.accessor.row() ||
2733 (accessor.row() == other.accessor.row() &&
2734 accessor.index() < other.accessor.index()));
2735 }
2736
2737
2738 template <bool Constness>
2739 inline bool
2740 Iterator<Constness>::operator>(const Iterator<Constness> &other) const
2741 {
2742 return (other < *this);
2743 }
2744
2745 } // namespace SparseMatrixIterators
2746
2747
2748
2750 SparseMatrix::begin() const
2751 {
2752 return begin(0);
2753 }
2754
2755
2756
2758 SparseMatrix::end() const
2759 {
2760 return const_iterator(this, m(), 0);
2761 }
2762
2763
2764
2766 SparseMatrix::begin(const size_type r) const
2767 {
2768 AssertIndexRange(r, m());
2769 if (in_local_range(r) && (row_length(r) > 0))
2770 return const_iterator(this, r, 0);
2771 else
2772 return end(r);
2773 }
2774
2775
2776
2778 SparseMatrix::end(const size_type r) const
2779 {
2780 AssertIndexRange(r, m());
2781
2782 // place the iterator on the first entry
2783 // past this line, or at the end of the
2784 // matrix
2785 for (size_type i = r + 1; i < m(); ++i)
2786 if (in_local_range(i) && (row_length(i) > 0))
2787 return const_iterator(this, i, 0);
2788
2789 // if there is no such line, then take the
2790 // end iterator of the matrix
2791 return end();
2792 }
2793
2794
2795
2798 {
2799 return begin(0);
2800 }
2801
2802
2803
2806 {
2807 return iterator(this, m(), 0);
2808 }
2809
2810
2811
2814 {
2815 AssertIndexRange(r, m());
2816 if (in_local_range(r) && (row_length(r) > 0))
2817 return iterator(this, r, 0);
2818 else
2819 return end(r);
2820 }
2821
2822
2823
2826 {
2827 AssertIndexRange(r, m());
2828
2829 // place the iterator on the first entry
2830 // past this line, or at the end of the
2831 // matrix
2832 for (size_type i = r + 1; i < m(); ++i)
2833 if (in_local_range(i) && (row_length(i) > 0))
2834 return iterator(this, i, 0);
2835
2836 // if there is no such line, then take the
2837 // end iterator of the matrix
2838 return end();
2839 }
2840
2841
2842
2843 inline bool
2844 SparseMatrix::in_local_range(const size_type index) const
2845 {
2847# ifndef DEAL_II_WITH_64BIT_INDICES
2848 begin = matrix->RowMap().MinMyGID();
2849 end = matrix->RowMap().MaxMyGID() + 1;
2850# else
2851 begin = matrix->RowMap().MinMyGID64();
2852 end = matrix->RowMap().MaxMyGID64() + 1;
2853# endif
2854
2855 return ((index >= static_cast<size_type>(begin)) &&
2856 (index < static_cast<size_type>(end)));
2857 }
2858
2859
2860
2861 inline bool
2863 {
2864 return compressed;
2865 }
2866
2867
2868
2869 // Inline the set() and add() functions, since they will be called
2870 // frequently, and the compiler can optimize away some unnecessary loops
2871 // when the sizes are given at compile time.
2872 template <>
2873 void
2874 SparseMatrix::set<TrilinosScalar>(const size_type row,
2875 const size_type n_cols,
2876 const size_type * col_indices,
2877 const TrilinosScalar *values,
2878 const bool elide_zero_values);
2879
2880
2881
2882 template <typename Number>
2883 void
2884 SparseMatrix::set(const size_type row,
2885 const size_type n_cols,
2886 const size_type *col_indices,
2887 const Number * values,
2888 const bool elide_zero_values)
2889 {
2890 std::vector<TrilinosScalar> trilinos_values(n_cols);
2891 std::copy(values, values + n_cols, trilinos_values.begin());
2892 this->set(
2893 row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2894 }
2895
2896
2897
2898 inline void
2900 const size_type j,
2901 const TrilinosScalar value)
2902 {
2903 AssertIsFinite(value);
2904
2905 set(i, 1, &j, &value, false);
2906 }
2907
2908
2909
2910 inline void
2911 SparseMatrix::set(const std::vector<size_type> & indices,
2913 const bool elide_zero_values)
2914 {
2915 Assert(indices.size() == values.m(),
2916 ExcDimensionMismatch(indices.size(), values.m()));
2917 Assert(values.m() == values.n(), ExcNotQuadratic());
2918
2919 for (size_type i = 0; i < indices.size(); ++i)
2920 set(indices[i],
2921 indices.size(),
2922 indices.data(),
2923 &values(i, 0),
2924 elide_zero_values);
2925 }
2926
2927
2928
2929 inline void
2931 const size_type j,
2932 const TrilinosScalar value)
2933 {
2934 AssertIsFinite(value);
2935
2936 if (value == 0)
2937 {
2938 // we have to check after Insert/Add in any case to be consistent
2939 // with the MPI communication model, but we can save some
2940 // work if the addend is zero. However, these actions are done in case
2941 // we pass on to the other function.
2942
2943 // TODO: fix this (do not run compress here, but fail)
2944 if (last_action == Insert)
2945 {
2946 const int ierr = matrix->GlobalAssemble(*column_space_map,
2947 matrix->RowMap(),
2948 false);
2949
2950 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2951 }
2952
2953 last_action = Add;
2954
2955 return;
2956 }
2957 else
2958 add(i, 1, &j, &value, false);
2959 }
2960
2961
2962
2963 // inline "simple" functions that are called frequently and do only involve
2964 // a call to some Trilinos function.
2966 SparseMatrix::m() const
2967 {
2968# ifndef DEAL_II_WITH_64BIT_INDICES
2969 return matrix->NumGlobalRows();
2970# else
2971 return matrix->NumGlobalRows64();
2972# endif
2973 }
2974
2975
2976
2978 SparseMatrix::n() const
2979 {
2980 // If the matrix structure has not been fixed (i.e., we did not have a
2981 // sparsity pattern), it does not know about the number of columns so we
2982 // must always take this from the additional column space map
2983 Assert(column_space_map.get() != nullptr, ExcInternalError());
2984 return n_global_elements(*column_space_map);
2985 }
2986
2987
2988
2989 inline unsigned int
2991 {
2992 return matrix->NumMyRows();
2993 }
2994
2995
2996
2997 inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
2999 {
3001# ifndef DEAL_II_WITH_64BIT_INDICES
3002 begin = matrix->RowMap().MinMyGID();
3003 end = matrix->RowMap().MaxMyGID() + 1;
3004# else
3005 begin = matrix->RowMap().MinMyGID64();
3006 end = matrix->RowMap().MaxMyGID64() + 1;
3007# endif
3008
3009 return std::make_pair(begin, end);
3010 }
3011
3012
3013
3014 inline std::uint64_t
3016 {
3017 // Trilinos uses 64bit functions internally for attribute access, which
3018 // return `long long`. They also offer 32bit variants that return `int`,
3019 // however those call the 64bit version and convert the values to 32bit.
3020 // There is no necessity in using the 32bit versions at all.
3021 return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
3022 }
3023
3024
3025
3026 template <typename SparsityPatternType>
3027 inline void
3028 SparseMatrix::reinit(const IndexSet & parallel_partitioning,
3029 const SparsityPatternType &sparsity_pattern,
3030 const MPI_Comm & communicator,
3031 const bool exchange_data)
3032 {
3033 reinit(parallel_partitioning,
3034 parallel_partitioning,
3035 sparsity_pattern,
3036 communicator,
3037 exchange_data);
3038 }
3039
3040
3041
3042 template <typename number>
3043 inline void
3044 SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3045 const ::SparseMatrix<number> &sparse_matrix,
3046 const MPI_Comm & communicator,
3047 const double drop_tolerance,
3048 const bool copy_values,
3049 const ::SparsityPattern * use_this_sparsity)
3050 {
3051 Epetra_Map map =
3052 parallel_partitioning.make_trilinos_map(communicator, false);
3053 reinit(parallel_partitioning,
3054 parallel_partitioning,
3055 sparse_matrix,
3056 drop_tolerance,
3057 copy_values,
3058 use_this_sparsity);
3059 }
3060
3061
3062
3063 inline const Epetra_CrsMatrix &
3065 {
3066 return static_cast<const Epetra_CrsMatrix &>(*matrix);
3067 }
3068
3069
3070
3071 inline const Epetra_CrsGraph &
3073 {
3074 return matrix->Graph();
3075 }
3076
3077
3078
3079 inline IndexSet
3081 {
3082 return IndexSet(matrix->DomainMap());
3083 }
3084
3085
3086
3087 inline IndexSet
3089 {
3090 return IndexSet(matrix->RangeMap());
3091 }
3092
3093
3094
3095 inline void
3097 {
3098 // nothing to do here
3099 }
3100
3101
3102
3103 inline void
3105 {
3106 // nothing to do here
3107 }
3108
3109
3110 namespace internal
3111 {
3112 namespace LinearOperatorImplementation
3113 {
3114 template <typename EpetraOpType>
3115 TrilinosPayload::TrilinosPayload(
3116 EpetraOpType & op,
3117 const bool supports_inverse_operations,
3118 const bool use_transpose,
3119 const MPI_Comm &mpi_communicator,
3120 const IndexSet &locally_owned_domain_indices,
3121 const IndexSet &locally_owned_range_indices)
3122 : use_transpose(use_transpose)
3123 , communicator(mpi_communicator)
3124 , domain_map(
3125 locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3126 , range_map(
3127 locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3128 {
3129 vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3130 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3131 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3132 Assert(&tril_src != &tril_dst,
3135 tril_src,
3136 tril_dst,
3137 op.UseTranspose());
3138
3139 const int ierr = op.Apply(tril_src, tril_dst);
3140 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3141 };
3142
3143 Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3144 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3145 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3146 Assert(&tril_src != &tril_dst,
3149 tril_src,
3150 tril_dst,
3151 !op.UseTranspose());
3152
3153 op.SetUseTranspose(!op.UseTranspose());
3154 const int ierr = op.Apply(tril_src, tril_dst);
3155 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3156 op.SetUseTranspose(!op.UseTranspose());
3157 };
3158
3159 if (supports_inverse_operations)
3160 {
3161 inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3162 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3163 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3164 Assert(
3165 &tril_src != &tril_dst,
3168 tril_src,
3169 tril_dst,
3170 !op.UseTranspose());
3171
3172 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3173 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3174 };
3175
3176 inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
3177 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3178 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3179 Assert(
3180 &tril_src != &tril_dst,
3183 tril_src,
3184 tril_dst,
3185 op.UseTranspose());
3186
3187 op.SetUseTranspose(!op.UseTranspose());
3188 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3189 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3190 op.SetUseTranspose(!op.UseTranspose());
3191 };
3192 }
3193 else
3194 {
3195 inv_vmult = [](Domain &, const Range &) {
3196 Assert(false,
3197 ExcMessage(
3198 "Uninitialized TrilinosPayload::inv_vmult called. "
3199 "The operator does not support inverse operations."));
3200 };
3201
3202 inv_Tvmult = [](Range &, const Domain &) {
3203 Assert(false,
3204 ExcMessage(
3205 "Uninitialized TrilinosPayload::inv_Tvmult called. "
3206 "The operator does not support inverse operations."));
3207 };
3208 }
3209 }
3210
3211
3212 template <typename Solver, typename Preconditioner>
3213 typename std::enable_if<
3214 std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3215 std::is_base_of<TrilinosWrappers::PreconditionBase,
3216 Preconditioner>::value,
3217 TrilinosPayload>::type
3219 Solver & solver,
3220 const Preconditioner &preconditioner) const
3221 {
3222 const auto &payload = *this;
3223
3224 TrilinosPayload return_op(payload);
3225
3226 // Capture by copy so the payloads are always valid
3227
3228 return_op.inv_vmult = [payload, &solver, &preconditioner](
3229 TrilinosPayload::Domain & tril_dst,
3230 const TrilinosPayload::Range &tril_src) {
3231 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3232 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3233 Assert(&tril_src != &tril_dst,
3236 tril_src,
3237 tril_dst,
3238 !payload.UseTranspose());
3239 solver.solve(payload, tril_dst, tril_src, preconditioner);
3240 };
3241
3242 return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3243 TrilinosPayload::Range & tril_dst,
3244 const TrilinosPayload::Domain &tril_src) {
3245 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3246 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3247 Assert(&tril_src != &tril_dst,
3250 tril_src,
3251 tril_dst,
3252 payload.UseTranspose());
3253
3254 const_cast<TrilinosPayload &>(payload).transpose();
3255 solver.solve(payload, tril_dst, tril_src, preconditioner);
3256 const_cast<TrilinosPayload &>(payload).transpose();
3257 };
3258
3259 // If the input operator is already setup for transpose operations, then
3260 // we must do similar with its inverse.
3261 if (return_op.UseTranspose() == true)
3262 std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3263
3264 return return_op;
3265 }
3266
3267 template <typename Solver, typename Preconditioner>
3268 typename std::enable_if<
3269 !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3270 std::is_base_of<TrilinosWrappers::PreconditionBase,
3271 Preconditioner>::value),
3272 TrilinosPayload>::type
3273 TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3274 {
3275 TrilinosPayload return_op(*this);
3276
3277 return_op.inv_vmult = [](TrilinosPayload::Domain &,
3278 const TrilinosPayload::Range &) {
3279 AssertThrow(false,
3280 ExcMessage("Payload inv_vmult disabled because of "
3281 "incompatible solver/preconditioner choice."));
3282 };
3283
3284 return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3285 const TrilinosPayload::Domain &) {
3286 AssertThrow(false,
3287 ExcMessage("Payload inv_vmult disabled because of "
3288 "incompatible solver/preconditioner choice."));
3289 };
3290
3291 return return_op;
3292 }
3293 } // namespace LinearOperatorImplementation
3294 } // namespace internal
3295
3296 template <>
3297 void
3298 SparseMatrix::set<TrilinosScalar>(const size_type row,
3299 const size_type n_cols,
3300 const size_type * col_indices,
3301 const TrilinosScalar *values,
3302 const bool elide_zero_values);
3303# endif // DOXYGEN
3304
3305} /* namespace TrilinosWrappers */
3306
3307
3309
3310
3311# endif // DEAL_II_WITH_TRILINOS
3312
3313
3314/*----------------------- trilinos_sparse_matrix.h --------------------*/
3315
3316#endif
3317/*----------------------- trilinos_sparse_matrix.h --------------------*/
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition: index_set.cc:601
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
bool operator>(const Iterator< Constness > &) const
typename Accessor< Constness >::MatrixType MatrixType
Iterator(const Iterator< Other > &other)
bool operator==(const Iterator< Constness > &) const
bool operator!=(const Iterator< Constness > &) const
Iterator(MatrixType *matrix, const size_type row, const size_type index)
bool operator<(const Iterator< Constness > &) const
void set(const size_type i, const size_type j, const TrilinosScalar value)
std::enable_if<!std::is_same< typenameVectorType::value_type, TrilinosScalar >::value >::type vmult(VectorType &dst, const VectorType &src) const
std::unique_ptr< Epetra_Map > column_space_map
std::enable_if< std::is_same< typenameVectorType::value_type, TrilinosScalar >::value >::type vmult(VectorType &dst, const VectorType &src) 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
std::unique_ptr< Epetra_FECrsMatrix > matrix
std::enable_if<!std::is_same< typenameVectorType::value_type, TrilinosScalar >::value >::type Tvmult(VectorType &dst, const VectorType &src) const
::types::global_dof_index size_type
const Epetra_CrsMatrix & trilinos_matrix() const
TrilinosScalar matrix_norm_square(const MPI::Vector &v) 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)
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
void clear_row(const size_type row, 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
void compress(::VectorOperation::values operation)
SparseMatrix & operator/=(const TrilinosScalar factor)
void Tvmult_add(VectorType &dst, const VectorType &src) const
void clear_rows(const std::vector< size_type > &rows, const TrilinosScalar new_diag_value=0)
TrilinosScalar el(const size_type i, const size_type j) const
IndexSet locally_owned_domain_indices() const
TrilinosScalar residual(MPI::Vector &dst, const MPI::Vector &x, const MPI::Vector &b) 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
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)
void reinit(const IndexSet &parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm &communicator=MPI_COMM_WORLD, const bool exchange_data=false)
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)
std::enable_if< std::is_same< typenameVectorType::value_type, TrilinosScalar >::value >::type Tvmult(VectorType &dst, const VectorType &src) const
TrilinosScalar operator()(const size_type i, const size_type j) const
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
std::unique_ptr< Epetra_CrsMatrix > nonlocal_matrix
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
std::enable_if< std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value, TrilinosPayload >::type inverse_payload(Solver &, const Preconditioner &) const
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)
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
virtual int Apply(const VectorType &X, VectorType &Y) const override
std::function< void(VectorType &, const VectorType &)> inv_vmult
std::enable_if<!(std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value), TrilinosPayload >::type inverse_payload(Solver &, const Preconditioner &) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
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:464
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
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)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcAccessToNonPresentElement(size_type arg1, size_type arg2)
#define AssertIsFinite(number)
Definition: exceptions.h:1758
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
static ::ExceptionBase & ExcBeyondEndOfMatrix()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
Expression operator>(const Expression &lhs, const Expression &rhs)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition: exceptions.h:555
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
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)
Definition: exceptions.h:1583
@ matrix
Contents is actually a matrix.
static const char V
SparseMatrix< double > SparseMatrix
types::global_dof_index size_type
Definition: cuda_kernels.h:45
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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 copy(const T *begin, const T *end, U *dest)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
STL namespace.
unsigned int global_dof_index
Definition: types.h:76
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::value_type value_type
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::difference_type difference_type
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
double TrilinosScalar
Definition: types.h:163