From: Pasquale Africa Date: Fri, 12 Apr 2024 14:02:16 +0000 (+0000) Subject: Implement EpetraWrappers::internal::VectorReference X-Git-Tag: v9.6.0-rc1~372^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3b252666a1189251a35f1281be1e452fd0f810e4;p=dealii.git Implement EpetraWrappers::internal::VectorReference --- diff --git a/include/deal.II/lac/trilinos_epetra_vector.h b/include/deal.II/lac/trilinos_epetra_vector.h index 53fa509a3c..ad2f901d92 100644 --- a/include/deal.II/lac/trilinos_epetra_vector.h +++ b/include/deal.II/lac/trilinos_epetra_vector.h @@ -50,6 +50,166 @@ namespace LinearAlgebra */ namespace EpetraWrappers { + /** + * @cond internal + */ + + /** + * A namespace for internal implementation details of the TrilinosWrapper + * members. + * + * @ingroup TrilinosWrappers + */ + namespace internal + { + /** + * Declare type for container size. + */ + using size_type = dealii::types::global_dof_index; + + /** + * Declare type for container value type. + */ + using value_type = double; + + /** + * This class implements a wrapper for accessing the Trilinos Epetra + * vector in the same way as we access deal.II objects: it is initialized + * with a vector and an element within it, and has a conversion operator + * to extract the scalar value of this element. It also has a variety of + * assignment operator for writing to this one element. + * + * @ingroup TrilinosWrappers + */ + class VectorReference + { + private: + /** + * Constructor. It is made private so as to only allow the actual vector + * class to create it. + */ + VectorReference(Vector &vector, const size_type index); + + public: + /** + * Copy constructor. + */ + VectorReference(const VectorReference &) = default; + + /** + * This looks like a copy operator, but does something different than + * usual. In particular, it does not copy the member variables of this + * reference. Rather, it handles the situation where we have two vectors + * @p v and @p w, and assign elements like in v(i)=w(i). Here, + * both left and right hand side of the assignment have data type + * VectorReference, but what we really mean is to assign the vector + * elements represented by the two references. This operator implements + * this operation. Note also that this allows us to make the assignment + * operator const. + */ + const VectorReference & + operator=(const VectorReference &r) const; + + /** + * Same as above but for non-const reference objects. + */ + VectorReference & + operator=(const VectorReference &r); + + /** + * Set the referenced element of the vector to s. + */ + const VectorReference & + operator=(const value_type &s) const; + + /** + * Add s to the referenced element of the vector-> + */ + const VectorReference & + operator+=(const value_type &s) const; + + /** + * Subtract s from the referenced element of the vector-> + */ + const VectorReference & + operator-=(const value_type &s) const; + + /** + * Multiply the referenced element of the vector by s. + */ + const VectorReference & + operator*=(const value_type &s) const; + + /** + * Divide the referenced element of the vector by s. + */ + const VectorReference & + operator/=(const value_type &s) const; + + /** + * Convert the reference to an actual value, i.e. return the value of + * the referenced element of the vector. + */ + operator value_type() const; + + /** + * Exception + */ + DeclException1(ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occurred while calling a Trilinos function"); + + /* + * Access to a an element that is not (locally-)owned. + * + * @ingroup Exceptions + */ + DeclException4( + ExcAccessToNonLocalElement, + size_type, + size_type, + size_type, + size_type, + << "You are trying to access element " << arg1 + << " of a distributed vector, but this element is not stored " + << "on the current processor. Note: There are " << arg2 + << " elements stored " + << "on the current processor from within the range [" << arg3 << ',' + << arg4 << "] but Trilinos vectors need not store contiguous " + << "ranges on each processor, and not every element in " + << "this range may in fact be stored locally." + << "\n\n" + << "A common source for this kind of problem is that you " + << "are passing a 'fully distributed' vector into a function " + << "that needs read access to vector elements that correspond " + << "to degrees of freedom on ghost cells (or at least to " + << "'locally active' degrees of freedom that are not also " + << "'locally owned'). You need to pass a vector that has these " + << "elements as ghost entries."); + + private: + /** + * Point to the vector we are referencing. + */ + Vector &vector; + + /** + * Index of the referenced element of the vector. + */ + const size_type index; + + // Make the vector class a friend, so that it can create objects of the + // present type. + friend class ::dealii::LinearAlgebra::EpetraWrappers::Vector; + }; // class VectorReference + + } // namespace internal + /** + * @endcond + */ + + /** * This class implements a wrapper to the Trilinos distributed vector * class Epetra_FEVector. This class requires Trilinos to be compiled @@ -58,12 +218,18 @@ namespace LinearAlgebra * @ingroup TrilinosWrappers * @ingroup Vectors */ - class Vector : public ReadVector, public Subscriptor + class Vector : public ReadVector, public Subscriptor { public: - using value_type = double; - using size_type = types::global_dof_index; - using real_type = double; + using value_type = internal::value_type; + using size_type = types::global_dof_index; + using reference = internal::VectorReference; + using const_reference = const internal::VectorReference; + + /** + * @name 1: Basic Object-handling + */ + /** @{ */ /** * Constructor. Create a vector of dimension zero. @@ -155,6 +321,57 @@ namespace LinearAlgebra import_elements(V, operation, communication_pattern); } + /** @} */ + + + /** + * @name 2: Data-Access + */ + /** @{ */ + + /** + * Provide access to a given element, both read and write. + * + * When using a vector distributed with MPI, this operation only makes + * sense for elements that are actually present on the calling processor. + * Otherwise, an exception is thrown. + */ + reference + operator()(const size_type index); + + /** + * Provide read-only access to an element. + * + * When using a vector distributed with MPI, this operation only makes + * sense for elements that are actually present on the calling processor. + * Otherwise, an exception is thrown. + */ + value_type + operator()(const size_type index) const; + + /** + * Provide access to a given element, both read and write. + * + * Exactly the same as operator(). + */ + reference + operator[](const size_type index); + + /** + * Provide read-only access to an element. + * + * Exactly the same as operator(). + */ + value_type + operator[](const size_type index) const; + + /** @} */ + + /** + * @name 3: Modification of vectors + */ + /** @{ */ + /** * Multiply the entire vector by a fixed factor. */ @@ -234,6 +451,14 @@ namespace LinearAlgebra bool all_zero() const; + /** @} */ + + + /** + * @name 4: Scalar products, norms and related operations + */ + /** @{ */ + /** * Return the mean value of the element of this vector. */ @@ -285,6 +510,14 @@ namespace LinearAlgebra */ double add_and_dot(const double a, const Vector &V, const Vector &W); + + /** @} */ + + /** + * @name 5: Scalar products, norms and related operations + */ + /** @{ */ + /** * This function always returns false and is present only for backward * compatibility. @@ -326,6 +559,14 @@ namespace LinearAlgebra ::dealii::IndexSet locally_owned_elements() const; + /** @} */ + + + /** + * @name 6: Mixed stuff + */ + /** @{ */ + /** * Compress the underlying representation of the Trilinos object, i.e. * flush the buffers of the vector object if it has any. This function is @@ -368,6 +609,8 @@ namespace LinearAlgebra std::size_t memory_consumption() const; + /** @} */ + /** * The vectors have different partitioning, i.e. their IndexSet objects * don't represent the same indices. @@ -416,8 +659,123 @@ namespace LinearAlgebra * source_stored_elements IndexSet and the current vector. */ std::shared_ptr epetra_comm_pattern; + + // Make the reference class a friend. + friend class internal::VectorReference; }; +# ifndef DOXYGEN + + // VectorReference + namespace internal + { + inline VectorReference::VectorReference(Vector &vector, + const size_type index) + : vector(vector) + , index(index) + {} + + + + inline const VectorReference & + VectorReference::operator=(const VectorReference &r) const + { + // as explained in the class + // documentation, this is not the copy + // operator. so simply pass on to the + // "correct" assignment operator + *this = static_cast(r); + + return *this; + } + + + + inline VectorReference & + VectorReference::operator=(const VectorReference &r) + { + // as above + *this = static_cast(r); + + return *this; + } + + + + inline const VectorReference & + VectorReference::operator=(const value_type &value) const + { + (*vector.vector)[0][index] = value; + + return *this; + } + + + + inline const VectorReference & + VectorReference::operator+=(const value_type &value) const + { + value_type new_value = static_cast(*this) + value; + (*vector.vector)[0][index] = new_value; + + return *this; + } + + + + inline const VectorReference & + VectorReference::operator-=(const value_type &value) const + { + value_type new_value = static_cast(*this) - value; + (*vector.vector)[0][index] = new_value; + + return *this; + } + + + + inline const VectorReference & + VectorReference::operator*=(const value_type &value) const + { + value_type new_value = static_cast(*this) * value; + (*vector.vector)[0][index] = new_value; + + return *this; + } + + + + inline const VectorReference & + VectorReference::operator/=(const value_type &value) const + { + value_type new_value = static_cast(*this) / value; + (*vector.vector)[0][index] = new_value; + + return *this; + } + } // namespace internal + +# endif /* DOXYGEN */ + + + inline internal::VectorReference + Vector::operator()(const size_type index) + { + return internal::VectorReference(*this, index); + } + + inline internal::VectorReference + Vector::operator[](const size_type index) + { + return operator()(index); + } + + inline Vector::value_type + Vector::operator[](const size_type index) const + { + return operator()(index); + } + inline bool Vector::has_ghost_elements() const diff --git a/source/lac/trilinos_epetra_vector.cc b/source/lac/trilinos_epetra_vector.cc index 09766ec35c..54556de525 100644 --- a/source/lac/trilinos_epetra_vector.cc +++ b/source/lac/trilinos_epetra_vector.cc @@ -37,6 +37,33 @@ namespace LinearAlgebra { namespace EpetraWrappers { +# ifndef DOXYGEN + namespace internal + { + VectorReference::operator value_type() const + { + AssertIndexRange(index, vector.size()); + + // Trilinos allows for vectors to be referenced by the [] or () + // operators but only () checks index bounds. We check these bounds by + // ourselves, so we can use []. Note that we can only get local values. + + const TrilinosWrappers::types::int_type local_index = + vector.vector->Map().LID( + static_cast(index)); + + Assert(local_index >= 0, + ExcAccessToNonLocalElement(index, + vector.vector->Map().NumMyElements(), + vector.vector->Map().MinMyGID(), + vector.vector->Map().MaxMyGID())); + + return (*(vector.vector))[0][local_index]; + } + } // namespace internal +# endif + + // Check that the class we declare here satisfies the // vector-space-vector concept. If we catch it here, // any mistake in the vector class declaration would diff --git a/tests/trilinos/epetra_vector_reference_01.cc b/tests/trilinos/epetra_vector_reference_01.cc new file mode 100644 index 0000000000..c0905d2e0f --- /dev/null +++ b/tests/trilinos/epetra_vector_reference_01.cc @@ -0,0 +1,73 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2017 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// LinearAlgebra::EpetraWrappers::internal::VectorReference had its non-const +// copy operator return a const reference. That was non-intuitive but, because +// of the particular semantics of copying vector reference objects, made no +// difference. Either way, we now check these semantics. + +#include + +#include + +#include + +#include "../tests.h" + + +void +test() +{ + LinearAlgebra::EpetraWrappers::Vector v; + v.reinit(complete_index_set(3), MPI_COMM_WORLD); + v(0) = 0; + v(1) = 1; + v(2) = 2; + + LinearAlgebra::EpetraWrappers::internal::VectorReference a(v(0)); + LinearAlgebra::EpetraWrappers::internal::VectorReference b(v(1)); + LinearAlgebra::EpetraWrappers::internal::VectorReference c(v(2)); + + // Copy the VectorReference objects. Note that operator= for these + // objects does not copy the *content* of the reference object, but + // rather assigns the value of the vector element referenced on the + // right to the vector element referenced on the left. + // So the applied operations result in the following actions: + // 1. We first assign c to point to where a points (c points to entry 0) + // 2. We next assign c to point to where b points (c points to entry 1) + // 3. Lastly, we assign b to point to where a points (b points to entry 0) + (c = a) = b; + b = a; + + deallog << static_cast(a) + << std::endl; // should point to v(0) + deallog << static_cast(b) + << std::endl; // should point to v(0) + deallog << static_cast(c) + << std::endl; // should point to v(1) +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + initlog(); + + test(); +} diff --git a/tests/trilinos/epetra_vector_reference_01.output b/tests/trilinos/epetra_vector_reference_01.output new file mode 100644 index 0000000000..440392e8db --- /dev/null +++ b/tests/trilinos/epetra_vector_reference_01.output @@ -0,0 +1,4 @@ + +DEAL::0.00000 +DEAL::0.00000 +DEAL::1.00000