NonlinearSolveAlg: give inner NonlinearFunction W's structure as jac_prototype#3819
NonlinearSolveAlg: give inner NonlinearFunction W's structure as jac_prototype#3819ChrisRackauckas-Claude wants to merge 2 commits into
Conversation
…prototype Under W-reuse the inner NonlinearSolve builds its Jacobian buffer without a jac_prototype, and NonlinearSolveBase allocates a dense matrix for analytic-jac functions. A sparse/structured W then degrades to dense LU, and sparse-only linear solvers piped in from the ODE algorithm crash outright (KLUFactorization: `FieldError: Matrix has no field nzval`). Pass `jac_prototype = similar(W_for_reuse)` so the inner Jacobian mirrors W's type and sparsity; `copyto!(J_out, W)` in the jac closure is then structure-preserving and sparse factorizations work. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
|
CI triage: the failing lanes here are master-baseline, not from this change — StochasticDiffEq QA/downgrade + GPU lanes fail on all concurrent PRs including the 2-line #3816; the OrdinaryDiffEqNonlinearSolve Core lane failure is the pre-existing |
The jac_prototype fix surfaced a latent resize hazard: the build-time jac
closure captured the original W array, and the build-time jac_prototype kept
its original size, so after `resize!(integrator, n)` the length-mismatch
rebuild in initialize! produced an inner Jacobian/linear cache at the stale
size (CI Cache Tests: `DimensionMismatch: B has leading dimension 2, but
needs 1`; on master the same stale closure was masked because the dense inner
J was sized from u instead of the prototype).
Replace the closure with a `WReuseJac` callable struct holding a
`Ref{typeof(W)}`. Its type depends only on W's type, not on a captured
binding, so when resize_J_W! replaces the W array the length-mismatch branch
just swaps the Ref target and remakes the NonlinearFunction with a same-typed
`similar(cache.W)` prototype — cache.prob's concrete type is preserved.
Adds a growing-cell resize test (mirroring ode_cache_tests) asserting the
resize completes without dimension errors for TRBDF2/KenCarp4/ImplicitEuler
and a Success retcode for TRBDF2. Retcodes for the other methods on this
marginally-stable scenario are method-sensitive on master too
(ImplicitEuler+NLNewton is Unstable on the deterministic variant), so they
are deliberately not asserted.
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
|
The Integrators_II |
2 participants
Bug
Under W-reuse (#3693/#3753) the inner
NonlinearFunctiongets the analytic-jac closurecopyto!(J_out, W)but nojac_prototype. NonlinearSolveBase allocates a dense J for analytic-jac functions without a prototype (jacobian.jl:safe_similar(fu, n, n)), so for sparse/structured W:_nlalg_with_linsolve): hard crash —FBDF(linsolve = KLUFactorization(), nlsolve = NonlinearSolveAlg(NewtonRaphson()))on a sparse-jac_prototypebrusselator dies withFieldError: Matrix has no field nzval;Fix
Pass
jac_prototype = similar(W_for_reuse)so the inner Jacobian mirrors W's type and structure; thecopyto!closure is then structure-preserving and sparse factorizations work.Verified on the 512-unknown sparse brusselator (before → after): inner J
Matrix{Float64}→SparseMatrixCSC(nnz=3072), KLU crash → clean solve. Combined with the J/W-update split and SciML/NonlinearSolve.jl#1003,FBDF+NSA+KLUreaches 1.00× ofFBDF+NLNewton+KLUwall time (220.3 vs 219.5 ms) — from a crash/3.9× before.Tests
New
nsa_sparse_tests.jl: asserts the inner NonlinearSolve J isSparseMatrixCSCunder W-reuse, and that FBDF/TRBDF2 + NSA +KLUFactorizationsolve a sparse-jac_prototypeproblem to the reference (6/6 pass locally). Local suite runs on this branch (Julia 1.12): Newton 6/6, Linear Nonlinear 36/36, W-Operator Prototype, NSA Sparse 6/6, Sparse Algebraic 12/12, Sparse DAE Init 27/27, Linear Solver 113/113, Split ODE 10/10 — no failures (full Core group exceeds a 1h local window; remaining suites left to CI).Part of a series: #3816 (nsolve stat), #3818 (recompute_jacobian on W===nothing paths), SciML/NonlinearSolve.jl#1003 (refactorization reuse), and a J/W-update-split PR to follow.
🤖 Generated with Claude Code