diff --git a/ext/LinearSolveForwardDiffExt.jl b/ext/LinearSolveForwardDiffExt.jl
index 64b87adef..e8fb3e5ef 100644
--- a/ext/LinearSolveForwardDiffExt.jl
+++ b/ext/LinearSolveForwardDiffExt.jl
@@ -763,4 +763,15 @@ function update_partials_list!(partial_matrix::SparseMatrixCSC, list_cache)
     return list_cache
 end
 
+
+# The MC64 matching in the vendored SupernodalLU solver runs its combinatorial
+# side (the max-product assignment) in Float64, so it needs a real magnitude
+# for a `Dual` entry.  Only the *ordering* of candidate pivots is decided from
+# these numbers; the factorization itself stays in Dual arithmetic, so taking
+# the primal here loses nothing.  Without this, `snlu` on a Dual matrix with a
+# missing or weak structural diagonal - exactly when `matching = :auto`
+# engages - throws `MethodError: no method matching Float64(::Dual)`.
+@inline LinearSolve.SupernodalLU._costabs(x::Dual) =
+    LinearSolve.SupernodalLU._costabs(ForwardDiff.value(x))
+
 end
diff --git a/test/Core/forwarddiff_overloads.jl b/test/Core/forwarddiff_overloads.jl
index 1aaa75462..262986107 100644
--- a/test/Core/forwarddiff_overloads.jl
+++ b/test/Core/forwarddiff_overloads.jl
@@ -750,3 +750,28 @@ struct ReinterpretTestTag end
     cache = init(LinearProblem(A, b), nothing)
     @test solve!(cache).u ≈ sol.u
 end
+
+@testset "SupernodalLU matching accepts Dual entries" begin
+    # MC64 matching decides its assignment in Float64, so it needs a real
+    # magnitude per entry.  An anti-diagonal matrix has no structural diagonal
+    # at all, so `matching = :auto` engages and the factorization goes through
+    # `_costabs` — which threw on Duals until the ForwardDiff extension
+    # supplied the primal-extracting overload its docstring promised.
+    SNLU = LinearSolve.SupernodalLU
+    n = 40
+    vals = [ForwardDiff.Dual{Nothing}(2.0 + i / n, 1.0) for i in 1:n]
+    P = sparse(collect(n:-1:1), collect(1:n), vals, n, n)
+    @test SNLU.needs_matching(P)
+    F = SNLU.snlu(P)
+    @test F.matched
+
+    # Solving must be right in both the value and the derivative.  With
+    # A(t) = diag-reversed(2 + i/n + t) and b fixed, each x_i = b_j / a_i, so
+    # dx_i/dt = -b_j / a_i^2.
+    b = randn(n)
+    x = similar(b, eltype(vals))
+    SNLU.solve!(x, F, ForwardDiff.Dual{Nothing}.(b, 0.0))
+    @test ForwardDiff.value.(x) ≈ [b[n + 1 - i] / (2.0 + i / n) for i in 1:n] rtol = 1.0e-10
+    @test ForwardDiff.partials.(x, 1) ≈
+        [-b[n + 1 - i] / (2.0 + i / n)^2 for i in 1:n] rtol = 1.0e-10
+end