From 31e05ab59948d0180129800baa2b9f15b33b9edb Mon Sep 17 00:00:00 2001
From: deXtoRious <dextorious@gmail.com>
Date: Wed, 8 Nov 2017 01:24:31 +0200
Subject: [PATCH 1/2] DiffEqDiffTools API changes

---
 src/integrators/rosenbrock_integrators.jl | 42 +++++++++++++++--------
 src/integrators/sdirk_integrators.jl      | 36 ++++++++++++-------
 2 files changed, 52 insertions(+), 26 deletions(-)

diff --git a/src/integrators/rosenbrock_integrators.jl b/src/integrators/rosenbrock_integrators.jl
index 49005b139d..f6d73c66db 100644
--- a/src/integrators/rosenbrock_integrators.jl
+++ b/src/integrators/rosenbrock_integrators.jl
@@ -32,7 +32,8 @@ end
       if alg_autodiff(integrator.alg)
         ForwardDiff.derivative!(dT, tf, vec(du2), t)
       else
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, vec(du2), Val{:DiffEqDerivativeWrapper})
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
       end
     end
   end
@@ -55,7 +56,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
       for j in 1:length(u), i in 1:length(u)
@@ -149,7 +151,8 @@ end
       if alg_autodiff(integrator.alg)
         ForwardDiff.derivative!(dT, tf, vec(du2), t)
       else
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, vec(du2), Val{:DiffEqDerivativeWrapper})
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
       end
     end
   end
@@ -171,7 +174,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
       for j in 1:length(u), i in 1:length(u)
@@ -456,7 +460,8 @@ end
       if alg_autodiff(integrator.alg)
         ForwardDiff.derivative!(dT, tf, vec(du2), t)
       else
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, vec(du2), Val{:DiffEqDerivativeWrapper})
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
       end
     end
   end
@@ -479,7 +484,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
       for j in 1:length(u), i in 1:length(u)
@@ -664,7 +670,8 @@ end
       if alg_autodiff(integrator.alg)
         ForwardDiff.derivative!(dT, tf, vec(du2), t)
       else
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, vec(du2), Val{:DiffEqDerivativeWrapper})
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
       end
     end
   end
@@ -686,7 +693,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
       for j in 1:length(u), i in 1:length(u)
@@ -896,7 +904,8 @@ end
       if alg_autodiff(integrator.alg)
         ForwardDiff.derivative!(dT, tf, vec(du2), t)
       else
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, vec(du2), Val{:DiffEqDerivativeWrapper})
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
       end
     end
   end
@@ -919,7 +928,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
       for j in 1:length(u), i in 1:length(u)
@@ -1152,7 +1162,8 @@ end
       if alg_autodiff(integrator.alg)
         ForwardDiff.derivative!(dT, tf, vec(du2), t)
       else
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, vec(du2), Val{:DiffEqDerivativeWrapper})
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
       end
     end
   end
@@ -1177,7 +1188,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
       for j in 1:length(u), i in 1:length(u)
@@ -1537,7 +1549,8 @@ end
       if alg_autodiff(integrator.alg)
         ForwardDiff.derivative!(dT, tf, vec(du2), t)
       else
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, vec(du2), Val{:DiffEqDerivativeWrapper})
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
       end
     end
   end
@@ -1562,7 +1575,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
       for j in 1:length(u), i in 1:length(u)
diff --git a/src/integrators/sdirk_integrators.jl b/src/integrators/sdirk_integrators.jl
index a0b63cc827..4c496a2395 100644
--- a/src/integrators/sdirk_integrators.jl
+++ b/src/integrators/sdirk_integrators.jl
@@ -135,7 +135,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -364,7 +365,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -911,7 +913,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -1234,7 +1237,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -1557,7 +1561,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -1907,7 +1912,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -2420,7 +2426,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -3028,7 +3035,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -3603,7 +3611,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -4160,7 +4169,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -4815,7 +4825,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end
@@ -5568,7 +5579,8 @@ end
         if alg_autodiff(integrator.alg)
           ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
         else
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), Val{integrator.alg.diff_type}, vec(du1), Val{:JacobianWrapper})
+          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
         end
       end
     end

From f6e7a6b71a24d6dcbd9e0ad9b10a029574fbb47f Mon Sep 17 00:00:00 2001
From: dextorious <dextorious@gmail.com>
Date: Wed, 8 Nov 2017 21:30:53 +0200
Subject: [PATCH 2/2] Convenience API for derivatives and Jacobians.

---
 src/derivative_wrappers.jl                | 20 +++++
 src/integrators/rosenbrock_integrators.jl | 98 ++++-------------------
 src/integrators/sdirk_integrators.jl      | 84 +++----------------
 3 files changed, 46 insertions(+), 156 deletions(-)

diff --git a/src/derivative_wrappers.jl b/src/derivative_wrappers.jl
index 4a35ea0a56..d508a244ae 100644
--- a/src/derivative_wrappers.jl
+++ b/src/derivative_wrappers.jl
@@ -46,3 +46,23 @@ mutable struct UDerivativeWrapper{F,tType} <: Function
   t::tType
 end
 (p::UDerivativeWrapper)(u) = p.f(p.t,u)
+
+function derivative!(df::AbstractArray{<:Number}, f, x::Union{Number,AbstractArray{<:Number}}, fx::AbstractArray{<:Number}, integrator::DEIntegrator)
+    if alg_autodiff(integrator.alg)
+        ForwardDiff.derivative!(df, f, fx, x)
+    else
+        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+        DiffEqDiffTools.finite_difference!(df, f, x, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, fx)
+    end
+    nothing
+end
+
+function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, fx::AbstractArray{<:Number}, integrator::DEIntegrator, jac_config)
+    if alg_autodiff(integrator.alg)
+      ForwardDiff.jacobian!(J, f, fx, x, jac_config)
+    else
+      RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
+      DiffEqDiffTools.finite_difference_jacobian!(J, f, x, integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, fx)
+    end
+    nothing
+end
diff --git a/src/integrators/rosenbrock_integrators.jl b/src/integrators/rosenbrock_integrators.jl
index f6d73c66db..487ebb4776 100644
--- a/src/integrators/rosenbrock_integrators.jl
+++ b/src/integrators/rosenbrock_integrators.jl
@@ -29,12 +29,7 @@ end
     else
       tf.vf.sizeu = sizeu
       tf.uprev = uprev
-      if alg_autodiff(integrator.alg)
-        ForwardDiff.derivative!(dT, tf, vec(du2), t)
-      else
-        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
-      end
+      derivative!(dT, tf, t, vec(du2), integrator)
     end
   end
 
@@ -53,12 +48,7 @@ end
       else
         uf.vfr.sizeu = sizeu
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
       for j in 1:length(u), i in 1:length(u)
           @inbounds W[i,j] = mass_matrix[i,j] - γ*J[i,j]
@@ -148,12 +138,7 @@ end
     else
       tf.vf.sizeu = sizeu
       tf.uprev = uprev
-      if alg_autodiff(integrator.alg)
-        ForwardDiff.derivative!(dT, tf, vec(du2), t)
-      else
-        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
-      end
+      derivative!(dT, tf, t, vec(du2), integrator)
     end
   end
 
@@ -171,12 +156,7 @@ end
       else
         uf.vfr.sizeu = sizeu
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
       for j in 1:length(u), i in 1:length(u)
         @inbounds W[i,j] = mass_matrix[i,j] - γ*J[i,j]
@@ -457,12 +437,7 @@ end
     else
       tf.vf.sizeu = sizeu
       tf.uprev = uprev
-      if alg_autodiff(integrator.alg)
-        ForwardDiff.derivative!(dT, tf, vec(du2), t)
-      else
-        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
-      end
+      derivative!(dT, tf, t, vec(du2), integrator)
     end
   end
 
@@ -481,12 +456,7 @@ end
       else
         uf.vfr.sizeu = sizeu
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
       for j in 1:length(u), i in 1:length(u)
           @inbounds W[i,j] = mass_matrix[i,j]/dtgamma - J[i,j]
@@ -667,12 +637,7 @@ end
     else
       tf.vf.sizeu = sizeu
       tf.uprev = uprev
-      if alg_autodiff(integrator.alg)
-        ForwardDiff.derivative!(dT, tf, vec(du2), t)
-      else
-        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
-      end
+      derivative!(dT, tf, t, vec(du2), integrator)
     end
   end
 
@@ -690,12 +655,7 @@ end
       else
         uf.vfr.sizeu = sizeu
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
       for j in 1:length(u), i in 1:length(u)
           @inbounds W[i,j] = mass_matrix[i,j]/dtgamma - J[i,j]
@@ -901,12 +861,7 @@ end
     else
       tf.vf.sizeu = sizeu
       tf.uprev = uprev
-      if alg_autodiff(integrator.alg)
-        ForwardDiff.derivative!(dT, tf, vec(du2), t)
-      else
-        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
-      end
+      derivative!(dT, tf, t, vec(du2), integrator)
     end
   end
 
@@ -925,12 +880,7 @@ end
       else
         uf.vfr.sizeu = sizeu
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
       for j in 1:length(u), i in 1:length(u)
           @inbounds W[i,j] = mass_matrix[i,j]/dtgamma - J[i,j]
@@ -1159,12 +1109,7 @@ end
     else
       tf.vf.sizeu = sizeu
       tf.uprev = uprev
-      if alg_autodiff(integrator.alg)
-        ForwardDiff.derivative!(dT, tf, vec(du2), t)
-      else
-        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
-      end
+      derivative!(dT, tf, t, vec(du2), integrator)
     end
   end
 
@@ -1185,12 +1130,7 @@ end
       else
         uf.vfr.sizeu = sizeu
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
       for j in 1:length(u), i in 1:length(u)
           @inbounds W[i,j] = mass_matrix[i,j]/dtgamma - J[i,j]
@@ -1546,12 +1486,7 @@ end
     else
       tf.vf.sizeu = sizeu
       tf.uprev = uprev
-      if alg_autodiff(integrator.alg)
-        ForwardDiff.derivative!(dT, tf, vec(du2), t)
-      else
-        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-        DiffEqDiffTools.finite_difference!(dT, tf, t, integrator.alg.diff_type, RealOrComplex, Val{:DiffEqDerivativeWrapper}, vec(du2))
-      end
+      derivative!(dT, tf, t, vec(du2), integrator)
     end
   end
 
@@ -1572,12 +1507,7 @@ end
       else
         uf.vfr.sizeu = sizeu
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J, uf, vec(du1), vec(uprev), jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
       for j in 1:length(u), i in 1:length(u)
           @inbounds W[i,j] = mass_matrix[i,j]/dtgamma - J[i,j]
diff --git a/src/integrators/sdirk_integrators.jl b/src/integrators/sdirk_integrators.jl
index 4c496a2395..d651a8d35e 100644
--- a/src/integrators/sdirk_integrators.jl
+++ b/src/integrators/sdirk_integrators.jl
@@ -132,12 +132,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -362,12 +357,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -910,12 +900,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -1234,12 +1219,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -1558,12 +1538,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -1909,12 +1884,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -2423,12 +2393,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -3032,12 +2997,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -3608,12 +3568,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -4166,12 +4121,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -4822,12 +4772,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated
@@ -5576,12 +5521,7 @@ end
         f(Val{:jac},t,uprev,J)
       else
         uf.t = t
-        if alg_autodiff(integrator.alg)
-          ForwardDiff.jacobian!(J,uf,vec(du1),vec(uprev),jac_config)
-        else
-          RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
-          DiffEqDiffTools.finite_difference_jacobian!(J, uf, vec(uprev), integrator.alg.diff_type, RealOrComplex, Val{:JacobianWrapper}, vec(du1))
-        end
+        jacobian!(J, uf, vec(uprev), vec(du1), integrator, jac_config)
       end
     end
     # skip calculation of W if step is repeated