Fix singular handling in the StaticArrays square fast path with SVD rescue for the default (#1085)

* Fix silent Success with non-finite u for singular static square matrices

The StaticArrays fast path in `solve` returned `retcode = Success` with
`u = [Inf NaN; ...]` (N <= 3, direct inverse formulas) or threw a
`SingularException` (N >= 4, checked `lu`) when the square matrix was
singular (#1084).

Now the static square default detects LU failure via
`lu(A, check = false)` + `issuccess` and rescues with an SVD
least-squares solve, returning the finite min-norm pseudo-solution with
`Success` — mirroring the dense default's LU -> pivoted-QR
`safetyfallback`. SVD is used instead of QR because
`qr(::SMatrix) \ b` least-squares is not defined in StaticArrays (the
same reason `defaultalg(::SMatrix)` uses `SVDFactorization()` for
non-square).

Explicit `LUFactorization()` on singular static input now returns
`retcode = Failure` with the zero-initialized `u` instead of throwing,
matching the dense `LUFactorization` behavior.

Nonsingular solves stay bit-identical to `\` and the fast path remains
allocation-free (AllocCheck passes); non-square static problems are
unchanged.

Fixes #1084

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

* Let the small-N static solve inline (review)

Drops the @noinline bit-identity barrier: StaticArrays small-size
formulas have inlining-context-dependent FMA contraction, so
nonsingular N <= 3 results may differ from a bare A \ b in the last
bit; the test asserts a 4-ulp match instead. Measured 2x2 default
solve: 13.2 ns inlined vs 11.8 ns with the barrier — the singularity
check dominates, not the call, so the barrier bought nothing.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

* Static GESVFactorization dispatch: one-shot solve, Failure on singular (review)

GESVFactorization previously fell through to the generic init path,
whose dense-only strided requirement rejects SMatrix. The static
dispatch mirrors the dense gesv driver semantics: a one-shot
factorize-and-solve whose singular failures are reported through the
return code — no SVD rescue (that is the default's behavior, which
gesv semantics do not include). N <= 3 uses the direct small-size
formulas with the (measured-free) finiteness/issuccess check; larger N
uses lu(check = false), never touching the throwing triangular solve
on failure. Non-square input gets an informative ArgumentError.

Measured 2x2 API: gesv 16.4 ns, default 12.8 ns, LUFactorization
21.7 ns.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

* docs: StaticArrays tutorial + fix pre-existing docs build failure

Add a StaticArrays tutorial page documenting the non-caching static
fast path (0-allocation direct solve, immutable so no init!/solve!),
the algorithm tier (DirectLdiv! / GESVFactorization / LUFactorization /
default) with their singular-input behavior and measured overhead, and
non-square SVD routing. Note the static GESV dispatch in the
GESVFactorization docstring.

Also raise the HTML size_threshold to 500 KiB: solvers.md (a long
auto-listed solver catalog) exceeds the default 200 KiB limit, which
has been failing the Documentation CI on main independently of this PR.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

---------

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

Author: 3 people

docs/Project.toml

@@ -4,15 +4,16 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 LinearSolveAutotune = "67398393-80e8-4254-b7e4-1b9a36a3c5b6"
 LinearSolvePyAMG = "7a56c47d-7ab1-4e99-b0e3-2952e463d64a"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[sources]
+LinearSolve = {path = ".."}
+LinearSolveAutotune = {path = "../lib/LinearSolveAutotune"}
+LinearSolvePyAMG = {path = "../lib/LinearSolvePyAMG"}
 
 [compat]
 Documenter = "1"
 LinearSolve = "4"
 LinearSolveAutotune = "1.1"
 LinearSolvePyAMG = "1"
 SciMLOperators = "1"
-
-[sources]
-LinearSolve = {path = ".."}
-LinearSolveAutotune = {path = "../lib/LinearSolveAutotune"}
-LinearSolvePyAMG = {path = "../lib/LinearSolvePyAMG"}

docs/make.jl

@@ -30,7 +30,11 @@ makedocs(
     ],
     format = Documenter.HTML(
         assets = ["assets/favicon.ico"],
-        canonical = "https://docs.sciml.ai/LinearSolve/stable/"
+        canonical = "https://docs.sciml.ai/LinearSolve/stable/",
+        # solvers.md is a long auto-listed solver catalog that exceeds the
+        # default 200 KiB per-page HTML limit; it is reference material meant
+        # to be read in one page.
+        size_threshold = 500 * 1024
     ),
     pages = pages
 )

docs/pages.jl

@@ -9,6 +9,7 @@ pages = [
         "tutorials/partitionedsolvers.md",
         "tutorials/petsc_mpi.md",
         "tutorials/accelerating_choices.md",
+        "tutorials/static_arrays.md",
         "tutorials/gpu.md",
         "tutorials/autotune.md",
         "tutorials/eigenvalue.md",

docs/src/tutorials/static_arrays.md

@@ -0,0 +1,75 @@
+# Solving with StaticArrays
+
+`LinearProblem`s whose `A` and `b` are
+[StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) (`SMatrix`,
+`SVector`) take a dedicated non-caching fast path: the solve happens directly
+on the immutable arrays with **zero heap allocations**, and the returned
+solution's `u` is a static array of the matching shape.
+
+```@example static
+using LinearSolve, StaticArrays
+
+A = SA[2.0 1.0; 1.0 3.0]
+b = SA[1.0, 2.0]
+prob = LinearProblem(A, b)
+
+sol = solve(prob)
+sol.u
+```
+
+Because static problems are immutable, there is no `init!`/`solve!` caching
+interface for them — each `solve` is a self-contained direct solve. (There is
+also nothing to cache: no factorization object is retained.)
+
+## Algorithm choices and singular-input behavior
+
+Static square problems support a tier of algorithms trading safety against
+overhead. All of them are allocation-free; timings below are for a 2×2
+`Float64` system on a prebuilt `LinearProblem` (bare `A \ b` is 4.4 ns on the
+same machine) and scale with the usual method costs at larger sizes.
+
+| algorithm | ~2×2 cost | on singular `A` |
+|---|---|---|
+| `DirectLdiv!()` | 5.6 ns | unchecked: whatever `\` gives — non-finite values for `N ≤ 3`, a thrown `SingularException` for larger `N` |
+| `GESVFactorization()` | 8–9 ns | `ReturnCode.Failure` with zeroed `u`; no rescue (one-shot factorize-and-solve, the static analog of LAPACK's `gesv` driver) |
+| `LUFactorization()` | 22 ns | `ReturnCode.Failure` with zeroed `u`, matching the dense `LUFactorization` behavior |
+| default (`solve(prob)`) | 9.2 ns | rescued: an SVD least-squares min-norm pseudo-solution with `ReturnCode.Success`, mirroring the dense default's singular-LU → pivoted-QR safety fallback |
+
+Guidance:
+
+  - Use the **default** unless you have a reason not to: it is nearly as fast
+    as the unchecked path (the finiteness check fuses with the solve) and a
+    singular matrix degrades gracefully instead of silently producing `Inf`/
+    `NaN`.
+  - Use **`DirectLdiv!`** when the matrix is known nonsingular and every
+    nanosecond counts: it is a bare `A \ b` behind one algorithm-type branch.
+  - Use **`GESVFactorization`** when you want failure *reported* (a checkable
+    return code) but not *repaired* — e.g. inside an iteration that has its own
+    recovery logic. It uses the direct small-size kernels, so it is faster than
+    `LUFactorization` at small sizes.
+  - `LUFactorization` on static arrays exists for algorithm-genericity (code
+    that passes an algorithm through to both dense and static problems); the
+    static `lu` costs about 3× the direct inverse formulas at `N ≤ 3`.
+
+`QRFactorization`, `CholeskyFactorization`, `NormalCholeskyFactorization`, and
+`SVDFactorization` also have direct static dispatches with their usual
+applicability requirements.
+
+Nonsingular results from the checked algorithms match `A \ b` up to the last
+bit or ulp-level differences from inlining-dependent FMA contraction in
+StaticArrays' small-size kernels.
+
+## Non-square static problems
+
+Non-square static problems route to an SVD least-squares solve by default
+(static QR least-squares division is not available in StaticArrays), returning
+the minimum-norm solution:
+
+```@example static
+Ans = SA[1.0 2.0; 3.0 4.0; 5.0 6.0]
+bns = SA[1.0, 2.0, 3.0]
+solve(LinearProblem(Ans, bns)).u
+```
+
+`GESVFactorization` is square-only and throws an `ArgumentError` for non-square
+static input, matching the LAPACK driver it mirrors.

src/common.jl

@@ -695,11 +695,68 @@ function SciMLBase.solve(prob::StaticLinearProblem, args...; kwargs...)
     return SciMLBase.solve(prob, nothing, args...; kwargs...)
 end
 
+"""
+    __static_default_ldiv(A, b)
+
+Solve for the static-array default algorithm. For square `A`, singular input is
+rescued with an SVD least-squares solve (min-norm pseudo-solution), mirroring
+the dense default's LU -> pivoted-QR `safetyfallback`. SVD is used instead of
+QR because `qr(::SMatrix) \\ b` least-squares is not defined in StaticArrays
+(see `defaultalg(::SMatrix)`).
+"""
+function __static_default_ldiv(A::SMatrix{N, N}, b) where {N}
+    # StaticArrays' square `\` uses direct inverse formulas for N <= 3 (never
+    # throws; singular input silently yields Inf/NaN) and `lu(A) \ b` with
+    # `check = true` (throws SingularException) for larger sizes. Keep `A \ b`
+    # for N <= 3 so nonsingular results stay bit-identical to `\`; for larger
+    # sizes `lu(A, check = false) \ b` is the same computation as `A \ b`
+    # without the throw.
+    if N <= 3
+        # Fully inlined: StaticArrays' small-size formulas have
+        # inlining-context-dependent FMA contraction, so results may differ
+        # from a bare `A \ b` in the last bit.
+        u = A \ b
+        # Only rescue on factorization failure: a nonsingular `A` with
+        # non-finite `u` (e.g. non-finite `b`) returns `u` as-is, matching the
+        # dense LU behavior.
+        if all(isfinite, u) || LinearAlgebra.issuccess(lu(A, check = false))
+            return u
+        end
+    else
+        F = lu(A, check = false)
+        if LinearAlgebra.issuccess(F)
+            return F \ b
+        end
+    end
+    return svd(A) \ b
+end
+__static_default_ldiv(A, b) = A \ b
+
+function __static_gesv_ldiv(A::SMatrix{N, N}, b) where {N}
+    if N <= 3
+        u = A \ b
+        ok = all(isfinite, u) || LinearAlgebra.issuccess(lu(A, check = false))
+        return u, ok
+    else
+        F = lu(A, check = false)
+        ok = LinearAlgebra.issuccess(F)
+        # StaticArrays' triangular solve throws on a zero pivot, so the failed
+        # branch must not touch `F \ b`. The placeholder mirrors `F \ b`'s
+        # promoted eltype to keep the return type-stable; the caller replaces
+        # `u` on failure.
+        u = ok ? F \ b : zero(b) / oneunit(eltype(A))
+        return u, ok
+    end
+end
+function __static_gesv_ldiv(A, b)
+    throw(ArgumentError("GESVFactorization requires a square matrix, got size $(size(A))"))
+end
+
 function SciMLBase.solve(
         prob::StaticLinearProblem,
         alg::Nothing, args...; kwargs...
     )
-    u = prob.A \ prob.b
+    u = __static_default_ldiv(prob.A, prob.b)
     return SciMLBase.build_linear_solution(
         alg, u, nothing, prob; retcode = ReturnCode.Success
     )
@@ -712,7 +769,30 @@ function SciMLBase.solve(
     if alg === nothing || alg isa DirectLdiv!
         u = prob.A \ prob.b
     elseif alg isa LUFactorization
-        u = lu(prob.A) \ prob.b
+        F = lu(prob.A, check = false)
+        if !LinearAlgebra.issuccess(F)
+            # Match dense `LUFactorization` on singular input: report Failure
+            # with the (zero-)initialized `u` instead of throwing.
+            u = prob.u0 !== nothing ? prob.u0 : __init_u0_from_Ab(prob.A, prob.b)
+            return SciMLBase.build_linear_solution(
+                alg, u, nothing, prob; retcode = ReturnCode.Failure
+            )
+        end
+        u = F \ prob.b
+    elseif alg isa GESVFactorization
+        # The static analog of the dense `gesv` driver: a one-shot
+        # factorize-and-solve whose singular failures are reported through the
+        # return code — no SVD rescue (that is the default algorithm's
+        # behavior, which gesv semantics do not include). Uses the direct
+        # small-size formulas, so this is the fastest explicit static
+        # algorithm at small N.
+        u, gesv_ok = __static_gesv_ldiv(prob.A, prob.b)
+        if !gesv_ok
+            u = prob.u0 !== nothing ? prob.u0 : __init_u0_from_Ab(prob.A, prob.b)
+            return SciMLBase.build_linear_solution(
+                alg, u, nothing, prob; retcode = ReturnCode.Failure
+            )
+        end
     elseif alg isa QRFactorization
         u = qr(prob.A) \ prob.b
     elseif alg isa CholeskyFactorization

src/factorization.jl

@@ -539,7 +539,11 @@ the user's matrix is left untouched.
 
 Only dense strided BLAS floating-point matrices
 (`StridedMatrix{<:Union{Float32, Float64, ComplexF32, ComplexF64}}`) are
-supported; other operator types throw an informative error at `init`.
+supported through the caching interface; other operator types throw an
+informative error at `init`. Square StaticArrays problems have a dedicated
+direct dispatch with the same semantics (one-shot solve, singular input
+reported as `ReturnCode.Failure`, no factorization retained) — see the
+StaticArrays tutorial page.
 """
 struct GESVFactorization <: AbstractDenseFactorization end
 

test/Core/retcodes.jl

@@ -94,3 +94,89 @@ staticarrayalgs = (
 
     @test_broken sol = solve(prob1, QRFactorization()) # Needs StaticArrays `qr` fix
 end
+
+@testset "StaticArray Singular" begin
+    A2 = SA[0.0 0.0; 1.0 -1.0]
+    A4 = SMatrix{4, 4}(
+        [
+            1.0 2.0 3.0 4.0
+            2.0 4.0 6.0 8.0
+            0.0 1.0 0.0 1.0
+            1.0 0.0 1.0 0.0
+        ]
+    )
+    for (A, b, B) in (
+            (A2, SA[1.0, 2.0], SMatrix{2, 2}(1.0I)),
+            (A4, SA[1.0, 2.0, 3.0, 4.0], SMatrix{4, 4}(1.0I)),
+        )
+        for rhs in (b, B)
+            # Default alg rescues singular LU with an SVD least-squares solve,
+            # returning the finite min-norm pseudo-solution (like the dense
+            # default's pivoted-QR rescue).
+            sol = solve(LinearProblem(A, rhs))
+            @test SciMLBase.successful_retcode(sol.retcode)
+            @test all(isfinite, sol.u)
+            @test sol.u isa (rhs isa SVector ? SVector : SMatrix)
+            @test sol.u ≈ pinv(Matrix(A)) * (rhs isa SVector ? Vector(rhs) : Matrix(rhs))
+
+            # Explicit LUFactorization reports Failure without throwing,
+            # matching the dense LUFactorization behavior on singular input.
+            sole = solve(LinearProblem(A, rhs), LUFactorization())
+            @test sole.retcode == ReturnCode.Failure
+            @test all(isfinite, sole.u)
+        end
+    end
+
+    # Nonsingular static solves match `\` (up to FMA contraction differences
+    # in StaticArrays' inlined small-size formulas, per review the check is
+    # allowed to inline rather than preserving bit-identity)
+    A2n = SA[2.0 1.0; 1.0 3.0]
+    A4n = SMatrix{4, 4}(
+        [
+            4.0 1.0 0.0 2.0
+            1.0 5.0 1.0 0.0
+            0.0 1.0 6.0 1.0
+            2.0 0.0 1.0 7.0
+        ]
+    )
+    for (A, b, B) in (
+            (A2n, SA[1.0, 2.0], SA[1.0 0.5; 2.0 1.5]),
+            (
+                A4n, SA[1.0, 2.0, 3.0, 4.0],
+                SMatrix{4, 2}([1.0 0.5; 2.0 1.5; 3.0 2.5; 4.0 3.5]),
+            ),
+        )
+        for rhs in (b, B)
+            @test solve(LinearProblem(A, rhs)).u ≈ A \ rhs rtol = 4 * eps()
+            @test solve(LinearProblem(A, rhs), LUFactorization()).u === lu(A) \ rhs
+        end
+    end
+
+    # GESVFactorization on static problems: the static analog of the dense
+    # gesv driver — one-shot solve, Failure through the retcode on singular,
+    # no SVD rescue.
+    for (A, b) in ((A2n, SA[1.0, 2.0]), (A4n, SVector{4}(1.0:4.0)))
+        sol = solve(LinearProblem(A, b), GESVFactorization())
+        @test sol.retcode == ReturnCode.Success
+        @test sol.u ≈ A \ b rtol = 4 * eps()
+    end
+    B2 = SMatrix{2, 2}(1.0I)
+    @test solve(LinearProblem(A2n, B2), GESVFactorization()).u ≈ A2n \ B2
+    As2 = SA[0.0 0.0; 1.0 -1.0]
+    As4 = SMatrix{4, 4}([1.0 2.0 3.0 4.0; 2.0 4.0 6.0 8.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0])
+    for (As, b) in ((As2, SA[1.0, 2.0]), (As4, SVector{4}(1.0:4.0)))
+        sol = solve(LinearProblem(As, b), GESVFactorization())
+        @test sol.retcode == ReturnCode.Failure
+        @test all(iszero, sol.u)
+    end
+    @test_throws ArgumentError solve(
+        LinearProblem(SMatrix{3, 2}(ones(3, 2)), SA[1.0, 2.0, 3.0]), GESVFactorization()
+    )
+
+    # Non-square static problems are unchanged by the singular handling
+    Ans = SMatrix{3, 2}([1.0 2.0; 3.0 4.0; 5.0 6.0])
+    bns = SA[1.0, 2.0, 3.0]
+    solns = solve(LinearProblem(Ans, bns))
+    @test SciMLBase.successful_retcode(solns.retcode)
+    @test solns.u ≈ Matrix(Ans) \ Vector(bns)
+end