Wrap Jacobian with both dense and sparse FunctionWrapper signatures

When an ODEFunction has a sparse sparsity pattern but no jac_prototype,
the solver's build_J_W may create a SparseMatrixCSC J from the sparsity
pattern (when using AutoSparse). The previous code only wrapped the
Jacobian with a dense Matrix{Float64} signature, causing "No matching
function wrapper was found!" at solve time.

Now detects non-dense sparsity patterns and creates a FunctionWrappersWrapper
with both dense and sparse matrix signatures so the Jacobian call works
regardless of which matrix type the solver allocates.

Uses Base.promote_op to compute the sparse matrix type without allocating.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

src/solve.jl

@@ -810,12 +810,33 @@ function promote_f(
         if f.tgrad !== nothing && !(f.tgrad isa FunctionWrappersWrappers.FunctionWrappersWrapper)
             f = @set f.tgrad = wrapfun_jac_iip(f.tgrad, (u0, u0, p, t))
         end
-        # Wrap the Jacobian if present, so its type is also erased
+        # Wrap the Jacobian if present, so its type is also erased.
+        # Include both dense and sparse matrix signatures when the function
+        # has a sparsity pattern, since the solver may use either depending on
+        # the autodiff configuration (AutoSparse creates sparse J from sparsity).
         if f.jac !== nothing && !(f.jac isa FunctionWrappersWrappers.FunctionWrappersWrapper)
             n = length(u0)
-            J_proto = f.jac_prototype !== nothing ? similar(f.jac_prototype, uElType) :
-                zeros(uElType, n, n)
-            f = @set f.jac = wrapfun_jac_iip(f.jac, (J_proto, u0, p, t))
+            if f.jac_prototype !== nothing
+                J_proto = similar(f.jac_prototype, uElType)
+                f = @set f.jac = wrapfun_jac_iip(f.jac, (J_proto, u0, p, t))
+            elseif isdefined(f, :sparsity) && f.sparsity isa AbstractMatrix &&
+                    !(f.sparsity isa Matrix)
+                # The sparsity pattern is a non-dense matrix (e.g. SparseMatrixCSC).
+                # The solver may call the Jacobian with either a dense or sparse matrix
+                # depending on the autodiff config, so wrap for both signatures.
+                # Compute types without allocating the matrices.
+                dense_sig = Tuple{Matrix{uElType}, typeof(u0), typeof(p), typeof(t)}
+                sparse_J_T = Base.promote_op(similar, typeof(f.sparsity), Type{uElType})
+                sparse_sig = Tuple{sparse_J_T, typeof(u0), typeof(p), typeof(t)}
+                f = @set f.jac = FunctionWrappersWrappers.FunctionWrappersWrapper(
+                    Void(f.jac),
+                    (dense_sig, sparse_sig),
+                    (Nothing, Nothing)
+                )
+            else
+                J_proto = zeros(uElType, n, n)
+                f = @set f.jac = wrapfun_jac_iip(f.jac, (J_proto, u0, p, t))
+            end
         end
         return unwrapped_f(f, wrapfun_iip(f.f, (u0, u0, p, t), Val(CS)))
     else