From: Sebastian Kinnewig <sebastian@kinnewig.org>
Date: Wed, 20 Mar 2024 18:54:30 +0000 (+0100)
Subject: Add local_range() and in_local_range() to LA::TpetraWrappers::Vector
X-Git-Tag: v9.6.0-rc1~456^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16781%2Fhead;p=dealii.git

Add local_range() and in_local_range() to LA::TpetraWrappers::Vector
---

diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h
index 9016045001..2990a86469 100644
--- a/include/deal.II/lac/trilinos_tpetra_vector.h
+++ b/include/deal.II/lac/trilinos_tpetra_vector.h
@@ -786,6 +786,41 @@ namespace LinearAlgebra
       size_type
       locally_owned_size() const;
 
+      /**
+       * Return a pair of indices indicating which elements of this vector are
+       * stored locally. The first number is the index of the first element
+       * stored, the second the index of the one past the last one that is
+       * stored locally. If this is a sequential vector, then the result will be
+       * the pair <code>(0,N)</code>, otherwise it will be a pair
+       * <code>(i,i+n)</code>, where <code>n</code> is the number of elements
+       * stored on this processor and and <code>i</code> is the first element of
+       * the vector stored on this processor, corresponding to the half open
+       * interval $[i,i+n)$
+       *
+       * @note The description above is true most of the time, but not always.
+       * In particular, Trilinos vectors need not store contiguous ranges of
+       * elements such as $[i,i+n)$. Rather, it can store vectors where the
+       * elements are distributed in an arbitrary way across all processors and
+       * each processor simply stores a particular subset, not necessarily
+       * contiguous. In this case, this function clearly makes no sense since it
+       * could, at best, return a range that includes all elements that are
+       * stored locally. Thus, the function only succeeds if the locally stored
+       * range is indeed contiguous. It will trigger an assertion if the local
+       * portion of the vector is not contiguous.
+       */
+      std::pair<size_type, size_type>
+      local_range() const;
+
+      /**
+       * Return whether @p index is in the local range or not, see also
+       * local_range().
+       *
+       * @note The same limitation for the applicability of this function
+       * applies as listed in the documentation of local_range().
+       */
+      bool
+      in_local_range(const size_type index) const;
+
       /**
        * Return the state of the vector, i.e., whether compress() needs to be
        * called after an operation requiring data exchange. A call to compress()
diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
index 2a2352c8be..2da046ba3d 100644
--- a/include/deal.II/lac/trilinos_tpetra_vector.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
@@ -1024,6 +1024,41 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    std::pair<typename Vector<Number, MemorySpace>::size_type,
+              typename Vector<Number, MemorySpace>::size_type>
+    Vector<Number, MemorySpace>::local_range() const
+    {
+      const size_type begin = vector->getMap()->getMinGlobalIndex();
+      const size_type end   = vector->getMap()->getMaxGlobalIndex() + 1;
+
+      Assert(
+#  if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0)
+        end - begin == vector->getMap()->getLocalNumElements(),
+#  else
+        end - begin == vector->getMap()->getNodeNumElements(),
+#  endif
+        ExcMessage(
+          "This function only makes sense if the elements that this "
+          "vector stores on the current processor form a contiguous range. "
+          "This does not appear to be the case for the current vector."));
+
+      return std::make_pair(begin, end);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    bool
+    Vector<Number, MemorySpace>::in_local_range(const size_type index) const
+    {
+      std::pair<size_type, size_type> range = local_range();
+
+      return ((index >= range.first) && (index < range.second));
+    }
+
+
+
     template <typename Number, typename MemorySpace>
     MPI_Comm
     Vector<Number, MemorySpace>::get_mpi_communicator() const