Add componentwise function application and componentwise product for

vectors

dune/matrix-vector/CMakeLists.txt

@@ -11,6 +11,7 @@ install(FILES
   axy.hh
   blockmatrixview.hh
   componentwisematrixmap.hh
+  componentwisevectormaps.hh
   concepts.hh
   crossproduct.hh
   genericvectortools.hh

dune/matrix-vector/componentwisevectormaps.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2
+#ifndef DUNE_MATRIX_VECTOR_COMPONENTWISE_VECTOR_MAPS
+#define DUNE_MATRIX_VECTOR_COMPONENTWISE_VECTOR_MAPS
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/common/iteratorfacades.hh>
+
+namespace Dune {
+namespace MatrixVector {
+
+  /** Apply a function f componentwise to a vector v and store it in u, i.e.
+   *  u_i = f(v_i)
+   */
+  template<typename F, typename V>
+  void applyComponentwise(F&& f, V& u, const V& v) {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<std::decay_t<V>>(), [&](auto&& id) {
+        u= f(id(v));},
+        [&](auto id) {
+        H::forEach(H::integralRange(H::size(id(v))), [&](auto i) {
+            applyComponentwise(f,H::elementAt(u,i), H::elementAt(v,i));
+            });
+        });
+  }
+
+  /** Multiply the i-th component of u by the i-th component of v, i.e.
+   *  u_i <- u_i*v_i
+   */
+  template<typename V>
+  void multiplyComponentwise(V& u, const V& v) {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<std::decay_t<V>>(), [&](auto&& id) {
+        u*=id(v);},
+        [&](auto id) {
+        H::forEach(H::integralRange(H::size(id(v))), [&](auto i) {
+            multiplyComponentwise(H::elementAt(u,i), H::elementAt(v,i));
+            });
+        });
+  }
+
+}
+}
+#endif

dune/matrix-vector/test/CMakeLists.txt

 dune_add_test(SOURCES arithmetictest.cc)
+dune_add_test(SOURCES componentwisevectortest.cc)
 dune_add_test(SOURCES resizetest.cc)
 dune_add_test(SOURCES staticmatrixtoolstest.cc)
 dune_add_test(SOURCES triangularsolvetest.cc)

dune/matrix-vector/test/componentwisevectortest.cc

+#include <config.h>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/common/indices.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/vbvector.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+#include <dune/matrix-vector/componentwisevectormaps.hh>
+
+#include <cmath>
+#include <string>
+
+template<typename V>
+auto test_vector(V v, std::string type) {
+  using namespace Dune;
+  auto suite = TestSuite("test componentwise on vector");
+
+  auto u = v;
+
+  // test function application
+  auto f = [](const auto& in) {
+    return in*2.;
+  };
+
+  MatrixVector::applyComponentwise(f, u, v);
+
+  // now, u_i should be f(v_i) = 2*v_i.
+  // we test this by computing the sum
+  // over both vectors.
+  // It should hold sum(u) == 2.*sum(v) (modulo floating point arithmetic)
+  auto ones =v;
+  ones = 1.;
+
+  auto v_sum = ones*v;
+
+  suite.check(std::abs(ones*u - 2.*v_sum) < 1e-15, "Test applyComponentwise") << "Failed for type " << type;
+
+  // test componentwise multiply
+  u = 2.; // reset u's state
+
+  MatrixVector::multiplyComponentwise(u, v); // u <- u.*v
+
+  // test again sums
+  suite.check(std::abs(ones*u - 2.*v_sum) < 1e-15, "Test multiplyComponentwise") << "Failed for type " << type;
+
+  return suite;
+
+}
+
+int main(int argc, char** argv) {
+
+  auto& helper = Dune::MPIHelper::instance(argc, argv);
+
+  auto suite = Dune::TestSuite("Test componentwise vector maps");
+
+  const auto n = 10;
+
+  // Test BlockVector<FieldVector<double, 1>>
+  {
+    auto vector = Dune::BlockVector<Dune::FieldVector<double, 1>>(n);
+    vector = 0.5;
+    suite.subTest(test_vector(std::move(vector), "BlockVector<FieldVector<double, 1>>"));
+  }
+
+  // Test BlockVector<FieldVector<double, 2>>
+  {
+    auto vector = Dune::BlockVector<Dune::FieldVector<double, 2>>(n);
+    vector = 0.5;
+    suite.subTest(test_vector(std::move(vector), "BlockVector<FieldVector<double, 2>>"));
+  }
+
+  // Test BlockVector<double>>
+  {
+    auto vector = Dune::BlockVector<double>(n);
+    vector = 0.5;
+    suite.subTest(test_vector(std::move(vector), "BlockVector<double>"));
+  }
+
+  // Test VariableBlockVector<double>
+  {
+    auto vector = Dune::VariableBlockVector<double>(n);
+    int i = 1;
+    for(auto it = vector.createbegin(); it != vector.createend(); ++it) {
+      *it = i;
+      i = i%2 +1; // non homogenous block structure
+    }
+    vector = 0.5;
+    suite.subTest(test_vector(std::move(vector), "VariableBlockVector<double>"));
+  }
+
+  // Test MultiTypeBlockVector
+  {
+    using V1 = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+    using V2 = Dune::BlockVector<Dune::FieldVector<double, 2>>;
+
+    using MT = Dune::MultiTypeBlockVector<V1, V2>;
+    MT mt{};
+    using namespace Dune::Indices;
+    mt[_0].resize(n);
+    mt[_1].resize(n);
+    mt = 0.5;
+    suite.subTest(test_vector(std::move(mt), "MultiTypeBlockVector"));
+  }
+
+  return suite.exit();
+}