DEAL_II_NAMESPACE_OPEN
+# ifndef DOXYGEN
namespace SUNDIALS
{
namespace internal
{
+ template <typename VectorType>
+ 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 <typename VectorType>
+ NVectorView<VectorType>
+ 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
* @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 <typename VectorType>
VectorType *
* @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 <typename VectorType>
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 <typename VectorType>
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.
std::unique_ptr<_generic_N_Vector, std::function<void(N_Vector)>>
vector_ptr;
};
-
} // namespace internal
} // 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.
*/
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);
* 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();
{
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();
}
+template <typename VectorType>
+SUNDIALS::internal::NVectorView<VectorType>
+SUNDIALS::internal::make_nvector_view(VectorType &vector)
+{
+ return NVectorView<VectorType>(vector);
+}
+
+
+
template <typename VectorType>
VectorType *
SUNDIALS::internal::unwrap_nvector(N_Vector v)
template <typename VectorType>
SUNDIALS::internal::NVectorView<VectorType>::NVectorView(VectorType &vector)
- : vector_ptr(create_nvector(new NVectorContent<VectorType>(&vector)),
- [](N_Vector v) { N_VDestroy(v); })
-{}
-
-
-
-template <typename VectorType>
-SUNDIALS::internal::NVectorView<VectorType>::NVectorView(
- const VectorType &vector)
- : vector_ptr(create_nvector(new NVectorContent<VectorType>(&vector)),
- [](N_Vector v) { N_VDestroy(v); })
+ : vector_ptr(
+ create_nvector(
+ new NVectorContent<typename std::remove_const<VectorType>::type>(
+ &vector)),
+ [](N_Vector v) { N_VDestroy(v); })
{}
N_Vector v = create_empty_nvector<VectorType>();
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);
{
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;
}
// the corresponding delete is called in destroy()
auto cloned = new NVectorContent<VectorType>();
auto *w_dealii = unwrap_nvector_const<VectorType>(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;
}
void
SUNDIALS::internal::NVectorOperations::destroy(N_Vector v)
{
+ // support destroying a nullptr because SUNDIALS vectors also do it
if (v == nullptr)
return;
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;
# endif
#endif
-#endif
\ No newline at end of file
+#endif
// serial Vector and BlockVector
for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
{
- template class SUNDIALS::internal::NVectorView<V<S>>;
+ template SUNDIALS::internal::NVectorView<V<S>>
+ SUNDIALS::internal::make_nvector_view<>(V<S> &);
+ template SUNDIALS::internal::NVectorView<const V<S>>
+ SUNDIALS::internal::make_nvector_view<>(const V<S> &);
template V<S> * SUNDIALS::internal::unwrap_nvector<V<S>>(N_Vector);
template const V<S> *SUNDIALS::internal::unwrap_nvector_const<V<S>>(
N_Vector);
// parallel Vector and BlockVector
for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
{
- template class SUNDIALS::internal::NVectorView<
- LinearAlgebra::distributed::V<S>>;
+ template SUNDIALS::internal::NVectorView<LinearAlgebra::distributed::V<S>>
+ SUNDIALS::internal::make_nvector_view<>(LinearAlgebra::distributed::V<S> &);
+ template SUNDIALS::internal::NVectorView<
+ const LinearAlgebra::distributed::V<S>>
+ SUNDIALS::internal::make_nvector_view<>(
+ const LinearAlgebra::distributed::V<S> &);
template LinearAlgebra::distributed::V<S>
*SUNDIALS::internal::unwrap_nvector<LinearAlgebra::distributed::V<S>>(
N_Vector);
// TrilinosWrappers Vector and BlockVector
for (V : DEAL_II_VEC_TEMPLATES)
{
- template class SUNDIALS::internal::NVectorView<TrilinosWrappers::MPI::V>;
+ template SUNDIALS::internal::NVectorView<TrilinosWrappers::MPI::V>
+ SUNDIALS::internal::make_nvector_view<>(TrilinosWrappers::MPI::V &);
+ template SUNDIALS::internal::NVectorView<const TrilinosWrappers::MPI::V>
+ SUNDIALS::internal::make_nvector_view<>(const TrilinosWrappers::MPI::V &);
template TrilinosWrappers::MPI::V
*SUNDIALS::internal::unwrap_nvector<TrilinosWrappers::MPI::V>(N_Vector);
template const TrilinosWrappers::MPI::V
//
//-----------------------------------------------------------
+// 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 <deal.II/base/logstream.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_vector.h>
+#include "../../include/deal.II/sundials/n_vector.templates.h"
#include <deal.II/sundials/n_vector.h>
#include "../tests.h"
void
test_nvector_view_unwrap()
{
- auto vector = create_test_vector<VectorType>();
- const auto const_vector = create_test_vector<VectorType>();
- NVectorView<VectorType> n_vector(vector);
- NVectorView<VectorType> const_n_vector(const_vector);
+ auto vector = create_test_vector<VectorType>();
+ const auto const_vector = create_test_vector<VectorType>();
+ 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());
Assert(const_n_vector != nullptr, NVectorTestError());
Assert((const_n_vector)->content != nullptr, NVectorTestError());
+ // unwrap non-const as non-const
auto *vector_unwrapped = unwrap_nvector<VectorType>(n_vector);
- // here we check for pointer equality
Assert(vector_unwrapped == &vector, NVectorTestError());
// unwrap non-const as const
void
test_get_vector_id()
{
- auto vector = create_test_vector<VectorType>();
- NVectorView<VectorType> n_vector(vector);
- auto id = N_VGetVectorID(n_vector);
+ auto vector = create_test_vector<VectorType>();
+ 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;
}
void
test_clone()
{
- auto vector = create_test_vector<VectorType>();
- NVectorView<VectorType> n_vector(vector);
- auto cloned = N_VClone(n_vector);
+ auto vector = create_test_vector<VectorType>();
+ auto n_vector = make_nvector_view(vector);
+ auto cloned = N_VClone(n_vector);
Assert(cloned != nullptr, NVectorTestError());
AssertDimension(unwrap_nvector<VectorType>(cloned)->size(), vector.size());
{
GrowingVectorMemory<VectorType> mem;
typename GrowingVectorMemory<VectorType>::Pointer vector(mem);
- NVectorView<VectorType> n_vector(*vector);
+ auto n_vector = make_nvector_view(*vector);
Assert(n_vector != nullptr, NVectorTestError());
auto cloned = N_VClone(n_vector);
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;
void
test_get_communicator()
{
- auto vector = create_test_vector<VectorType>();
- NVectorView<VectorType> n_vector(vector);
+ auto vector = create_test_vector<VectorType>();
+ auto n_vector = make_nvector_view(vector);
Assert(N_VGetCommunicator(n_vector) == nullptr, NVectorTestError());
deallog << "test_get_communicator OK" << std::endl;
void
test_get_communicator()
{
- auto vector = create_test_vector<VectorType>();
- NVectorView<VectorType> n_vector(vector);
+ auto vector = create_test_vector<VectorType>();
+ auto n_vector = make_nvector_view(vector);
Assert(N_VGetCommunicator(n_vector) == MPI_COMM_WORLD, NVectorTestError());
deallog << "test_get_communicator OK" << std::endl;
void
test_length()
{
- auto vector = create_test_vector<VectorType>();
- NVectorView<VectorType> n_vector(vector);
+ auto vector = create_test_vector<VectorType>();
+ auto n_vector = make_nvector_view(vector);
Assert(N_VGetLength(n_vector) == static_cast<int>(vector.size()),
NVectorTestError());
auto vc = create_test_vector<VectorType>();
auto expected = create_test_vector<VectorType>();
- NVectorView<VectorType> nv_a(va);
- NVectorView<VectorType> nv_b(vb);
- NVectorView<VectorType> 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;
const auto vb = create_test_vector<VectorType>(3.0);
const auto size = va.size();
- NVectorView<VectorType> nv_a(va);
- NVectorView<VectorType> 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;
auto expected = create_test_vector<VectorType>();
expected = 1.0;
- NVectorView<VectorType> n_vector(vector);
+ auto n_vector = make_nvector_view(vector);
N_VConst(1.0, n_vector);
Assert(vector_equal(vector, expected), NVectorTestError());
auto result = create_test_vector<VectorType>(-1.0);
auto expected = create_test_vector<VectorType>(3.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> 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());
auto vector_result = create_test_vector<VectorType>();
auto expected = create_test_vector<VectorType>(6.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> nv_b(vector_b);
- NVectorView<VectorType> 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);
auto vector_c = create_test_vector<VectorType>(2.0);
- NVectorView<VectorType> nv_c(vector_c);
+ auto nv_c = make_nvector_view(vector_c);
// c .*= b
N_VProd(nv_c, nv_b, nv_c);
auto vector_result = create_test_vector<VectorType>();
auto expected = create_test_vector<VectorType>(2.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> nv_b(vector_b);
- NVectorView<VectorType> 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());
auto vector_result = create_test_vector<VectorType>();
auto expected = create_test_vector<VectorType>(1.0 / 6.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> 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());
auto vector_result = create_test_vector<VectorType>();
auto expected = create_test_vector<VectorType>(2.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> 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());
const auto vector_a = create_test_vector<VectorType>(2.0);
const auto vector_b = create_test_vector<VectorType>(3.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> 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());
const auto vector_a = create_test_vector<VectorType>(2.0);
const auto vector_b = create_test_vector<VectorType>(-3.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> 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());
auto vector_b = create_test_vector<VectorType>(2.0);
const auto expected = create_test_vector<VectorType>(4.0);
- NVectorView<VectorType> nv_a(vector_a);
- NVectorView<VectorType> 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);
projects / dealii.git / commitdiff