Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
trilinos_sparse_matrix.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2023 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
27
36
37# include <Epetra_Comm.h>
38# include <Epetra_CrsGraph.h>
39# include <Epetra_Export.h>
40# include <Epetra_FECrsMatrix.h>
41# include <Epetra_Map.h>
42# include <Epetra_MpiComm.h>
43# include <Epetra_MultiVector.h>
44# include <Epetra_Operator.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 template <bool OtherConstness>
419 bool
421
425 template <bool OtherConstness>
426 bool
428
434 template <bool OtherConstness>
435 bool
436 operator<(const Iterator<OtherConstness> &) const;
437
441 template <bool OtherConstness>
442 bool
444
449 size_type,
450 size_type,
451 << "Attempt to access element " << arg2 << " of row "
452 << arg1 << " which doesn't have that many elements.");
453
454 private:
459
460 template <bool Other>
461 friend class Iterator;
462 };
463
464 } // namespace SparseMatrixIterators
465} // namespace TrilinosWrappers
466
468
469namespace std
470{
471 template <bool Constness>
474 {
475 using iterator_category = forward_iterator_tag;
477 typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
478 Constness>::value_type;
480 typename ::TrilinosWrappers::SparseMatrixIterators::Iterator<
481 Constness>::difference_type;
482 };
483} // namespace std
484
486
487
488namespace TrilinosWrappers
489{
551 {
552 public:
557
562 std::size_t,
563 << "You tried to access row " << arg1
564 << " of a non-contiguous locally owned row set."
565 << " The row " << arg1
566 << " is not stored locally and can't be accessed.");
567
575 struct Traits
576 {
581 static const bool zero_addition_can_be_elided = true;
582 };
583
588
593
598
606 SparseMatrix();
607
616 const size_type n,
617 const unsigned int n_max_entries_per_row);
618
627 const size_type n,
628 const std::vector<unsigned int> &n_entries_per_row);
629
633 SparseMatrix(const SparsityPattern &InputSparsityPattern);
634
639 SparseMatrix(SparseMatrix &&other) noexcept;
640
644 SparseMatrix(const SparseMatrix &) = delete;
645
650 operator=(const SparseMatrix &) = delete;
651
655 virtual ~SparseMatrix() override = default;
656
672 template <typename SparsityPatternType>
673 void
674 reinit(const SparsityPatternType &sparsity_pattern);
675
688 void
689 reinit(const SparsityPattern &sparsity_pattern);
690
699 void
700 reinit(const SparseMatrix &sparse_matrix);
701
722 template <typename number>
723 void
724 reinit(const ::SparseMatrix<number> &dealii_sparse_matrix,
725 const double drop_tolerance = 1e-13,
726 const bool copy_values = true,
727 const ::SparsityPattern * use_this_sparsity = nullptr);
728
734 void
735 reinit(const Epetra_CrsMatrix &input_matrix, const bool copy_values = true);
754 SparseMatrix(const IndexSet & parallel_partitioning,
755 const MPI_Comm communicator = MPI_COMM_WORLD,
756 const unsigned int n_max_entries_per_row = 0);
757
765 SparseMatrix(const IndexSet & parallel_partitioning,
766 const MPI_Comm communicator,
767 const std::vector<unsigned int> &n_entries_per_row);
768
783 SparseMatrix(const IndexSet &row_parallel_partitioning,
784 const IndexSet &col_parallel_partitioning,
785 const MPI_Comm communicator = MPI_COMM_WORLD,
786 const size_type n_max_entries_per_row = 0);
787
802 SparseMatrix(const IndexSet & row_parallel_partitioning,
803 const IndexSet & col_parallel_partitioning,
804 const MPI_Comm communicator,
805 const std::vector<unsigned int> &n_entries_per_row);
806
827 template <typename SparsityPatternType>
828 void
829 reinit(const IndexSet & parallel_partitioning,
830 const SparsityPatternType &sparsity_pattern,
831 const MPI_Comm communicator = MPI_COMM_WORLD,
832 const bool exchange_data = false);
833
846 template <typename SparsityPatternType>
847 std::enable_if_t<
848 !std::is_same<SparsityPatternType, ::SparseMatrix<double>>::value>
849 reinit(const IndexSet & row_parallel_partitioning,
850 const IndexSet & col_parallel_partitioning,
851 const SparsityPatternType &sparsity_pattern,
852 const MPI_Comm communicator = MPI_COMM_WORLD,
853 const bool exchange_data = false);
854
871 template <typename number>
872 void
873 reinit(const IndexSet & parallel_partitioning,
874 const ::SparseMatrix<number> &dealii_sparse_matrix,
875 const MPI_Comm communicator = MPI_COMM_WORLD,
876 const double drop_tolerance = 1e-13,
877 const bool copy_values = true,
878 const ::SparsityPattern * use_this_sparsity = nullptr);
879
893 template <typename number>
894 void
895 reinit(const IndexSet & row_parallel_partitioning,
896 const IndexSet & col_parallel_partitioning,
897 const ::SparseMatrix<number> &dealii_sparse_matrix,
898 const MPI_Comm communicator = MPI_COMM_WORLD,
899 const double drop_tolerance = 1e-13,
900 const bool copy_values = true,
901 const ::SparsityPattern * use_this_sparsity = nullptr);
912 m() const;
913
918 n() const;
919
928 unsigned int
929 local_size() const;
930
939 std::pair<size_type, size_type>
940 local_range() const;
941
946 bool
947 in_local_range(const size_type index) const;
948
953 std::uint64_t
955
959 unsigned int
960 row_length(const size_type row) const;
961
968 bool
970
977 memory_consumption() const;
978
983 get_mpi_communicator() const;
984
1000 SparseMatrix &
1001 operator=(const double d);
1002
1010 void
1011 clear();
1012
1040 void
1042
1064 void
1065 set(const size_type i, const size_type j, const TrilinosScalar value);
1066
1099 void
1100 set(const std::vector<size_type> & indices,
1101 const FullMatrix<TrilinosScalar> &full_matrix,
1102 const bool elide_zero_values = false);
1103
1109 void
1110 set(const std::vector<size_type> & row_indices,
1111 const std::vector<size_type> & col_indices,
1112 const FullMatrix<TrilinosScalar> &full_matrix,
1113 const bool elide_zero_values = false);
1114
1142 void
1143 set(const size_type row,
1144 const std::vector<size_type> & col_indices,
1145 const std::vector<TrilinosScalar> &values,
1146 const bool elide_zero_values = false);
1147
1175 template <typename Number>
1176 void
1177 set(const size_type row,
1178 const size_type n_cols,
1179 const size_type *col_indices,
1180 const Number * values,
1181 const bool elide_zero_values = false);
1182
1192 void
1193 add(const size_type i, const size_type j, const TrilinosScalar value);
1194
1213 void
1214 add(const std::vector<size_type> & indices,
1215 const FullMatrix<TrilinosScalar> &full_matrix,
1216 const bool elide_zero_values = true);
1217
1223 void
1224 add(const std::vector<size_type> & row_indices,
1225 const std::vector<size_type> & col_indices,
1226 const FullMatrix<TrilinosScalar> &full_matrix,
1227 const bool elide_zero_values = true);
1228
1242 void
1243 add(const size_type row,
1244 const std::vector<size_type> & col_indices,
1245 const std::vector<TrilinosScalar> &values,
1246 const bool elide_zero_values = true);
1247
1261 void
1262 add(const size_type row,
1263 const size_type n_cols,
1264 const size_type * col_indices,
1265 const TrilinosScalar *values,
1266 const bool elide_zero_values = true,
1267 const bool col_indices_are_sorted = false);
1268
1272 SparseMatrix &
1273 operator*=(const TrilinosScalar factor);
1274
1278 SparseMatrix &
1279 operator/=(const TrilinosScalar factor);
1280
1284 void
1285 copy_from(const SparseMatrix &source);
1286
1294 void
1295 add(const TrilinosScalar factor, const SparseMatrix &matrix);
1296
1323 void
1324 clear_row(const size_type row, const TrilinosScalar new_diag_value = 0);
1325
1346 void
1347 clear_rows(const std::vector<size_type> &rows,
1348 const TrilinosScalar new_diag_value = 0);
1349
1359 void
1360 transpose();
1361
1377 operator()(const size_type i, const size_type j) const;
1378
1396 el(const size_type i, const size_type j) const;
1397
1405 diag_element(const size_type i) const;
1406
1443 template <typename VectorType>
1444 std::enable_if_t<
1445 std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1446 vmult(VectorType &dst, const VectorType &src) const;
1447
1454 template <typename VectorType>
1455 std::enable_if_t<
1456 !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1457 vmult(VectorType &dst, const VectorType &src) const;
1458
1473 template <typename VectorType>
1474 std::enable_if_t<
1475 std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1476 Tvmult(VectorType &dst, const VectorType &src) const;
1477
1484 template <typename VectorType>
1485 std::enable_if_t<
1486 !std::is_same<typename VectorType::value_type, TrilinosScalar>::value>
1487 Tvmult(VectorType &dst, const VectorType &src) const;
1488
1498 template <typename VectorType>
1499 void
1500 vmult_add(VectorType &dst, const VectorType &src) const;
1501
1512 template <typename VectorType>
1513 void
1514 Tvmult_add(VectorType &dst, const VectorType &src) const;
1515
1538 matrix_norm_square(const MPI::Vector &v) const;
1539
1560 matrix_scalar_product(const MPI::Vector &u, const MPI::Vector &v) const;
1561
1579 residual(MPI::Vector & dst,
1580 const MPI::Vector &x,
1581 const MPI::Vector &b) const;
1582
1597 void
1598 mmult(SparseMatrix & C,
1599 const SparseMatrix &B,
1600 const MPI::Vector & V = MPI::Vector()) const;
1601
1602
1619 void
1620 Tmmult(SparseMatrix & C,
1621 const SparseMatrix &B,
1622 const MPI::Vector & V = MPI::Vector()) const;
1623
1638 l1_norm() const;
1639
1649 linfty_norm() const;
1650
1656 frobenius_norm() const;
1657
1668 const Epetra_CrsMatrix &
1670
1675 const Epetra_CrsGraph &
1677
1689 IndexSet
1691
1697 IndexSet
1699
1726 begin() const;
1727
1731 iterator
1733
1739 end() const;
1740
1744 iterator
1746
1776 begin(const size_type r) const;
1777
1781 iterator
1783
1794 end(const size_type r) const;
1795
1799 iterator
1800 end(const size_type r);
1801
1813 void
1814 write_ascii();
1815
1823 void
1824 print(std::ostream &out,
1825 const bool write_extended_trilinos_info = false) const;
1826
1837 int,
1838 << "An error with error number " << arg1
1839 << " occurred while calling a Trilinos function. "
1840 "\n\n"
1841 "For historical reasons, many Trilinos functions express "
1842 "errors by returning specific integer values to indicate "
1843 "certain errors. Unfortunately, different Trilinos functions "
1844 "often use the same integer values for different kinds of "
1845 "errors, and in most cases it is also not documented what "
1846 "each error code actually means. As a consequence, it is often "
1847 "difficult to say what a particular error (in this case, "
1848 "the error with integer code '"
1849 << arg1
1850 << "') represents and how one should fix a code to avoid it. "
1851 "The best one can often do is to look up the call stack to "
1852 "see which deal.II function generated the error, and which "
1853 "Trilinos function the error code had originated from; "
1854 "then look up the Trilinos source code of that function (for "
1855 "example on github) to see what code path set that error "
1856 "code. Short of going through all of that, the only other "
1857 "option is to guess the cause of the error from "
1858 "the context in which the error appeared.");
1859
1860
1865 size_type,
1866 size_type,
1867 << "The entry with index <" << arg1 << ',' << arg2
1868 << "> does not exist.");
1869
1874 "You are attempting an operation on two vectors that "
1875 "are the same object, but the operation requires that the "
1876 "two objects are in fact different.");
1877
1882
1887 size_type,
1888 size_type,
1889 size_type,
1890 size_type,
1891 << "You tried to access element (" << arg1 << '/' << arg2
1892 << ')'
1893 << " of a distributed matrix, but only rows in range ["
1894 << arg3 << ',' << arg4
1895 << "] are stored locally and can be accessed.");
1896
1901 size_type,
1902 size_type,
1903 << "You tried to access element (" << arg1 << '/' << arg2
1904 << ')' << " of a sparse matrix, but it appears to not"
1905 << " exist in the Trilinos sparsity pattern.");
1910 protected:
1921 void
1923
1929 void
1931
1932
1933
1934 private:
1939 std::unique_ptr<Epetra_Map> column_space_map;
1940
1946 std::unique_ptr<Epetra_FECrsMatrix> matrix;
1947
1953 std::unique_ptr<Epetra_CrsMatrix> nonlocal_matrix;
1954
1958 std::unique_ptr<Epetra_Export> nonlocal_matrix_exporter;
1959
1971 Epetra_CombineMode last_action;
1972
1978
1979 // To allow calling protected prepare_add() and prepare_set().
1980 friend class BlockMatrixBase<SparseMatrix>;
1981 };
1982
1983
1984
1985 // forwards declarations
1986 class SolverBase;
1987 class PreconditionBase;
1988
1989 namespace internal
1990 {
1991 inline void
1992 check_vector_map_equality(const Epetra_CrsMatrix & mtrx,
1993 const Epetra_MultiVector &src,
1994 const Epetra_MultiVector &dst,
1995 const bool transpose)
1996 {
1997 if (transpose == false)
1998 {
1999 Assert(src.Map().SameAs(mtrx.DomainMap()) == true,
2000 ExcMessage(
2001 "Column map of matrix does not fit with vector map!"));
2002 Assert(dst.Map().SameAs(mtrx.RangeMap()) == true,
2003 ExcMessage("Row map of matrix does not fit with vector map!"));
2004 }
2005 else
2006 {
2007 Assert(src.Map().SameAs(mtrx.RangeMap()) == true,
2008 ExcMessage(
2009 "Column map of matrix does not fit with vector map!"));
2010 Assert(dst.Map().SameAs(mtrx.DomainMap()) == true,
2011 ExcMessage("Row map of matrix does not fit with vector map!"));
2012 }
2013 (void)mtrx; // removes -Wunused-variable in optimized mode
2014 (void)src;
2015 (void)dst;
2016 }
2017
2018 inline void
2020 const Epetra_MultiVector &src,
2021 const Epetra_MultiVector &dst,
2022 const bool transpose)
2023 {
2024 if (transpose == false)
2025 {
2026 Assert(src.Map().SameAs(op.OperatorDomainMap()) == true,
2027 ExcMessage(
2028 "Column map of operator does not fit with vector map!"));
2029 Assert(dst.Map().SameAs(op.OperatorRangeMap()) == true,
2030 ExcMessage(
2031 "Row map of operator does not fit with vector map!"));
2032 }
2033 else
2034 {
2035 Assert(src.Map().SameAs(op.OperatorRangeMap()) == true,
2036 ExcMessage(
2037 "Column map of operator does not fit with vector map!"));
2038 Assert(dst.Map().SameAs(op.OperatorDomainMap()) == true,
2039 ExcMessage(
2040 "Row map of operator does not fit with vector map!"));
2041 }
2042 (void)op; // removes -Wunused-variable in optimized mode
2043 (void)src;
2044 (void)dst;
2045 }
2046
2047 namespace LinearOperatorImplementation
2048 {
2069 {
2070 public:
2074 using VectorType = Epetra_MultiVector;
2075
2080
2085
2099
2103 TrilinosPayload(const TrilinosWrappers::SparseMatrix &matrix_exemplar,
2104 const TrilinosWrappers::SparseMatrix &matrix);
2105
2109 TrilinosPayload(const TrilinosPayload & payload_exemplar,
2110 const TrilinosWrappers::SparseMatrix &matrix);
2111
2116 const TrilinosWrappers::SparseMatrix & matrix_exemplar,
2117 const TrilinosWrappers::PreconditionBase &preconditioner);
2118
2123 const TrilinosWrappers::PreconditionBase &preconditioner_exemplar,
2124 const TrilinosWrappers::PreconditionBase &preconditioner);
2125
2130 const TrilinosPayload & payload_exemplar,
2131 const TrilinosWrappers::PreconditionBase &preconditioner);
2132
2136 TrilinosPayload(const TrilinosPayload &payload);
2137
2145 TrilinosPayload(const TrilinosPayload &first_op,
2146 const TrilinosPayload &second_op);
2147
2151 virtual ~TrilinosPayload() override = default;
2152
2157 identity_payload() const;
2158
2163 null_payload() const;
2164
2169 transpose_payload() const;
2170
2187 template <typename Solver, typename Preconditioner>
2188 std::enable_if_t<
2189 std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2190 std::is_base_of<TrilinosWrappers::PreconditionBase,
2191 Preconditioner>::value,
2193 inverse_payload(Solver &, const Preconditioner &) const;
2194
2212 template <typename Solver, typename Preconditioner>
2213 std::enable_if_t<
2214 !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
2215 std::is_base_of<TrilinosWrappers::PreconditionBase,
2216 Preconditioner>::value),
2218 inverse_payload(Solver &, const Preconditioner &) const;
2219
2232 IndexSet
2234
2240 IndexSet
2242
2246 MPI_Comm
2247 get_mpi_communicator() const;
2248
2255 void
2256 transpose();
2257
2265 std::function<void(VectorType &, const VectorType &)> vmult;
2266
2274 std::function<void(VectorType &, const VectorType &)> Tvmult;
2275
2284 std::function<void(VectorType &, const VectorType &)> inv_vmult;
2285
2294 std::function<void(VectorType &, const VectorType &)> inv_Tvmult;
2295
2309 virtual bool
2310 UseTranspose() const override;
2311
2327 virtual int
2328 SetUseTranspose(bool UseTranspose) override;
2329
2341 virtual int
2342 Apply(const VectorType &X, VectorType &Y) const override;
2343
2362 virtual int
2363 ApplyInverse(const VectorType &Y, VectorType &X) const override;
2377 virtual const char *
2378 Label() const override;
2379
2387 virtual const Epetra_Comm &
2388 Comm() const override;
2389
2397 virtual const Epetra_Map &
2398 OperatorDomainMap() const override;
2399
2408 virtual const Epetra_Map &
2409 OperatorRangeMap() const override;
2412 private:
2422 template <typename EpetraOpType>
2423 TrilinosPayload(EpetraOpType & op,
2424 const bool supports_inverse_operations,
2425 const bool use_transpose,
2426 const MPI_Comm mpi_communicator,
2429
2435
2440 Epetra_MpiComm communicator;
2441
2446 Epetra_Map domain_map;
2447
2452 Epetra_Map range_map;
2453
2462 virtual bool
2463 HasNormInf() const override;
2464
2472 virtual double
2473 NormInf() const override;
2474 };
2475
2481 operator+(const TrilinosPayload &first_op,
2482 const TrilinosPayload &second_op);
2483
2489 operator*(const TrilinosPayload &first_op,
2490 const TrilinosPayload &second_op);
2491
2492 } // namespace LinearOperatorImplementation
2493 } /* namespace internal */
2494
2495
2496
2497 // ----------------------- inline and template functions --------------------
2498
2499# ifndef DOXYGEN
2500
2501 namespace SparseMatrixIterators
2502 {
2504 size_type row,
2505 size_type index)
2506 : matrix(matrix)
2507 , a_row(row)
2508 , a_index(index)
2509 {
2510 visit_present_row();
2511 }
2512
2513
2514 inline AccessorBase::size_type
2515 AccessorBase::row() const
2516 {
2517 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2518 return a_row;
2519 }
2520
2521
2522 inline AccessorBase::size_type
2523 AccessorBase::column() const
2524 {
2525 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2526 return (*colnum_cache)[a_index];
2527 }
2528
2529
2530 inline AccessorBase::size_type
2531 AccessorBase::index() const
2532 {
2533 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2534 return a_index;
2535 }
2536
2537
2538 inline Accessor<true>::Accessor(MatrixType * matrix,
2539 const size_type row,
2540 const size_type index)
2541 : AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2542 {}
2543
2544
2545 template <bool Other>
2546 inline Accessor<true>::Accessor(const Accessor<Other> &other)
2547 : AccessorBase(other)
2548 {}
2549
2550
2551 inline TrilinosScalar
2552 Accessor<true>::value() const
2553 {
2554 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2555 return (*value_cache)[a_index];
2556 }
2557
2558
2559 inline Accessor<false>::Reference::Reference(const Accessor<false> &acc)
2560 : accessor(const_cast<Accessor<false> &>(acc))
2561 {}
2562
2563
2564 inline Accessor<false>::Reference::operator TrilinosScalar() const
2565 {
2566 return (*accessor.value_cache)[accessor.a_index];
2567 }
2568
2569 inline const Accessor<false>::Reference &
2570 Accessor<false>::Reference::operator=(const TrilinosScalar n) const
2571 {
2572 (*accessor.value_cache)[accessor.a_index] = n;
2573 accessor.matrix->set(accessor.row(),
2574 accessor.column(),
2575 static_cast<TrilinosScalar>(*this));
2576 return *this;
2577 }
2578
2579
2580 inline const Accessor<false>::Reference &
2581 Accessor<false>::Reference::operator+=(const TrilinosScalar n) const
2582 {
2583 (*accessor.value_cache)[accessor.a_index] += n;
2584 accessor.matrix->set(accessor.row(),
2585 accessor.column(),
2586 static_cast<TrilinosScalar>(*this));
2587 return *this;
2588 }
2589
2590
2591 inline const Accessor<false>::Reference &
2592 Accessor<false>::Reference::operator-=(const TrilinosScalar n) const
2593 {
2594 (*accessor.value_cache)[accessor.a_index] -= n;
2595 accessor.matrix->set(accessor.row(),
2596 accessor.column(),
2597 static_cast<TrilinosScalar>(*this));
2598 return *this;
2599 }
2600
2601
2602 inline const Accessor<false>::Reference &
2603 Accessor<false>::Reference::operator*=(const TrilinosScalar n) const
2604 {
2605 (*accessor.value_cache)[accessor.a_index] *= n;
2606 accessor.matrix->set(accessor.row(),
2607 accessor.column(),
2608 static_cast<TrilinosScalar>(*this));
2609 return *this;
2610 }
2611
2612
2613 inline const Accessor<false>::Reference &
2614 Accessor<false>::Reference::operator/=(const TrilinosScalar n) const
2615 {
2616 (*accessor.value_cache)[accessor.a_index] /= n;
2617 accessor.matrix->set(accessor.row(),
2618 accessor.column(),
2619 static_cast<TrilinosScalar>(*this));
2620 return *this;
2621 }
2622
2623
2624 inline Accessor<false>::Accessor(MatrixType * matrix,
2625 const size_type row,
2626 const size_type index)
2627 : AccessorBase(matrix, row, index)
2628 {}
2629
2630
2631 inline Accessor<false>::Reference
2632 Accessor<false>::value() const
2633 {
2634 Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
2635 return {*this};
2636 }
2637
2638
2639
2640 template <bool Constness>
2641 inline Iterator<Constness>::Iterator(MatrixType * matrix,
2642 const size_type row,
2643 const size_type index)
2644 : accessor(matrix, row, index)
2645 {}
2646
2647
2648 template <bool Constness>
2649 template <bool Other>
2650 inline Iterator<Constness>::Iterator(const Iterator<Other> &other)
2651 : accessor(other.accessor)
2652 {}
2653
2654
2655 template <bool Constness>
2656 inline Iterator<Constness> &
2657 Iterator<Constness>::operator++()
2658 {
2659 Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
2660
2661 ++accessor.a_index;
2662
2663 // If at end of line: do one
2664 // step, then cycle until we
2665 // find a row with a nonzero
2666 // number of entries.
2667 if (accessor.a_index >= accessor.colnum_cache->size())
2668 {
2669 accessor.a_index = 0;
2670 ++accessor.a_row;
2671
2672 while ((accessor.a_row < accessor.matrix->m()) &&
2673 ((accessor.matrix->in_local_range(accessor.a_row) == false) ||
2674 (accessor.matrix->row_length(accessor.a_row) == 0)))
2675 ++accessor.a_row;
2676
2677 accessor.visit_present_row();
2678 }
2679 return *this;
2680 }
2681
2682
2683 template <bool Constness>
2684 inline Iterator<Constness>
2685 Iterator<Constness>::operator++(int)
2686 {
2687 const Iterator<Constness> old_state = *this;
2688 ++(*this);
2689 return old_state;
2690 }
2691
2692
2693
2694 template <bool Constness>
2695 inline const Accessor<Constness> &
2696 Iterator<Constness>::operator*() const
2697 {
2698 return accessor;
2699 }
2700
2701
2702
2703 template <bool Constness>
2704 inline const Accessor<Constness> *
2705 Iterator<Constness>::operator->() const
2706 {
2707 return &accessor;
2708 }
2709
2710
2711
2712 template <bool Constness>
2713 template <bool OtherConstness>
2714 inline bool
2715 Iterator<Constness>::operator==(const Iterator<OtherConstness> &other) const
2716 {
2717 return (accessor.a_row == other.accessor.a_row &&
2718 accessor.a_index == other.accessor.a_index);
2719 }
2720
2721
2722
2723 template <bool Constness>
2724 template <bool OtherConstness>
2725 inline bool
2726 Iterator<Constness>::operator!=(const Iterator<OtherConstness> &other) const
2727 {
2728 return !(*this == other);
2729 }
2730
2731
2732
2733 template <bool Constness>
2734 template <bool OtherConstness>
2735 inline bool
2736 Iterator<Constness>::operator<(const Iterator<OtherConstness> &other) const
2737 {
2738 return (accessor.row() < other.accessor.row() ||
2739 (accessor.row() == other.accessor.row() &&
2740 accessor.index() < other.accessor.index()));
2741 }
2742
2743
2744 template <bool Constness>
2745 template <bool OtherConstness>
2746 inline bool
2747 Iterator<Constness>::operator>(const Iterator<OtherConstness> &other) const
2748 {
2749 return (other < *this);
2750 }
2751
2752 } // namespace SparseMatrixIterators
2753
2754
2755
2757 SparseMatrix::begin() const
2758 {
2759 return begin(0);
2760 }
2761
2762
2763
2765 SparseMatrix::end() const
2766 {
2767 return const_iterator(this, m(), 0);
2768 }
2769
2770
2771
2773 SparseMatrix::begin(const size_type r) const
2774 {
2775 AssertIndexRange(r, m());
2776 if (in_local_range(r) && (row_length(r) > 0))
2777 return const_iterator(this, r, 0);
2778 else
2779 return end(r);
2780 }
2781
2782
2783
2785 SparseMatrix::end(const size_type r) const
2786 {
2787 AssertIndexRange(r, m());
2788
2789 // place the iterator on the first entry
2790 // past this line, or at the end of the
2791 // matrix
2792 for (size_type i = r + 1; i < m(); ++i)
2793 if (in_local_range(i) && (row_length(i) > 0))
2794 return const_iterator(this, i, 0);
2795
2796 // if there is no such line, then take the
2797 // end iterator of the matrix
2798 return end();
2799 }
2800
2801
2802
2805 {
2806 return begin(0);
2807 }
2808
2809
2810
2813 {
2814 return iterator(this, m(), 0);
2815 }
2816
2817
2818
2820 SparseMatrix::begin(const size_type r)
2821 {
2822 AssertIndexRange(r, m());
2823 if (in_local_range(r) && (row_length(r) > 0))
2824 return iterator(this, r, 0);
2825 else
2826 return end(r);
2827 }
2828
2829
2830
2832 SparseMatrix::end(const size_type r)
2833 {
2834 AssertIndexRange(r, m());
2835
2836 // place the iterator on the first entry
2837 // past this line, or at the end of the
2838 // matrix
2839 for (size_type i = r + 1; i < m(); ++i)
2840 if (in_local_range(i) && (row_length(i) > 0))
2841 return iterator(this, i, 0);
2842
2843 // if there is no such line, then take the
2844 // end iterator of the matrix
2845 return end();
2846 }
2847
2848
2849
2850 inline bool
2851 SparseMatrix::in_local_range(const size_type index) const
2852 {
2854# ifndef DEAL_II_WITH_64BIT_INDICES
2855 begin = matrix->RowMap().MinMyGID();
2856 end = matrix->RowMap().MaxMyGID() + 1;
2857# else
2858 begin = matrix->RowMap().MinMyGID64();
2859 end = matrix->RowMap().MaxMyGID64() + 1;
2860# endif
2861
2862 return ((index >= static_cast<size_type>(begin)) &&
2863 (index < static_cast<size_type>(end)));
2864 }
2865
2866
2867
2868 inline bool
2870 {
2871 return compressed;
2872 }
2873
2874
2875
2876 // Inline the set() and add() functions, since they will be called
2877 // frequently, and the compiler can optimize away some unnecessary loops
2878 // when the sizes are given at compile time.
2879 template <>
2880 void
2881 SparseMatrix::set<TrilinosScalar>(const size_type row,
2882 const size_type n_cols,
2883 const size_type * col_indices,
2884 const TrilinosScalar *values,
2885 const bool elide_zero_values);
2886
2887
2888
2889 template <typename Number>
2890 void
2891 SparseMatrix::set(const size_type row,
2892 const size_type n_cols,
2893 const size_type *col_indices,
2894 const Number * values,
2895 const bool elide_zero_values)
2896 {
2897 std::vector<TrilinosScalar> trilinos_values(n_cols);
2898 std::copy(values, values + n_cols, trilinos_values.begin());
2899 this->set(
2900 row, n_cols, col_indices, trilinos_values.data(), elide_zero_values);
2901 }
2902
2903
2904
2905 inline void
2906 SparseMatrix::set(const size_type i,
2907 const size_type j,
2908 const TrilinosScalar value)
2909 {
2910 AssertIsFinite(value);
2911
2912 set(i, 1, &j, &value, false);
2913 }
2914
2915
2916
2917 inline void
2918 SparseMatrix::set(const std::vector<size_type> & indices,
2919 const FullMatrix<TrilinosScalar> &values,
2920 const bool elide_zero_values)
2921 {
2922 Assert(indices.size() == values.m(),
2923 ExcDimensionMismatch(indices.size(), values.m()));
2924 Assert(values.m() == values.n(), ExcNotQuadratic());
2925
2926 for (size_type i = 0; i < indices.size(); ++i)
2927 set(indices[i],
2928 indices.size(),
2929 indices.data(),
2930 &values(i, 0),
2931 elide_zero_values);
2932 }
2933
2934
2935
2936 inline void
2937 SparseMatrix::add(const size_type i,
2938 const size_type j,
2939 const TrilinosScalar value)
2940 {
2941 AssertIsFinite(value);
2942
2943 if (value == 0)
2944 {
2945 // we have to check after Insert/Add in any case to be consistent
2946 // with the MPI communication model, but we can save some
2947 // work if the addend is zero. However, these actions are done in case
2948 // we pass on to the other function.
2949
2950 // TODO: fix this (do not run compress here, but fail)
2951 if (last_action == Insert)
2952 {
2953 const int ierr = matrix->GlobalAssemble(*column_space_map,
2954 matrix->RowMap(),
2955 false);
2956
2957 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
2958 }
2959
2960 last_action = Add;
2961
2962 return;
2963 }
2964 else
2965 add(i, 1, &j, &value, false);
2966 }
2967
2968
2969
2970 // inline "simple" functions that are called frequently and do only involve
2971 // a call to some Trilinos function.
2973 SparseMatrix::m() const
2974 {
2975# ifndef DEAL_II_WITH_64BIT_INDICES
2976 return matrix->NumGlobalRows();
2977# else
2978 return matrix->NumGlobalRows64();
2979# endif
2980 }
2981
2982
2983
2985 SparseMatrix::n() const
2986 {
2987 // If the matrix structure has not been fixed (i.e., we did not have a
2988 // sparsity pattern), it does not know about the number of columns so we
2989 // must always take this from the additional column space map
2990 Assert(column_space_map.get() != nullptr, ExcInternalError());
2991 return n_global_elements(*column_space_map);
2992 }
2993
2994
2995
2996 inline unsigned int
2998 {
2999 return matrix->NumMyRows();
3000 }
3001
3002
3003
3004 inline std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
3006 {
3008# ifndef DEAL_II_WITH_64BIT_INDICES
3009 begin = matrix->RowMap().MinMyGID();
3010 end = matrix->RowMap().MaxMyGID() + 1;
3011# else
3012 begin = matrix->RowMap().MinMyGID64();
3013 end = matrix->RowMap().MaxMyGID64() + 1;
3014# endif
3015
3016 return std::make_pair(begin, end);
3017 }
3018
3019
3020
3021 inline std::uint64_t
3023 {
3024 // Trilinos uses 64bit functions internally for attribute access, which
3025 // return `long long`. They also offer 32bit variants that return `int`,
3026 // however those call the 64bit version and convert the values to 32bit.
3027 // There is no necessity in using the 32bit versions at all.
3028 return static_cast<std::uint64_t>(matrix->NumGlobalNonzeros64());
3029 }
3030
3031
3032
3033 template <typename SparsityPatternType>
3034 inline void
3035 SparseMatrix::reinit(const IndexSet & parallel_partitioning,
3036 const SparsityPatternType &sparsity_pattern,
3037 const MPI_Comm communicator,
3038 const bool exchange_data)
3039 {
3040 reinit(parallel_partitioning,
3041 parallel_partitioning,
3042 sparsity_pattern,
3043 communicator,
3044 exchange_data);
3045 }
3046
3047
3048
3049 template <typename number>
3050 inline void
3051 SparseMatrix::reinit(const IndexSet &parallel_partitioning,
3052 const ::SparseMatrix<number> &sparse_matrix,
3053 const MPI_Comm communicator,
3054 const double drop_tolerance,
3055 const bool copy_values,
3056 const ::SparsityPattern * use_this_sparsity)
3057 {
3058 Epetra_Map map =
3059 parallel_partitioning.make_trilinos_map(communicator, false);
3060 reinit(parallel_partitioning,
3061 parallel_partitioning,
3062 sparse_matrix,
3063 drop_tolerance,
3064 copy_values,
3065 use_this_sparsity);
3066 }
3067
3068
3069
3070 inline const Epetra_CrsMatrix &
3072 {
3073 return static_cast<const Epetra_CrsMatrix &>(*matrix);
3074 }
3075
3076
3077
3078 inline const Epetra_CrsGraph &
3080 {
3081 return matrix->Graph();
3082 }
3083
3084
3085
3086 inline IndexSet
3088 {
3089 return IndexSet(matrix->DomainMap());
3090 }
3091
3092
3093
3094 inline IndexSet
3096 {
3097 return IndexSet(matrix->RangeMap());
3098 }
3099
3100
3101
3102 inline void
3104 {
3105 // nothing to do here
3106 }
3107
3108
3109
3110 inline void
3112 {
3113 // nothing to do here
3114 }
3115
3116
3117 namespace internal
3118 {
3119 namespace LinearOperatorImplementation
3120 {
3121 template <typename EpetraOpType>
3122 TrilinosPayload::TrilinosPayload(
3123 EpetraOpType & op,
3124 const bool supports_inverse_operations,
3125 const bool use_transpose,
3126 const MPI_Comm mpi_communicator,
3127 const IndexSet &locally_owned_domain_indices,
3128 const IndexSet &locally_owned_range_indices)
3129 : use_transpose(use_transpose)
3130 , communicator(mpi_communicator)
3131 , domain_map(
3132 locally_owned_domain_indices.make_trilinos_map(communicator.Comm()))
3133 , range_map(
3134 locally_owned_range_indices.make_trilinos_map(communicator.Comm()))
3135 {
3136 vmult = [&op](Range &tril_dst, const Domain &tril_src) {
3137 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3138 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3139 Assert(&tril_src != &tril_dst,
3142 tril_src,
3143 tril_dst,
3144 op.UseTranspose());
3145
3146 const int ierr = op.Apply(tril_src, tril_dst);
3147 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3148 };
3149
3150 Tvmult = [&op](Domain &tril_dst, const Range &tril_src) {
3151 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3152 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3153 Assert(&tril_src != &tril_dst,
3156 tril_src,
3157 tril_dst,
3158 !op.UseTranspose());
3159
3160 op.SetUseTranspose(!op.UseTranspose());
3161 const int ierr = op.Apply(tril_src, tril_dst);
3162 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3163 op.SetUseTranspose(!op.UseTranspose());
3164 };
3165
3166 if (supports_inverse_operations)
3167 {
3168 inv_vmult = [&op](Domain &tril_dst, const Range &tril_src) {
3169 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3170 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3171 Assert(
3172 &tril_src != &tril_dst,
3175 tril_src,
3176 tril_dst,
3177 !op.UseTranspose());
3178
3179 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3180 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3181 };
3182
3183 inv_Tvmult = [&op](Range &tril_dst, const Domain &tril_src) {
3184 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3185 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3186 Assert(
3187 &tril_src != &tril_dst,
3190 tril_src,
3191 tril_dst,
3192 op.UseTranspose());
3193
3194 op.SetUseTranspose(!op.UseTranspose());
3195 const int ierr = op.ApplyInverse(tril_src, tril_dst);
3196 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
3197 op.SetUseTranspose(!op.UseTranspose());
3198 };
3199 }
3200 else
3201 {
3202 inv_vmult = [](Domain &, const Range &) {
3203 Assert(false,
3204 ExcMessage(
3205 "Uninitialized TrilinosPayload::inv_vmult called. "
3206 "The operator does not support inverse operations."));
3207 };
3208
3209 inv_Tvmult = [](Range &, const Domain &) {
3210 Assert(false,
3211 ExcMessage(
3212 "Uninitialized TrilinosPayload::inv_Tvmult called. "
3213 "The operator does not support inverse operations."));
3214 };
3215 }
3216 }
3217
3218
3219 template <typename Solver, typename Preconditioner>
3220 std::enable_if_t<
3221 std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3222 std::is_base_of<TrilinosWrappers::PreconditionBase,
3223 Preconditioner>::value,
3224 TrilinosPayload>
3226 Solver & solver,
3227 const Preconditioner &preconditioner) const
3228 {
3229 const auto &payload = *this;
3230
3231 TrilinosPayload return_op(payload);
3232
3233 // Capture by copy so the payloads are always valid
3234
3235 return_op.inv_vmult = [payload, &solver, &preconditioner](
3236 TrilinosPayload::Domain & tril_dst,
3237 const TrilinosPayload::Range &tril_src) {
3238 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3239 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3240 Assert(&tril_src != &tril_dst,
3243 tril_src,
3244 tril_dst,
3245 !payload.UseTranspose());
3246 solver.solve(payload, tril_dst, tril_src, preconditioner);
3247 };
3248
3249 return_op.inv_Tvmult = [payload, &solver, &preconditioner](
3250 TrilinosPayload::Range & tril_dst,
3251 const TrilinosPayload::Domain &tril_src) {
3252 // Duplicated from TrilinosWrappers::PreconditionBase::vmult
3253 // as well as from TrilinosWrappers::SparseMatrix::Tvmult
3254 Assert(&tril_src != &tril_dst,
3257 tril_src,
3258 tril_dst,
3259 payload.UseTranspose());
3260
3261 const_cast<TrilinosPayload &>(payload).transpose();
3262 solver.solve(payload, tril_dst, tril_src, preconditioner);
3263 const_cast<TrilinosPayload &>(payload).transpose();
3264 };
3265
3266 // If the input operator is already setup for transpose operations, then
3267 // we must do similar with its inverse.
3268 if (return_op.UseTranspose() == true)
3269 std::swap(return_op.inv_vmult, return_op.inv_Tvmult);
3270
3271 return return_op;
3272 }
3273
3274 template <typename Solver, typename Preconditioner>
3275 std::enable_if_t<
3276 !(std::is_base_of<TrilinosWrappers::SolverBase, Solver>::value &&
3277 std::is_base_of<TrilinosWrappers::PreconditionBase,
3278 Preconditioner>::value),
3279 TrilinosPayload>
3280 TrilinosPayload::inverse_payload(Solver &, const Preconditioner &) const
3281 {
3282 TrilinosPayload return_op(*this);
3283
3284 return_op.inv_vmult = [](TrilinosPayload::Domain &,
3285 const TrilinosPayload::Range &) {
3286 AssertThrow(false,
3287 ExcMessage("Payload inv_vmult disabled because of "
3288 "incompatible solver/preconditioner choice."));
3289 };
3290
3291 return_op.inv_Tvmult = [](TrilinosPayload::Range &,
3292 const TrilinosPayload::Domain &) {
3293 AssertThrow(false,
3294 ExcMessage("Payload inv_vmult disabled because of "
3295 "incompatible solver/preconditioner choice."));
3296 };
3297
3298 return return_op;
3299 }
3300 } // namespace LinearOperatorImplementation
3301 } // namespace internal
3302
3303 template <>
3304 void
3305 SparseMatrix::set<TrilinosScalar>(const size_type row,
3306 const size_type n_cols,
3307 const size_type * col_indices,
3308 const TrilinosScalar *values,
3309 const bool elide_zero_values);
3310# endif // DOXYGEN
3311
3312} /* namespace TrilinosWrappers */
3313
3314
3316
3317
3318# endif // DEAL_II_WITH_TRILINOS
3319
3320
3321/*----------------------- trilinos_sparse_matrix.h --------------------*/
3322
3323#endif
3324/*----------------------- trilinos_sparse_matrix.h --------------------*/
Epetra_Map make_trilinos_map(const MPI_Comm communicator=MPI_COMM_WORLD, const bool overlapping=false) const
Definition index_set.cc:789
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
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
::types::global_dof_index size_type
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
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
const_iterator begin() const
void reinit(const SparsityPatternType &sparsity_pattern)
std::enable_if_t< std::is_same< typename VectorType::value_type, TrilinosScalar >::value > Tvmult(VectorType &dst, const VectorType &src) const
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
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)
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
std::enable_if_t< !std::is_same< SparsityPatternType, ::SparseMatrix< double > >::value > 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)
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
std::enable_if_t< std::is_same< typename VectorType::value_type, TrilinosScalar >::value > vmult(VectorType &dst, const VectorType &src) 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::enable_if_t< !(std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value), TrilinosPayload > inverse_payload(Solver &, const Preconditioner &) const
std::function< void(VectorType &, const VectorType &)> inv_Tvmult
virtual int ApplyInverse(const VectorType &Y, VectorType &X) const override
virtual int Apply(const VectorType &X, VectorType &Y) const override
std::enable_if_t< std::is_base_of< TrilinosWrappers::SolverBase, Solver >::value &&std::is_base_of< TrilinosWrappers::PreconditionBase, Preconditioner >::value, TrilinosPayload > 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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
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:465
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:579
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:533
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcSourceEqualsDestination()
#define DeclException3(Exception3, type1, type2, type3, outsequence)
Definition exceptions.h:556
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
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
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:82
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::value_type value_type
typename ::TrilinosWrappers::SparseMatrixIterators::Iterator< Constness >::difference_type difference_type
double TrilinosScalar
Definition types.h:175