Skip to content
Snippets Groups Projects
Commit 64fdd4e7 authored by Max Kahnt's avatar Max Kahnt
Browse files

Use hybrid algorithms for addProduct operation.

parent b18ba9cc
No related branches found
No related tags found
1 merge request!3Feature/multitype arithmetic
......@@ -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);
});
});
}
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment