From dea47047bc5adfc5a031c3a88d56130edf413808 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Mon, 23 Feb 2026 07:44:46 -0500
Subject: [PATCH] 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.

All type computations use Base.promote_op to avoid allocating prototype
matrices just for their types.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
---
 src/solve.jl | 32 +++++++++++++++++++++++++++-----
 1 file changed, 27 insertions(+), 5 deletions(-)

diff --git a/src/solve.jl b/src/solve.jl
index 6e3512040..1d2107331 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -810,12 +810,34 @@ 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_T = Base.promote_op(similar, typeof(f.jac_prototype), Type{uElType})
+                sig = Tuple{J_T, typeof(u0), typeof(p), typeof(t)}
+                f = @set f.jac = FunctionWrappersWrappers.FunctionWrappersWrapper(
+                    Void(f.jac), (sig,), (Nothing,))
+            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.
+                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
+                sig = Tuple{Matrix{uElType}, typeof(u0), typeof(p), typeof(t)}
+                f = @set f.jac = FunctionWrappersWrappers.FunctionWrappersWrapper(
+                    Void(f.jac), (sig,), (Nothing,))
+            end
         end
         return unwrapped_f(f, wrapfun_iip(f.f, (u0, u0, p, t), Val(CS)))
     else