diff --git a/dune/solvers/iterationsteps/blockgsstep.cc b/dune/solvers/iterationsteps/blockgsstep.cc index 1cb6959a7d807117af993b2ba03448b519cb4b7e..284351e051ae919ca2775d9164f715bf7ffaf6c4 100644 --- a/dune/solvers/iterationsteps/blockgsstep.cc +++ b/dune/solvers/iterationsteps/blockgsstep.cc @@ -3,25 +3,19 @@ #include <cassert> +#include <dune/solvers/common/arithmetic.hh> + template<class MatrixType, class DiscFuncType, class BitVectorType> inline void BlockGSStep<MatrixType, DiscFuncType, BitVectorType>:: residual(int index, VectorBlock& r) const { - const MatrixType& mat = *this->mat_; - const auto& row = mat[index]; - r = (*this->rhs_)[index]; - /* The following loop subtracts - * \f[ sum_i = \sum_j A_{ij}w_j \f] - */ - auto cIt = row.begin(); - auto cEndIt = row.end(); - - for (; cIt!=cEndIt; ++cIt) { + const auto& row = (*this->mat_)[index]; + for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) { // r_i -= A_ij x_j - cIt->mmv((*this->x_)[cIt.index()], r); + Arithmetic::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]); } }