From eaf7fafc83f04bfc791a93e0d71377cc80c81190 Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Tue, 12 Jan 2021 10:47:58 +0100 Subject: [PATCH] make_nvector_view and cleanup --- include/deal.II/sundials/n_vector.h | 118 ++++++++++++++---- include/deal.II/sundials/n_vector.templates.h | 75 +++++------ source/sundials/n_vector.inst.in | 18 ++- tests/sundials/n_vector.cc | 96 +++++++------- 4 files changed, 193 insertions(+), 114 deletions(-) diff --git a/include/deal.II/sundials/n_vector.h b/include/deal.II/sundials/n_vector.h index d583924824..0b3e6ae34d 100644 --- a/include/deal.II/sundials/n_vector.h +++ b/include/deal.II/sundials/n_vector.h @@ -25,10 +25,41 @@ DEAL_II_NAMESPACE_OPEN +# ifndef DOXYGEN namespace SUNDIALS { namespace internal { + template + class NVectorView; + } +} // namespace SUNDIALS +# endif + +namespace SUNDIALS +{ + namespace internal + { + /** + * Create a NVectorView of the given @p vector. + * + * This call is intended to be used as + * + * @code + * auto view = make_nvector_view(vector); + * @endcode + * + * @tparam VectorType Type of the viewed vector. This parameter can be + * deduced automatically and will respect a potential const-qualifier. + * @param vector The vector to view. + * @return A NVectorView of the @p vector. + * + * @related NVectorView + */ + template + NVectorView + make_nvector_view(VectorType &vector); + /** * Retrieve the underlying vector attached to N_Vector @p v. This call will * only succeed if the underlying vector is not const. Use @@ -37,9 +68,9 @@ namespace SUNDIALS * @note Users must ensure that they ask for the correct VectorType when * calling this function and there are no type-safety checks in place. * - * @tparam VectorType type of the vector that is stored in @p v - * @param v vector to unwrap - * @return the vector that is stored inside @p v + * @tparam VectorType Type of the vector that is stored in @p v + * @param v Vector to unwrap + * @return The vector that is stored inside @p v */ template VectorType * @@ -52,56 +83,94 @@ namespace SUNDIALS * @note Users must ensure that they ask for the correct VectorType when * calling this function and there are no type-safety checks in place. * - * @tparam VectorType type of the vector that is stored in @p v - * @param v vector to unwrap - * @return the vector that is stored inside @p v + * @tparam VectorType Type of the vector that is stored in @p v + * @param v Vector to unwrap + * @return The vector that is stored inside @p v */ template const VectorType * unwrap_nvector_const(N_Vector v); - /** - * A view to an N_Vector which can be used whenever a N_Vector is required. - * An object of this class will automatically clean up all internal data for - * the view when it is destroyed. This doesn't mean that the actual vector - * is deleted though, which completely depends on how the N_Vector has been - * created. + * A view to a vector which can be used whenever a N_Vector is required. + * + * Objects of this class should preferrably be created by + * make_nvector_view() as + * + * @code + * auto view = make_nvector_view(vector); + * @endcode + * + * The resulting N_Vector is a view of the actual vector and not owning + * memory. Also, N_VDestroy() cannot be called on the resulting N_Vector + * since this would lead to a double delete in the destructor. + * + * @note SUNDIALS will never call N_VDestroy() on a vector it didn't create + * itself and thus the above constraint is not limiting the user. * - * @note N_VDestroy() should not be called on the resulting N_Vector since - * this would lead to a double delete. Let the destructor do this work. + * @tparam VectorType Type of the vector that is stored in @p v */ template class NVectorView { public: /** - * Constructor. This overload is chosen to view a non-const @p vector. + * Default constructor. + * + * @note This constructor exists for compilers that are not eliding the + * copy in a statement like: + * @code + * auto view = make_nvector_view(vector); + * @endcode + */ + NVectorView() = default; + + /** + * Constructor. Create view of @p vector. */ NVectorView(VectorType &vector); /** - * Constructor. This overload is chosen to view a const @p vector. + * Move assignment. */ - NVectorView(const VectorType &vector); + NVectorView(NVectorView &&) noexcept = default; /** - * Implicit conversion to N_Vector. This makes the NVectorView look like - * an actual N_Vector and it can be used directly as an argument. + * Move constructor. */ - operator N_Vector() const; + NVectorView & + operator=(NVectorView &&) noexcept = default; /** - * Access the N_Vector that is viewed by this object. + * Explicitly delete copy ctor. This class is move-only. */ - N_Vector operator->() const; + NVectorView(const NVectorView &) = delete; /** - * Destructor. Automatically release all memory that was only allocated - * for the view. + * Explicitly delete copy assignment. This class is move-only. + */ + NVectorView & + operator=(const NVectorView &) = delete; + + /** + * Destructor. + * + * @note This will not destroy the viewed vector. */ ~NVectorView() = default; + /** + * Implicit conversion to N_Vector. This operator makes the NVectorView + * look like an actual N_Vector and it can be used directly as an + * argument in many SUNDIALS functions. + */ + operator N_Vector() const; + + /** + * Access the N_Vector that is viewed by this object. + */ + N_Vector operator->() const; + private: /** * Actual pointer to a vector viewed by this class. @@ -109,7 +178,6 @@ namespace SUNDIALS std::unique_ptr<_generic_N_Vector, std::function> vector_ptr; }; - } // namespace internal } // namespace SUNDIALS diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index 3d58334f9c..fd693ea597 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -44,7 +44,7 @@ namespace SUNDIALS * field in the SUNDIALS N_Vector module. In addition, this class has a * flag to store whether the stored vector should be treated as const. When * get() is called on a non-const object of this class the flag is checked - * and an exception is thrown if the vector is actually const. Thus we + * and an exception is thrown if the vector is actually const. Thus, we * retain a kind of "runtime const correctness" although the static const * correctness is lost because SUNDIALS N_Vector does not support constness. */ @@ -54,14 +54,14 @@ namespace SUNDIALS public: /** * Create a non-owning content with an existing @p vector. - * @param vector + * @param vector The underlying vector to wrap in thi object. */ NVectorContent(VectorType *vector); /** * Create a non-owning content with an existing const @p vector. If this * constructor is used, access is only allowed via the get() const method. - * @param vector + * @param vector The underlying vector to wrap in thi object. */ NVectorContent(const VectorType *vector); @@ -69,7 +69,7 @@ namespace SUNDIALS * Allocate a new (non-const) vector wrapped in a new content object. The * vector will be deallocated automatically when this object is destroyed. * - * This constructor is called by the N_VClone call of SUNDIALS. + * @note This constructor is intended for the N_VClone() call of SUNDIALS. */ NVectorContent(); @@ -281,12 +281,11 @@ SUNDIALS::internal::NVectorContent::get() { AssertThrow( !is_const, - ExcMessage( - "Tried to access a constant vector content as non-const." - " This most likely happened because a vector that is passed to an" - " create_nvector() call should not be const.\n" - "Alternatively, if you tried to access the vector, use" - " unwrap_nvector_const() for this vector.")); + ExcMessage("Tried to access a constant vector content as non-const." + " This most likely happened because a vector that is passed to a" + " NVectorView() call should not be const.\n" + "Alternatively, if you tried to access the vector, use" + " unwrap_nvector_const() for this vector.")); return vector.get(); } @@ -301,6 +300,15 @@ SUNDIALS::internal::NVectorContent::get() const +template +SUNDIALS::internal::NVectorView +SUNDIALS::internal::make_nvector_view(VectorType &vector) +{ + return NVectorView(vector); +} + + + template VectorType * SUNDIALS::internal::unwrap_nvector(N_Vector v) @@ -328,17 +336,11 @@ SUNDIALS::internal::unwrap_nvector_const(N_Vector v) template SUNDIALS::internal::NVectorView::NVectorView(VectorType &vector) - : vector_ptr(create_nvector(new NVectorContent(&vector)), - [](N_Vector v) { N_VDestroy(v); }) -{} - - - -template -SUNDIALS::internal::NVectorView::NVectorView( - const VectorType &vector) - : vector_ptr(create_nvector(new NVectorContent(&vector)), - [](N_Vector v) { N_VDestroy(v); }) + : vector_ptr( + create_nvector( + new NVectorContent::type>( + &vector)), + [](N_Vector v) { N_VDestroy(v); }) {} @@ -367,7 +369,6 @@ SUNDIALS::internal::create_nvector(NVectorContent *content) N_Vector v = create_empty_nvector(); Assert(v != nullptr, ExcInternalError()); - // Create a non-owning vector content using a pointer and attach to N_Vector v->content = content; Assert(v->content != nullptr, ExcInternalError()); return (v); @@ -387,17 +388,13 @@ SUNDIALS::internal::NVectorOperations::clone_empty(N_Vector w) { Assert(w != nullptr, ExcInternalError()); N_Vector v = N_VNewEmpty(); - if (v == nullptr) - return (nullptr); + Assert(v != nullptr, ExcInternalError()); - // Copy operations - if (N_VCopyOps(w, v)) - { - N_VDestroy(v); - return (nullptr); - } + int status = N_VCopyOps(w, v); + Assert(status == 0, ExcInternalError()); + (void)status; - return (v); + return v; } @@ -411,16 +408,12 @@ SUNDIALS::internal::NVectorOperations::clone(N_Vector w) // the corresponding delete is called in destroy() auto cloned = new NVectorContent(); auto *w_dealii = unwrap_nvector_const(w); - cloned->get()->reinit(*w_dealii); + // reinit the cloned vector based on the layout of the source vector + cloned->get()->reinit(*w_dealii); v->content = cloned; - if (v->content == nullptr) - { - N_VDestroy(v); - return nullptr; - } - return (v); + return v; } @@ -429,6 +422,7 @@ template void SUNDIALS::internal::NVectorOperations::destroy(N_Vector v) { + // support destroying a nullptr because SUNDIALS vectors also do it if (v == nullptr) return; @@ -754,8 +748,7 @@ N_Vector SUNDIALS::internal::create_empty_nvector() { N_Vector v = N_VNewEmpty(); - if (v == nullptr) - return (nullptr); + Assert(v != nullptr, ExcInternalError()); /* constructors, destructors, and utility operations */ v->ops->nvgetvectorid = NVectorOperations::get_vector_id; @@ -806,4 +799,4 @@ DEAL_II_NAMESPACE_CLOSE # endif #endif -#endif \ No newline at end of file +#endif diff --git a/source/sundials/n_vector.inst.in b/source/sundials/n_vector.inst.in index f8a6b6e36a..fa61e4714a 100644 --- a/source/sundials/n_vector.inst.in +++ b/source/sundials/n_vector.inst.in @@ -16,7 +16,10 @@ // serial Vector and BlockVector for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) { - template class SUNDIALS::internal::NVectorView>; + template SUNDIALS::internal::NVectorView> + SUNDIALS::internal::make_nvector_view<>(V &); + template SUNDIALS::internal::NVectorView> + SUNDIALS::internal::make_nvector_view<>(const V &); template V * SUNDIALS::internal::unwrap_nvector>(N_Vector); template const V *SUNDIALS::internal::unwrap_nvector_const>( N_Vector); @@ -25,8 +28,12 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) // parallel Vector and BlockVector for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) { - template class SUNDIALS::internal::NVectorView< - LinearAlgebra::distributed::V>; + template SUNDIALS::internal::NVectorView> + SUNDIALS::internal::make_nvector_view<>(LinearAlgebra::distributed::V &); + template SUNDIALS::internal::NVectorView< + const LinearAlgebra::distributed::V> + SUNDIALS::internal::make_nvector_view<>( + const LinearAlgebra::distributed::V &); template LinearAlgebra::distributed::V *SUNDIALS::internal::unwrap_nvector>( N_Vector); @@ -37,7 +44,10 @@ for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES) // TrilinosWrappers Vector and BlockVector for (V : DEAL_II_VEC_TEMPLATES) { - template class SUNDIALS::internal::NVectorView; + template SUNDIALS::internal::NVectorView + SUNDIALS::internal::make_nvector_view<>(TrilinosWrappers::MPI::V &); + template SUNDIALS::internal::NVectorView + SUNDIALS::internal::make_nvector_view<>(const TrilinosWrappers::MPI::V &); template TrilinosWrappers::MPI::V *SUNDIALS::internal::unwrap_nvector(N_Vector); template const TrilinosWrappers::MPI::V diff --git a/tests/sundials/n_vector.cc b/tests/sundials/n_vector.cc index 06cb51267b..1649bc5343 100644 --- a/tests/sundials/n_vector.cc +++ b/tests/sundials/n_vector.cc @@ -13,6 +13,9 @@ // //----------------------------------------------------------- +// Test SUNDIALS' vector operations on N_Vector implementation. The N_Vectors +// are created by calling NVectorView on one of the internal vector types. + #include "../../include/deal.II/sundials/n_vector.h" #include @@ -23,6 +26,7 @@ #include #include +#include "../../include/deal.II/sundials/n_vector.templates.h" #include #include "../tests.h" @@ -154,10 +158,10 @@ template void test_nvector_view_unwrap() { - auto vector = create_test_vector(); - const auto const_vector = create_test_vector(); - NVectorView n_vector(vector); - NVectorView const_n_vector(const_vector); + auto vector = create_test_vector(); + const auto const_vector = create_test_vector(); + auto n_vector = make_nvector_view(vector); + auto const_n_vector = make_nvector_view(const_vector); Assert(n_vector != nullptr, NVectorTestError()); Assert((n_vector)->content != nullptr, NVectorTestError()); @@ -165,8 +169,8 @@ test_nvector_view_unwrap() Assert(const_n_vector != nullptr, NVectorTestError()); Assert((const_n_vector)->content != nullptr, NVectorTestError()); + // unwrap non-const as non-const auto *vector_unwrapped = unwrap_nvector(n_vector); - // here we check for pointer equality Assert(vector_unwrapped == &vector, NVectorTestError()); // unwrap non-const as const @@ -202,9 +206,9 @@ template void test_get_vector_id() { - auto vector = create_test_vector(); - NVectorView n_vector(vector); - auto id = N_VGetVectorID(n_vector); + auto vector = create_test_vector(); + auto n_vector = make_nvector_view(vector); + auto id = N_VGetVectorID(n_vector); Assert(id == SUNDIALS_NVEC_CUSTOM, NVectorTestError()); deallog << "test_get_vector_id OK" << std::endl; } @@ -215,9 +219,9 @@ template void test_clone() { - auto vector = create_test_vector(); - NVectorView n_vector(vector); - auto cloned = N_VClone(n_vector); + auto vector = create_test_vector(); + auto n_vector = make_nvector_view(vector); + auto cloned = N_VClone(n_vector); Assert(cloned != nullptr, NVectorTestError()); AssertDimension(unwrap_nvector(cloned)->size(), vector.size()); @@ -234,7 +238,7 @@ test_destroy() { GrowingVectorMemory mem; typename GrowingVectorMemory::Pointer vector(mem); - NVectorView n_vector(*vector); + auto n_vector = make_nvector_view(*vector); Assert(n_vector != nullptr, NVectorTestError()); auto cloned = N_VClone(n_vector); @@ -248,9 +252,13 @@ test_destroy() N_VDestroy(cloned); // NOTE: destroying a view vector is not possible and would lead to a double - // delete + // delete at the end of scope when the destructor of NVectorView calls destroy + // again // N_VDestroy(n_vector); + // destroying a nullptr does nothing but is possible + N_VDestroy(nullptr); + // the destructor of NVectorView will do a destroy, so this is also // checked here deallog << "test_destroy OK" << std::endl; @@ -263,8 +271,8 @@ template (); - NVectorView n_vector(vector); + auto vector = create_test_vector(); + auto n_vector = make_nvector_view(vector); Assert(N_VGetCommunicator(n_vector) == nullptr, NVectorTestError()); deallog << "test_get_communicator OK" << std::endl; @@ -277,8 +285,8 @@ template (); - NVectorView n_vector(vector); + auto vector = create_test_vector(); + auto n_vector = make_nvector_view(vector); Assert(N_VGetCommunicator(n_vector) == MPI_COMM_WORLD, NVectorTestError()); deallog << "test_get_communicator OK" << std::endl; @@ -290,8 +298,8 @@ template void test_length() { - auto vector = create_test_vector(); - NVectorView n_vector(vector); + auto vector = create_test_vector(); + auto n_vector = make_nvector_view(vector); Assert(N_VGetLength(n_vector) == static_cast(vector.size()), NVectorTestError()); @@ -309,9 +317,9 @@ test_linear_sum() auto vc = create_test_vector(); auto expected = create_test_vector(); - NVectorView nv_a(va); - NVectorView nv_b(vb); - NVectorView nv_c(vc); + auto nv_a = make_nvector_view(va); + auto nv_b = make_nvector_view(vb); + auto nv_c = make_nvector_view(vc); va = 1.0; vb = 2.0; @@ -346,8 +354,8 @@ test_dot_product() const auto vb = create_test_vector(3.0); const auto size = va.size(); - NVectorView nv_a(va); - NVectorView nv_b(vb); + auto nv_a = make_nvector_view(va); + auto nv_b = make_nvector_view(vb); auto result = N_VDotProd(nv_a, nv_b); auto expected = 6.0 * size; @@ -371,7 +379,7 @@ test_set_constant() auto expected = create_test_vector(); expected = 1.0; - NVectorView n_vector(vector); + auto n_vector = make_nvector_view(vector); N_VConst(1.0, n_vector); Assert(vector_equal(vector, expected), NVectorTestError()); @@ -393,8 +401,8 @@ test_add_constant() auto result = create_test_vector(-1.0); auto expected = create_test_vector(3.0); - NVectorView nv_a(vector_a); - NVectorView nv_result(result); + auto nv_a = make_nvector_view(vector_a); + auto nv_result = make_nvector_view(result); N_VAddConst(nv_a, 1.0, nv_result); Assert(vector_equal(result, expected), NVectorTestError()); @@ -418,9 +426,9 @@ test_elementwise_product() auto vector_result = create_test_vector(); auto expected = create_test_vector(6.0); - NVectorView nv_a(vector_a); - NVectorView nv_b(vector_b); - NVectorView nv_result(vector_result); + auto nv_a = make_nvector_view(vector_a); + auto nv_b = make_nvector_view(vector_b); + auto nv_result = make_nvector_view(vector_result); // result = a.*b N_VProd(nv_a, nv_b, nv_result); @@ -428,7 +436,7 @@ test_elementwise_product() auto vector_c = create_test_vector(2.0); - NVectorView nv_c(vector_c); + auto nv_c = make_nvector_view(vector_c); // c .*= b N_VProd(nv_c, nv_b, nv_c); @@ -457,9 +465,9 @@ test_elementwise_div() auto vector_result = create_test_vector(); auto expected = create_test_vector(2.0); - NVectorView nv_a(vector_a); - NVectorView nv_b(vector_b); - NVectorView nv_result(vector_result); + auto nv_a = make_nvector_view(vector_a); + auto nv_b = make_nvector_view(vector_b); + auto nv_result = make_nvector_view(vector_result); N_VDiv(nv_a, nv_b, nv_result); Assert(vector_equal(vector_result, expected), NVectorTestError()); @@ -477,8 +485,8 @@ test_elementwise_inv() auto vector_result = create_test_vector(); auto expected = create_test_vector(1.0 / 6.0); - NVectorView nv_a(vector_a); - NVectorView nv_result(vector_result); + auto nv_a = make_nvector_view(vector_a); + auto nv_result = make_nvector_view(vector_result); N_VInv(nv_a, nv_result); Assert(vector_equal(vector_result, expected), NVectorTestError()); @@ -501,8 +509,8 @@ test_elementwise_abs() auto vector_result = create_test_vector(); auto expected = create_test_vector(2.0); - NVectorView nv_a(vector_a); - NVectorView nv_result(vector_result); + auto nv_a = make_nvector_view(vector_a); + auto nv_result = make_nvector_view(vector_result); N_VAbs(nv_a, nv_result); Assert(vector_equal(vector_result, expected), NVectorTestError()); @@ -524,8 +532,8 @@ test_weighted_rms_norm() const auto vector_a = create_test_vector(2.0); const auto vector_b = create_test_vector(3.0); - NVectorView nv_a(vector_a); - NVectorView nv_b(vector_b); + auto nv_a = make_nvector_view(vector_a); + auto nv_b = make_nvector_view(vector_b); const auto result = N_VWrmsNorm(nv_a, nv_b); Assert(std::fabs(result - 6.0) < 1e-12, NVectorTestError()); @@ -542,8 +550,8 @@ test_max_norm() const auto vector_a = create_test_vector(2.0); const auto vector_b = create_test_vector(-3.0); - NVectorView nv_a(vector_a); - NVectorView nv_b(vector_b); + auto nv_a = make_nvector_view(vector_a); + auto nv_b = make_nvector_view(vector_b); auto result = N_VMaxNorm(nv_a); Assert(std::fabs(result - 2.0) < 1e-12, NVectorTestError()); @@ -564,8 +572,8 @@ test_scale() auto vector_b = create_test_vector(2.0); const auto expected = create_test_vector(4.0); - NVectorView nv_a(vector_a); - NVectorView nv_b(vector_b); + auto nv_a = make_nvector_view(vector_a); + auto nv_b = make_nvector_view(vector_b); // b *= 2 N_VScale(2.0, nv_b, nv_b); -- 2.39.5