diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 25fd514582a9873eb26cf59f50897b1dbe59cdc9..c0f4900b8fff2f77b7cbefd593214cfc8d3af87d 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -94,6 +94,7 @@ class PerftoolFormOptionsArray(ImmutableRecord):
     vectorization_target = PerftoolOption(_type=float, helpstr="The cost function target for the 'target' cost model. Only needed to verify the cost model itself, do not use light-heartedly!!!")
     simplify = PerftoolOption(default=False, helpstr="Whether to simplify expressions using sympy")
     generate_jacobians = PerftoolOption(default=True, helpstr="Whether jacobian_* methods should be generated. This is set to false automatically, when numerical_jacobian is set to true.")
+    generate_jacobian_apply = PerftoolOption(default=False, helpstr="Wether jacobian_allpy_* methods should be generated.")
     generate_residuals = PerftoolOption(default=True, helpstr="Whether alpha_* methods should be generated.")
     unroll_dimension_loops = PerftoolOption(default=False, helpstr="whether loops over the geometric dimension should be unrolled")
     precompute_quadrature_info = PerftoolOption(default=True, helpstr="compute quadrature points and weights in the constructor of the local operator")
@@ -172,7 +173,7 @@ def process_form_options(opt, form):
         opt = opt.copy(unroll_dimension_loops=True)
 
     if opt.numerical_jacobian:
-        opt = opt.copy(generate_jacobians=False)
+        opt = opt.copy(generate_jacobians=False, generate_jacobian_apply=False)
 
     if opt.form is None:
         opt = opt.copy(form=form)
@@ -189,6 +190,10 @@ def process_form_options(opt, form):
                        generate_jacobians=False,
                        matrix_free=True,
                        )
+
+    if opt.matrix_free:
+        opt = opt.copy(generate_jacobian_apply=True)
+
     return opt
 
 
diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py
index 3877165c574ac8ef7b6d7197b56f97d071f38f42..591c022ddedd8656102f1594d5148714bd7561de 100644
--- a/python/dune/perftool/pdelab/driver/__init__.py
+++ b/python/dune/perftool/pdelab/driver/__init__.py
@@ -253,7 +253,7 @@ def generate_driver():
         from dune.perftool.loopy.target import type_floatingpoint
         pre_include("#define HP_TIMER_OPCOUNTER {}".format(type_floatingpoint()), filetag="driver")
         evaluate_residual_timer()
-        if get_form_option("matrix_free"):
+        if get_form_option("generate_jacobian_apply"):
             apply_jacobian_timer()
     elif is_stationary():
         from dune.perftool.pdelab.driver.solve import dune_solve
diff --git a/python/dune/perftool/pdelab/driver/solve.py b/python/dune/perftool/pdelab/driver/solve.py
index cdf41b82ddd3a4f33e3beb0b8ef29cfab6d012e5..25747850ba9188f6141b6f848584daed3dd09921 100644
--- a/python/dune/perftool/pdelab/driver/solve.py
+++ b/python/dune/perftool/pdelab/driver/solve.py
@@ -52,8 +52,10 @@ def dune_solve():
         snp = name_stationarynonlinearproblemsolver(go_type, go)
         solve = "{}.apply();".format(snp)
 
-    print_residual()
-    print_matrix()
+    if get_form_option("generate_residuals"):
+        print_residual()
+    if get_form_option("generate_jacobians"):
+        print_matrix()
 
     if get_option('instrumentation_level') >= 2:
         from dune.perftool.pdelab.driver.timings import setup_timer, name_timing_stream, name_timing_identifier
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 4daeb45a7627f7bda1d83c178e83a20a5362607c..1f808ef1c90f4b358b60d6f4744a2322a726e055 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -797,7 +797,7 @@ def generate_residual_kernels(form, original_form):
                                      )
 
                     # In the case of matrix free operator evaluation we need jacobian apply methods
-                    if get_form_option("matrix_free"):
+                    if get_form_option("generate_jacobian_apply"):
                         from dune.perftool.pdelab.driver import is_linear
                         if is_linear(original_form):
                             # Numeical jacobian apply base class
@@ -840,35 +840,7 @@ def generate_jacobian_kernels(form, original_form):
         jacform = offdiagonal_block_jacobian(jacform)
 
     operator_kernels = {}
-    with global_context(form_type="jacobian"):
-        if get_form_option("generate_jacobians"):
-            if get_form_option("sumfact"):
-                was_sumfact = True
-                if get_form_option("sumfact_regular_jacobians"):
-                    set_form_option("sumfact", False)
-            for measure in set(i.integral_type() for i in jacform.integrals()):
-                logger.info("generate_jacobian_kernels: measure {}".format(measure))
-                with global_context(integral_type=measure):
-                    from dune.perftool.pdelab.signatures import assembler_routine_name
-                    with global_context(kernel=assembler_routine_name()):
-                        kernel = [k for k in get_backend(interface="generate_kernels_per_integral")(jacform.integrals_by_type(measure))]
-                operator_kernels[(measure, 'jacobian')] = kernel
-
-            # Generate dummy functions for those kernels, that vanished in the differentiation process
-            # We *could* solve this problem by using lambda_* terms but we do not really want that, so
-            # we use empty jacobian assembly methods instead
-            alpha_measures = set(i.integral_type() for i in form.integrals())
-            jacobian_measures = set(i.integral_type() for i in jacform.integrals())
-            for it in alpha_measures - jacobian_measures:
-                with global_context(integral_type=it):
-                    from dune.perftool.pdelab.signatures import assembly_routine_signature
-                    operator_kernels[(it, 'jacobian')] = [LoopyKernelMethod(assembly_routine_signature(), kernel=None)]
-            if get_form_option("sumfact_regular_jacobians"):
-                if was_sumfact:
-                    set_form_option("sumfact", True)
-
-    # Jacobian apply methods for matrix-free computations
-    if get_form_option("matrix_free"):
+    if get_form_option("generate_jacobian_apply"):
         # The apply vector has reserved index 1 so we directly use Coefficient class from ufl
         from ufl import Coefficient
         apply_coefficient = Coefficient(form.coefficients()[0].ufl_element(), 1)
@@ -895,6 +867,33 @@ def generate_jacobian_kernels(form, original_form):
                     with global_context(integral_type=it):
                         from dune.perftool.pdelab.signatures import assembly_routine_signature
                         operator_kernels[(it, 'jacobian_apply')] = [LoopyKernelMethod(assembly_routine_signature(), kernel=None)]
+    if get_form_option("generate_jacobians"):
+        with global_context(form_type="jacobian"):
+            if get_form_option("generate_jacobians"):
+                if get_form_option("sumfact"):
+                    was_sumfact = True
+                    if get_form_option("sumfact_regular_jacobians"):
+                        set_form_option("sumfact", False)
+                for measure in set(i.integral_type() for i in jacform.integrals()):
+                    logger.info("generate_jacobian_kernels: measure {}".format(measure))
+                    with global_context(integral_type=measure):
+                        from dune.perftool.pdelab.signatures import assembler_routine_name
+                        with global_context(kernel=assembler_routine_name()):
+                            kernel = [k for k in get_backend(interface="generate_kernels_per_integral")(jacform.integrals_by_type(measure))]
+                    operator_kernels[(measure, 'jacobian')] = kernel
+
+                # Generate dummy functions for those kernels, that vanished in the differentiation process
+                # We *could* solve this problem by using lambda_* terms but we do not really want that, so
+                # we use empty jacobian assembly methods instead
+                alpha_measures = set(i.integral_type() for i in form.integrals())
+                jacobian_measures = set(i.integral_type() for i in jacform.integrals())
+                for it in alpha_measures - jacobian_measures:
+                    with global_context(integral_type=it):
+                        from dune.perftool.pdelab.signatures import assembly_routine_signature
+                        operator_kernels[(it, 'jacobian')] = [LoopyKernelMethod(assembly_routine_signature(), kernel=None)]
+                if get_form_option("sumfact_regular_jacobians"):
+                    if was_sumfact:
+                        set_form_option("sumfact", True)
 
     return operator_kernels