diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh index 707f6aed8e4892749cc5c6271d957484f5c2dfab..cb4bf43af9999bfb6d9885991a003a1e59236331 100644 --- a/dune/tectonic/myblockproblem.hh +++ b/dune/tectonic/myblockproblem.hh @@ -104,7 +104,7 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem { v /= vnorm; // Rescale for numerical stability VectorType tmp = problem.f; - problem.A.mmv(u, tmp); + Arithmetic::addProduct(tmp, -1.0, problem.A, u); double const localb = tmp * v; problem.A.mv(v, tmp); @@ -187,9 +187,9 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem { // b should be a descent direction { VectorType const direction = linearization.b; - VectorType tmp = linearization.b; // b - linearization.A.mmv(u, tmp); // b-Au - double const localA = tmp * direction; // <b-Au,v> + VectorType tmp = linearization.b; // b + Arithmetic::addProduct(tmp, -1.0, linearization.A, u); // b-Au + double const localA = tmp * direction; // <b-Au,v> linearization.A.mv(direction, tmp); // Av double const localb = tmp * direction; // <Av,v> @@ -316,7 +316,7 @@ class MyBlockProblem<MyConvexProblemTypeTEMPLATE>::IterateObject { if (j == m) localA = &(*it); // localA = A[m][m] else - it->mmv(u[j], localb); // localb -= A[m][j] * u[j] + Arithmetic::addProduct(localb, -1.0, *it, u[j]); // b-Au } assert(localA != nullptr); diff --git a/src/timestepping.cc b/src/timestepping.cc index 1922e0f2dd5a77c59169b1342b212e0f00fe53f9..07339e5ca074e707c3cf347a88b10bad97db0690 100644 --- a/src/timestepping.cc +++ b/src/timestepping.cc @@ -33,7 +33,7 @@ void ImplicitEuler<VectorType, MatrixType, FunctionType, dim>::setup( tau = _tau; problem_rhs = ell; - A.mmv(u_old, problem_rhs); + Arithmetic::addProduct(problem_rhs, -1.0, A, u_old); // For fixed tau, we'd only really have to do this once problem_A = A;