Supply the _costabs Dual overload SupernodalLU's matching documents (#1118)

`src/SupernodalLU/matching.jl` says of `_costabs`:

    Value types with no exact Float64 conversion (e.g. ForwardDiff.Dual) get
    an overload in the package extension that extracts the primal value.

No such overload exists. It was left behind when the solver was vendored from
PurePardiso.jl, so the comment documents behaviour the package does not have.

MC64 matching decides its max-product assignment in Float64, so `_costabs`
needs a real magnitude per entry. Without the overload, factoring a `Dual`
matrix whose diagonal is missing or weak - exactly when `matching = :auto`
engages - throws:

    MethodError: no method matching Float64(::ForwardDiff.Dual{Nothing, Float64, 1})
      _costabs at matching.jl:206 [inlined]
      mc64_matching(::SparseMatrixCSC{Dual{Nothing, Float64, 1}, Int64})

Only the *ordering* of candidate pivots comes from these numbers; the
factorization itself stays in Dual arithmetic, so taking the primal loses
nothing.

Reachable today through `LinearSolve.__init`. The public path is unaffected:
`__dual_init` strips duals before any algorithm sees them, so `solve` and
`init` on a Dual sparse problem factor in Float64 and never reach this code.

`_scalarval` needs no such fix - returning `NaN` for a Dual is intended and
documented there (`anorm` is informational; thresholds stay in native
arithmetic).

Test: an anti-diagonal matrix, which has no structural diagonal at all so
matching always engages, checked in both value and derivative against the
closed form. It fails on the unfixed code at the `snlu` call.

GROUP=Core passes.



Claude-Session: https://claude.ai/code/session_017rr42T1vFRRT8D7DMAmJzf

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

Author: 3 people

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

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