From acdb69a9f1bcd2763d178609899ece9d6f739dbc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Carsten=20Gr=C3=A4ser?= <graeser@dune-project.org>
Date: Tue, 11 Jul 2017 16:55:48 +0200
Subject: [PATCH] Implement axy using sparseRangeFor()
This way it works for multi-type matrices again.
This helps to fix the dune-tnnmg test which currently
breaks.
Notice that this may be slightly less efficient, because
we use a new temporary for each loop instance. This is
necessary because we may need a different type each time.
This may be fixed by adding a specialization for the
case that the corresponding vector is a classic dynamic
container.
---
dune/matrix-vector/axy.hh | 64 ++++++++++++++++++++-------------------
1 file changed, 33 insertions(+), 31 deletions(-)
diff --git a/dune/matrix-vector/axy.hh b/dune/matrix-vector/axy.hh
index c378c89..c626b12 100644
--- a/dune/matrix-vector/axy.hh
+++ b/dune/matrix-vector/axy.hh
@@ -5,6 +5,7 @@
#include "axpy.hh"
#include "matrixtraits.hh"
+#include "algorithm.hh"
namespace Dune {
namespace MatrixVector {
@@ -40,22 +41,23 @@ namespace MatrixVector {
static typename VectorType::field_type Axy(const MatrixType &A,
const VectorType &x,
const VectorType2 &y) {
- assert(x.N() == A.M());
- assert(y.N() == A.N());
+ assert(x.size() == A.M());
+ assert(y.size() == A.N());
- typename VectorType::field_type outer = 0.0;
- typename VectorType2::block_type inner;
- typename MatrixType::ConstIterator endi = A.end();
- for (typename MatrixType::ConstIterator i = A.begin(); i != endi; ++i) {
- inner = 0.0;
- const size_t iindex = i.index();
- typename MatrixType::row_type::ConstIterator endj = i->end();
- for (typename MatrixType::row_type::ConstIterator j = i->begin();
- j != endj; ++j)
- addProduct(inner, *j, x[j.index()]);
- outer += inner * y[iindex];
- }
- return outer;
+ using Result = std::decay_t<decltype(y*y)>;
+
+ Result result = 0;
+ sparseRangeFor(A, [&](auto&& Ai, auto i) {
+ // ToDo: Provide a specialization using a single temporary
+ // for the case that each y[i] has the same type.
+ auto Aix = y[i];
+ Aix = 0;
+ sparseRangeFor(Ai, [&](auto&& Aij, auto j) {
+ addProduct(Aix, Aij, x[j]);
+ });
+ result += Aix * y[i];
+ });
+ return result;
}
template <class VectorType, class VectorType2>
@@ -63,23 +65,23 @@ namespace MatrixVector {
const VectorType2 &b,
const VectorType &x,
const VectorType2 &y) {
- assert(x.N() == A.M());
- assert(y.N() == A.N());
- assert(y.N() == b.N());
+ assert(x.size() == A.M());
+ assert(y.size() == A.N());
+ assert(y.size() == b.size());
+
+ using Result = std::decay_t<decltype(y*y)>;
- typename VectorType::field_type outer = 0.0;
- typename VectorType2::block_type inner;
- typename MatrixType::ConstIterator endi = A.end();
- for (typename MatrixType::ConstIterator i = A.begin(); i != endi; ++i) {
- const size_t iindex = i.index();
- inner = b[iindex];
- typename MatrixType::row_type::ConstIterator endj = i->end();
- for (typename MatrixType::row_type::ConstIterator j = i->begin();
- j != endj; ++j)
- subtractProduct(inner, *j, x[j.index()]);
- outer += inner * y[iindex];
- }
- return outer;
+ Result result = 0;
+ sparseRangeFor(A, [&](auto&& Ai, auto i) {
+ // ToDo: Provide a specialization using a single temporary
+ // for the case that each b[i] has the same type.
+ auto Aix = b[i];
+ sparseRangeFor(Ai, [&](auto&& Aij, auto j) {
+ subtractProduct(Aix, Aij, x[j]);
+ });
+ result += Aix * y[i];
+ });
+ return result;
}
};
--
GitLab