]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce SUNDIALS N_Vector module
authorSebastian Proell <proell@lnm.mw.tum.de>
Fri, 18 Dec 2020 12:11:37 +0000 (13:11 +0100)
committerSebastian Proell <proell@lnm.mw.tum.de>
Thu, 7 Jan 2021 08:54:28 +0000 (09:54 +0100)
include/deal.II/sundials/n_vector.h [new file with mode: 0644]
include/deal.II/sundials/n_vector.templates.h [new file with mode: 0644]
source/sundials/CMakeLists.txt
source/sundials/n_vector.cc [new file with mode: 0644]
source/sundials/n_vector.inst.in [new file with mode: 0644]
tests/sundials/n_vector.cc [new file with mode: 0644]
tests/sundials/n_vector.with_sundials.geq.5.0.0.output [new file with mode: 0644]
tests/sundials/n_vector.with_sundials.geq.5.0.0.with_mpi=true.with_trilinos=true.mpirun=2.output [new file with mode: 0644]

diff --git a/include/deal.II/sundials/n_vector.h b/include/deal.II/sundials/n_vector.h
new file mode 100644 (file)
index 0000000..d583924
--- /dev/null
@@ -0,0 +1,120 @@
+//-----------------------------------------------------------
+//
+//    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
diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h
new file mode 100644 (file)
index 0000000..3d58334
--- /dev/null
@@ -0,0 +1,809 @@
+//-----------------------------------------------------------
+//
+//    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
index e7dca9c723cb898dbb4a0aa6bca7cdfbe0515521..d01e607689f630a4d7d6a53cc9d42b0b7178b872 100644 (file)
@@ -20,9 +20,11 @@ SET(_src
   ida.cc
   copy.cc
   kinsol.cc
+  n_vector.cc
   )
 
 SET(_inst
+  n_vector.inst.in
   )
 
 FILE(GLOB _header
diff --git a/source/sundials/n_vector.cc b/source/sundials/n_vector.cc
new file mode 100644 (file)
index 0000000..cdcb320
--- /dev/null
@@ -0,0 +1,32 @@
+//-----------------------------------------------------------
+//
+//    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
diff --git a/source/sundials/n_vector.inst.in b/source/sundials/n_vector.inst.in
new file mode 100644 (file)
index 0000000..f8a6b6e
--- /dev/null
@@ -0,0 +1,46 @@
+// ---------------------------------------------------------------------
+//
+// 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);
+  }
diff --git a/tests/sundials/n_vector.cc b/tests/sundials/n_vector.cc
new file mode 100644 (file)
index 0000000..06cb512
--- /dev/null
@@ -0,0 +1,625 @@
+//-----------------------------------------------------------
+//
+//    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");
+}
diff --git a/tests/sundials/n_vector.with_sundials.geq.5.0.0.output b/tests/sundials/n_vector.with_sundials.geq.5.0.0.output
new file mode 100644 (file)
index 0000000..0d53a6a
--- /dev/null
@@ -0,0 +1,86 @@
+
+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
diff --git a/tests/sundials/n_vector.with_sundials.geq.5.0.0.with_mpi=true.with_trilinos=true.mpirun=2.output b/tests/sundials/n_vector.with_sundials.geq.5.0.0.with_mpi=true.with_trilinos=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..f57dfff
--- /dev/null
@@ -0,0 +1,173 @@
+
+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
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.