Eltype generic linearization

src/linearization.jl

@@ -96,39 +96,46 @@ function linearization_function(sys::AbstractSystem, inputs,
             return du
         end
     end
+    T = promote_type(eltype(u0), typeof(t0))
+    ct0 = DI.Constant(T(t0))
+    u0T = T.(prob.u0)
+    cu0T = DI.Constant(u0T)
+    cp = DI.Constant(p)
+
     if u0 === nothing
         uf_jac = h_jac = pf_jac = nothing
-        T = p isa MTKParameters ? eltype(p.tunable) : eltype(p)
+        Tp = promote_type(p isa MTKParameters ? eltype(p.tunable) : eltype(p), typeof(t0))
         hp_jac = PreparedJacobian{true}(
-            hp_fun, zeros(T, size(outputs)), autodiff, inputvals,
-            DI.Constant(prob.u0), DI.Constant(p), DI.Constant(t0))
+            hp_fun, zeros(Tp, size(outputs)), autodiff, inputvals,
+            cu0T, cp, DI.Constant(t0))
     else
         uf_fun = let fun = prob.f
             function uff(du, u, p, t)
                 SciMLBase.UJacobianWrapper(fun, t, p)(du, u)
             end
         end
+
         uf_jac = PreparedJacobian{true}(
-            uf_fun, similar(prob.u0), autodiff, prob.u0, DI.Constant(p), DI.Constant(t0))
+            uf_fun, similar(prob.u0, T), autodiff, u0T, cp, ct0)
         # observed function is a `GeneratedFunctionWrapper` with iip component
-        h_jac = PreparedJacobian{true}(h, similar(prob.u0, size(outputs)), autodiff,
-            prob.u0, DI.Constant(p), DI.Constant(t0))
+        h_jac = PreparedJacobian{true}(h, similar(prob.u0, T, size(outputs)), autodiff,
+            u0T, cp, ct0)
         pf_fun = let fun = prob.f, setter = setp_oop(sys, inputs)
             function pff(du, input, u, p, t)
                 p = setter(p, input)
                 SciMLBase.ParamJacobianWrapper(fun, t, u)(du, p)
             end
         end
-        pf_jac = PreparedJacobian{true}(pf_fun, similar(prob.u0), autodiff, inputvals,
-            DI.Constant(prob.u0), DI.Constant(p), DI.Constant(t0))
+        pf_jac = PreparedJacobian{true}(pf_fun, similar(prob.u0, T), autodiff, inputvals,
+            cu0T, cp, ct0)
         hp_jac = PreparedJacobian{true}(
-            hp_fun, similar(prob.u0, size(outputs)), autodiff, inputvals,
-            DI.Constant(prob.u0), DI.Constant(p), DI.Constant(t0))
+            hp_fun, similar(prob.u0, T, size(outputs)), autodiff, inputvals,
+            cu0T, cp, ct0)
     end
 
     lin_fun = LinearizationFunction(
         diff_idxs, alge_idxs, inputs, length(unknowns(sys)),
-        prob, h, u0 === nothing ? nothing : similar(u0), uf_jac, h_jac, pf_jac,
+        prob, h, u0 === nothing ? nothing : similar(u0, T), uf_jac, h_jac, pf_jac,
         hp_jac, initializealg, initialization_kwargs)
     return lin_fun, sys
 end

test/linearize.jl

@@ -37,6 +37,8 @@ eqs = [u ~ kp * (r - y)
 @named sys = System(eqs, t)
 
 lsys, ssys, extras = linearize(sys, [r], [y])
+lsysbig, ssysbig, extrasbig = linearize(sys, [r], [y], t=big(0.0))
+@test lsys.A[] == Float64.(lsysbig.A[])
 lprob = LinearizationProblem(sys, [r], [y])
 lsys2, extras2 = solve(lprob)
 lsys3, _ = linearize(sys, [r], [y]; autodiff = AutoFiniteDiff())