--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#ifndef dealii_sundials_n_vector_h
+#define dealii_sundials_n_vector_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
+
+# include <sundials/sundials_nvector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+ namespace internal
+ {
+ /**
+ * Retrieve the underlying vector attached to N_Vector @p v. This call will
+ * only succeed if the underlying vector is not const. Use
+ * unwrap_nvector_const() for this case.
+ *
+ * @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
+ */
+ template <typename VectorType>
+ VectorType *
+ unwrap_nvector(N_Vector v);
+
+ /**
+ * Retrieve the underlying vector attached to N_Vector @p v as a constant
+ * pointer.
+ *
+ * @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
+ */
+ 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.
+ *
+ * @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.
+ */
+ template <typename VectorType>
+ class NVectorView
+ {
+ public:
+ /**
+ * Constructor. This overload is chosen to view a non-const @p vector.
+ */
+ NVectorView(VectorType &vector);
+
+ /**
+ * Constructor. This overload is chosen to view a const @p vector.
+ */
+ NVectorView(const VectorType &vector);
+
+ /**
+ * Implicit conversion to N_Vector. This makes the NVectorView look like
+ * an actual N_Vector and it can be used directly as an argument.
+ */
+ operator N_Vector() const;
+
+ /**
+ * Access the N_Vector that is viewed by this object.
+ */
+ N_Vector operator->() const;
+
+ /**
+ * Destructor. Automatically release all memory that was only allocated
+ * for the view.
+ */
+ ~NVectorView() = default;
+
+ 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
+
+DEAL_II_NAMESPACE_CLOSE
+
+# endif
+#endif
+#endif
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#ifndef dealii_sundials_n_vector_templates_h
+#define dealii_sundials_n_vector_templates_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/sundials/n_vector.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
+
+# include <deal.II/base/exceptions.h>
+
+# include <deal.II/lac/block_vector.h>
+# include <deal.II/lac/la_parallel_block_vector.h>
+# include <deal.II/lac/la_parallel_vector.h>
+# include <deal.II/lac/la_vector.h>
+# include <deal.II/lac/trilinos_parallel_block_vector.h>
+# include <deal.II/lac/trilinos_vector.h>
+# include <deal.II/lac/vector_memory.h>
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+ namespace internal
+ {
+ /**
+ * An internal class that stores a pointer to a vector and manages the
+ * memory if necessary. Objects of this class are used as the `content`
+ * 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
+ * retain a kind of "runtime const correctness" although the static const
+ * correctness is lost because SUNDIALS N_Vector does not support constness.
+ */
+ template <typename VectorType>
+ class NVectorContent
+ {
+ public:
+ /**
+ * Create a non-owning content with an existing @p vector.
+ * @param vector
+ */
+ 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
+ */
+ 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.
+ */
+ NVectorContent();
+
+ /**
+ * Non-const access to the stored vector. Only allowed if a constructor
+ * different than NVectorContent(const VectorType *vector) was used.
+ * @return
+ */
+ VectorType *
+ get();
+
+ /**
+ * Const access to the stored vector. Always allowed.
+ */
+ const VectorType *
+ get() const;
+
+ private:
+ using PointerType =
+ std::unique_ptr<VectorType, std::function<void(VectorType *)>>;
+
+ /**
+ * Vector memory which might be used to allocate storage if
+ * NVectorContent() is called.
+ */
+ GrowingVectorMemory<VectorType> mem;
+
+ /**
+ * Actually stored vector content.
+ */
+ PointerType vector;
+
+ /**
+ * Flag storing whether the stored pointer is to be treated as const. If
+ * the pointer passed in the constructor was indeed const, it is cast away
+ * but this flag will be set to true. Access to the pointer must then
+ * check that this flag is set correctly.
+ */
+ const bool is_const;
+ };
+
+ /**
+ * Helper to create a vector with all operations and the given @p content.
+ *
+ * @param content The vector content to attach to the N_Vector.
+ * @return A new N_Vector
+ */
+ template <typename VectorType>
+ N_Vector
+ create_nvector(NVectorContent<VectorType> *content);
+
+ /**
+ * Helper to create an empty vector with all operations but no content.
+ * @return A new N_Vector
+ */
+ template <typename VectorType>
+ N_Vector
+ create_empty_nvector();
+
+ /**
+ * Collection of all operations specified by SUNDIALS N_Vector
+ * documentation. These functions are attached to the generic N_Vector
+ * structure.
+ */
+ namespace NVectorOperations
+ {
+ N_Vector_ID
+ get_vector_id(N_Vector v);
+
+ N_Vector
+ clone_empty(N_Vector w);
+
+ template <typename VectorType>
+ N_Vector
+ clone(N_Vector w);
+
+ template <typename VectorType>
+ void
+ destroy(N_Vector v);
+
+ template <typename VectorType>
+ sunindextype
+ get_global_length(N_Vector v);
+
+ template <typename VectorType>
+ void
+ linear_sum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z);
+
+ template <typename VectorType>
+ realtype
+ dot_product(N_Vector x, N_Vector y);
+
+ template <typename VectorType>
+ void
+ elementwise_product(N_Vector x, N_Vector y, N_Vector z);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<!IsBlockVector<VectorType>::value, int> = 0>
+ void
+ elementwise_div(N_Vector x, N_Vector y, N_Vector z);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<IsBlockVector<VectorType>::value, int> = 0>
+ void
+ elementwise_div(N_Vector x, N_Vector y, N_Vector z);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<!IsBlockVector<VectorType>::value, int> = 0>
+ void
+ elementwise_inv(N_Vector x, N_Vector z);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<IsBlockVector<VectorType>::value, int> = 0>
+ void
+ elementwise_inv(N_Vector x, N_Vector z);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<!IsBlockVector<VectorType>::value, int> = 0>
+ void
+ elementwise_abs(N_Vector x, N_Vector z);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<IsBlockVector<VectorType>::value, int> = 0>
+ void
+ elementwise_abs(N_Vector x, N_Vector z);
+
+ template <typename VectorType>
+ realtype
+ weighted_rms_norm(N_Vector x, N_Vector w);
+
+ template <typename VectorType>
+ realtype
+ max_norm(N_Vector x);
+
+ template <typename VectorType>
+ void
+ scale(realtype c, N_Vector x, N_Vector z);
+
+ template <typename VectorType>
+ void
+ set_constant(realtype c, N_Vector v);
+
+ template <typename VectorType>
+ void
+ add_constant(N_Vector x, realtype b, N_Vector z);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<is_serial_vector<VectorType>::value, int> = 0>
+ void *get_communicator(N_Vector);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<!is_serial_vector<VectorType>::value &&
+ !IsBlockVector<VectorType>::value,
+ int> = 0>
+ void *
+ get_communicator(N_Vector v);
+
+ template <
+ typename VectorType,
+ typename std::enable_if_t<!is_serial_vector<VectorType>::value &&
+ IsBlockVector<VectorType>::value,
+ int> = 0>
+ void *
+ get_communicator(N_Vector v);
+ } // namespace NVectorOperations
+ } // namespace internal
+} // namespace SUNDIALS
+
+
+
+template <typename VectorType>
+SUNDIALS::internal::NVectorContent<VectorType>::NVectorContent()
+ : vector(typename VectorMemory<VectorType>::Pointer(mem))
+ , is_const(false)
+{}
+
+
+
+template <typename VectorType>
+SUNDIALS::internal::NVectorContent<VectorType>::NVectorContent(
+ VectorType *vector)
+ : vector(vector, [](VectorType *) { /* not owning memory -> don't free*/ })
+ , is_const(false)
+{}
+
+
+
+template <typename VectorType>
+SUNDIALS::internal::NVectorContent<VectorType>::NVectorContent(
+ const VectorType *vector)
+ : vector(const_cast<VectorType *>(vector),
+ [](VectorType *) { /* not owning memory -> don't free*/ })
+ , is_const(true)
+{}
+
+
+
+template <typename VectorType>
+VectorType *
+SUNDIALS::internal::NVectorContent<VectorType>::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."));
+ return vector.get();
+}
+
+
+
+template <typename VectorType>
+const VectorType *
+SUNDIALS::internal::NVectorContent<VectorType>::get() const
+{
+ return vector.get();
+}
+
+
+
+template <typename VectorType>
+VectorType *
+SUNDIALS::internal::unwrap_nvector(N_Vector v)
+{
+ Assert(v != nullptr, ExcInternalError());
+ Assert(v->content != nullptr, ExcInternalError());
+ auto *pContent = reinterpret_cast<NVectorContent<VectorType> *>(v->content);
+ return pContent->get();
+}
+
+
+
+template <typename VectorType>
+const VectorType *
+SUNDIALS::internal::unwrap_nvector_const(N_Vector v)
+{
+ Assert(v != nullptr, ExcInternalError());
+ Assert(v->content != nullptr, ExcInternalError());
+ auto *pContent =
+ reinterpret_cast<const NVectorContent<VectorType> *>(v->content);
+ return pContent->get();
+}
+
+
+
+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); })
+{}
+
+
+
+template <typename VectorType>
+SUNDIALS::internal::NVectorView<VectorType>::operator N_Vector() const
+{
+ return vector_ptr.get();
+}
+
+
+
+template <typename VectorType>
+N_Vector SUNDIALS::internal::NVectorView<VectorType>::operator->() const
+{
+ return vector_ptr.get();
+}
+
+
+
+template <typename VectorType>
+N_Vector
+SUNDIALS::internal::create_nvector(NVectorContent<VectorType> *content)
+{
+ // Create an N_Vector with operators attached and empty content
+ 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);
+}
+
+
+
+N_Vector_ID SUNDIALS::internal::NVectorOperations::get_vector_id(N_Vector)
+{
+ return SUNDIALS_NVEC_CUSTOM;
+}
+
+
+
+N_Vector
+SUNDIALS::internal::NVectorOperations::clone_empty(N_Vector w)
+{
+ Assert(w != nullptr, ExcInternalError());
+ N_Vector v = N_VNewEmpty();
+ if (v == nullptr)
+ return (nullptr);
+
+ // Copy operations
+ if (N_VCopyOps(w, v))
+ {
+ N_VDestroy(v);
+ return (nullptr);
+ }
+
+ return (v);
+}
+
+
+
+template <typename VectorType>
+N_Vector
+SUNDIALS::internal::NVectorOperations::clone(N_Vector w)
+{
+ N_Vector v = clone_empty(w);
+
+ // 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);
+
+ v->content = cloned;
+ if (v->content == nullptr)
+ {
+ N_VDestroy(v);
+ return nullptr;
+ }
+
+ return (v);
+}
+
+
+
+template <typename VectorType>
+void
+SUNDIALS::internal::NVectorOperations::destroy(N_Vector v)
+{
+ if (v == nullptr)
+ return;
+
+ if (v->content != nullptr)
+ {
+ auto *content =
+ reinterpret_cast<NVectorContent<VectorType> *>(v->content);
+ // the NVectorContent knows if it owns the memory and will free
+ // correctly
+ delete content;
+ v->content = nullptr;
+ }
+
+ /* free ops and vector */
+ if (v->ops != nullptr)
+ {
+ free(v->ops);
+ v->ops = nullptr;
+ }
+ free(v);
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<is_serial_vector<VectorType>::value, int>>
+void *SUNDIALS::internal::NVectorOperations::get_communicator(N_Vector)
+{
+ return nullptr;
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<!is_serial_vector<VectorType>::value &&
+ IsBlockVector<VectorType>::value,
+ int>>
+void *
+SUNDIALS::internal::NVectorOperations::get_communicator(N_Vector v)
+{
+ return unwrap_nvector_const<VectorType>(v)->block(0).get_mpi_communicator();
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<!is_serial_vector<VectorType>::value &&
+ !IsBlockVector<VectorType>::value,
+ int>>
+void *
+SUNDIALS::internal::NVectorOperations::get_communicator(N_Vector v)
+{
+ return unwrap_nvector_const<VectorType>(v)->get_mpi_communicator();
+}
+
+
+
+template <typename VectorType>
+sunindextype
+SUNDIALS::internal::NVectorOperations::get_global_length(N_Vector v)
+{
+ return unwrap_nvector_const<VectorType>(v)->size();
+}
+
+
+
+template <typename VectorType>
+void
+SUNDIALS::internal::NVectorOperations::linear_sum(realtype a,
+ N_Vector x,
+ realtype b,
+ N_Vector y,
+ N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *y_dealii = unwrap_nvector_const<VectorType>(y);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ if (z_dealii == x_dealii)
+ z_dealii->sadd(a, b, *y_dealii);
+ else if (z_dealii == y_dealii)
+ z_dealii->sadd(b, a, *x_dealii);
+ else
+ {
+ *z_dealii = 0;
+ z_dealii->add(a, *x_dealii, b, *y_dealii);
+ }
+}
+
+
+
+template <typename VectorType>
+realtype
+SUNDIALS::internal::NVectorOperations::dot_product(N_Vector x, N_Vector y)
+{
+ return *unwrap_nvector_const<VectorType>(x) *
+ *unwrap_nvector_const<VectorType>(y);
+}
+
+
+
+template <typename VectorType>
+void
+SUNDIALS::internal::NVectorOperations::set_constant(realtype c, N_Vector v)
+{
+ auto *v_dealii = unwrap_nvector<VectorType>(v);
+ *v_dealii = c;
+}
+
+
+
+template <typename VectorType>
+void
+SUNDIALS::internal::NVectorOperations::add_constant(N_Vector x,
+ realtype b,
+ N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ if (x_dealii == z_dealii)
+ z_dealii->add(b);
+ else
+ {
+ *z_dealii = *x_dealii;
+ z_dealii->add(b);
+ }
+}
+
+
+
+template <typename VectorType>
+void
+SUNDIALS::internal::NVectorOperations::elementwise_product(N_Vector x,
+ N_Vector y,
+ N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *y_dealii = unwrap_nvector_const<VectorType>(y);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+ if (z_dealii == x_dealii)
+ z_dealii->scale(*y_dealii);
+ else if (z_dealii == y_dealii)
+ z_dealii->scale(*x_dealii);
+ else
+ {
+ *z_dealii = *y_dealii;
+ z_dealii->scale(*x_dealii);
+ }
+}
+
+
+
+template <typename VectorType>
+realtype
+SUNDIALS::internal::NVectorOperations::weighted_rms_norm(N_Vector x, N_Vector w)
+{
+ // TODO copy can be avoided by a custom kernel
+ VectorType tmp = *unwrap_nvector_const<VectorType>(x);
+ auto * w_dealii = unwrap_nvector_const<VectorType>(w);
+ const auto n = tmp.size();
+ tmp.scale(*w_dealii);
+ return tmp.l2_norm() / std::sqrt(n);
+}
+
+
+
+template <typename VectorType>
+realtype
+SUNDIALS::internal::NVectorOperations::max_norm(N_Vector x)
+{
+ return unwrap_nvector_const<VectorType>(x)->linfty_norm();
+}
+
+
+
+template <typename VectorType>
+void
+SUNDIALS::internal::NVectorOperations::scale(realtype c, N_Vector x, N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ if (x_dealii == z_dealii)
+ (*z_dealii) *= c;
+ else
+ z_dealii->sadd(0.0, c, *x_dealii);
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<!IsBlockVector<VectorType>::value, int>>
+void
+SUNDIALS::internal::NVectorOperations::elementwise_div(N_Vector x,
+ N_Vector y,
+ N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *y_dealii = unwrap_nvector_const<VectorType>(y);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ AssertDimension(x_dealii->size(), z_dealii->size());
+ AssertDimension(x_dealii->size(), y_dealii->size());
+
+ std::transform(x_dealii->begin(),
+ x_dealii->end(),
+ y_dealii->begin(),
+ z_dealii->begin(),
+ std::divides<>{});
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<IsBlockVector<VectorType>::value, int>>
+void
+SUNDIALS::internal::NVectorOperations::elementwise_div(N_Vector x,
+ N_Vector y,
+ N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *y_dealii = unwrap_nvector_const<VectorType>(y);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ AssertDimension(x_dealii->size(), z_dealii->size());
+ AssertDimension(x_dealii->size(), y_dealii->size());
+ AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks());
+ AssertDimension(x_dealii->n_blocks(), y_dealii->n_blocks());
+
+ for (unsigned i = 0; i < x_dealii->n_blocks(); ++i)
+ std::transform(x_dealii->block(i).begin(),
+ x_dealii->block(i).end(),
+ y_dealii->block(i).begin(),
+ z_dealii->block(i).begin(),
+ std::divides<>{});
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<!IsBlockVector<VectorType>::value, int>>
+void
+SUNDIALS::internal::NVectorOperations::elementwise_inv(N_Vector x, N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ AssertDimension(x_dealii->size(), z_dealii->size());
+
+ std::transform(x_dealii->begin(),
+ x_dealii->end(),
+ z_dealii->begin(),
+ [](typename VectorType::value_type v) { return 1.0 / v; });
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<IsBlockVector<VectorType>::value, int>>
+void
+SUNDIALS::internal::NVectorOperations::elementwise_inv(N_Vector x, N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ AssertDimension(x_dealii->size(), z_dealii->size());
+ AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks());
+
+ for (unsigned i = 0; i < x_dealii->n_blocks(); ++i)
+ std::transform(x_dealii->block(i).begin(),
+ x_dealii->block(i).end(),
+ z_dealii->block(i).begin(),
+ [](typename VectorType::value_type v) { return 1.0 / v; });
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<!IsBlockVector<VectorType>::value, int>>
+void
+SUNDIALS::internal::NVectorOperations::elementwise_abs(N_Vector x, N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ AssertDimension(x_dealii->size(), z_dealii->size());
+
+ std::transform(x_dealii->begin(),
+ x_dealii->end(),
+ z_dealii->begin(),
+ [](typename VectorType::value_type v) {
+ return std::fabs(v);
+ });
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<IsBlockVector<VectorType>::value, int>>
+void
+SUNDIALS::internal::NVectorOperations::elementwise_abs(N_Vector x, N_Vector z)
+{
+ auto *x_dealii = unwrap_nvector_const<VectorType>(x);
+ auto *z_dealii = unwrap_nvector<VectorType>(z);
+
+ AssertDimension(x_dealii->size(), z_dealii->size());
+ AssertDimension(x_dealii->n_blocks(), z_dealii->n_blocks());
+
+ for (unsigned i = 0; i < x_dealii->n_blocks(); ++i)
+ std::transform(x_dealii->block(i).begin(),
+ x_dealii->block(i).end(),
+ z_dealii->block(i).begin(),
+ [](typename VectorType::value_type v) {
+ return std::fabs(v);
+ });
+}
+
+
+
+template <typename VectorType>
+N_Vector
+SUNDIALS::internal::create_empty_nvector()
+{
+ N_Vector v = N_VNewEmpty();
+ if (v == nullptr)
+ return (nullptr);
+
+ /* constructors, destructors, and utility operations */
+ v->ops->nvgetvectorid = NVectorOperations::get_vector_id;
+ v->ops->nvclone = NVectorOperations::clone<VectorType>;
+ v->ops->nvcloneempty = NVectorOperations::clone_empty;
+ v->ops->nvdestroy = NVectorOperations::destroy<VectorType>;
+ // v->ops->nvspace = undef;
+ v->ops->nvgetcommunicator = NVectorOperations::get_communicator<VectorType>;
+ v->ops->nvgetlength = NVectorOperations::get_global_length<VectorType>;
+
+ /* standard vector operations */
+ v->ops->nvlinearsum = NVectorOperations::linear_sum<VectorType>;
+ v->ops->nvconst = NVectorOperations::set_constant<VectorType>;
+ v->ops->nvprod = NVectorOperations::elementwise_product<VectorType>;
+ v->ops->nvdiv = NVectorOperations::elementwise_div<VectorType>;
+ v->ops->nvscale = NVectorOperations::scale<VectorType>;
+ v->ops->nvabs = NVectorOperations::elementwise_abs<VectorType>;
+ v->ops->nvinv = NVectorOperations::elementwise_inv<VectorType>;
+ v->ops->nvaddconst = NVectorOperations::add_constant<VectorType>;
+ v->ops->nvdotprod = NVectorOperations::dot_product<VectorType>;
+ v->ops->nvmaxnorm = NVectorOperations::max_norm<VectorType>;
+ v->ops->nvwrmsnorm = NVectorOperations::weighted_rms_norm<VectorType>;
+ // v->ops->nvwrmsnormmask = undef;
+ // v->ops->nvmin = undef;
+ // v->ops->nvwl2norm = undef;
+ // v->ops->nvl1norm = undef;
+ // v->ops->nvcompare = undef;
+ // v->ops->nvinvtest = undef;
+ // v->ops->nvconstrmask = undef;
+ // v->ops->nvminquotient = undef;
+
+ /* fused and vector array operations are disabled (NULL) by default */
+
+ /* local reduction operations */
+ // v->ops->nvdotprodlocal = undef;
+ // v->ops->nvmaxnormlocal = undef;
+ // v->ops->nvminlocal = undef;
+ // v->ops->nvl1normlocal = undef;
+ // v->ops->nvinvtestlocal = undef;
+ // v->ops->nvconstrmasklocal = undef;
+ // v->ops->nvminquotientlocal = undef;
+ // v->ops->nvwsqrsumlocal = undef;
+ // v->ops->nvwsqrsummasklocal = undef;
+ return (v);
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+# endif
+#endif
+#endif
\ No newline at end of file
ida.cc
copy.cc
kinsol.cc
+ n_vector.cc
)
SET(_inst
+ n_vector.inst.in
)
FILE(GLOB _header
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/sundials/n_vector.h>
+#include <deal.II/sundials/n_vector.templates.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+# if DEAL_II_SUNDIALS_VERSION_GTE(5, 0, 0)
+
+DEAL_II_NAMESPACE_OPEN
+
+# include "n_vector.inst"
+
+DEAL_II_NAMESPACE_CLOSE
+
+# endif
+#endif
\ No newline at end of file
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// serial Vector and BlockVector
+for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)
+ {
+ template class SUNDIALS::internal::NVectorView<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 LinearAlgebra::distributed::V<S>
+ *SUNDIALS::internal::unwrap_nvector<LinearAlgebra::distributed::V<S>>(
+ N_Vector);
+ template const LinearAlgebra::distributed::V<S> *SUNDIALS::internal::
+ unwrap_nvector_const<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 TrilinosWrappers::MPI::V
+ *SUNDIALS::internal::unwrap_nvector<TrilinosWrappers::MPI::V>(N_Vector);
+ template const TrilinosWrappers::MPI::V
+ *SUNDIALS::internal::unwrap_nvector_const<TrilinosWrappers::MPI::V>(
+ N_Vector);
+ }
--- /dev/null
+//-----------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+
+#include "../../include/deal.II/sundials/n_vector.h"
+
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/la_parallel_block_vector.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+#include <deal.II/sundials/n_vector.h>
+
+#include "../tests.h"
+
+using namespace SUNDIALS::internal;
+
+// anonymous namespace groups helper functions for testing
+namespace
+{
+ DeclExceptionMsg(NVectorTestError,
+ "The internal N_Vector implementation didn't pass a test.");
+
+
+ IndexSet
+ create_parallel_index_set()
+ {
+ const unsigned int current_process =
+ Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+ const unsigned int n_processes =
+ Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ // give the zeroth process 2 dofs, the first 4, etc
+ const types::global_dof_index n_local_dofs = 2 * (1 + current_process);
+ const types::global_dof_index begin_dof =
+ (current_process * (current_process + 1));
+ const types::global_dof_index end_dof = begin_dof + n_local_dofs;
+
+ const types::global_dof_index n_global_dofs =
+ ((n_processes - 1) * (n_processes)) + 2 * n_processes;
+ IndexSet local_dofs(n_global_dofs);
+ local_dofs.add_range(begin_dof, end_dof);
+ local_dofs.compress();
+ AssertDimension(local_dofs.n_elements(), n_local_dofs);
+ return local_dofs;
+ }
+
+ /**
+ * Create a vector for testing and intialize all entries to @p value.
+ */
+ template <typename VectorType>
+ VectorType
+ create_test_vector(double value = 0.0);
+
+
+
+ template <>
+ Vector<double>
+ create_test_vector(double value)
+ {
+ Vector<double> vector(3 /*size*/);
+ vector = value;
+ return vector;
+ }
+
+
+
+ template <>
+ BlockVector<double>
+ create_test_vector(double value)
+ {
+ const int num_blocks = 2;
+ const int size_block = 3;
+ BlockVector<double> vector(num_blocks, size_block);
+ vector = value;
+ return vector;
+ }
+
+ template <>
+ LinearAlgebra::distributed::Vector<double>
+ create_test_vector(double value)
+ {
+ IndexSet local_dofs = create_parallel_index_set();
+ LinearAlgebra::distributed::Vector<double> vector(local_dofs,
+ MPI_COMM_WORLD);
+ vector = value;
+ return vector;
+ }
+
+ template <>
+ LinearAlgebra::distributed::BlockVector<double>
+ create_test_vector(double value)
+ {
+ const unsigned n_processes =
+ Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ const unsigned block_size =
+ ((n_processes - 1) * (n_processes)) / 2 + n_processes;
+ const auto partitioning =
+ create_parallel_index_set().split_by_block({block_size, block_size});
+ LinearAlgebra::distributed::BlockVector<double> vector(partitioning,
+ MPI_COMM_WORLD);
+ vector = value;
+ return vector;
+ }
+
+
+ template <>
+ TrilinosWrappers::MPI::Vector
+ create_test_vector(double value)
+ {
+ IndexSet local_dofs = create_parallel_index_set();
+
+ TrilinosWrappers::MPI::Vector vector(local_dofs, MPI_COMM_WORLD);
+ vector = value;
+ return vector;
+ }
+
+ template <typename VectorType>
+ bool
+ vector_equal(const VectorType &a, const VectorType &b)
+ {
+ auto elements_a = a.locally_owned_elements();
+ auto elements_b = b.locally_owned_elements();
+ if (elements_a.size() != elements_b.size())
+ return false;
+
+ const auto is_equal = [&](const IndexSet::size_type idxa,
+ const IndexSet::size_type idxb) {
+ return std::fabs(a[idxa] - b[idxb]) < 1e-15;
+ };
+
+ return std::equal(elements_a.begin(),
+ elements_a.end(),
+ elements_b.begin(),
+ is_equal);
+ }
+} // namespace
+
+
+template <typename VectorType>
+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);
+
+ Assert(n_vector != nullptr, NVectorTestError());
+ Assert((n_vector)->content != nullptr, NVectorTestError());
+
+ Assert(const_n_vector != nullptr, NVectorTestError());
+ Assert((const_n_vector)->content != nullptr, NVectorTestError());
+
+ auto *vector_unwrapped = unwrap_nvector<VectorType>(n_vector);
+ // here we check for pointer equality
+ Assert(vector_unwrapped == &vector, NVectorTestError());
+
+ // unwrap non-const as const
+ const auto *const_vector_unwrapped =
+ unwrap_nvector_const<VectorType>(n_vector);
+ Assert(const_vector_unwrapped == &vector, NVectorTestError());
+
+ // unwrap const vector as const
+ const_vector_unwrapped = unwrap_nvector_const<VectorType>(const_n_vector);
+ Assert(const_vector_unwrapped == &const_vector, NVectorTestError());
+
+ // unwrap const vector as non-const throws
+ try
+ {
+ unwrap_nvector<VectorType>(const_n_vector);
+ Assert(false, NVectorTestError());
+ }
+ catch (ExcMessage e)
+ {
+ const std::string msg(e.what());
+ Assert(msg.find(
+ "Tried to access a constant vector content as non-const.") !=
+ std::string::npos,
+ NVectorTestError());
+ }
+
+ deallog << "test_nvector_view_unwrap OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_get_vector_id()
+{
+ auto vector = create_test_vector<VectorType>();
+ NVectorView<VectorType> n_vector(vector);
+ auto id = N_VGetVectorID(n_vector);
+ Assert(id == SUNDIALS_NVEC_CUSTOM, NVectorTestError());
+ deallog << "test_get_vector_id OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_clone()
+{
+ auto vector = create_test_vector<VectorType>();
+ NVectorView<VectorType> n_vector(vector);
+ auto cloned = N_VClone(n_vector);
+
+ Assert(cloned != nullptr, NVectorTestError());
+ AssertDimension(unwrap_nvector<VectorType>(cloned)->size(), vector.size());
+
+ N_VDestroy(cloned);
+ deallog << "test_clone OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_destroy()
+{
+ GrowingVectorMemory<VectorType> mem;
+ typename GrowingVectorMemory<VectorType>::Pointer vector(mem);
+ NVectorView<VectorType> n_vector(*vector);
+
+ Assert(n_vector != nullptr, NVectorTestError());
+ auto cloned = N_VClone(n_vector);
+ Assert(cloned != nullptr, NVectorTestError());
+
+ // destroy memory-owning vector
+ N_VDestroy(cloned);
+
+ // clone empty(=without an actual vector attached) and destroy
+ cloned = N_VCloneEmpty(n_vector);
+ N_VDestroy(cloned);
+
+ // NOTE: destroying a view vector is not possible and would lead to a double
+ // delete
+ // N_VDestroy(n_vector);
+
+ // the destructor of NVectorView will do a destroy, so this is also
+ // checked here
+ deallog << "test_destroy OK" << std::endl;
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<is_serial_vector<VectorType>::value, int> = 0>
+void
+test_get_communicator()
+{
+ auto vector = create_test_vector<VectorType>();
+ NVectorView<VectorType> n_vector(vector);
+ Assert(N_VGetCommunicator(n_vector) == nullptr, NVectorTestError());
+
+ deallog << "test_get_communicator OK" << std::endl;
+}
+
+
+
+template <typename VectorType,
+ std::enable_if_t<!is_serial_vector<VectorType>::value, int> = 0>
+void
+test_get_communicator()
+{
+ auto vector = create_test_vector<VectorType>();
+ NVectorView<VectorType> n_vector(vector);
+ Assert(N_VGetCommunicator(n_vector) == MPI_COMM_WORLD, NVectorTestError());
+
+ deallog << "test_get_communicator OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_length()
+{
+ auto vector = create_test_vector<VectorType>();
+ NVectorView<VectorType> n_vector(vector);
+ Assert(N_VGetLength(n_vector) == static_cast<int>(vector.size()),
+ NVectorTestError());
+
+ deallog << "test_length OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_linear_sum()
+{
+ auto va = create_test_vector<VectorType>();
+ auto vb = create_test_vector<VectorType>();
+ 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);
+
+ va = 1.0;
+ vb = 2.0;
+
+ expected = -3.0;
+
+ // test sum into third vector
+ N_VLinearSum(1.0, nv_a, -2.0, nv_b, nv_c);
+ Assert(vector_equal(vc, expected), NVectorTestError());
+ // repeat to test that sum overwrites initial content
+ N_VLinearSum(1.0, nv_a, -2.0, nv_b, nv_c);
+ Assert(vector_equal(vc, expected), NVectorTestError());
+
+ // test store sum into one of the summands
+ N_VLinearSum(1.0, nv_a, -2.0, nv_b, nv_a);
+ Assert(vector_equal(va, expected), NVectorTestError());
+ va = 1.0;
+
+ N_VLinearSum(1.0, nv_a, -2.0, nv_b, nv_b);
+ Assert(vector_equal(vb, expected), NVectorTestError());
+
+ deallog << "test_linear_sum OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_dot_product()
+{
+ const auto va = create_test_vector<VectorType>(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 result = N_VDotProd(nv_a, nv_b);
+ auto expected = 6.0 * size;
+ Assert(std::fabs(result - expected) < 1e-12, NVectorTestError());
+
+ // test same vector
+ result = N_VDotProd(nv_a, nv_a);
+ expected = 4.0 * size;
+ Assert(std::fabs(result - expected) < 1e-12, NVectorTestError());
+
+ deallog << "test_dot_product OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_set_constant()
+{
+ auto vector = create_test_vector<VectorType>();
+ auto expected = create_test_vector<VectorType>();
+ expected = 1.0;
+
+ NVectorView<VectorType> n_vector(vector);
+
+ N_VConst(1.0, n_vector);
+ Assert(vector_equal(vector, expected), NVectorTestError());
+
+ N_VConst(2.0, n_vector);
+ expected = 2.0;
+ Assert(vector_equal(vector, expected), NVectorTestError());
+
+ deallog << "test_set_constant OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_add_constant()
+{
+ const auto vector_a = create_test_vector<VectorType>(2.0);
+ 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);
+
+ N_VAddConst(nv_a, 1.0, nv_result);
+ Assert(vector_equal(result, expected), NVectorTestError());
+
+ // test overwrite self
+ result = 2.0;
+ N_VAddConst(nv_result, 1.0, nv_result);
+ Assert(vector_equal(result, expected), NVectorTestError());
+
+ deallog << "test_add_constant OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_elementwise_product()
+{
+ const auto vector_a = create_test_vector<VectorType>(2.0);
+ const auto vector_b = create_test_vector<VectorType>(3.0);
+ 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);
+
+ // result = a.*b
+ N_VProd(nv_a, nv_b, nv_result);
+ Assert(vector_equal(vector_result, expected), NVectorTestError());
+
+ auto vector_c = create_test_vector<VectorType>(2.0);
+
+ NVectorView<VectorType> nv_c(vector_c);
+
+ // c .*= b
+ N_VProd(nv_c, nv_b, nv_c);
+ Assert(vector_equal(vector_c, expected), NVectorTestError());
+
+ // same as above with flipped arguments
+ vector_c = 2.0;
+ N_VProd(nv_b, nv_c, nv_c);
+ Assert(vector_equal(vector_c, expected), NVectorTestError());
+
+ // all operands the same
+ vector_c = 2.0;
+ expected = 4.0;
+ N_VProd(nv_c, nv_c, nv_c);
+ Assert(vector_equal(vector_c, expected), NVectorTestError());
+
+ deallog << "test_elementwise_product OK" << std::endl;
+}
+
+template <typename VectorType>
+void
+test_elementwise_div()
+{
+ const auto vector_a = create_test_vector<VectorType>(6.0);
+ const auto vector_b = create_test_vector<VectorType>(3.0);
+ 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);
+
+ N_VDiv(nv_a, nv_b, nv_result);
+ Assert(vector_equal(vector_result, expected), NVectorTestError());
+
+ deallog << "test_elementwise_div OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_elementwise_inv()
+{
+ const auto vector_a = create_test_vector<VectorType>(6.0);
+ 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);
+
+ N_VInv(nv_a, nv_result);
+ Assert(vector_equal(vector_result, expected), NVectorTestError());
+
+ // test overwrite self
+ vector_result = 6.0;
+ N_VInv(nv_result, nv_result);
+ Assert(vector_equal(vector_result, expected), NVectorTestError());
+
+ deallog << "test_elementwise_inv OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_elementwise_abs()
+{
+ const auto vector_a = create_test_vector<VectorType>(-2.0);
+ 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);
+
+ N_VAbs(nv_a, nv_result);
+ Assert(vector_equal(vector_result, expected), NVectorTestError());
+
+ // test overwrite self
+ vector_result = -2.0;
+ N_VAbs(nv_result, nv_result);
+ Assert(vector_equal(vector_result, expected), NVectorTestError());
+
+ deallog << "test_elementwise_abs OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_weighted_rms_norm()
+{
+ 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);
+
+ const auto result = N_VWrmsNorm(nv_a, nv_b);
+ Assert(std::fabs(result - 6.0) < 1e-12, NVectorTestError());
+
+ deallog << "test_weighted_rms_norm OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_max_norm()
+{
+ 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 result = N_VMaxNorm(nv_a);
+ Assert(std::fabs(result - 2.0) < 1e-12, NVectorTestError());
+
+ result = N_VMaxNorm(nv_b);
+ Assert(std::fabs(result - 3.0) < 1e-12, NVectorTestError());
+
+ deallog << "test_max_norm OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+test_scale()
+{
+ const auto vector_a = create_test_vector<VectorType>(2.0);
+ 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);
+
+ // b *= 2
+ N_VScale(2.0, nv_b, nv_b);
+ Assert(vector_equal(vector_b, expected), NVectorTestError());
+
+ // b = 2*a
+ N_VScale(2.0, nv_a, nv_b);
+ Assert(vector_equal(vector_b, expected), NVectorTestError());
+
+ deallog << "test_scale OK" << std::endl;
+}
+
+
+
+template <typename VectorType>
+void
+run_all_tests(const std::string &prefix)
+{
+ LogStream::Prefix p(prefix);
+ // test conversion between vectors
+ test_nvector_view_unwrap<VectorType>();
+ test_get_vector_id<VectorType>();
+
+ // test vector operations
+ test_clone<VectorType>();
+ test_destroy<VectorType>();
+ test_get_communicator<VectorType>();
+ test_length<VectorType>();
+ test_linear_sum<VectorType>();
+ test_dot_product<VectorType>();
+ test_set_constant<VectorType>();
+ test_add_constant<VectorType>();
+ test_elementwise_product<VectorType>();
+ test_elementwise_div<VectorType>();
+ test_elementwise_inv<VectorType>();
+ test_elementwise_abs<VectorType>();
+ test_weighted_rms_norm<VectorType>();
+ test_max_norm<VectorType>();
+ test_scale<VectorType>();
+}
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ MPILogInitAll log_all;
+
+ using VectorType = Vector<double>;
+
+ run_all_tests<Vector<double>>("Vector<double>");
+ run_all_tests<BlockVector<double>>("BlockVector<double>");
+ run_all_tests<LinearAlgebra::distributed::Vector<double>>(
+ "LinearAlgebra::distributed::Vector<double>");
+ run_all_tests<LinearAlgebra::distributed::BlockVector<double>>(
+ "LinearAlgebra::distributed::BlockVector<double>");
+ run_all_tests<TrilinosWrappers::MPI::Vector>("TrilinosWrappers::MPI::Vector");
+}
--- /dev/null
+
+DEAL:0:Vector<double>::test_nvector_view_unwrap OK
+DEAL:0:Vector<double>::test_get_vector_id OK
+DEAL:0:Vector<double>::test_clone OK
+DEAL:0:Vector<double>::test_destroy OK
+DEAL:0:Vector<double>::test_get_communicator OK
+DEAL:0:Vector<double>::test_length OK
+DEAL:0:Vector<double>::test_linear_sum OK
+DEAL:0:Vector<double>::test_dot_product OK
+DEAL:0:Vector<double>::test_set_constant OK
+DEAL:0:Vector<double>::test_add_constant OK
+DEAL:0:Vector<double>::test_elementwise_product OK
+DEAL:0:Vector<double>::test_elementwise_div OK
+DEAL:0:Vector<double>::test_elementwise_inv OK
+DEAL:0:Vector<double>::test_elementwise_abs OK
+DEAL:0:Vector<double>::test_weighted_rms_norm OK
+DEAL:0:Vector<double>::test_max_norm OK
+DEAL:0:Vector<double>::test_scale OK
+DEAL:0:BlockVector<double>::test_nvector_view_unwrap OK
+DEAL:0:BlockVector<double>::test_get_vector_id OK
+DEAL:0:BlockVector<double>::test_clone OK
+DEAL:0:BlockVector<double>::test_destroy OK
+DEAL:0:BlockVector<double>::test_get_communicator OK
+DEAL:0:BlockVector<double>::test_length OK
+DEAL:0:BlockVector<double>::test_linear_sum OK
+DEAL:0:BlockVector<double>::test_dot_product OK
+DEAL:0:BlockVector<double>::test_set_constant OK
+DEAL:0:BlockVector<double>::test_add_constant OK
+DEAL:0:BlockVector<double>::test_elementwise_product OK
+DEAL:0:BlockVector<double>::test_elementwise_div OK
+DEAL:0:BlockVector<double>::test_elementwise_inv OK
+DEAL:0:BlockVector<double>::test_elementwise_abs OK
+DEAL:0:BlockVector<double>::test_weighted_rms_norm OK
+DEAL:0:BlockVector<double>::test_max_norm OK
+DEAL:0:BlockVector<double>::test_scale OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_nvector_view_unwrap OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_get_vector_id OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_clone OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_destroy OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_get_communicator OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_length OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_linear_sum OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_dot_product OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_set_constant OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_add_constant OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_product OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_div OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_inv OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_abs OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_weighted_rms_norm OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_max_norm OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_scale OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_nvector_view_unwrap OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_get_vector_id OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_clone OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_destroy OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_get_communicator OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_length OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_linear_sum OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_dot_product OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_set_constant OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_add_constant OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_product OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_div OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_inv OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_abs OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_weighted_rms_norm OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_max_norm OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_scale OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_nvector_view_unwrap OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_get_vector_id OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_clone OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_destroy OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_get_communicator OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_length OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_linear_sum OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_dot_product OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_set_constant OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_add_constant OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_product OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_div OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_inv OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_abs OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_weighted_rms_norm OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_max_norm OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_scale OK
--- /dev/null
+
+DEAL:0:Vector<double>::test_nvector_view_unwrap OK
+DEAL:0:Vector<double>::test_get_vector_id OK
+DEAL:0:Vector<double>::test_clone OK
+DEAL:0:Vector<double>::test_destroy OK
+DEAL:0:Vector<double>::test_get_communicator OK
+DEAL:0:Vector<double>::test_length OK
+DEAL:0:Vector<double>::test_linear_sum OK
+DEAL:0:Vector<double>::test_dot_product OK
+DEAL:0:Vector<double>::test_set_constant OK
+DEAL:0:Vector<double>::test_add_constant OK
+DEAL:0:Vector<double>::test_elementwise_product OK
+DEAL:0:Vector<double>::test_elementwise_div OK
+DEAL:0:Vector<double>::test_elementwise_inv OK
+DEAL:0:Vector<double>::test_elementwise_abs OK
+DEAL:0:Vector<double>::test_weighted_rms_norm OK
+DEAL:0:Vector<double>::test_max_norm OK
+DEAL:0:Vector<double>::test_scale OK
+DEAL:0:BlockVector<double>::test_nvector_view_unwrap OK
+DEAL:0:BlockVector<double>::test_get_vector_id OK
+DEAL:0:BlockVector<double>::test_clone OK
+DEAL:0:BlockVector<double>::test_destroy OK
+DEAL:0:BlockVector<double>::test_get_communicator OK
+DEAL:0:BlockVector<double>::test_length OK
+DEAL:0:BlockVector<double>::test_linear_sum OK
+DEAL:0:BlockVector<double>::test_dot_product OK
+DEAL:0:BlockVector<double>::test_set_constant OK
+DEAL:0:BlockVector<double>::test_add_constant OK
+DEAL:0:BlockVector<double>::test_elementwise_product OK
+DEAL:0:BlockVector<double>::test_elementwise_div OK
+DEAL:0:BlockVector<double>::test_elementwise_inv OK
+DEAL:0:BlockVector<double>::test_elementwise_abs OK
+DEAL:0:BlockVector<double>::test_weighted_rms_norm OK
+DEAL:0:BlockVector<double>::test_max_norm OK
+DEAL:0:BlockVector<double>::test_scale OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_nvector_view_unwrap OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_get_vector_id OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_clone OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_destroy OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_get_communicator OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_length OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_linear_sum OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_dot_product OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_set_constant OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_add_constant OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_product OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_div OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_inv OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_elementwise_abs OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_weighted_rms_norm OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_max_norm OK
+DEAL:0:LinearAlgebra::distributed::Vector<double>::test_scale OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_nvector_view_unwrap OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_get_vector_id OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_clone OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_destroy OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_get_communicator OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_length OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_linear_sum OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_dot_product OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_set_constant OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_add_constant OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_product OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_div OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_inv OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_abs OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_weighted_rms_norm OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_max_norm OK
+DEAL:0:LinearAlgebra::distributed::BlockVector<double>::test_scale OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_nvector_view_unwrap OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_get_vector_id OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_clone OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_destroy OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_get_communicator OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_length OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_linear_sum OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_dot_product OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_set_constant OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_add_constant OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_product OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_div OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_inv OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_elementwise_abs OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_weighted_rms_norm OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_max_norm OK
+DEAL:0:TrilinosWrappers::MPI::Vector::test_scale OK
+
+DEAL:1:Vector<double>::test_nvector_view_unwrap OK
+DEAL:1:Vector<double>::test_get_vector_id OK
+DEAL:1:Vector<double>::test_clone OK
+DEAL:1:Vector<double>::test_destroy OK
+DEAL:1:Vector<double>::test_get_communicator OK
+DEAL:1:Vector<double>::test_length OK
+DEAL:1:Vector<double>::test_linear_sum OK
+DEAL:1:Vector<double>::test_dot_product OK
+DEAL:1:Vector<double>::test_set_constant OK
+DEAL:1:Vector<double>::test_add_constant OK
+DEAL:1:Vector<double>::test_elementwise_product OK
+DEAL:1:Vector<double>::test_elementwise_div OK
+DEAL:1:Vector<double>::test_elementwise_inv OK
+DEAL:1:Vector<double>::test_elementwise_abs OK
+DEAL:1:Vector<double>::test_weighted_rms_norm OK
+DEAL:1:Vector<double>::test_max_norm OK
+DEAL:1:Vector<double>::test_scale OK
+DEAL:1:BlockVector<double>::test_nvector_view_unwrap OK
+DEAL:1:BlockVector<double>::test_get_vector_id OK
+DEAL:1:BlockVector<double>::test_clone OK
+DEAL:1:BlockVector<double>::test_destroy OK
+DEAL:1:BlockVector<double>::test_get_communicator OK
+DEAL:1:BlockVector<double>::test_length OK
+DEAL:1:BlockVector<double>::test_linear_sum OK
+DEAL:1:BlockVector<double>::test_dot_product OK
+DEAL:1:BlockVector<double>::test_set_constant OK
+DEAL:1:BlockVector<double>::test_add_constant OK
+DEAL:1:BlockVector<double>::test_elementwise_product OK
+DEAL:1:BlockVector<double>::test_elementwise_div OK
+DEAL:1:BlockVector<double>::test_elementwise_inv OK
+DEAL:1:BlockVector<double>::test_elementwise_abs OK
+DEAL:1:BlockVector<double>::test_weighted_rms_norm OK
+DEAL:1:BlockVector<double>::test_max_norm OK
+DEAL:1:BlockVector<double>::test_scale OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_nvector_view_unwrap OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_get_vector_id OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_clone OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_destroy OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_get_communicator OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_length OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_linear_sum OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_dot_product OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_set_constant OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_add_constant OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_elementwise_product OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_elementwise_div OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_elementwise_inv OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_elementwise_abs OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_weighted_rms_norm OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_max_norm OK
+DEAL:1:LinearAlgebra::distributed::Vector<double>::test_scale OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_nvector_view_unwrap OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_get_vector_id OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_clone OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_destroy OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_get_communicator OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_length OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_linear_sum OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_dot_product OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_set_constant OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_add_constant OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_product OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_div OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_inv OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_elementwise_abs OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_weighted_rms_norm OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_max_norm OK
+DEAL:1:LinearAlgebra::distributed::BlockVector<double>::test_scale OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_nvector_view_unwrap OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_get_vector_id OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_clone OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_destroy OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_get_communicator OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_length OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_linear_sum OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_dot_product OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_set_constant OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_add_constant OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_elementwise_product OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_elementwise_div OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_elementwise_inv OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_elementwise_abs OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_weighted_rms_norm OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_max_norm OK
+DEAL:1:TrilinosWrappers::MPI::Vector::test_scale OK
+
projects / dealii.git / commitdiff