diff --git a/include/numerics/eigen_sparse_vector.h b/include/numerics/eigen_sparse_vector.h index 7415b6b06c..ec1fc6b183 100644 --- a/include/numerics/eigen_sparse_vector.h +++ b/include/numerics/eigen_sparse_vector.h @@ -231,10 +231,20 @@ class EigenSparseVector final : public NumericVector virtual void pointwise_divide (const NumericVector & vec1, const NumericVector & vec2) override; + virtual void create_subvector(NumericVector & subvector, + const std::vector & rows, + bool supplying_global_rows = true) const override; + virtual void swap (NumericVector & v) override; virtual std::size_t max_allowed_id() const override; + virtual std::unique_ptr> + get_subvector(const std::vector & rows) override; + + virtual void restore_subvector(std::unique_ptr> subvector, + const std::vector & rows) override; + /** * References to the underlying Eigen data types. * diff --git a/include/numerics/numeric_vector.h b/include/numerics/numeric_vector.h index 9db70c6649..ede24c0159 100644 --- a/include/numerics/numeric_vector.h +++ b/include/numerics/numeric_vector.h @@ -752,39 +752,48 @@ class NumericVector : public ReferenceCountedObject>, /** * 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 & , - const std::vector &, - 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 & /* subvector */, + const std::vector & /* 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> - get_subvector(const std::vector &) + get_subvector(const std::vector & /* 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>, - const std::vector &) + virtual void restore_subvector(std::unique_ptr> /* subvector */, + const std::vector & /* rows */) { libmesh_not_implemented(); } diff --git a/src/numerics/eigen_sparse_vector.C b/src/numerics/eigen_sparse_vector.C index 8743df6f4f..54a40beaab 100644 --- a/src/numerics/eigen_sparse_vector.C +++ b/src/numerics/eigen_sparse_vector.C @@ -426,6 +426,57 @@ void EigenSparseVector::pointwise_divide (const NumericVector & /*vec1*/, } + +template +void EigenSparseVector::create_subvector(NumericVector & subvector, + const std::vector & rows, + const bool /* supplying_global_rows */) const +{ + // Make sure the passed in subvector is really an EigenSparseVector + EigenSparseVector * eigen_subvector = cast_ptr *>(&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 +std::unique_ptr> +EigenSparseVector::get_subvector(const std::vector & rows) +{ + auto returnval = std::make_unique>(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 +void +EigenSparseVector::restore_subvector(std::unique_ptr> subvector, + const std::vector & rows) +{ + auto * const eigen_subvector = cast_ptr *>(subvector.get()); + + this->vec()(rows) = eigen_subvector->vec(); + this->_is_closed = true; +} + + + template Real EigenSparseVector::max() const { diff --git a/tests/numerics/eigen_sparse_vector_test.C b/tests/numerics/eigen_sparse_vector_test.C index 6648617dec..ebe4e12fab 100644 --- a/tests/numerics/eigen_sparse_vector_test.C +++ b/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(); }; diff --git a/tests/numerics/numeric_vector_test.h b/tests/numerics/numeric_vector_test.h index fd5239b3fa..7e5d164b29 100644 --- a/tests/numerics/numeric_vector_test.h +++ b/tests/numerics/numeric_vector_test.h @@ -7,6 +7,7 @@ // libMesh includes #include #include +#include #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 NumericVectorTest : public CppUnit::TestCase { @@ -150,6 +155,98 @@ class NumericVectorTest : public CppUnit::TestCase { libMesh::TOLERANCE*libMesh::TOLERANCE); } + + template + void Subvectors() + { + auto v_ptr = std::make_unique(*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(n+1)); + v.close(); + + std::unique_ptr s_ptr = + std::make_unique(*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 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 * subvec_ptr = + v.get_subvector(rows).release(); + s_ptr.reset(dynamic_cast(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 void WriteAndRead() { @@ -406,6 +503,20 @@ class NumericVectorTest : public CppUnit::TestCase { Operations,DerivedClass>(); } + void testSubvectors() + { + LOG_UNIT_TEST; + + Subvectors(); + } + + void testSubvectorsBase() + { + LOG_UNIT_TEST; + + Subvectors,DerivedClass>(); + } + void testWriteAndRead() { LOG_UNIT_TEST; diff --git a/tests/numerics/petsc_vector_test.C b/tests/numerics/petsc_vector_test.C index 04a191db3f..33fe70aec8 100644 --- a/tests/numerics/petsc_vector_test.C +++ b/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();