15#ifndef dealii_trilinos_tpetra_vector_h
16#define dealii_trilinos_tpetra_vector_h
25#ifdef DEAL_II_TRILINOS_WITH_TPETRA
36# include <Teuchos_Comm.hpp>
37# include <Teuchos_OrdinalTraits.hpp>
38# include <Tpetra_Core.hpp>
39# include <Tpetra_Vector.hpp>
40# include <Tpetra_Version.hpp>
52template <
typename Number>
56# ifdef HAVE_TPETRA_INST_FLOAT
62# ifdef HAVE_TPETRA_INST_DOUBLE
68# ifdef DEAL_II_WITH_COMPLEX_VALUES
69# ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
75# ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
86 template <
typename Number>
87 class ReadWriteVector;
104 namespace TpetraWrappers
143 template <
typename Number,
145 class VectorReference
153 const size_type index);
159 VectorReference(
const VectorReference &) =
default;
172 const VectorReference &
173 operator=(
const VectorReference &r)
const;
179 operator=(
const VectorReference &r);
184 const VectorReference &
185 operator=(
const Number &s)
const;
190 const VectorReference &
191 operator+=(
const Number &s)
const;
196 const VectorReference &
197 operator-=(
const Number &s)
const;
202 const VectorReference &
203 operator*=(
const Number &s)
const;
208 const VectorReference &
209 operator/=(
const Number &s)
const;
215 operator Number()
const;
222 <<
"An error with error number " << arg1
223 <<
" occurred while calling a Trilinos function");
231 ExcAccessToNonLocalElement,
236 <<
"You are trying to access element " << arg1
237 <<
" of a distributed vector, but this element is not stored "
238 <<
"on the current processor. Note: There are " << arg2
239 <<
" elements stored "
240 <<
"on the current processor from within the range [" << arg3 <<
','
241 << arg4 <<
"] but Trilinos vectors need not store contiguous "
242 <<
"ranges on each processor, and not every element in "
243 <<
"this range may in fact be stored locally."
245 <<
"A common source for this kind of problem is that you "
246 <<
"are passing a 'fully distributed' vector into a function "
247 <<
"that needs read access to vector elements that correspond "
248 <<
"to degrees of freedom on ghost cells (or at least to "
249 <<
"'locally active' degrees of freedom that are not also "
250 <<
"'locally owned'). You need to pass a vector that has these "
251 <<
"elements as ghost entries.");
262 const size_type index;
286 template <
typename Number,
typename MemorySpace = ::MemorySpace::Host>
296 using reference = internal::VectorReference<Number, MemorySpace>;
298 const internal::VectorReference<Number, MemorySpace>;
461 template <
typename OtherNumber>
484 const Teuchos::RCP<const Utilities::MPI::CommunicationPatternBase>
495 const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
515 std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
632 add(
const std::vector<size_type> &indices,
633 const std::vector<Number> &values);
644 const Number *values);
664 const Number *values);
833 std::pair<size_type, size_type>
917 Teuchos::RCP<const TpetraTypes::VectorType<Number, MemorySpace>>
924 Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>>
932 const unsigned int precision = 3,
933 const bool scientific =
true,
934 const bool across =
true)
const;
977 <<
"You are trying to access element " <<
arg1
978 <<
" of a distributed vector, but this element is not stored "
979 <<
"on the current processor. Note: There are " <<
arg2
980 <<
" elements stored "
981 <<
"on the current processor from within the range [" <<
arg3 <<
','
982 <<
arg4 <<
"] but Trilinos vectors need not store contiguous "
983 <<
"ranges on each processor, and not every element in "
984 <<
"this range may in fact be stored locally."
986 <<
"A common source for this kind of problem is that you "
987 <<
"are passing a 'fully distributed' vector into a function "
988 <<
"that needs read access to vector elements that correspond "
989 <<
"to degrees of freedom on ghost cells (or at least to "
990 <<
"'locally active' degrees of freedom that are not also "
991 <<
"'locally owned'). You need to pass a vector that has these "
992 <<
"elements as ghost entries.");
1000 "To compress a vector, a locally_relevant_dofs "
1001 "index set, and a locally_owned_dofs index set "
1002 "must be provided. These index sets must be "
1003 "provided either when the vector is initialized "
1004 "or when compress is called. See the documentation "
1005 "of compress() for more information.");
1014 <<
"An error with error number " <<
arg1
1015 <<
" occurred while calling a Trilinos function");
1046 Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>>
vector;
1053 Teuchos::RCP<TpetraTypes::VectorType<Number, MemorySpace>>
1065 Teuchos::RCP<const TpetraWrappers::CommunicationPattern<MemorySpace>>
1069 friend class internal::VectorReference<Number,
MemorySpace>;
1075 template <
typename Number,
typename MemorySpace>
1084 template <
typename Number,
typename MemorySpace>
1093 template <
typename Number,
typename MemorySpace>
1100 template <
typename Number,
typename MemorySpace>
1104 vector.swap(v.vector);
1108 template <
typename Number,
typename MemorySpace>
1111 const std::vector<Number> &values)
1117 add(indices.size(), indices.data(), values.data());
1122 template <
typename Number,
typename MemorySpace>
1126 const Number *values)
1132 auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
1133 Tpetra::Access::ReadWrite);
1141 auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0);
1142 using ViewType1d =
decltype(vector_1d_local);
1143 std::optional<ViewType1d> vector_1d_nonlocal;
1145# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1148 vector->template modify<Kokkos::HostSpace>();
1151 for (
size_type i = 0; i < n_elements; ++i)
1159 vector->getMap()->getLocalElement(row);
1160 local_row != Teuchos::OrdinalTraits<int>::invalid())
1162 vector_1d_local(local_row) += values[i];
1167 if (nonlocal_vector.get() !=
nullptr)
1177 nonlocal_vector->getMap()->getLocalElement(row);
1179# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1180 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1181 ExcAccessToNonLocalElement(
1183 vector->getMap()->getLocalNumElements(),
1184 vector->getMap()->getMinLocalIndex(),
1185 vector->getMap()->getMaxLocalIndex()));
1187 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1188 ExcAccessToNonLocalElement(
1190 vector->getMap()->getNodeNumElements(),
1191 vector->getMap()->getMinLocalIndex(),
1192 vector->getMap()->getMaxLocalIndex()));
1200 if (!vector_1d_nonlocal)
1202# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1203 auto vector_2d_nonlocal =
1204 nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
1205 Tpetra::Access::ReadWrite);
1207 auto vector_2d_nonlocal =
1208 nonlocal_vector->template getLocalView<Kokkos::HostSpace>();
1211 vector_1d_nonlocal =
1212 Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0);
1214# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1216 nonlocal_vector->template modify<Kokkos::HostSpace>();
1219 (*vector_1d_nonlocal)(nonlocal_row) += values[i];
1224# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1225 vector->template sync<
1226 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1227 device_type::memory_space>();
1231 if (vector_1d_nonlocal)
1232 nonlocal_vector->template sync<
1233 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1234 device_type::memory_space>();
1240 template <
typename Number,
typename MemorySpace>
1244 const Number *values)
1250# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1251 auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
1252 Tpetra::Access::ReadWrite);
1254 vector->template sync<Kokkos::HostSpace>();
1256 auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>();
1265 auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0);
1266 using ViewType1d =
decltype(vector_1d_local);
1267 std::optional<ViewType1d> vector_1d_nonlocal;
1269# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1272 vector->template modify<Kokkos::HostSpace>();
1275 for (
size_type i = 0; i < n_elements; ++i)
1283 vector->getMap()->getLocalElement(row);
1284 local_row != Teuchos::OrdinalTraits<int>::invalid())
1286 vector_1d_local(local_row) = values[i];
1291 if (nonlocal_vector.get() !=
nullptr)
1301 nonlocal_vector->getMap()->getLocalElement(row);
1303# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1304 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1305 ExcAccessToNonLocalElement(
1307 vector->getMap()->getLocalNumElements(),
1308 vector->getMap()->getMinLocalIndex(),
1309 vector->getMap()->getMaxLocalIndex()));
1311 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1312 ExcAccessToNonLocalElement(
1314 vector->getMap()->getNodeNumElements(),
1315 vector->getMap()->getMinLocalIndex(),
1316 vector->getMap()->getMaxLocalIndex()));
1324 if (!vector_1d_nonlocal)
1326# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1327 auto vector_2d_nonlocal =
1328 nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
1329 Tpetra::Access::ReadWrite);
1331 auto vector_2d_nonlocal =
1332 nonlocal_vector->template getLocalView<Kokkos::HostSpace>();
1335 vector_1d_nonlocal =
1336 Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0);
1338# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1340 nonlocal_vector->template modify<Kokkos::HostSpace>();
1343 (*vector_1d_nonlocal)(nonlocal_row) = values[i];
1348# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1349 vector->template sync<
1350 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1351 device_type::memory_space>();
1355 if (vector_1d_nonlocal)
1356 nonlocal_vector->template sync<
1357 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1358 device_type::memory_space>();
1364 template <
typename Number,
typename MemorySpace>
1365 inline internal::VectorReference<Number, MemorySpace>
1368 return internal::VectorReference(*
this, index);
1371 template <
typename Number,
typename MemorySpace>
1372 inline internal::VectorReference<Number, MemorySpace>
1375 return operator()(index);
1378 template <
typename Number,
typename MemorySpace>
1382 return operator()(index);
1390 template <
typename Number,
typename MemorySpace>
1391 inline VectorReference<Number, MemorySpace>::VectorReference(
1393 const size_type index)
1400 template <
typename Number,
typename MemorySpace>
1401 inline const VectorReference<Number, MemorySpace> &
1402 VectorReference<Number, MemorySpace>::operator=(
1403 const VectorReference<Number, MemorySpace> &r)
const
1409 *
this =
static_cast<Number
>(r);
1416 template <
typename Number,
typename MemorySpace>
1417 inline VectorReference<Number, MemorySpace> &
1418 VectorReference<Number, MemorySpace>::operator=(
1419 const VectorReference<Number, MemorySpace> &r)
1422 *
this =
static_cast<Number
>(r);
1429 template <
typename Number,
typename MemorySpace>
1430 inline const VectorReference<Number, MemorySpace> &
1431 VectorReference<Number, MemorySpace>::operator=(
const Number &value)
const
1433 vector.set(1, &index, &value);
1439 template <
typename Number,
typename MemorySpace>
1440 inline const VectorReference<Number, MemorySpace> &
1441 VectorReference<Number, MemorySpace>::operator+=(
1442 const Number &value)
const
1444 vector.add(1, &index, &value);
1450 template <
typename Number,
typename MemorySpace>
1451 inline const VectorReference<Number, MemorySpace> &
1452 VectorReference<Number, MemorySpace>::operator-=(
1453 const Number &value)
const
1455 Number new_value = -
value;
1456 vector.add(1, &index, &new_value);
1462 template <
typename Number,
typename MemorySpace>
1463 inline const VectorReference<Number, MemorySpace> &
1464 VectorReference<Number, MemorySpace>::operator*=(
1465 const Number &value)
const
1467 Number new_value =
static_cast<Number
>(*this) *
value;
1468 vector.set(1, &index, &new_value);
1474 template <
typename Number,
typename MemorySpace>
1475 inline const VectorReference<Number, MemorySpace> &
1476 VectorReference<Number, MemorySpace>::operator/=(
1477 const Number &value)
const
1479 Number new_value =
static_cast<Number
>(*this) /
value;
1480 vector.set(1, &index, &new_value);
1496 namespace LinearOperatorImplementation
1505 template <
typename Number,
typename MemorySpace>
1510 template <
typename Matrix>
1513 const Matrix &matrix,
1515 bool omit_zeroing_entries)
1517 v.
reinit(matrix.locally_owned_range_indices(),
1518 matrix.get_mpi_communicator(),
1519 omit_zeroing_entries);
1522 template <
typename Matrix>
1525 const Matrix &matrix,
1527 bool omit_zeroing_entries)
1529 v.
reinit(matrix.locally_owned_domain_indices(),
1530 matrix.get_mpi_communicator(),
1531 omit_zeroing_entries);
1541template <
typename Number,
typename MemorySpace>
1543 LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace>> : std::false_type
632 add(
const std::vector<size_type> &indices, {
…}
void reinit(const Vector< Number, MemorySpace > &V, const bool omit_zeroing_entries=false)
void equ(const Number a, const Vector< Number, MemorySpace > &V)
Teuchos::RCP< const TpetraWrappers::CommunicationPattern< MemorySpace > > tpetra_comm_pattern
void reinit(const IndexSet &locally_owned_entries, const IndexSet &locally_relevant_or_ghost_entries, const MPI_Comm communicator=MPI_COMM_WORLD, const bool vector_writable=false)
void add(const Number a, const Vector< Number, MemorySpace > &V, const Number b, const Vector< Number, MemorySpace > &W)
bool is_non_negative() const
Number add_and_dot(const Number a, const Vector< Number, MemorySpace > &V, const Vector< Number, MemorySpace > &W)
TpetraTypes::VectorType< Number, MemorySpace > & trilinos_vector()
void reinit(const IndexSet ¶llel_partitioner, const MPI_Comm communicator=MPI_COMM_WORLD, const bool omit_zeroing_entries=false)
MPI_Comm get_mpi_communicator() const
Teuchos::RCP< TpetraTypes::VectorType< Number, MemorySpace > > vector
Number operator[](const size_type index) const
bool operator==(const Vector< Number, MemorySpace > &v) const
Vector(const IndexSet ¶llel_partitioner, const MPI_Comm communicator)
void scale(const Vector< Number, MemorySpace > &scaling_factors)
void set(const size_type n_elements, const size_type *indices, const Number *values)
void import_elements(const ReadWriteVector< Number > &V, VectorOperation::values operation, const std::shared_ptr< const Utilities::MPI::CommunicationPatternBase > &communication_pattern)
MPI_Comm mpi_comm() const
Vector(const Teuchos::RCP< TpetraTypes::VectorType< Number, MemorySpace > > V)
internal::VectorReference< Number, MemorySpace > reference
size_type locally_owned_size() const
reference operator()(const size_type index)
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, const ArrayView< Number > &elements) const override
const internal::VectorReference< Number, MemorySpace > const_reference
void compress(const VectorOperation::values operation)
bool is_compressed() const
bool operator!=(const Vector< Number, MemorySpace > &v) const
Number operator()(const size_type index) const
Vector & operator/=(const Number factor)
std::pair< size_type, size_type > local_range() const
::IndexSet source_stored_elements
void sadd(const Number s, const Number a, const Vector< Number, MemorySpace > &V)
reference operator[](const size_type index)
void add(const Number a, const Vector< Number, MemorySpace > &V)
Vector & operator*=(const Number factor)
void add(const size_type n_elements, const size_type *indices, const Number *values)
void import_elements(const ReadWriteVector< Number > &V, VectorOperation::values operation, const Teuchos::RCP< const Utilities::MPI::CommunicationPatternBase > &communication_pattern)
real_type l2_norm() const
void add(const std::vector< size_type > &indices, const std::vector< Number > &values)
typename numbers::NumberTraits< Number >::real_type real_type
Vector & operator-=(const Vector< Number, MemorySpace > &V)
Vector & operator=(const ::Vector< OtherNumber > &V)
::IndexSet locally_owned_elements() const
Vector & operator=(const Vector &V)
Vector & operator=(const Number s)
Vector & operator+=(const Vector< Number, MemorySpace > &V)
const TpetraTypes::VectorType< Number, MemorySpace > & trilinos_vector() const
real_type linfty_norm() const
real_type norm_sqr() const
virtual size_type size() const override
void import_elements(const ReadWriteVector< Number > &V, VectorOperation::values operation)
void create_tpetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm)
real_type l1_norm() const
Teuchos::RCP< const TpetraTypes::VectorType< Number, MemorySpace > > trilinos_rcp() const
Number operator*(const Vector< Number, MemorySpace > &V) const
Number mean_value() const
std::size_t memory_consumption() const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
bool has_ghost_elements() const
bool in_local_range(const size_type index) const
Vector(const IndexSet &locally_owned_entries, const IndexSet &ghost_entries, const MPI_Comm communicator, const bool vector_writable=false)
Teuchos::RCP< TpetraTypes::VectorType< Number, MemorySpace > > nonlocal_vector
Teuchos::RCP< TpetraTypes::VectorType< Number, MemorySpace > > trilinos_rcp()
virtual void swap(Vector &v) noexcept
static void reinit_range_vector(const Matrix &matrix, LinearAlgebra::TpetraWrappers::Vector< Number, MemorySpace > &v, bool omit_zeroing_entries)
static void reinit_domain_vector(const Matrix &matrix, LinearAlgebra::TpetraWrappers::Vector< Number, MemorySpace > &v, bool omit_zeroing_entries)
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define DeclException0(Exception0)
static ::ExceptionBase & ExcGhostsPresent()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
static ::ExceptionBase & ExcVectorTypeNotCompatible()
static ::ExceptionBase & ExcAccessToNonLocalElement(size_type arg1, size_type arg2, size_type arg3, size_type arg4)
#define Assert(cond, exc)
static ::ExceptionBase & ExcMissingIndexSet()
#define AssertDimension(dim1, dim2)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcDifferentParallelPartitioning()
static ::ExceptionBase & ExcTrilinosError(int arg1)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
void swap(Vector< Number, MemorySpace > &u, Vector< Number, MemorySpace > &v) noexcept
unsigned int global_dof_index