ENH: Multithread the state loop in the inner-maximization sweeps

Motivation

The per-state inner maximization (bellman_operator!, compute_greedy!, and the evaluation sweep in evaluate!) is embarrassingly parallel: states are independent, writes go to disjoint entries of Tv/X, and the coefficient vector is read-only during a sweep. After #93/#99/#100/#101 the per-state work is allocation-free, so there is no GC pressure to spoil scaling — near-linear speedup in the number of cores is expected for the sweep, which now dominates VFI solve time.

The groundwork is already in place: CDPWorkspace, FunEvalCache, and DerivFunEvalCache are documented as one-per-thread (#95, #98), and per-state independence means threaded results are bitwise identical to single-threaded ones, independent of scheduling.

Design sketch

• Threading is opt-in (e.g. solve(cdp, method; nthreads=1) or threaded=false by default): the user's f, g, x_lb, x_ub must be thread-safe (pure functions of their arguments), which the package cannot verify.
• CDPWorkspace holds one evaluation context per parallel task, constructed once per solve:

struct EvalContext{TF<:FunEvalCache,TD}
    fec::TF
    dfecs::TD  # NTuple of DerivFunEvalCache, or nothing
end

• The sweep is chunked, with contexts indexed by chunk — not by Threads.threadid(), which is unsafe under task migration:

function _s_wise_max_foc_sweep_threaded!(cdp, C, Tv, X, contexts)
    foreach(ctx -> foreach(d -> set_coefs!(d, C), ctx.dfecs), contexts)
    n = size(cdp.interp.S, 1)
    chunks = collect(Iterators.partition(1:n, cld(n, length(contexts))))
    @sync for (c, idxs) in enumerate(chunks)
        Threads.@spawn begin
            ctx = contexts[c]
            for i in idxs
                s = _row(cdp.interp.S, i)
                Tv[i], X[i] = try
                    _s_wise_max_foc!(cdp, s, C, ctx.fec, ctx.dfecs, X[i])
                catch err
                    err isa InterruptException && rethrow()
                    _s_wise_max!(cdp, s, C, ctx.fec)
                end
            end
        end
    end
    return Tv, X
end

• The Brent-based sweeps (s_wise_max!, compute_greedy!, and the evaluate! sweep) get the same treatment; evaluate_policy!'s system assembly stays single-threaded for now (it is ~0.2 ms after #100; threading it would need per-chunk triplet buffers concatenated before sparse).
• Warm starts are unaffected: each state reads and writes only its own X[i].
• Load balancing: a static equal split should be adequate since per-state cost is roughly uniform; if corner/fallback states cluster (their FOC attempt + Brent fallback is costlier), smaller chunks with round-robin assignment are the first tuning knob.

Testing and benchmarking

• Correctness: threaded and single-threaded sweeps must agree exactly (==), not approximately — per-state independence makes this a well-defined requirement.
• CI: one test job with julia -t auto (runners have few cores; the point is exercising the code path, not measuring speedup).
• Benchmarks: PkgBenchmark comparisons must hold JULIA_NUM_THREADS fixed between baseline and target; a separate local scaling measurement (1/2/4/8 threads) belongs in the PR body rather than the suite.

Open questions

• API: keyword name and placement (solve kwarg vs. CDPWorkspace option), and whether nthreads should cap at Threads.nthreads() silently or error.
• Whether the threaded path should also cover the evaluate! sweep from the start or in a follow-up.

🤖 Drafted with Claude Code (Claude Fable 5)

Start by tracing solve, CDPWorkspace, bellman_operator!, compute_greedy!, and evaluate! to understand the existing sweeps and evaluation contexts. Review the stated threading design and correctness requirements, then resolve the open API and evaluate! scope questions. Done means opt-in threaded sweeps, exact agreement with single-threaded results, CI coverage with Julia threads, and documented benchmark comparisons.