Fix type inference for ForwardDiff Dual problems

The `init` function for `DualAbstractLinearProblem` with `DefaultLinearSolver`
previously had a runtime conditional that checked whether the auto-selected
algorithm was `GenericLUFactorization`. This caused the return type to be
`Union{LinearCache, DualLinearCache}`, breaking type inference.

This fix:
1. Removes the runtime conditional from `init`, ensuring it always returns
   `DualLinearCache` for type stability
2. Moves the `GenericLUFactorization` optimization to `solve!` instead,
   where it checks at solve-time and uses direct dual solving if applicable
3. Adds `_use_direct_dual_solve` helper functions to detect when direct
   dual solving should be used
4. Adds `_solve_direct_dual!` function that solves the dual system directly
   without separating primal/partials

This preserves the performance optimization for `GenericLUFactorization`
(which can work directly with any number type including Duals) while
ensuring `init` has a concrete return type.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Opus 4.5 <[email protected]>

Author: ChrisRackauckas

ext/LinearSolveForwardDiffExt.jl

@@ -265,22 +265,19 @@ function SciMLBase.init(prob::DualAbstractLinearProblem, alg::SciMLLinearSolveAl
     return __dual_init(prob, alg, args...; kwargs...)
 end
 
-# Opt out for GenericLUFactorization
-function SciMLBase.init(prob::DualAbstractLinearProblem, alg::GenericLUFactorization, args...; kwargs...)
-    return __init(prob, alg, args...; kwargs...)
-end
+# NOTE: Removed GenericLUFactorization opt-out from init to fix type inference.
+# The special handling for GenericLUFactorization is now done in solve! instead.
+# This ensures init always returns DualLinearCache for type stability.
 
-# Opt out for SparspakFactorization
+# Opt out for SparspakFactorization (sparse solvers can't handle Duals in the same way)
 function SciMLBase.init(prob::DualAbstractLinearProblem, alg::SparspakFactorization, args...; kwargs...)
     return __init(prob, alg, args...; kwargs...)
 end
 
+# NOTE: Removed the runtime conditional for DefaultLinearSolver that checked for
+# GenericLUFactorization. Now always use __dual_init for type stability.
 function SciMLBase.init(prob::DualAbstractLinearProblem, alg::DefaultLinearSolver, args...; kwargs...)
-    if alg.alg === DefaultAlgorithmChoice.GenericLUFactorization
-        return __init(prob, alg, args...; kwargs...)
-    else
-        return __dual_init(prob, alg, args...; kwargs...)
-    end
+    return __dual_init(prob, alg, args...; kwargs...)
 end
 
 function SciMLBase.init(prob::DualAbstractLinearProblem, alg::Nothing,
@@ -376,9 +373,25 @@ function SciMLBase.solve!(cache::DualLinearCache, args...; kwargs...)
     solve!(cache, getfield(cache, :linear_cache).alg, args...; kwargs...)
 end
 
+# Check if the algorithm should use the direct dual solve path
+# (algorithms that can work directly with Dual numbers without the primal/partials separation)
+function _use_direct_dual_solve(alg)
+    return alg isa GenericLUFactorization
+end
+
+function _use_direct_dual_solve(alg::DefaultLinearSolver)
+    return alg.alg === DefaultAlgorithmChoice.GenericLUFactorization
+end
+
 function SciMLBase.solve!(
         cache::DualLinearCache{DT}, alg::SciMLLinearSolveAlgorithm, args...; kwargs...) where {DT <:
                                                                                                ForwardDiff.Dual}
+    # Check if this algorithm can work directly with Duals (e.g., GenericLUFactorization)
+    # In that case, we solve the dual problem directly without separating primal/partials
+    if _use_direct_dual_solve(getfield(cache, :linear_cache).alg)
+        return _solve_direct_dual!(cache, alg, args...; kwargs...)
+    end
+
     primal_sol = linearsolve_forwarddiff_solve!(
         cache::DualLinearCache, getfield(cache, :linear_cache).alg, args...; kwargs...)
 
@@ -396,6 +409,37 @@ function SciMLBase.solve!(
     )
 end
 
+# Direct solve path for algorithms that can work with Dual numbers directly
+# This avoids the primal/partials separation overhead
+function _solve_direct_dual!(
+        cache::DualLinearCache{DT}, alg, args...; kwargs...) where {DT <: ForwardDiff.Dual}
+    # Reconstruct the dual A and b
+    dual_A = getfield(cache, :dual_A)
+    dual_b = getfield(cache, :dual_b)
+
+    # Solve directly with Duals using the generic LU path
+    # This works because GenericLUFactorization doesn't use BLAS and can handle any number type
+    dual_u = dual_A \ dual_b
+
+    # Update the cache
+    if getfield(cache, :dual_u) isa AbstractArray
+        getfield(cache, :dual_u) .= dual_u
+    else
+        setfield!(cache, :dual_u, dual_u)
+    end
+
+    # Also update the primal cache for consistency
+    primal_u = nodual_value.(dual_u)
+    if getfield(cache, :linear_cache).u isa AbstractArray
+        getfield(cache, :linear_cache).u .= primal_u
+    end
+
+    return SciMLBase.build_linear_solution(
+        getfield(cache, :linear_cache).alg, getfield(cache, :dual_u), nothing, cache;
+        retcode = ReturnCode.Success, iters = 1, stats = nothing
+    )
+end
+
 function setA!(dc::DualLinearCache, A)
     # Put the Dual-stripped versions in the LinearCache
     prop = nodual_value!(getproperty(dc.linear_cache, :A), A) # Update in-place