diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 9d4ccb3f5da38ac1d83bd5d7b63f7c7a1059c99d..504aba445a925f60bf1d95f6dec1edf05f9890d0 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -93,6 +93,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")
@@ -171,7 +172,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)
@@ -188,6 +189,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 061a14395867903e47774f31963ec93d0a3040a4..e0a90caeb8393e4ed354f0a272cd586f583670a8 100644
--- a/python/dune/perftool/pdelab/driver/solve.py
+++ b/python/dune/perftool/pdelab/driver/solve.py
@@ -50,8 +50,9 @@ def dune_solve():
         snp = name_stationarynonlinearproblemsolver(go_type, go)
         solve = "{}.apply();".format(snp)
 
-    print_residual()
-    if not matrix_free:
+    if get_form_option("generate_residuals"):
+        print_residual()
+    if get_form_option("generate_jacobians"):
         print_matrix()
 
     if get_option('instrumentation_level') >= 2:
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 3f12af75d7f5a7bc4dfafcc5df23d15fe76f6158..0a3636ccc73fcb44ef2ca24815de2a291b8c4b73 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -783,7 +783,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
@@ -826,10 +826,7 @@ def generate_jacobian_kernels(form, original_form):
         jacform = offdiagonal_block_jacobian(jacform)
 
     operator_kernels = {}
-    # Decide if the jacobian apply or jacobian method is needed
-    if get_form_option("matrix_free"):
-        # Matrix-free computations need the jacobian apply method
-
+    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)
@@ -856,8 +853,7 @@ 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)]
-    else:
-        # Non matrix-free computations need the jacobian method
+    if get_form_option("generate_jacobians"):
         with global_context(form_type="jacobian"):
             if get_form_option("generate_jacobians"):
                 if get_form_option("sumfact"):