Make the operator really matrix-free

Closes: #2

src/linearsolver.hh

@@ -54,36 +54,37 @@ namespace PPS {
   };
 
   template<class GridOperator, class X, class Y>
-  class FakeLinearizedMatrixFreeOperator :
+  class LinearizedMatrixFreeOperator :
     public Dune::LinearOperator<X, Y>,
     public LinearizedOperatorMixin<X>
   {
     using field_type = typename Y::field_type;
-    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<field_type, 1, 1> >;
 
     const GridOperator *go_;
-    Matrix matrix_;
+    X linearizationPoint_;
 
   public:
-    FakeLinearizedMatrixFreeOperator(const GridOperator * go) :
-      go_(go), matrix_(go->pattern().template make<Matrix>())
+    LinearizedMatrixFreeOperator(const GridOperator * go) :
+      go_(go)
     { }
 
     void linearizeAt(const X &linearizationPoint) override {
-      matrix_ = 0;
-      go_->jacobian(linearizationPoint, matrix_);
+      linearizationPoint_ = linearizationPoint;
     }
 
     //! apply operator to x:  \f$ y = A(x) \f$
     void apply (const X& x, Y& y) const override
     {
-      matrix_.mv(x,y);
+      go_->nonlinear_jacobian_apply(linearizationPoint_, x, y);
     }
 
     //! 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
     {
-      matrix_.usmv(alpha,x,y);
+      Y y2(y.size());
+      y2 = 0;
+      go_->nonlinear_jacobian_apply(linearizationPoint_, x, y2);
+      y.axpy(alpha, y2);
     }
 
     //! Category of the solver (see SolverCategory::Category)
@@ -106,7 +107,7 @@ namespace PPS {
       return std::make_unique<Operator>(go);
     }
     if(type == "on-the-fly" || type == "matrix-free") {
-      using Operator = FakeLinearizedMatrixFreeOperator<GridOperator, X, Y>;
+      using Operator = LinearizedMatrixFreeOperator<GridOperator, X, Y>;
       return std::make_unique<Operator>(go);
     }
     DUNE_THROW(Dune::NotImplemented, "operator=" << type);