From 64fdd4e7943a1be4b24b9944816f0a696dddf2dd Mon Sep 17 00:00:00 2001 From: Max Kahnt <max.kahnt@fu-berlin.de> Date: Fri, 12 May 2017 14:46:20 +0200 Subject: [PATCH] Use hybrid algorithms for addProduct operation. --- dune/matrix-vector/axpy.hh | 67 ++++++++++++++------------------------ 1 file changed, 25 insertions(+), 42 deletions(-) diff --git a/dune/matrix-vector/axpy.hh b/dune/matrix-vector/axpy.hh index 792dc03..3a5e130 100644 --- a/dune/matrix-vector/axpy.hh +++ b/dune/matrix-vector/axpy.hh @@ -8,6 +8,7 @@ #include <dune/istl/bcrsmatrix.hh> #include <dune/istl/scaledidmatrix.hh> +#include "algorithm.hh" #include "matrixtraits.hh" #include "scalartraits.hh" @@ -30,18 +31,13 @@ namespace MatrixVector { class ADummy = A, std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0> static void addProduct(A& a, const B& b, const C& c) { - typename B::ConstRowIterator bi = b.begin(); - typename B::ConstRowIterator bEnd = b.end(); - for (; bi != bEnd; ++bi) { - typename B::ConstColIterator bik = bi->begin(); - typename B::ConstColIterator biEnd = bi->end(); - for (; bik != biEnd; ++bik) { - typename C::ConstColIterator ckj = c[bik.index()].begin(); - typename C::ConstColIterator ckEnd = c[bik.index()].end(); - for (; ckj != ckEnd; ++ckj) - a[bi.index()][ckj.index()] += (*bik) * (*ckj); - } - } + rangeForEach(b, [&](auto&& bi, auto&& i) { + rangeForEach(bi, [&](auto&& bik, auto&& k) { + rangeForEach(c[k], [&](auto&& ckj, auto&& j) { + Dune::MatrixVector::addProduct(a[i][j], bik, ckj); + }); + }); + }); } }; @@ -61,18 +57,13 @@ namespace MatrixVector { class ADummy = A, std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0> static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) { - typename B::ConstRowIterator bi = b.begin(); - typename B::ConstRowIterator bEnd = b.end(); - for (; bi != bEnd; ++bi) { - typename B::ConstColIterator bik = b[bi.index()].begin(); - typename B::ConstColIterator biEnd = b[bi.index()].end(); - for (; bik != biEnd; ++bik) { - typename C::ConstColIterator ckj = c[bik.index()].begin(); - typename C::ConstColIterator ckEnd = c[bik.index()].end(); - for (; ckj != ckEnd; ++ckj) - a[bi.index()][ckj.index()] += scalar * (*bik) * (*ckj); - } - } + rangeForEach(b, [&](auto&& bi, auto&& i) { + rangeForEach(bi, [&](auto&& bik, auto&& k) { + rangeForEach(c[k], [&](auto&& ckj, auto&& j) { + Dune::MatrixVector::addProduct(a[i][j], scalar, bik, ckj); + }); + }); + }); } }; @@ -257,15 +248,11 @@ namespace MatrixVector { class ADummy = A, std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0> static void addProduct(A& a, const B& b, const C& c) { - typename C::ConstRowIterator ci = c.begin(); - typename C::ConstRowIterator cEnd = c.end(); - for (; ci != cEnd; ++ci) { - typename C::ConstColIterator cik = c[ci.index()].begin(); - typename C::ConstColIterator ciEnd = c[ci.index()].end(); - for (; cik != ciEnd; ++cik) { - a[ci.index()][cik.index()] += b * (*cik); - } - } + rangeForEach(c, [&](auto&& ci, auto && i) { + rangeForEach(ci, [&](auto&& cij, auto && j) { + Dune::MatrixVector::addProduct(a[i][j], b, cij); + }); + }); } }; @@ -286,15 +273,11 @@ namespace MatrixVector { class ADummy = A, std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0> static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) { - typename C::ConstRowIterator ci = c.begin(); - typename C::ConstRowIterator cEnd = c.end(); - for (; ci != cEnd; ++ci) { - typename C::ConstColIterator cik = c[ci.index()].begin(); - typename C::ConstColIterator ciEnd = c[ci.index()].end(); - for (; cik != ciEnd; ++cik) { - a[ci.index()][cik.index()] += scalar * b * (*cik); - } - } + rangeForEach(c, [&](auto&& ci, auto&& i) { + rangeForEach(ci, [&](auto&& cij, auto&& j) { + Dune::MatrixVector::addProduct(a[i][j], scalar, b, cij); + }); + }); } }; -- GitLab