Fix singular handling in the StaticArrays square fast path with SVD rescue for the default

[ChrisRackauckas-Claude]

Fixes #1084

Note: This PR should be ignored until reviewed by @ChrisRackauckas.

Problem

The StaticArrays fast path in solve had no singular handling for square matrices:

using LinearSolve, StaticArrays, LinearAlgebra
As = SA[0.0 0.0; 1.0 -1.0]
sol = solve(LinearProblem(As, SMatrix{2,2}(1.0I)))
sol.retcode  # Success   ← wrong
sol.u        # [Inf NaN; Inf NaN]

For N ≤ 3, StaticArrays' \ uses direct inverse formulas that silently produce Inf/NaN on singular input, so the fast path returned Success with a non-finite solution. For N ≥ 4, \ is lu(A) \ b with check = true, so the default path threw SingularException instead. Explicit LUFactorization() (lu(prob.A) \ prob.b) also threw. (The init/solve! cache route was already fine — the generic LU machinery reports Failure there.)

Fix

Two parts, both confined to the static fast path in src/common.jl:

1. Explicit static LUFactorization reports failure. The LUFactorization branch now factors with lu(A, check = false) and checks LinearAlgebra.issuccess (which StaticArrays defines for its LU by checking the U diagonal — the same criterion as LAPACK's info). On failure it returns retcode = ReturnCode.Failure with the zero-initialized u, exactly matching the dense LUFactorization behavior on singular input (verified: dense returns Failure with u = [0.0, 0.0]), instead of throwing. issuccess is used rather than a finiteness check on u because it tests the factorization itself, matching dense semantics (a nonsingular A with non-finite b still reports Success with the non-finite u, on both paths), and it is also cheaper (N diagonal checks, and no wasted triangular solves on failure — StaticArrays' triangular \ throws on a zero diagonal, so solving first isn't even an option here).

2. Static square default gets an SVD rescue. solve(::StaticLinearProblem, ::Nothing) — the static default — now detects LU failure the same way and rescues with svd(A) \ b, returning the finite min-norm least-squares pseudo-solution with retcode = Success. This mirrors the dense default's singular-LU → pivoted-QR safetyfallback, which likewise returns the rescue solve's retcode (Success, verified on the dense reproducer).

Why SVD and not QR

qr(::SMatrix) \ b least-squares is not defined in StaticArrays — the same reason defaultalg(::SMatrix) already uses SVDFactorization() for the non-square branch ("QR(...) \ b is not defined currently"). StaticArrays' svd \ is tolerance-thresholded, so it gives the min-norm pseudo-solution: the rescue result matches pinv(Matrix(A)) * b to the last ulp, consistent with the intent of the dense pivoted-QR rescue.

Design choice

The check lives directly in the static fast-path methods rather than in a new static default-algorithm struct or the dense DefaultLinearSolver machinery. Default provenance is already explicit in the fast path's dispatch: solve(prob, nothing) is the static default, and solve(prob, alg::SciMLLinearSolveAlgorithm) is the explicit-algorithm path, so no provenance flag is needed. The DefaultLinearSolver cacheval union machinery is built around mutable caches and in-place factorization reuse and is not static-friendly; forcing static arrays through it would sacrifice the zero-allocation fast path for no benefit.

Two implementation subtleties, both commented in the source:

• For N ≤ 3 the nonsingular solve stays A \ b (direct inverse formulas) so results remain bit-identical to \; for N ≥ 4 it becomes lu(A, check = false) \ b, which is the identical computation to A \ b minus the throw.
• The plain \ kernel sits behind a @noinline barrier: StaticArrays' small-size formulas contain muladds whose FMA contraction is inlining-context dependent, and inlining them next to the new finiteness check changed the result in the last bit. The barrier restores bit-identity (verified === for N = 1…5, vector and matrix RHS, default and explicit LU).

Measured behavior

Before (main, c4ebafe):

case	result
singular 2×2, default	Success, u = [Inf NaN; Inf NaN]
singular 4×4, default	throws SingularException
singular static, LUFactorization()	throws SingularException

After:

case	result
singular 2×2/4×4, default, SVector & SMatrix RHS	Success, finite u ≈ pinv(Matrix(A)) * b (min-norm)
singular static, LUFactorization()	Failure, zero u, no throw
nonsingular, N = 1…5, both RHS kinds	bit-identical (===) to A \ b / lu(A) \ b
non-square static	unchanged

Performance (BenchmarkTools minimum, @inferred passes, AllocCheck @check_allocs passes — the whole path including the rescue branch is statically allocation-free):

case	main	this PR
default 2×2	5.6 ns, 0 allocs	11.8 ns, 0 allocs (@noinline barrier)
default 5×5	205 ns, 0 allocs	248 ns, 0 allocs
LUFactorization() 2×2 / 5×5	20.2 / 252 ns, 0 allocs	20.4 / 250 ns, 0 allocs
singular 2×2 rescue	— (was wrong/threw)	1.47 µs, 0 allocs

Tests

New "StaticArray Singular" testset in test/Core/retcodes.jl: singular 2×2 and 4×4 with SVector and SMatrix RHS (default → finite min-norm ≈ pinv, correct retcode; explicit LU → Failure, finite u, no throw), nonsingular bit-identity (===) for both sizes and both RHS kinds against \ and lu(A) \, and a non-square problem unchanged. Ran locally: GROUP=Core Pkg.test() and test/AD/static_arrays.jl (AllocCheck/@inferred static tests, 23 pass / 1 pre-existing broken).

Notes

• Interaction with LinearSolve 5.0: lightweight solutions — stop populating LinearSolution.cache in solve! returns #1083 (LinearSolve 5.0 lightweight solution return): this change is deliberately confined to the three static fast-path methods in src/common.jl (plus tests), touching only how u/retcode are computed before build_linear_solution, so whichever of the two merges second rebases trivially.
• Downstream: Quasi-Newton J⁻¹ init as a linsolve workspace: factorization + identity solve through the LinearSolve default NonlinearSolve.jl#1039 keeps a pinv special case for SMatrix precisely because of this bug; once this merges, that case can route through the LinearSolve static path (follow-up, not done here).

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Done — the barrier is removed and the small-N path inlines (4d5001f2 → new head). The nonsingular test now asserts a 4-ulp isapprox against A \ b instead of bit-identity (the explicit-LU assertion stays exact — it is lu(A) \ b either way). Honest perf note: 2×2 default measured 13.2 ns inlined vs 11.8 ns with the barrier vs 5.6 ns on main — the singularity check dominates, not the function call, so the barrier was buying nothing. Retcodes testset re-run post-change: 61 pass / 1 pre-existing broken.

[ChrisRackauckas-Claude]

Correction on the perf numbers + answers to the review questions, with a fair same-session A/B (the earlier "5.6 ns on main" was measured in a different session/machine state — my error propagating it):

	2×2 default solve (API)
main	12.8 ns
this PR (inlined, 7f7fdee6)	13.2 ns

True cost of the singular handling: ~0.4 ns — the LinearProblem/solution wrapper (~9 ns) dominates at this size regardless.

Component profile (raw kernels, N=2/4/8): all(isfinite, u) is free (fuses with the solve: 4.0 vs 4.4 ns bare); issuccess on an already-computed static LU is free (13.6 vs 13.7); the rescue branch costs ~2 ns at N=2 at kernel level and vanishes by N=8 — invisible through the API.

So: no check-skipping LUFactorization dispatch is warranted — the check costs nothing measurable, and explicit static LUFactorization is slower than the default for a different reason entirely (static lu() costs 3× the direct inverse formulas at N=2: 13.7 vs 4.4 ns raw). Users wanting maximal small-N speed should use the default. Also tried @inline on the kernel (worse: 21.2 ns — drags the cold rescue into the wrapper) and outlining the rescue tail (no change: 13.2) — the committed fully-inline version stands.

[ChrisRackauckas-Claude]

Added the GESVFactorization static dispatch per review. It previously fell through to the generic init path (dense-strided-only → rejected SMatrix outright). Semantics mirror the dense gesv driver exactly: one-shot factorize-and-solve, singular input reported as Failure through the retcode with zeroed u, no SVD rescue (that's the default algorithm's behavior, which gesv semantics don't include). N ≤ 3 uses the direct small-size formulas + the measured-free checks; larger N uses lu(check = false) with the failure branch never touching StaticArrays' throwing triangular solve. Non-square → informative ArgumentError. @inferred-stable on both branches.

Measured 2×2 through the API: gesv 16.4 ns vs default 12.8 vs explicit LU 21.7 — faster than LU (skips the static lu() at small N), slightly above the default (elseif-chain + tuple plumbing); its value is the semantics, not raw speed. Retcodes testset: 71 pass / 1 pre-existing broken (10 new assertions: nonsingular 2×2/4×4 vector+matrix RHS, singular → Failure both size regimes, non-square throw).

[ChrisRackauckas-Claude]

Overhead decomposition for the static path (2×2, prebuilt LinearProblem, all 0 allocations):

algorithm	ns	semantics on singular
bare A \ b	4.4	Inf/NaN (N≤3) or throw (N≥4)
DirectLdiv!()	5.6	same as \ — the unchecked zero-overhead path (~1 ns = one alg-type branch)
GESVFactorization()	~8-9	Failure retcode, no rescue
default (nothing)	9.2	Failure-safe + SVD min-norm rescue

(LinearProblem construction itself costs 4.4 ns — earlier quoted numbers included it; per-solve marginal cost on a prebuilt problem is the column above.) So the tiering is: DirectLdiv! when you know your matrix is fine and want raw \ speed, gesv for failure reporting, default for full safety — all allocation-free.

[ChrisRackauckas-Claude]

Documented the StaticArrays behavior (1e13cb1c): new docs/src/tutorials/static_arrays.md covering the non-caching static fast path (0-allocation direct solve; immutable ⇒ no init!/solve! interface), the algorithm tier with singular-input behavior and measured overhead —

algorithm	~2×2	singular A
DirectLdiv!	5.6 ns	unchecked (Inf/NaN or throw)
GESVFactorization	8–9 ns	Failure, no rescue
LUFactorization	22 ns	Failure
default	9.2 ns	SVD min-norm rescue, Success

— plus non-square SVD routing, and a note in the GESVFactorization docstring about its static dispatch. The @example blocks execute during the build (they render correctly).

While there I also fixed a pre-existing Documentation CI failure unrelated to this PR: solvers/solvers.md (a ~970-line auto-listed solver catalog) exceeds Documenter's default 200 KiB per-page HTML limit, which has been failing the docs build on main (confirmed on c4ebafe and 571d0a8). Raised size_threshold to 500 KiB — full docs build is now green locally (exit 0, no cross-ref/threshold errors).