Vector(const Teuchos::RCP<VectorType> V);
/**
+ * TODO: This is not used
* This constructor takes an IndexSet that defines how to distribute the
* individual components among the MPI processors. Since it also
* includes information about the size of the vector, this is all we
explicit Vector(const IndexSet ¶llel_partitioner,
const MPI_Comm communicator);
+ /**
+ * In addition to just specifying one index set as in all the other
+ * methods above, this method allows to supply an additional set of
+ * ghost entries.
+ *
+ * Depending on whether the @p locally_relevant_or_ghost_entries argument uniquely
+ * subdivides elements among processors or not, the resulting vector may
+ * or may not have ghost elements. See the general documentation of this
+ * class for more information.
+ *
+ * @see
+ * @ref GlossGhostedVector "vectors with ghost elements"
+
+ */
+ explicit Vector(const IndexSet &locally_owned_entries,
+ const IndexSet &ghost_entries,
+ const MPI_Comm communicator,
+ const bool vector_writable = false);
+
/**
* Reinit functionality. This function destroys the old vector content
* and generates a new one based on the input partitioning. The flag
const MPI_Comm communicator,
const bool omit_zeroing_entries = false);
-
/**
* Reinit functionality. This function destroys the old vector content
* and generates a new one based on the input partitioning. In addition
void
reinit(const IndexSet &locally_owned_entries,
const IndexSet &locally_relevant_or_ghost_entries,
- const MPI_Comm communicator = MPI_COMM_WORLD);
+ const MPI_Comm communicator = MPI_COMM_WORLD,
+ const bool vector_writable = false);
/**
* Change the dimension to that of the vector @p V. The elements of @p V are not
create_tpetra_comm_pattern(const IndexSet &source_index_set,
const MPI_Comm mpi_comm);
+ /**
+ * A boolean variable to hold information on whether the vector is
+ * compressed or not.
+ */
+ bool compressed;
+
/**
* Store whether the vector has ghost elements or not.
*
Teuchos::RCP<VectorType> vector;
/**
- * Pointer to the Tpetra::Map storing the owned elemenets
- * per rank. This map has a one-to-one mapping.
- */
- Teuchos::RCP<MapType> locally_owned_map;
-
- /**
- * Pointer to the Tpetra::Map storing the relevant elements
- * per rank. There are entries in this map, that belong to multiple
- * ranks.
+ * A vector object in Trilinos to be used for collecting the non-local
+ * elements if the vector was constructed with an additional IndexSet
+ * describing ghost elements.
*/
- Teuchos::RCP<MapType> locally_relevant_map;
+ Teuchos::RCP<VectorType> nonlocal_vector;
/**
* IndexSet of the elements of the last imported vector.
inline bool
Vector<Number>::is_compressed() const
{
- return !has_ghost;
+ return compressed;
}
const Number *values)
{
# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
- auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>(
+ auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
Tpetra::Access::ReadWrite);
+ auto vector_2d_nonlocal =
+ nonlocal_vector->template getLocalView<Kokkos::HostSpace>(
+ Tpetra::Access::ReadWrite);
# else
vector->template sync<Kokkos::HostSpace>();
- auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>();
+ auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>();
+ auto vector_2d_nonlocal =
+ nonlocal_vector->template getLocalView<Kokkos::HostSpace>();
# endif
- auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+
+ auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0);
+ auto vector_1d_nonlocal =
+ Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0);
# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
vector->template modify<Kokkos::HostSpace>();
+ nonlocal_vector->template modify<Kokkos::HostSpace>();
# endif
for (size_type i = 0; i < n_elements; ++i)
TrilinosWrappers::types::int_type local_row =
vector->getMap()->getLocalElement(row);
+ // check if the index is in the non local set
+ bool nonlocal = false;
+ if (local_row == Teuchos::OrdinalTraits<int>::invalid())
+ {
+ local_row = nonlocal_vector->getMap()->getLocalElement(row);
+ nonlocal = true;
+ compressed = false;
+ }
+
# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
Assert(
local_row != Teuchos::OrdinalTraits<int>::invalid(),
# endif
if (local_row != Teuchos::OrdinalTraits<int>::invalid())
- vector_1d(local_row) += values[i];
+ if (nonlocal)
+ vector_1d_nonlocal(local_row) += values[i];
+ else
+ vector_1d_local(local_row) += values[i];
}
# if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
vector->template sync<
typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
device_type::memory_space>();
+ nonlocal_vector->template sync<
+ typename Tpetra::Vector<Number, int, types::signed_global_dof_index>::
+ device_type::memory_space>();
# endif
}
const size_type *indices,
const Number *values)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>(
Tpetra::Access::ReadWrite);
template <typename Number>
Vector<Number>::Vector()
: Subscriptor()
+ , compressed(true)
, has_ghost(false)
, vector(Utilities::Trilinos::internal::make_rcp<VectorType>(
Utilities::Trilinos::internal::make_rcp<MapType>(
template <typename Number>
Vector<Number>::Vector(const Vector<Number> &V)
: Subscriptor()
+ , compressed(V.compressed)
, has_ghost(V.has_ghost)
- , vector(Utilities::Trilinos::internal::make_rcp<VectorType>(
- V.trilinos_vector(),
- Teuchos::Copy))
- , locally_owned_map(V.locally_owned_map)
- , locally_relevant_map(V.locally_relevant_map)
+ , vector(V.vector)
+ , nonlocal_vector(V.nonlocal_vector)
{}
template <typename Number>
Vector<Number>::Vector(const Teuchos::RCP<VectorType> V)
: Subscriptor()
- , has_ghost(false)
+ , compressed(true)
+ , has_ghost(V->getMap()->isOneToOne() == false)
, vector(V)
- {
- if (vector->getMap()->isOneToOne())
- {
- locally_owned_map =
- Teuchos::rcp_const_cast<MapType>(vector->getMap());
- }
- else
- {
- locally_relevant_map =
- Teuchos::rcp_const_cast<MapType>(vector->getMap());
- }
- }
+ {}
Vector<Number>::Vector(const IndexSet ¶llel_partitioner,
const MPI_Comm communicator)
: Subscriptor()
+ , compressed(true)
, has_ghost(false)
- , locally_owned_map(
- parallel_partitioner.make_tpetra_map_rcp(communicator, false))
+ , vector(Utilities::Trilinos::internal::make_rcp<VectorType>(
+ parallel_partitioner.make_tpetra_map_rcp(communicator, true)))
+ {}
+
+
+
+ template <typename Number>
+ Vector<Number>::Vector(const IndexSet &locally_owned_entries,
+ const IndexSet &ghost_entries,
+ const MPI_Comm communicator,
+ const bool vector_writable)
+ : Subscriptor()
{
- vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
- locally_owned_map);
+ if (!vector_writable)
+ {
+ IndexSet parallel_partitioner = locally_owned_entries;
+ parallel_partitioner.add_indices(ghost_entries);
+
+ vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
+ parallel_partitioner.make_tpetra_map_rcp(communicator, true));
+
+ compressed = true;
+ }
+ else
+ {
+ Teuchos::RCP<MapType> map =
+ locally_owned_entries.make_tpetra_map_rcp(communicator, false);
+ vector = Utilities::Trilinos::internal::make_rcp<VectorType>(map);
+
+ IndexSet nonlocal_entries(ghost_entries);
+ nonlocal_entries.subtract_set(locally_owned_entries);
+ nonlocal_vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
+ nonlocal_entries.make_tpetra_map_rcp(communicator, true));
+
+ compressed = false;
+ }
+
+ has_ghost = (vector->getMap()->isOneToOne() == false);
+
+# ifdef DEBUG
+ MPI_Comm comm = Utilities::Trilinos::teuchos_comm_to_mpi_comm(
+ vector->getMap()->getComm());
+ const size_type n_elements_global =
+ Utilities::MPI::sum(vector->getLocalLength(), comm);
+ Assert(has_ghost || n_elements_global == size(), ExcInternalError());
+# endif
}
void
Vector<Number>::reinit(const IndexSet ¶llel_partitioner,
const MPI_Comm communicator,
- const bool omit_zeroing_entries)
+ const bool /*omit_zeroing_entries*/)
{
- // release memory before reallocation
vector.reset();
- locally_owned_map.reset();
- locally_relevant_map.reset();
- source_stored_elements.clear();
-
- // Here we assume, that the index set is non-overlapping.
- // If the index set is overlapping the function
- // make_tpetra_map_rcp() will throw an assert.
- has_ghost = false;
- locally_owned_map =
- parallel_partitioner.make_tpetra_map_rcp(communicator, false);
-
- if (!vector->getMap()->isSameAs(*locally_owned_map))
- vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
- locally_owned_map);
- if (!omit_zeroing_entries)
- {
- vector->putScalar(0.);
- }
+ nonlocal_vector.reset();
+
+ compressed = true;
+ has_ghost = false;
+ vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
+ parallel_partitioner.make_tpetra_map_rcp(communicator, true));
}
void
Vector<Number>::reinit(const IndexSet &locally_owned_entries,
const IndexSet &ghost_entries,
- const MPI_Comm communicator)
+ const MPI_Comm communicator,
+ const bool vector_writable)
{
// release memory before reallocation
- vector.reset();
- locally_owned_map.reset();
- locally_relevant_map.reset();
+ nonlocal_vector.reset();
+
+ if (!vector_writable)
+ {
+ vector.reset();
+
+ IndexSet parallel_partitioner = locally_owned_entries;
+ parallel_partitioner.add_indices(ghost_entries);
- locally_owned_map =
- locally_owned_entries.make_tpetra_map_rcp(communicator, false);
+ vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
+ parallel_partitioner.make_tpetra_map_rcp(communicator, true));
- IndexSet parallel_partitioner = locally_owned_entries;
- parallel_partitioner.add_indices(ghost_entries);
- locally_relevant_map =
- parallel_partitioner.make_tpetra_map_rcp(communicator, true);
+ compressed = true;
+ }
+ else
+ {
+ Teuchos::RCP<MapType> map =
+ locally_owned_entries.make_tpetra_map_rcp(communicator, false);
+
+ if (!vector->getMap()->isSameAs(*map))
+ {
+ vector.reset();
+ vector = Utilities::Trilinos::internal::make_rcp<VectorType>(map);
+ }
+ else
+ vector->putScalar(0);
- vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
- locally_relevant_map);
+ IndexSet nonlocal_entries(ghost_entries);
+ nonlocal_entries.subtract_set(locally_owned_entries);
- has_ghost = false;
+ nonlocal_vector = Utilities::Trilinos::internal::make_rcp<VectorType>(
+ nonlocal_entries.make_tpetra_map_rcp(communicator, true));
+
+ compressed = false;
+ }
+
+ has_ghost = (vector->getMap()->isOneToOne() == false);
}
V.vector->getMap());
Tpetra::deep_copy(*vector, *V.vector);
+ compressed = V.compressed;
has_ghost = V.has_ghost;
- locally_owned_map = V.locally_owned_map;
- locally_relevant_map = V.locally_relevant_map;
source_stored_elements = V.source_stored_elements;
tpetra_comm_pattern = V.tpetra_comm_pattern;
}
Vector<Number> &
Vector<Number>::operator=(const Number s)
{
- Assert(s == Number(0.),
+ (void)s;
+ Assert(s == Number(0.0),
ExcMessage("Only 0 can be assigned to a vector."));
- vector->putScalar(s);
+ // As checked above, we are only allowed to use d==0.0, so pass
+ // a constant zero (instead of a run-time value 'd' that *happens* to
+ // have a zero value) to the underlying class in hopes that the compiler
+ // can optimize this somehow.
+ vector->putScalar(/*s=*/0.0);
+
+ if (nonlocal_vector.is_null())
+ nonlocal_vector->putScalar(/*s=*/0.0);
return *this;
}
ExcDimensionMismatch(this->size(), V.size()));
Assert(vector->getMap()->isSameAs(*V.trilinos_vector().getMap()),
ExcDifferentParallelPartitioning());
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
return vector->dot(V.trilinos_vector());
}
// If the element is not present on the current processor, we can't
// continue. This is the main difference to the el() function.
- if (local_index == -1)
+ if (local_index == Teuchos::OrdinalTraits<int>::invalid())
{
# if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
Assert(
{
AssertIsFinite(a);
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
# if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>(
Tpetra::Access::ReadWrite);
Vector<Number>::add(const Number a, const Vector<Number> &V)
{
AssertIsFinite(a);
+
Assert(vector->getMap()->isSameAs(*(V.trilinos_vector().getMap())),
ExcDifferentParallelPartitioning());
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
vector->update(a, V.trilinos_vector(), 1.);
}
const Number b,
const Vector<Number> &W)
{
+ AssertIsFinite(a);
+ AssertIsFinite(b);
+
Assert(vector->getMap()->isSameAs(*(V.trilinos_vector().getMap())),
ExcDifferentParallelPartitioning());
Assert(vector->getMap()->isSameAs(*(W.trilinos_vector().getMap())),
ExcDifferentParallelPartitioning());
- AssertIsFinite(a);
- AssertIsFinite(b);
+
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
vector->update(a, V.trilinos_vector(), b, W.trilinos_vector(), 1.);
}
const Number a,
const Vector<Number> &V)
{
+ AssertIsFinite(s);
+ AssertIsFinite(a);
+
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
*this *= s;
Vector<Number> tmp(V);
*(scaling_factors.trilinos_vector().getMap())),
ExcDifferentParallelPartitioning());
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
vector->elementWiseMultiply(1., *scaling_factors.vector, *vector, 0.);
}
void
Vector<Number>::equ(const Number a, const Vector<Number> &V)
{
+ AssertIsFinite(a);
+
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
// If we don't have the same map, copy.
if (vector->getMap()->isSameAs(*V.trilinos_vector().getMap()) == false)
this->sadd(0., a, V);
Number
Vector<Number>::mean_value() const
{
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
return vector->meanValue();
}
typename Vector<Number>::real_type
Vector<Number>::l1_norm() const
{
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
return vector->norm1();
}
typename Vector<Number>::real_type
Vector<Number>::l2_norm() const
{
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
return vector->norm2();
}
typename Vector<Number>::real_type
Vector<Number>::linfty_norm() const
{
+ Assert(!has_ghost_elements(), ExcGhostsPresent());
+
return vector->normInf();
}
const Vector<Number> &V,
const Vector<Number> &W)
{
+ AssertIsFinite(a);
+
this->add(a, V);
return *this * W;
Vector<Number>::compress(const VectorOperation::values operation)
{
Assert(has_ghost == false,
- ExcMessage("Calling compress() is only useful if a vector "
- "has been written into, but this is a vector with ghost "
- "elements and consequently is read-only. It does "
- "not make sense to call compress() for such "
- "vectors."));
-
- if (!vector->getMap()->isOneToOne())
+ ExcMessage(
+ "Calling compress() is only useful if a vector "
+ "has been written into, but this is a vector with ghost "
+ "elements and consequently is read-only. It does "
+ "not make sense to call compress() for such "
+ "vectors."));
+
+ if (!compressed)
{
Tpetra::CombineMode tpetra_operation = Tpetra::ZERO;
if (operation == VectorOperation::insert)
else
Assert(false, ExcNotImplemented());
- if (locally_relevant_map.is_null())
- locally_relevant_map =
- Teuchos::rcp_const_cast<MapType>(vector->getMap());
-
- Assert(!locally_relevant_map.is_null(), ExcMissingIndexSet());
- Assert(!locally_owned_map.is_null(), ExcMissingIndexSet());
-
- VectorType dummy = VectorType(locally_owned_map.getConst());
-
Teuchos::RCP<const ExportType> exporter =
- Tpetra::createExport(locally_relevant_map.getConst(),
- locally_owned_map.getConst());
- dummy.doExport(*vector, *exporter, tpetra_operation);
+ Tpetra::createExport(nonlocal_vector->getMap(), vector->getMap());
+ vector->doExport(*nonlocal_vector, *exporter, tpetra_operation);
- *vector = dummy;
+ compressed = true;
}
}
projects / dealii.git / commitdiff