]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement EpetraWrappers::internal::VectorReference
authorPasquale Africa <pafrica@sissa.it>
Fri, 12 Apr 2024 14:02:16 +0000 (14:02 +0000)
committerPasquale Africa <pafrica@sissa.it>
Tue, 16 Apr 2024 07:33:20 +0000 (09:33 +0200)
include/deal.II/lac/trilinos_epetra_vector.h
source/lac/trilinos_epetra_vector.cc
tests/trilinos/epetra_vector_reference_01.cc [new file with mode: 0644]
tests/trilinos/epetra_vector_reference_01.output [new file with mode: 0644]

index 53fa509a3cb38a06e4f107d35e83563300b320b6..ad2f901d9266ff45c3b54dc0a0ffe51e8556811c 100644 (file)
@@ -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 <tt>v(i)=w(i)</tt>. 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 <tt>s</tt>.
+         */
+        const VectorReference &
+        operator=(const value_type &s) const;
+
+        /**
+         * Add <tt>s</tt> to the referenced element of the vector->
+         */
+        const VectorReference &
+        operator+=(const value_type &s) const;
+
+        /**
+         * Subtract <tt>s</tt> from the referenced element of the vector->
+         */
+        const VectorReference &
+        operator-=(const value_type &s) const;
+
+        /**
+         * Multiply the referenced element of the vector by <tt>s</tt>.
+         */
+        const VectorReference &
+        operator*=(const value_type &s) const;
+
+        /**
+         * Divide the referenced element of the vector by <tt>s</tt>.
+         */
+        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<double>, public Subscriptor
+    class Vector : public ReadVector<internal::value_type>, 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<const CommunicationPattern> 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<value_type>(r);
+
+        return *this;
+      }
+
+
+
+      inline VectorReference &
+      VectorReference::operator=(const VectorReference &r)
+      {
+        // as above
+        *this = static_cast<value_type>(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<value_type>(*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<value_type>(*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<value_type>(*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<value_type>(*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
index 09766ec35c6aa51b4673767eef00a88237b8f459..54556de525656cc40010e84007bcc55cc5bd2c0f 100644 (file)
@@ -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<TrilinosWrappers::types::int_type>(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 (file)
index 0000000..c0905d2
--- /dev/null
@@ -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 <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_epetra_vector.h>
+
+#include <iostream>
+
+#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<TrilinosScalar>(a)
+          << std::endl; // should point to v(0)
+  deallog << static_cast<TrilinosScalar>(b)
+          << std::endl; // should point to v(0)
+  deallog << static_cast<TrilinosScalar>(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 (file)
index 0000000..440392e
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::0.00000
+DEAL::0.00000
+DEAL::1.00000

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.