15#ifndef dealii_trilinos_tpetra_sparse_matrix_h
16#define dealii_trilinos_tpetra_sparse_matrix_h
20#ifdef DEAL_II_TRILINOS_WITH_TPETRA
33# include <Tpetra_Core.hpp>
34# include <Tpetra_CrsMatrix.hpp>
36# include <type_traits>
43template <
typename Number>
47template <
typename MatrixType>
52 namespace TpetraWrappers
54 template <
typename MemorySpace>
59 template <
typename Number,
typename MemorySpace,
bool Constness>
69 namespace TpetraWrappers
109 template <
typename Number,
typename MemorySpace = ::MemorySpace::Host>
123 <<
"You tried to access row " << arg1
124 <<
" of a non-contiguous locally owned row set."
125 <<
" The row " << arg1
126 <<
" is not stored locally and can't be accessed.");
165# if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0)
166 using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode<
167 typename MemorySpace::kokkos_space::execution_space,
168 typename MemorySpace::kokkos_space>;
170 using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode<
171 typename MemorySpace::kokkos_space::execution_space,
172 typename MemorySpace::kokkos_space>;
179 Tpetra::CrsMatrix<Number,
188 Tpetra::Map<int, ::types::signed_global_dof_index, NodeType>;
194 Tpetra::CrsGraph<int, ::types::signed_global_dof_index, NodeType>;
200 Vector<Number, int, ::types::signed_global_dof_index, NodeType>;
225 const unsigned int n_max_entries_per_row);
236 const std::vector<unsigned int> &n_entries_per_row);
282 template <
typename SparsityPatternType>
284 reinit(
const SparsityPatternType &sparsity_pattern);
316 const MPI_Comm communicator = MPI_COMM_WORLD,
317 const unsigned int n_max_entries_per_row = 0);
329 const std::vector<unsigned int> &n_entries_per_row);
346 const IndexSet &col_parallel_partitioning,
347 const MPI_Comm communicator = MPI_COMM_WORLD,
348 const size_type n_max_entries_per_row = 0);
359 const IndexSet &col_parallel_partitioning,
361 const std::vector<unsigned int> &n_entries_per_row);
382 template <
typename SparsityPatternType>
384 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
386 const SparsityPatternType &sparsity_pattern,
387 const MPI_Comm communicator = MPI_COMM_WORLD,
388 const bool exchange_data =
false);
402 template <
typename SparsityPatternType>
404 !std::is_same_v<SparsityPatternType, ::SparseMatrix<double>>>
406 const IndexSet &col_parallel_partitioning,
407 const SparsityPatternType &sparsity_pattern,
408 const MPI_Comm communicator = MPI_COMM_WORLD,
409 const bool exchange_data =
false);
429 const IndexSet &col_parallel_partitioning,
430 const ::SparseMatrix<Number> &dealii_sparse_matrix,
431 const MPI_Comm communicator = MPI_COMM_WORLD,
432 const double drop_tolerance = 1e-13,
433 const bool copy_values =
true,
434 const ::SparsityPattern *use_this_sparsity =
nullptr);
474 std::pair<size_type, size_type>
576 const Number *values,
577 const bool elide_zero_values =
true,
578 const bool col_indices_are_sorted =
false);
647 set(
const std::vector<size_type> &indices,
649 const bool elide_zero_values =
false);
657 set(
const std::vector<size_type> &row_indices,
658 const std::vector<size_type> &col_indices,
660 const bool elide_zero_values =
false);
691 const std::vector<size_type> &col_indices,
692 const std::vector<Number> &values,
693 const bool elide_zero_values =
false);
726 const Number *values,
727 const bool elide_zero_values =
false);
780 const Number new_diag_value = 0);
852 template <
typename InputVectorType>
854 vmult(InputVectorType &dst,
const InputVectorType &src)
const;
863 template <
typename InputVectorType>
865 Tvmult(InputVectorType &dst,
const InputVectorType &src)
const;
873 template <
typename InputVectorType>
875 vmult_add(InputVectorType &dst,
const InputVectorType &src)
const;
884 template <
typename InputVectorType>
886 Tvmult_add(InputVectorType &dst,
const InputVectorType &src)
const;
985 const bool print_detailed_trilinos_information =
false)
const;
1057 Teuchos::RCP<const MatrixType>
1069 Teuchos::RCP<MatrixType>
1211 "You are attempting an operation on two vectors that "
1212 "are the same object, but the operation requires that the "
1213 "two objects are in fact different.");
1219 "The column partitioning of a matrix does not match "
1220 "the partitioning of a vector you are trying to "
1221 "multiply it with. Are you multiplying the "
1222 "matrix with a vector that has ghost elements?");
1228 "The row partitioning of a matrix does not match "
1229 "the partitioning of a vector you are trying to "
1230 "put the result of a matrix-vector product in. "
1231 "Are you trying to put the product of the "
1232 "matrix with a vector into a vector that has "
1241 <<
"The entry with index <" << arg1 <<
',' << arg2
1242 <<
"> does not exist.");
1252 <<
"You tried to access element (" << arg1 <<
'/' << arg2
1254 <<
" of a distributed matrix, but only rows in range ["
1255 << arg3 <<
',' << arg4
1256 <<
"] are stored locally and can be accessed.");
1337 template <
typename Number,
typename MemorySpace>
1409 std::shared_ptr<std::vector<::types::signed_global_dof_index>>
1432 template <
typename Number,
typename MemorySpace,
bool Constness>
1451 template <
typename Number,
typename MemorySpace>
1479 template <
bool Other>
1492 template <
typename Number,
typename MemorySpace>
1507 operator Number()
const;
1513 operator=(
const Number n)
const;
1519 operator+=(
const Number n)
const;
1525 operator-=(
const Number n)
const;
1531 operator*=(
const Number n)
const;
1537 operator/=(
const Number n)
const;
1575 friend class Reference;
1591 template <
typename Number,
typename MemorySpace,
bool Constness>
1630 template <
bool Other>
1661 template <
bool OtherConstness>
1668 template <
bool OtherConstness>
1677 template <
bool OtherConstness>
1684 template <
bool OtherConstness>
1694 <<
"Attempt to access element " << arg2 <<
" of row "
1695 << arg1 <<
" which doesn't have that many elements.");
1717 template <
typename Number,
typename MemorySpace,
bool Constness>
1720 Iterator<Number, MemorySpace, Constness>>
1724 typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
1725 Iterator<Number, MemorySpace, Constness>::value_type;
1727 typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::
1728 Iterator<Number, MemorySpace, Constness>::difference_type;
1740 namespace TpetraWrappers
1742 template <
typename Number,
typename MemorySpace>
1748 set(i, 1, &j, &value,
false);
1753 template <
typename Number,
typename MemorySpace>
1759 add(i, 1, &j, &value,
false);
1764 template <
typename Number,
typename MemorySpace>
1780 template <
typename Number,
typename MemorySpace>
1784 return matrix->getRowMap()->getGlobalNumElements();
1789 template <
typename Number,
typename MemorySpace>
1797 return column_space_map->getGlobalNumElements();
1802 template <
typename Number,
typename MemorySpace>
1811 template <
typename Number,
typename MemorySpace>
1820 template <
typename Number,
typename MemorySpace>
1829 template <
typename Number,
typename MemorySpace>
1834 if (in_local_range(r) && (row_length(r) > 0))
1841 template <
typename Number,
typename MemorySpace>
1851 if (in_local_range(i) && (row_length(i) > 0))
1861 template <
typename Number,
typename MemorySpace>
1870 template <
typename Number,
typename MemorySpace>
1879 template <
typename Number,
typename MemorySpace>
1884 if (in_local_range(r) && (row_length(r) > 0))
1892 template <
typename Number,
typename MemorySpace>
1902 if (in_local_range(i) && (row_length(i) > 0))
1912 template <
typename Number,
typename MemorySpace>
1917 const size_type begin = matrix->getRowMap()->getMinGlobalIndex();
1918 const size_type end = matrix->getRowMap()->getMaxGlobalIndex() + 1;
1920 return ((index >= begin) && (index < end));
1925 template <
typename Number,
typename MemorySpace>
1929 auto n_entries = matrix->getNumEntriesInGlobalRow(row);
1931 Teuchos::OrdinalTraits<
decltype(n_entries)>::invalid(),
1932 ExcAccessToNonlocalRow(row));
1939 template <
typename Number,
typename MemorySpace>
1948 template <
typename Number,
typename MemorySpace>
1957 template <
typename Number,
typename MemorySpace>
1958 inline const Tpetra::CrsMatrix<
1970 template <
typename Number,
typename MemorySpace>
1971 inline Tpetra::CrsMatrix<
1983 template <
typename Number,
typename MemorySpace>
1984 inline Teuchos::RCP<
const Tpetra::CrsMatrix<
1991 return matrix.getConst();
1996 template <
typename Number,
typename MemorySpace>
1997 inline Teuchos::RCP<
1998 Tpetra::CrsMatrix<Number,
2009 template <
typename Number,
typename MemorySpace>
2013 return IndexSet(matrix->getDomainMap());
2018 template <
typename Number,
typename MemorySpace>
2022 return IndexSet(matrix->getRangeMap());
2028 template <
typename Number,
typename MemorySpace>
2041 template <
typename Number,
typename MemorySpace>
2050 template <
typename Number,
typename MemorySpace>
2055 return (*colnum_cache)[a_index];
2059 template <
typename Number,
typename MemorySpace>
2069 template <
typename Number,
typename MemorySpace>
2079 if ((this->a_row == matrix->m()) ||
2080 (matrix->in_local_range(this->a_row) ==
false))
2082 colnum_cache.reset();
2083 value_cache.reset();
2091 matrix->row_length(this->a_row);
2092 if (value_cache.get() ==
nullptr)
2094 value_cache = std::make_shared<std::vector<Number>>(colnums);
2095 colnum_cache = std::make_shared<
2096 std::vector<::types::signed_global_dof_index>>(colnums);
2100 value_cache->resize(colnums);
2101 colnum_cache->resize(colnums);
2104# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
2106 nonconst_global_inds_host_view_type col_indices(colnum_cache->data(),
2109 nonconst_values_host_view_type values(value_cache->data(), colnums);
2111 Teuchos::ArrayView<::types::signed_global_dof_index> col_indices(
2113 Teuchos::ArrayView<Number> values(*value_cache);
2115 matrix->trilinos_matrix().getGlobalRowCopy(this->a_row,
2132 template <
typename Number,
typename MemorySpace>
2145 template <
typename Number,
typename MemorySpace>
2146 template <
bool Other>
2154 template <
typename Number,
typename MemorySpace>
2167 template <
typename Number,
typename MemorySpace>
2175 template <
typename Number,
typename MemorySpace>
2179 return (*accessor.value_cache)[accessor.a_index];
2184 template <
typename Number,
typename MemorySpace>
2187 const Number n)
const
2189 (*accessor.value_cache)[accessor.a_index] = n;
2190 accessor.matrix->set(accessor.row(),
2192 static_cast<Number
>(*
this));
2198 template <
typename Number,
typename MemorySpace>
2201 const Number n)
const
2203 (*accessor.value_cache)[accessor.a_index] += n;
2204 accessor.matrix->set(accessor.row(),
2206 static_cast<Number
>(*
this));
2212 template <
typename Number,
typename MemorySpace>
2215 const Number n)
const
2217 (*accessor.value_cache)[accessor.a_index] -= n;
2218 accessor.matrix->set(accessor.row(),
2220 static_cast<Number
>(*
this));
2225 template <
typename Number,
typename MemorySpace>
2228 const Number n)
const
2230 (*accessor.value_cache)[accessor.a_index] *= n;
2231 accessor.matrix->set(accessor.row(),
2233 static_cast<Number
>(*
this));
2238 template <
typename Number,
typename MemorySpace>
2241 const Number n)
const
2243 (*accessor.value_cache)[accessor.a_index] /= n;
2244 accessor.matrix->set(accessor.row(),
2246 static_cast<Number
>(*
this));
2251 template <
typename Number,
typename MemorySpace>
2260 template <
typename Number,
typename MemorySpace>
2272 template <
typename Number,
typename MemorySpace,
bool Constness>
2277 : accessor(matrix, row, index)
2281 template <
typename Number,
typename MemorySpace,
bool Constness>
2282 template <
bool Other>
2285 : accessor(other.accessor)
2290 template <
typename Number,
typename MemorySpace,
bool Constness>
2302 if (accessor.a_index >= accessor.colnum_cache->size())
2304 accessor.a_index = 0;
2308 (accessor.a_row < accessor.matrix->m()) &&
2309 ((accessor.matrix->in_local_range(accessor.a_row) ==
false) ||
2310 (accessor.matrix->row_length(accessor.a_row) == 0)))
2313 accessor.visit_present_row();
2319 template <
typename Number,
typename MemorySpace,
bool Constness>
2330 template <
typename Number,
typename MemorySpace,
bool Constness>
2339 template <
typename Number,
typename MemorySpace,
bool Constness>
2348 template <
typename Number,
typename MemorySpace,
bool Constness>
2349 template <
bool OtherConstness>
2354 return (accessor.a_row == other.
accessor.a_row &&
2355 accessor.a_index == other.
accessor.a_index);
2360 template <
typename Number,
typename MemorySpace,
bool Constness>
2361 template <
bool OtherConstness>
2366 return !(*
this == other);
2371 template <
typename Number,
typename MemorySpace,
bool Constness>
2372 template <
bool OtherConstness>
2377 return (accessor.row() < other.
accessor.row() ||
2378 (accessor.row() == other.
accessor.row() &&
2379 accessor.index() < other.
accessor.index()));
2383 template <
typename Number,
typename MemorySpace,
bool Constness>
2384 template <
bool OtherConstness>
2389 return (other < *
this);
SparseMatrix< Number, MemorySpace > * matrix
std::shared_ptr< std::vector< Number > > value_cache
std::shared_ptr< std::vector<::types::signed_global_dof_index > > colnum_cache
AccessorBase(SparseMatrix< Number, MemorySpace > *matrix, const size_type row, const size_type index)
typename AccessorBase< Number, MemorySpace >::size_type size_type
typename AccessorBase< Number, MemorySpace >::size_type size_type
bool operator!=(const Iterator< Number, MemorySpace, OtherConstness > &) const
bool operator<(const Iterator< Number, MemorySpace, OtherConstness > &) const
bool operator>(const Iterator< Number, MemorySpace, OtherConstness > &) const
Iterator(MatrixType *matrix, const size_type row, const size_type index)
Accessor< Number, MemorySpace, Constness > accessor
const Accessor< Number, MemorySpace, Constness > * operator->() const
Iterator< Number, MemorySpace, Constness > & operator++()
typename Accessor< Number, MemorySpace, Constness >::MatrixType MatrixType
const Accessor< Number, MemorySpace, Constness > & operator*() const
bool operator==(const Iterator< Number, MemorySpace, OtherConstness > &) const
virtual ~SparseMatrix() override=default
Tpetra::CrsGraph< int, ::types::signed_global_dof_index, NodeType > GraphType
unsigned int row_length(const size_type row) const
SparseMatrix(const SparseMatrix< Number, MemorySpace > &)=delete
void clear_row(const size_type row, const Number new_diag_value=0)
SparseMatrix(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const MPI_Comm communicator=MPI_COMM_WORLD, const size_type n_max_entries_per_row=0)
void print(std::ostream &out, const bool print_detailed_trilinos_information=false) const
IndexSet locally_owned_domain_indices() const
Tpetra::Vector< Number, int, ::types::signed_global_dof_index, NodeType > VectorType
const MatrixType & trilinos_matrix() const
iterator begin(const size_type r)
SparseMatrix< Number, MemorySpace > & operator=(SparseMatrix< Number, MemorySpace > &&other) noexcept
void add(const size_type i, const size_type j, const Number value)
void set(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< Number > &full_matrix, const bool elide_zero_values=false)
SparseMatrix(SparseMatrix< Number, MemorySpace > &&other) noexcept
void clear_rows(const ArrayView< const size_type > &rows, const Number new_diag_value=0)
const_iterator end() const
void set(const std::vector< size_type > &indices, const FullMatrix< Number > &full_matrix, const bool elide_zero_values=false)
Number residual(Vector< Number, MemorySpace > &dst, const Vector< Number, MemorySpace > &x, const Vector< Number, MemorySpace > &b) const
Teuchos::RCP< MatrixType > matrix
std::pair< size_type, size_type > local_range() const
void set(const size_type row, const std::vector< size_type > &col_indices, const std::vector< Number > &values, const bool elide_zero_values=false)
Number frobenius_norm() const
void Tvmult(InputVectorType &dst, const InputVectorType &src) const
Tpetra::CrsMatrix< Number, int, ::types::signed_global_dof_index, NodeType > MatrixType
Number matrix_scalar_product(const Vector< Number, MemorySpace > &u, const Vector< Number, MemorySpace > &v) const
void reinit(const SparsityPatternType &sparsity_pattern)
SparseMatrix & operator=(const double d)
Number element(const size_type i, const size_type j, const bool no_error) const
SparseMatrix(const size_type m, const size_type n, const unsigned int n_max_entries_per_row)
SparseMatrix & operator/=(const Number factor)
Teuchos::RCP< const MatrixType > trilinos_rcp() const
SparseMatrix & operator*=(const Number factor)
bool in_local_range(const size_type index) const
MPI_Comm get_mpi_communicator() const
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const Number *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
const_iterator begin(const size_type r) const
const_iterator begin() const
SparseMatrix(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const MPI_Comm communicator, const std::vector< unsigned int > &n_entries_per_row)
SparseMatrix(const SparsityPattern< MemorySpace > &sparsity_pattern)
SparseMatrix(const IndexSet ¶llel_partitioning, const MPI_Comm communicator=MPI_COMM_WORLD, const unsigned int n_max_entries_per_row=0)
Number diag_element(const size_type i) const
const_iterator end(const size_type r) const
Number operator()(const size_type i, const size_type j) const
void vmult(InputVectorType &dst, const InputVectorType &src) const
SparseMatrix(const size_type m, const size_type n, const std::vector< unsigned int > &n_entries_per_row)
void set(const size_type i, const size_type j, const Number value)
void copy_from(const SparseMatrix< Number, MemorySpace > &source)
iterator end(const size_type r)
SparseMatrix(const IndexSet ¶llel_partitioning, const MPI_Comm communicator, const std::vector< unsigned int > &n_entries_per_row)
std::enable_if_t< !std::is_same_v< SparsityPatternType, ::SparseMatrix< double > > > reinit(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator=MPI_COMM_WORLD, const bool exchange_data=false)
SparseMatrix< Number, MemorySpace > & operator=(const SparseMatrix< Number, MemorySpace > &)=delete
IndexSet locally_owned_range_indices() const
Tpetra::Map< int, ::types::signed_global_dof_index, NodeType > MapType
size_t n_nonzero_elements() const
void Tvmult_add(InputVectorType &dst, const InputVectorType &src) const
void reinit(const SparsityPattern< MemorySpace > &sparsity_pattern)
Number el(const size_type i, const size_type j) const
void vmult_add(InputVectorType &dst, const InputVectorType &src) const
void compress(VectorOperation::values operation)
std::enable_if_t< !std::is_same_v< SparsityPatternType, ::SparseMatrix< double > > > reinit(const IndexSet ¶llel_partitioning, const SparsityPatternType &sparsity_pattern, const MPI_Comm communicator=MPI_COMM_WORLD, const bool exchange_data=false)
void add(const Number factor, const SparseMatrix< Number, MemorySpace > &matrix)
MatrixType & trilinos_matrix()
void reinit(const IndexSet &row_parallel_partitioning, const IndexSet &col_parallel_partitioning, const ::SparseMatrix< Number > &dealii_sparse_matrix, const MPI_Comm communicator=MPI_COMM_WORLD, const double drop_tolerance=1e-13, const bool copy_values=true, const ::SparsityPattern *use_this_sparsity=nullptr)
void set(const size_type row, const size_type n_cols, const size_type *col_indices, const Number *values, const bool elide_zero_values=false)
bool is_compressed() const
Number matrix_norm_square(const Vector< Number, MemorySpace > &v) const
unsigned int local_size() const
Teuchos::RCP< MapType > column_space_map
Teuchos::RCP< MatrixType > trilinos_rcp()
real_type l2_norm() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcDomainMapMismatch()
static ::ExceptionBase & ExcColMapMismatch()
#define DeclException0(Exception0)
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
#define Assert(cond, exc)
static ::ExceptionBase & ExcMatrixNotCompressed()
static ::ExceptionBase & ExcIteratorPastEnd()
static ::ExceptionBase & ExcAccessToNonlocalRow(std::size_t arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
static ::ExceptionBase & ExcBeyondEndOfMatrix()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcInvalidIndexWithinRow(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcSourceEqualsDestination()
int signed_global_dof_index
unsigned int global_dof_index
static const bool zero_addition_can_be_elided
forward_iterator_tag iterator_category
typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::Iterator< Number, MemorySpace, Constness >::difference_type difference_type
typename ::LinearAlgebra::TpetraWrappers::SparseMatrixIterators::Iterator< Number, MemorySpace, Constness >::value_type value_type