Skip to content
Snippets Groups Projects
Commit 58ce5c15 authored by Lasse Hinrichsen's avatar Lasse Hinrichsen
Browse files

Add inplaceBinary method

With this tools, one can write quite general code to apply
x = op(x,y)
entrywise for vectors x,y.
parent 1411c107
Branches
No related tags found
1 merge request!9[WIP] Add applyEntrywise method
Pipeline #30326 failed
...@@ -30,6 +30,16 @@ struct HasSize { ...@@ -30,6 +30,16 @@ struct HasSize {
auto require(C&& c) -> decltype(c.size()); auto require(C&& c) -> decltype(c.size());
}; };
struct IsLegalBinaryOperator {
template <class C, class D, class Op>
auto require(C&& c, D&& d, Op&& op) -> decltype(op(c,d));
};
template <class C, class D, class Op>
constexpr bool isLegalBinaryOperator(C&& c, D&& d, Op&& op) {
return models<IsLegalBinaryOperator, C, D, Op>();
}
} // end namespace Concept } // end namespace Concept
} // end namespace MatrixVector } // end namespace MatrixVector
} // end namespace Dune } // end namespace Dune
......
...@@ -13,6 +13,7 @@ ...@@ -13,6 +13,7 @@
#include <dune/common/hybridutilities.hh> #include <dune/common/hybridutilities.hh>
#include <dune/matrix-vector/algorithm.hh> #include <dune/matrix-vector/algorithm.hh>
#include <dune/matrix-vector/concepts.hh>
#include <dune/matrix-vector/traits/utilities.hh> #include <dune/matrix-vector/traits/utilities.hh>
//! \brief Various tools for working with istl vectors of arbitrary nesting depth //! \brief Various tools for working with istl vectors of arbitrary nesting depth
...@@ -44,6 +45,25 @@ void truncate(Vector& v, const BitVector& tr) { ...@@ -44,6 +45,25 @@ void truncate(Vector& v, const BitVector& tr) {
Impl::ScalarSwitch<Vector>::truncate(v, tr); Impl::ScalarSwitch<Vector>::truncate(v, tr);
} }
/**! For a given operator "op", compute x = op(x,y) ENTRYWISE.
*
* If the op is not defined for the given types,
* the operator is recursively applied to subblocks
* until it can be used.
*/
template<typename X, typename Y, typename BinaryOp>
void inplaceBinary(X&& x, Y&& y, BinaryOp&& op) {
if constexpr(Concept::isLegalBinaryOperator(x,y,op)) {
x = op(x,y);
}
else {
for (std::size_t i = 0; i < x.size(); ++i) {
inplaceBinary(x[i], y[i], op);
}
}
}
namespace Impl { namespace Impl {
//! recursion helper for scalars nested in iterables (vectors) //! recursion helper for scalars nested in iterables (vectors)
......
dune_add_test(SOURCES arithmetictest.cc) dune_add_test(SOURCES arithmetictest.cc)
dune_add_test(SOURCES genericvectortoolstest.cc)
dune_add_test(SOURCES resizetest.cc) dune_add_test(SOURCES resizetest.cc)
dune_add_test(SOURCES staticmatrixtoolstest.cc) dune_add_test(SOURCES staticmatrixtoolstest.cc)
dune_add_test(SOURCES triangularsolvetest.cc) dune_add_test(SOURCES triangularsolvetest.cc)
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <cstddef>
#include <vector>
#include <dune/common/bitsetvector.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <dune/matrix-vector/genericvectortools.hh>
using namespace Dune;
//template<typename X, typename Y>
//bool my_or(const X& x, const Y& y) {
//return x || y;
//}
TestSuite test_inplace_or() {
TestSuite suite;
using namespace Dune::MatrixVector::Generic;
// This would be the STL approach
//auto inplace_or = [](auto& x, const auto& y) {
//inplaceBinary(x,y, std::logical_or<void>());
//};
// For demonstration, let'ss cook something up.
// Note that we cant use "auto" for the types of
// xx and yy, because we explicitly check if the given input is
// bool (or something convertible) and in that case apply
// the operator
auto my_or = [](bool xx, bool yy) { return xx || yy; };
auto inplace_or = [&](auto& x, const auto& y) {
inplaceBinary(x,y, my_or);
};
// check bools (trivial case)
{
bool t = true;
bool f = false;
inplace_or(t,f); // t |= f
suite.check(t == true);
inplace_or(f,t); // f |= t
suite.check(f == true);
f=false;
inplace_or(f,f);
suite.check(f == false);
inplace_or(t,t);
suite.check(t == true);
}
// check BitSetVector
constexpr const int N = 3;
using BSV = BitSetVector<N>;
{
BSV x(2, false);
BSV y(2, true);
inplace_or(x,y);
for (std::size_t i = 0; i < x.size(); ++i) {
suite.check(x[i].all() == true);
}
}
// check vector<BSV>
{
auto x = std::vector<BSV>(4);
auto y = std::vector<BSV>(4);
for (auto&& xx : x)
xx = BSV(2, false);
for (auto&& yy : y)
yy = BSV(2, true);
inplace_or(x, y);
for (const auto& bsv : x) {
for (const auto& bs: bsv) {
suite.check(bs.all());
}
}
}
// mix types
{
auto x = std::vector<bool>(10, false);
auto y = std::vector<char>(10, true);
inplace_or(x,y);
for(const auto& xi : x)
suite.check(xi == true);
}
return suite;
}
TestSuite test_inplace_plus() {
TestSuite suite;
using namespace Dune::MatrixVector::Generic;
auto my_plus = std::plus<double>{};
auto inplace_plus = [&](auto& x, const auto& y) { inplaceBinary(x, y, my_plus); };
// test double
{
// both values should be exact in binary representation
auto x = 1.5;
auto y = 2.5;
inplace_plus(x,y); // x += y
suite.check(x== 4., "inplace_plus failed for double");
}
// test vector<double>
{
auto x = std::vector<double>(10, 1.5);
auto y = std::vector<double>(10, 2.5);
inplace_plus(x,y); // x += y
for (const auto& entry : x)
suite.check(entry == 4.0, "inplace_plus failed for vector<double> entry.");
}
// test BlockVector<FieldVector<double, 2>>
{
auto x = BlockVector<FieldVector<double, 2>>(10);
auto y = BlockVector<FieldVector<double, 2>>(10);
x = 1.5;
y = 2.5;
inplace_plus(x,y); // x += y
for (const auto& block : x)
for (const auto& entry : block)
suite.check(entry == 4.0, "inplace_plus failed for BlockVector<FieldVector<2>> entry.");
}
return suite;
}
TestSuite test_inplaceBinary() {
TestSuite suite;
// we test two example operators, the || and the + operator.
suite.subTest(test_inplace_or());
suite.subTest(test_inplace_plus());
return suite;
}
int main(int argc, char* argv[]) {
MPIHelper::instance(argc, argv);
TestSuite suite;
// TODO: There are more things in genericvectortools
// that need tests!
suite.subTest(test_inplaceBinary());
return suite.exit();;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment