Use addProduct() in place of .mmv()

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);
 

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;