15#ifndef dealii_trilinos_tpetra_vector_h
16#define dealii_trilinos_tpetra_vector_h
21#ifdef DEAL_II_TRILINOS_WITH_TPETRA
33# include <Teuchos_Comm.hpp>
34# include <Teuchos_OrdinalTraits.hpp>
35# include <Tpetra_Core.hpp>
36# include <Tpetra_Vector.hpp>
37# include <Tpetra_Version.hpp>
49template <
typename Number>
53# ifdef HAVE_TPETRA_INST_FLOAT
59# ifdef HAVE_TPETRA_INST_DOUBLE
65# ifdef DEAL_II_WITH_COMPLEX_VALUES
66# ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
72# ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
83 template <
typename Number>
84 class ReadWriteVector;
101 namespace TpetraWrappers
139 template <
typename Number,
141 class VectorReference
149 const size_type index);
155 VectorReference(
const VectorReference &) =
default;
168 const VectorReference &
169 operator=(
const VectorReference &r)
const;
175 operator=(
const VectorReference &r);
180 const VectorReference &
181 operator=(
const Number &s)
const;
186 const VectorReference &
187 operator+=(
const Number &s)
const;
192 const VectorReference &
193 operator-=(
const Number &s)
const;
198 const VectorReference &
199 operator*=(
const Number &s)
const;
204 const VectorReference &
205 operator/=(
const Number &s)
const;
211 operator Number()
const;
218 <<
"An error with error number " << arg1
219 <<
" occurred while calling a Trilinos function");
227 ExcAccessToNonLocalElement,
232 <<
"You are trying to access element " << arg1
233 <<
" of a distributed vector, but this element is not stored "
234 <<
"on the current processor. Note: There are " << arg2
235 <<
" elements stored "
236 <<
"on the current processor from within the range [" << arg3 <<
','
237 << arg4 <<
"] but Trilinos vectors need not store contiguous "
238 <<
"ranges on each processor, and not every element in "
239 <<
"this range may in fact be stored locally."
241 <<
"A common source for this kind of problem is that you "
242 <<
"are passing a 'fully distributed' vector into a function "
243 <<
"that needs read access to vector elements that correspond "
244 <<
"to degrees of freedom on ghost cells (or at least to "
245 <<
"'locally active' degrees of freedom that are not also "
246 <<
"'locally owned'). You need to pass a vector that has these "
247 <<
"elements as ghost entries.");
258 const size_type index;
282 template <
typename Number,
typename MemorySpace = ::MemorySpace::Host>
289# if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0)
290 using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode<
291 typename MemorySpace::kokkos_space::execution_space,
292 typename MemorySpace::kokkos_space>;
294 using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode<
295 typename MemorySpace::kokkos_space::execution_space,
296 typename MemorySpace::kokkos_space>;
301 using reference = internal::VectorReference<Number, MemorySpace>;
303 const internal::VectorReference<Number, MemorySpace>;
305 Tpetra::Map<int, ::types::signed_global_dof_index, NodeType>;
307 Vector<Number, int, ::types::signed_global_dof_index, NodeType>;
309 Tpetra::Export<int, ::types::signed_global_dof_index, NodeType>;
311 Tpetra::Import<int, ::types::signed_global_dof_index, NodeType>;
333 Vector(
const Teuchos::RCP<VectorType> V);
473 template <
typename OtherNumber>
496 const Teuchos::RCP<const Utilities::MPI::CommunicationPatternBase>
507 const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
527 std::shared_ptr<const Utilities::MPI::CommunicationPatternBase>
644 add(
const std::vector<size_type> &indices,
645 const std::vector<Number> &values);
656 const Number *values);
676 const Number *values);
845 std::pair<size_type, size_type>
915 const Tpetra::Vector<Number, int, types::signed_global_dof_index> &
922 Tpetra::Vector<Number, int, types::signed_global_dof_index> &
930 const Tpetra::Vector<Number, int, types::signed_global_dof_index>>
937 Teuchos::RCP<Tpetra::Vector<Number, int, types::signed_global_dof_index>>
945 const unsigned int precision = 3,
946 const bool scientific =
true,
947 const bool across =
true)
const;
990 <<
"You are trying to access element " <<
arg1
991 <<
" of a distributed vector, but this element is not stored "
992 <<
"on the current processor. Note: There are " <<
arg2
993 <<
" elements stored "
994 <<
"on the current processor from within the range [" <<
arg3 <<
','
995 <<
arg4 <<
"] but Trilinos vectors need not store contiguous "
996 <<
"ranges on each processor, and not every element in "
997 <<
"this range may in fact be stored locally."
999 <<
"A common source for this kind of problem is that you "
1000 <<
"are passing a 'fully distributed' vector into a function "
1001 <<
"that needs read access to vector elements that correspond "
1002 <<
"to degrees of freedom on ghost cells (or at least to "
1003 <<
"'locally active' degrees of freedom that are not also "
1004 <<
"'locally owned'). You need to pass a vector that has these "
1005 <<
"elements as ghost entries.");
1013 "To compress a vector, a locally_relevant_dofs "
1014 "index set, and a locally_owned_dofs index set "
1015 "must be provided. These index sets must be "
1016 "provided either when the vector is initialized "
1017 "or when compress is called. See the documentation "
1018 "of compress() for more information.");
1027 <<
"An error with error number " <<
arg1
1028 <<
" occurred while calling a Trilinos function");
1077 Teuchos::RCP<const TpetraWrappers::CommunicationPattern>
1081 friend class internal::VectorReference<Number,
MemorySpace>;
1087 template <
typename Number,
typename MemorySpace>
1095 template <
typename Number,
typename MemorySpace>
1104 template <
typename Number,
typename MemorySpace>
1111 template <
typename Number,
typename MemorySpace>
1119 template <
typename Number,
typename MemorySpace>
1122 const std::vector<Number> &values)
1128 add(indices.size(), indices.data(), values.data());
1133 template <
typename Number,
typename MemorySpace>
1137 const Number *values)
1143# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1144 auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
1145 Tpetra::Access::ReadWrite);
1147 vector->template sync<Kokkos::HostSpace>();
1149 auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>();
1158 auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0);
1159 using ViewType1d =
decltype(vector_1d_local);
1160 std::optional<ViewType1d> vector_1d_nonlocal;
1162# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1165 vector->template modify<Kokkos::HostSpace>();
1168 for (
size_type i = 0; i < n_elements; ++i)
1176 vector->getMap()->getLocalElement(row);
1177 local_row != Teuchos::OrdinalTraits<int>::invalid())
1179 vector_1d_local(local_row) += values[i];
1188 nonlocal_vector->getMap()->getLocalElement(row);
1190# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1191 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1192 ExcAccessToNonLocalElement(
1194 vector->getMap()->getLocalNumElements(),
1195 vector->getMap()->getMinLocalIndex(),
1196 vector->getMap()->getMaxLocalIndex()));
1198 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1199 ExcAccessToNonLocalElement(
1201 vector->getMap()->getNodeNumElements(),
1202 vector->getMap()->getMinLocalIndex(),
1203 vector->getMap()->getMaxLocalIndex()));
1211 if (!vector_1d_nonlocal)
1213# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1214 auto vector_2d_nonlocal =
1215 nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
1216 Tpetra::Access::ReadWrite);
1218 auto vector_2d_nonlocal =
1219 nonlocal_vector->template getLocalView<Kokkos::HostSpace>();
1222 vector_1d_nonlocal =
1223 Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0);
1225# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1227 nonlocal_vector->template modify<Kokkos::HostSpace>();
1230 (*vector_1d_nonlocal)(nonlocal_row) += values[i];
1235# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1236 vector->template sync<
1237 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1238 device_type::memory_space>();
1242 if (vector_1d_nonlocal)
1243 nonlocal_vector->template sync<
1244 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1245 device_type::memory_space>();
1251 template <
typename Number,
typename MemorySpace>
1255 const Number *values)
1261# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1262 auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
1263 Tpetra::Access::ReadWrite);
1265 vector->template sync<Kokkos::HostSpace>();
1267 auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>();
1276 auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0);
1277 using ViewType1d =
decltype(vector_1d_local);
1278 std::optional<ViewType1d> vector_1d_nonlocal;
1280# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1283 vector->template modify<Kokkos::HostSpace>();
1286 for (
size_type i = 0; i < n_elements; ++i)
1294 vector->getMap()->getLocalElement(row);
1295 local_row != Teuchos::OrdinalTraits<int>::invalid())
1297 vector_1d_local(local_row) = values[i];
1306 nonlocal_vector->getMap()->getLocalElement(row);
1308# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
1309 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1310 ExcAccessToNonLocalElement(
1312 vector->getMap()->getLocalNumElements(),
1313 vector->getMap()->getMinLocalIndex(),
1314 vector->getMap()->getMaxLocalIndex()));
1316 Assert(nonlocal_row != Teuchos::OrdinalTraits<int>::invalid(),
1317 ExcAccessToNonLocalElement(
1319 vector->getMap()->getNodeNumElements(),
1320 vector->getMap()->getMinLocalIndex(),
1321 vector->getMap()->getMaxLocalIndex()));
1329 if (!vector_1d_nonlocal)
1331# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1332 auto vector_2d_nonlocal =
1333 nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
1334 Tpetra::Access::ReadWrite);
1336 auto vector_2d_nonlocal =
1337 nonlocal_vector->template getLocalView<Kokkos::HostSpace>();
1340 vector_1d_nonlocal =
1341 Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0);
1343# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1345 nonlocal_vector->template modify<Kokkos::HostSpace>();
1348 (*vector_1d_nonlocal)(nonlocal_row) = values[i];
1353# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
1354 vector->template sync<
1355 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1356 device_type::memory_space>();
1360 if (vector_1d_nonlocal)
1361 nonlocal_vector->template sync<
1362 typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
1363 device_type::memory_space>();
1369 template <
typename Number,
typename MemorySpace>
1370 inline internal::VectorReference<Number, MemorySpace>
1373 return internal::VectorReference(*
this, index);
1376 template <
typename Number,
typename MemorySpace>
1377 inline internal::VectorReference<Number, MemorySpace>
1380 return operator()(index);
1383 template <
typename Number,
typename MemorySpace>
1387 return operator()(index);
1395 template <
typename Number,
typename MemorySpace>
1396 inline VectorReference<Number, MemorySpace>::VectorReference(
1398 const size_type index)
1405 template <
typename Number,
typename MemorySpace>
1406 inline const VectorReference<Number, MemorySpace> &
1407 VectorReference<Number, MemorySpace>::operator=(
1408 const VectorReference<Number, MemorySpace> &r)
const
1414 *
this =
static_cast<Number
>(r);
1421 template <
typename Number,
typename MemorySpace>
1422 inline VectorReference<Number, MemorySpace> &
1423 VectorReference<Number, MemorySpace>::operator=(
1424 const VectorReference<Number, MemorySpace> &r)
1427 *
this =
static_cast<Number
>(r);
1434 template <
typename Number,
typename MemorySpace>
1435 inline const VectorReference<Number, MemorySpace> &
1436 VectorReference<Number, MemorySpace>::operator=(
const Number &value)
const
1438 vector.set(1, &index, &value);
1444 template <
typename Number,
typename MemorySpace>
1445 inline const VectorReference<Number, MemorySpace> &
1446 VectorReference<Number, MemorySpace>::operator+=(
1447 const Number &value)
const
1449 vector.add(1, &index, &value);
1455 template <
typename Number,
typename MemorySpace>
1456 inline const VectorReference<Number, MemorySpace> &
1457 VectorReference<Number, MemorySpace>::operator-=(
1458 const Number &value)
const
1460 Number new_value = -
value;
1461 vector.add(1, &index, &new_value);
1467 template <
typename Number,
typename MemorySpace>
1468 inline const VectorReference<Number, MemorySpace> &
1469 VectorReference<Number, MemorySpace>::operator*=(
1470 const Number &value)
const
1472 Number new_value =
static_cast<Number
>(*this) *
value;
1473 vector.set(1, &index, &new_value);
1479 template <
typename Number,
typename MemorySpace>
1480 inline const VectorReference<Number, MemorySpace> &
1481 VectorReference<Number, MemorySpace>::operator/=(
1482 const Number &value)
const
1484 Number new_value =
static_cast<Number
>(*this) /
value;
1485 vector.set(1, &index, &new_value);
1501 namespace LinearOperatorImplementation
1510 template <
typename Number,
typename MemorySpace>
1515 template <
typename Matrix>
1518 const Matrix &matrix,
1520 bool omit_zeroing_entries)
1522 v.
reinit(matrix.locally_owned_range_indices(),
1523 matrix.get_mpi_communicator(),
1524 omit_zeroing_entries);
1527 template <
typename Matrix>
1530 const Matrix &matrix,
1532 bool omit_zeroing_entries)
1534 v.
reinit(matrix.locally_owned_domain_indices(),
1535 matrix.get_mpi_communicator(),
1536 omit_zeroing_entries);
1546template <
typename Number,
typename MemorySpace>
1548 LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace>> : std::false_type
::types::global_dof_index size_type
void reinit(const Vector< Number, MemorySpace > &V, const bool omit_zeroing_entries=false)
void equ(const Number a, const Vector< Number, MemorySpace > &V)
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)
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
Tpetra::Export< int, ::types::signed_global_dof_index, NodeType > ExportType
Number operator[](const size_type index) const
bool operator==(const Vector< Number, MemorySpace > &v) const
Vector(const IndexSet ¶llel_partitioner, const MPI_Comm communicator)
virtual void swap(Vector &v)
void scale(const Vector< Number, MemorySpace > &scaling_factors)
Teuchos::RCP< Tpetra::Vector< Number, int, types::signed_global_dof_index > > trilinos_rcp()
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
internal::VectorReference< Number, MemorySpace > reference
size_type locally_owned_size() const
Teuchos::RCP< VectorType > vector
reference operator()(const size_type index)
Tpetra::Import< int, ::types::signed_global_dof_index, NodeType > ImportType
Teuchos::RCP< const Tpetra::Vector< Number, int, types::signed_global_dof_index > > trilinos_rcp() const
const internal::VectorReference< Number, MemorySpace > const_reference
Teuchos::RCP< const TpetraWrappers::CommunicationPattern > tpetra_comm_pattern
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
virtual void extract_subvector_to(const ArrayView< const types::global_dof_index > &indices, ArrayView< Number > &elements) const override
Teuchos::RCP< VectorType > nonlocal_vector
::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)
const Tpetra::Vector< Number, int, types::signed_global_dof_index > & trilinos_vector() const
Vector & operator=(const ::Vector< OtherNumber > &V)
::IndexSet locally_owned_elements() const
Tpetra::Vector< Number, int, ::types::signed_global_dof_index, NodeType > VectorType
Vector & operator=(const Vector &V)
Vector & operator=(const Number s)
Vector & operator+=(const Vector< Number, MemorySpace > &V)
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)
Tpetra::Vector< Number, int, types::signed_global_dof_index > & trilinos_vector()
void create_tpetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm)
real_type l1_norm() const
Vector(const Teuchos::RCP< VectorType > V)
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
Tpetra::Map< int, ::types::signed_global_dof_index, NodeType > MapType
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)
Tpetra::KokkosCompat::KokkosDeviceWrapperNode< typename MemorySpace::kokkos_space::execution_space, typename MemorySpace::kokkos_space > NodeType
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)
void swap(Vector< Number, MemorySpace > &u, Vector< Number, MemorySpace > &v)
unsigned int global_dof_index