Eigen subvector support, unit tests

include/numerics/eigen_sparse_vector.h

@@ -231,10 +231,20 @@ class EigenSparseVector final : public NumericVector<T>
   virtual void pointwise_divide (const NumericVector<T> & vec1,
                                  const NumericVector<T> & vec2) override;
 
+  virtual void create_subvector(NumericVector<T> & subvector,
+                                const std::vector<numeric_index_type> & rows,
+                                bool supplying_global_rows = true) const override;
+
   virtual void swap (NumericVector<T> & v) override;
 
   virtual std::size_t max_allowed_id() const override;
 
+  virtual std::unique_ptr<NumericVector<T>>
+  get_subvector(const std::vector<numeric_index_type> & rows) override;
+
+  virtual void restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
+                                 const std::vector<numeric_index_type> & rows) override;
+
   /**
    * References to the underlying Eigen data types.
    *

include/numerics/numeric_vector.h

@@ -752,39 +752,48 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
 
   /**
    * Fills in \p subvector from this vector using the indices in \p rows.
-   * Similar to the \p create_submatrix() routine for the SparseMatrix class, it
-   * is currently only implemented for PetscVectors. The boolean parameter
-   * communicates whether the supplied vector of rows corresponds to all the
-   * rows that should be used in the subvector's index set, e.g. whether the
-   * rows correspond to the global collective. If the rows supplied are only the
-   * local indices, then the boolean parameter should be set to false
-   */
-  virtual void create_subvector(NumericVector<T> & ,
-                                const std::vector<numeric_index_type> &,
-                                bool = true) const
+   *
+   * \p supplying_global_rows communicates whether
+   * the supplied vector of rows corresponds to all the rows that
+   * should be used in the subvector's index set, e.g. whether the
+   * rows correspond to the global collective. If the rows supplied
+   * are only the local indices, then \p supplying_global_rows should
+   * be set to false
+   *
+   * This is currently only implemented for \p PetscVector and \p
+   * EigenSparseVector.
+   */
+  virtual void create_subvector(NumericVector<T> & /* subvector */,
+                                const std::vector<numeric_index_type> & /* rows */,
+                                bool /* supplying_global_rows */ = true) const
   {
     libmesh_not_implemented();
   }
 
   /**
    * Creates a view into this vector using the indices in \p
-   * rows.  Similar to the \p create_submatrix() routine for the
-   * SparseMatrix class, it is currently only implemented for
-   * PetscVectors.
+   * rows. Calls to this API should always be paired with a call to \p restore_subvector,
+   * else any changes made in the subvector will not be reflected in (\p this) original vector.
+   *
+   * This is currently only implemented for \p PetscVector and \p
+   * EigenSparseVector.
    */
   virtual std::unique_ptr<NumericVector<T>>
-  get_subvector(const std::vector<numeric_index_type> &)
+  get_subvector(const std::vector<numeric_index_type> & /* rows */)
   {
     libmesh_not_implemented();
   }
 
   /**
    * Restores a view into this vector using the indices in \p
-   * rows.  This  is currently only implemented for
-   * PetscVectors.
+   * rows.  The \p subvector should have been acquired from \p
+   * get_subvector() with the same rows.
+   *
+   * This is currently only implemented for \p PetscVector and \p
+   * EigenSparseVector.
    */
-  virtual void restore_subvector(std::unique_ptr<NumericVector<T>>,
-                                 const std::vector<numeric_index_type> &)
+  virtual void restore_subvector(std::unique_ptr<NumericVector<T>> /* subvector */,
+                                 const std::vector<numeric_index_type> & /* rows */)
   {
     libmesh_not_implemented();
   }

src/numerics/eigen_sparse_vector.C

@@ -426,6 +426,57 @@ void EigenSparseVector<T>::pointwise_divide (const NumericVector<T> & /*vec1*/,
 }
 
 
+
+template <typename T>
+void EigenSparseVector<T>::create_subvector(NumericVector<T> & subvector,
+                                            const std::vector<numeric_index_type> & rows,
+                                            const bool /* supplying_global_rows */) const
+{
+  // Make sure the passed in subvector is really an EigenSparseVector
+  EigenSparseVector<T> * eigen_subvector = cast_ptr<EigenSparseVector<T> *>(&subvector);
+
+  // If the eigen_subvector is already initialized, we assume that the
+  // user has already allocated the *correct* amount of space for it.
+  // If not, we do so.
+  if (!eigen_subvector->initialized())
+    eigen_subvector->init(rows.size());
+
+  eigen_subvector->vec() = this->vec()(rows);
+
+  eigen_subvector->_is_closed = true;
+}
+
+
+
+template <typename T>
+std::unique_ptr<NumericVector<T>>
+EigenSparseVector<T>::get_subvector(const std::vector<numeric_index_type> & rows)
+{
+  auto returnval = std::make_unique<EigenSparseVector<T>>(this->comm(), rows.size());
+
+  // We're just going to make a copy here, since we'll be calling a
+  // restore later anyway.  Someone with more familiarity with Eigen
+  // can see if there's a way to improve this later.
+  this->create_subvector(*returnval, rows);
+
+  this->_is_closed = false;
+
+  return returnval;
+}
+
+template <typename T>
+void
+EigenSparseVector<T>::restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
+                                        const std::vector<numeric_index_type> & rows)
+{
+  auto * const eigen_subvector = cast_ptr<EigenSparseVector<T> *>(subvector.get());
+
+  this->vec()(rows) = eigen_subvector->vec();
+  this->_is_closed = true;
+}
+
+
+
 template <typename T>
 Real EigenSparseVector<T>::max() const
 {

tests/numerics/eigen_sparse_vector_test.C

@@ -40,6 +40,8 @@ public:
   CPPUNIT_TEST_SUITE( EigenSparseVectorTest );
 
   NUMERICVECTORTEST
+  CPPUNIT_TEST( testSubvectors );
+  CPPUNIT_TEST( testSubvectorsBase );
 
   CPPUNIT_TEST_SUITE_END();
 };

tests/numerics/numeric_vector_test.h

@@ -7,6 +7,7 @@
 // libMesh includes
 #include <libmesh/parallel.h>
 #include <libmesh/fuzzy_equals.h>
+#include <libmesh/int_range.h>
 
 #include "libmesh_cppunit.h"
 
@@ -26,6 +27,10 @@
   CPPUNIT_TEST( testWriteAndRead );             \
   CPPUNIT_TEST( testWriteAndReadBase );
 
+// Tests to add manually, on subclasses supporting them
+//  CPPUNIT_TEST( testSubvectors );
+//  CPPUNIT_TEST( testSubvectorsBase );
+
 
 template <class DerivedClass>
 class NumericVectorTest : public CppUnit::TestCase {
@@ -150,6 +155,98 @@ class NumericVectorTest : public CppUnit::TestCase {
        libMesh::TOLERANCE*libMesh::TOLERANCE);
   }
 
+
+  template <class Base, class Derived>
+  void Subvectors()
+  {
+    auto v_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
+    Base & v = *v_ptr;
+
+    const libMesh::dof_id_type
+      first = v.first_local_index(),
+      last  = v.last_local_index();
+
+    const unsigned int sub_local_size = block_size/2;
+    const unsigned int sub_global_size = sub_local_size * my_comm->size();
+
+    for (auto supplying_global_rows : {false, true})
+      for (auto get_not_create : {false, true})
+        {
+          // We don't have such an option?
+          if (get_not_create && supplying_global_rows)
+            continue;
+
+          for (libMesh::dof_id_type n=first; n != last; n++)
+            v.set (n, static_cast<libMesh::Number>(n+1));
+          v.close();
+
+          std::unique_ptr<Base> s_ptr =
+            std::make_unique<Derived>(*my_comm, sub_global_size,
+                                      sub_local_size);
+
+          const libMesh::dof_id_type
+            sub_first = s_ptr->first_local_index(),
+            sub_last  = s_ptr->last_local_index();
+
+          CPPUNIT_ASSERT_EQUAL
+            (sub_first, libMesh::dof_id_type(sub_local_size *
+                                             my_comm->rank()));
+
+          CPPUNIT_ASSERT_EQUAL
+            (sub_last, libMesh::dof_id_type(sub_local_size *
+                                            (my_comm->rank()+1)));
+
+          std::vector<libMesh::dof_id_type> rows;
+          for (auto i : libMesh::make_range(sub_local_size))
+            rows.push_back(first+2*i);
+
+          if (supplying_global_rows)
+            my_comm->allgather(rows);
+
+          if (get_not_create)
+            {
+              libMesh::NumericVector<libMesh::Number> * subvec_ptr =
+                v.get_subvector(rows).release();
+              s_ptr.reset(dynamic_cast<Base *>(subvec_ptr));
+            }
+          else
+            v.create_subvector(*s_ptr, rows, supplying_global_rows);
+
+          Base & s = *s_ptr;
+
+          for (auto i : libMesh::make_range(sub_local_size))
+            LIBMESH_ASSERT_NUMBERS_EQUAL
+              (s(sub_first+i), libMesh::Real(first+2*i+1),
+               libMesh::TOLERANCE*libMesh::TOLERANCE);
+
+          for (libMesh::dof_id_type n=sub_first; n != sub_last; n++)
+            s.set(n, (my_comm->rank()+1) * n);
+          s.close();
+
+          // If we "get" a subvector then we're set up to "restore" it
+          // afterwards; if not then we're done.
+          if (!get_not_create)
+            continue;
+
+          v.restore_subvector(std::move(s_ptr), rows);
+
+          for (libMesh::dof_id_type n=first; n != last; n++)
+            {
+              const libMesh::dof_id_type offset = n - first;
+              if (offset < 2*sub_local_size && !(offset%2))
+                LIBMESH_ASSERT_NUMBERS_EQUAL
+                  (v(n),
+                   libMesh::Real((my_comm->rank()+1) *
+                                 ((offset/2) + sub_first)),
+                   libMesh::TOLERANCE*libMesh::TOLERANCE);
+              else
+                LIBMESH_ASSERT_NUMBERS_EQUAL
+                  (v(n), libMesh::Real(n+1),
+                   libMesh::TOLERANCE*libMesh::TOLERANCE);
+            }
+        }
+  }
+
   template <class Base, class Derived>
   void WriteAndRead()
   {
@@ -406,6 +503,20 @@ class NumericVectorTest : public CppUnit::TestCase {
     Operations<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
   }
 
+  void testSubvectors()
+  {
+    LOG_UNIT_TEST;
+
+    Subvectors<DerivedClass,DerivedClass >();
+  }
+
+  void testSubvectorsBase()
+  {
+    LOG_UNIT_TEST;
+
+    Subvectors<libMesh::NumericVector<libMesh::Number>,DerivedClass>();
+  }
+
   void testWriteAndRead()
   {
     LOG_UNIT_TEST;

tests/numerics/petsc_vector_test.C

@@ -22,8 +22,9 @@ public:
   NUMERICVECTORTEST
 
   CPPUNIT_TEST( testGetArray );
-
   CPPUNIT_TEST( testPetscOperations );
+  CPPUNIT_TEST( testSubvectors );
+  CPPUNIT_TEST( testSubvectorsBase );
 
   CPPUNIT_TEST_SUITE_END();