Allow matrix-backed operator Jacobians in LHL defaults

[ChrisRackauckas-Claude]

Draft — please ignore until reviewed by @ChrisRackauckas.

What changed and why

Now that OrdinaryDiffEq marks a WOperator stale when an operator-backed Jacobian moves in place, LinearSolve can remove the consumer-side exclusion added in #1215. The LHL default predicate now inspects the matrix backing a MatrixOperator, while matrix-free operators still fall through to Krylov.

DefaultLinearSolver also aliases WOperator inputs. The staleness flag belongs to the wrapper, so copying the wrapper at init disconnected the producer's mark_jacobian_updated! call from the wrapper observed by the LHL cache even though both wrappers shared the same matrix data.

The regression covers default selection, LHL cache initialization, wrapper identity, the first solve, and a marked in-place Jacobian update. The LHL docstring and shifted-systems tutorial describe the expanded default.

Failing before / passing after

Test-only state on unmodified main:

$ julia +1.10 --project -e 'using Pkg; Pkg.instantiate(); include("test/Core/lhl.jl")'
a dense operator Jacobian uses LHL by default: Test Failed
  LinearSolve._lhl_defaultable(W, assump)
  false

a dense operator Jacobian uses LHL by default: Test Failed
  defaultalg(...) == DefaultLinearSolver(LHLFactorization)
  KrylovJL_GMRES == LHLFactorization

a dense operator Jacobian uses LHL by default: Test Failed
  cache.cacheval.LHLFactorization !== nothing
  nothing !== nothing

Test Summary:                                 | Pass  Fail  Total
 a dense operator Jacobian uses LHL by default |    2     3      5
ERROR: Some tests did not pass: 2 passed, 3 failed

Widening the predicate alone exposed the second half of the bug: the default selected LHL, but the copied wrapper did not receive the staleness signal, so the updated-J assertion failed with the old solution (4 passed, 1 failed). After preserving WOperator identity:

Test Summary:                                 | Pass  Total  Time
a dense operator Jacobian uses LHL by default |    6      6  2.6s
Test Summary:     | Pass  Total     Time
LHL Factorization |  204    204  1m05.0s

Verification

$ GROUP=Core julia +1.10 --project -e 'using Pkg; Pkg.test()'
Test Summary:     | Pass  Total     Time
LHL Factorization |  204    204  1m05.0s
Testing LinearSolve tests passed

$ JULIA_PKG_PRECOMPILE_AUTO=0 GROUP=QA julia +1.10 --project -e 'using Pkg; Pkg.test()'
Test Summary: | Pass  Broken  Total
JET Tests     |   20      23     43
Test Summary: | Pass  Total
Allocation QA |   55     55
Test Summary:              | Pass  Total
SupernodalLU Allocation QA |    6      6
Test Summary:     | Pass  Total
Quality Assurance |   48     48
Testing LinearSolve tests passed

$ JULIA_CONDAPKG_BACKEND=Null JULIA_PYTHONCALL_EXE="$PWD/.docs-venv/bin/python" julia +1.12 --project=docs docs/make.jl
[ Info: Doctest: running doctests.
[ Info: CheckDocument: running document checks.
[ Info: RenderDocument: rendering document.
[ Info: HTMLWriter: rendering HTML pages.

$ julia -m Runic --check src/default.jl src/lhl.jl test/Core/lhl.jl
$ git diff -- src/default.jl src/lhl.jl test/Core/lhl.jl docs/src/tutorials/shifted_systems.md | typos -
$ git diff --check

The first standard QA invocation spent its full one-hour budget in automatic Reactant extension precompilation before any QA assertion ran. Disabling only automatic precompilation allowed the unchanged QA group to execute and pass. The docs environment likewise needed a temporary Python environment for its already-declared PyAMG dependency after the machine-local Pixi environment failed to load Python's encodings module; the actual Documenter build then completed.

CI

At commit f8276dd69e3a616abfc5fa2875fe698688ad1f18, 60 checks passed and one generated matrix entry was skipped. The sole failure was ModelingToolkit's downstream All suite:

MethodError: no method matching create_array(::Type{CasADi.MX}, ::Nothing, ::Val{1}, ...)
ERROR: Some tests did not pass: 114 passed, 0 failed, 12 errored, 1 broken.

The identical CasADi failure occurs on unmodified upstream main at 85b24a13f17cf6c0813a34c3f906f7cac92c527c, so this branch did not introduce it. A clean local checkout reproduced the same first Optimization error with LinearSolve 5.15.0, ModelingToolkit 11.40.1, CasADi 1.3.0, and SymbolicUtils 4.46.1 (1 passed, 1 errored).

Commit-level boundary testing identified CasADi 54f83b489d1d99afd53a7571d94e5e66d03bb1c7 as good and b041ed63a5beb49f24620cd97de0e35c514c3115 as the first bad commit. CasADi PR 42 already restores the removed SymbolicUtils create_array integration. With its merge commit a879b555aab7ea46ab24409c49025d7b99f15f35, the actual previously failing ModelingToolkit add_solve_constraints! path returned Opti and passed. A later full-file local run hit an unrelated pip-CasADi/Julia-Ipopt ABI symbol error, so the entire file was not verified locally with the fix.

No duplicate issue was opened. CasADi's draft 1.3.1 release PR is the existing delivery path and currently has eight passing checks.

• PR CI failure: https://github.com/SciML/LinearSolve.jl/actions/runs/33292439860/job/99206322545
• Matching upstream-main failure: https://github.com/SciML/LinearSolve.jl/actions/runs/33282762136/job/99180607405

Not verified

• Real CUDA, AMDGPU, or Metal execution locally; the CI GPU job passed.
• GROUP=Everything locally; CI ran the declared version, platform, sublibrary, and downstream matrix, with only the pre-existing ModelingToolkit failure documented above.
• An end-to-end OrdinaryDiffEq solve against this branch; the producer-side staleness regression is in Mark an operator Jacobian as updated when calc_W! moves it OrdinaryDiffEq.jl#4354.

Reviewer decision

This default is safe only when the producer honors the WOperator staleness contract. OrdinaryDiffEq master now does, but older OrdinaryDiffEq releases do not, and LinearSolve cannot express a compat bound on a downstream package. Release sequencing or an explicit capability signal may be preferable to immediately shipping the widened default.

Aliasing all WOperator inputs in DefaultLinearSolver is intentional because its eligible branches are operator algorithms (LHL or Krylov) that need the live operator. Please push back if copying a WOperator is expected public behavior for another default-solver path.

Links

• OrdinaryDiffEq issue: calc_W! never marks an operator Jacobian as updated, so a caching linear solver cannot see J move OrdinaryDiffEq.jl#4303
• Producer-side staleness fix: Mark an operator Jacobian as updated when calc_W! moves it OrdinaryDiffEq.jl#4354
• Original LHL default and consumer-side exclusion: Add LHLFactorization: O(n²) re-shift for repeated solves of I - γJ #1215
• CasADi regression commit: SciML/CasADi.jl@b041ed6
• CasADi fix PR: Restore SymbolicUtils code generation support CasADi.jl#42
• CasADi 1.3.1 release PR: Release v1.3.1 CasADi.jl#44

🤖 Generated with Codex CLI 0.151.0 (model: unknown; session: local session ID 01a0507d-8515-73c0-841a-d7231caf221a).