Merge branch 'applyscaleadd-no-extra-vector' into 'master'

Avoid allocating extra vector in matrix-free applyscaleadd()

Closes #4

See merge request jfahl_01/pacxx-projectseminar-2019!17

src/gridoperator.hh

@@ -133,10 +133,12 @@ namespace PPS {
      * \param lp Point where to linearize
      * \param x  Point to apply the Jacobian to
      * \param y  Where to add the result
+     * \param trafo A functor used to transform values before accumulation
+     *              into y (usually to apply a factor)
      */
-    template<class Domain, class Range>
+    template<class Domain, class Range, class Trafo>
     void nonlinear_jacobian_apply(const Domain &lp, const Domain & x,
-                                  Range & y) const
+                                  Range & y, Trafo trafo) const
     {
       std::vector<typename Domain::field_type> lpl, xl;
       std::vector<typename Range::field_type> yl;
@@ -181,7 +183,7 @@ namespace PPS {
 
         // Notify assembler engine about unbinds
         for(auto i : Dune::range(size))
-          y[iset.subIndex(element, i, element.dimension)] += yl[i];
+          y[iset.subIndex(element, i, element.dimension)] += trafo(yl[i]);
 
       } // it
 

src/linearsolver.hh

@@ -75,17 +75,16 @@ namespace PPS {
     //! apply operator to x:  \f$ y = A(x) \f$
     void apply (const X& x, Y& y) const override
     {
+      auto id = [](auto v){ return v; };
       y = 0;
-      go_->nonlinear_jacobian_apply(linearizationPoint_, x, y);
+      go_->nonlinear_jacobian_apply(linearizationPoint_, x, y, id);
     }
 
     //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
     void applyscaleadd (field_type alpha, const X& x, Y& y) const override
     {
-      Y y2(y.size());
-      y2 = 0;
-      go_->nonlinear_jacobian_apply(linearizationPoint_, x, y2);
-      y.axpy(alpha, y2);
+      auto scale = [=](auto v){ return alpha*v; };
+      go_->nonlinear_jacobian_apply(linearizationPoint_, x, y, scale);
     }
 
     //! Category of the solver (see SolverCategory::Category)