Fall back to pivoted QR when the default hits a rank-deficient least-squares problem

[AJ0070]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Fixes #531.

• defaultalg picks an unpivoted QRFactorization for overdetermined dense systems. It is not rank-revealing, so a rank-deficient A returns cache.u untouched (all zeros) with ReturnCode.Failure, while A \ b gives the correct least-squares solution.
• Worse variant not in the issue: with a nearly dependent column, no exact zero shows up on the R diagonal and the result is finite, so a wrong answer is returned with a success retcode. A[:, 1] .= 1e-14 .* A[:, 2] gives [4.33e28, -4.33e14, 0.020, 0.305] where A \ b gives [0.0, 0.809, 0.100, 0.286].
• Pivoting unconditionally would fix both but taxes every full-rank solve: qr!(A, ColumnNorm()) measured 1.5x to 2.8x qr!(A, NoPivot()) on typical tall-skinny shapes here (1000x100: 2.80 ms vs 7.82 ms).
• Instead, unpivoted QR stays on the fast path and now hooks into the safety-fallback machinery that already exists for LU.
• _qr_rank_deficient(F) scans the R diagonal in O(min(m, n)) against min(m, n) * eps, the threshold LAPACK's xGELSY (and so A \ b) uses to truncate the rank. It dispatches on the factorization's storage, so SPQR, GPU, and Dual factorizations keep their current behavior.
• The check is a heuristic, since unpivoted QR does not order diag(R) by magnitude, and can miss a Kahan-style deficiency. The docstring says so.
• _default_qr_solve_with_fallback re-solves with QRFactorization(ColumnNorm()) on failure, non-finite output, or rank deficiency. It reuses _do_qr_fallback, so A restoration and fell_back_to_qr cache reuse work as they do for LU. Its guards are decided by the cache's type, so the branch folds away at compile time.
• Both the exactly-singular and numerically-singular cases now reproduce A \ b, @inferred solve! still returns a concrete LinearSolution, and full-rank problems never enter the fallback.
• An explicit QRFactorization() or LUFactorization() still reports ReturnCode.Failure on rank-deficient input, matching the conclusion in the issue thread. The QRFactorization docstring and solvers.md now say so and point at ColumnNorm(), SVDFactorization, and the least-squares Krylov methods.
• Two assertions in test/Core/default_algs.jl encoded the old failure behavior and now assert the corrected result. The second one's purpose (the default alg choice must not bake in a rank check, since rank can change across cache reuses) is unaffected, and it now also covers the fallback firing when rank drops mid-cache.
• New coverage in test/Core/nonsquare.jl: tall with a zero column, tall with a duplicated column, tall numerically rank-deficient, wide rank-deficient, square singular under VeryIllConditioned, cache reuse across the fallback, and a full-rank case leaving fell_back_to_qr unset.
• Full Core group passes locally on Julia 1.12.6 (25 testsets, 0 failures).

AI Disclosure: Used Opus 5

[ChrisRackauckas]

Is this GPU safe? Easier is to just check the blas return info and go based on that.

[AJ0070]

Yes on GPU safety, but you are right that it was resting on something too indirect. It relied on ArrayInterface.fast_scalar_indexing, whose false for GPU arrays comes from a weak-dependency extension. Pushed 69b6480 to dispatch on the factorization type instead, mirroring the GPU method on _notsuccessful just above it. Both methods constrain the same element type so the GPU one is strictly more specific, and detect_ambiguities reports none. GPU also never reaches the check anyway: the caller is guarded on cache.A isa DenseMatrix and _qr_fallback_pivot(cache.A) isa ColumnNorm, and that helper returns NoPivot() for GPU arrays.

On the BLAS info: there isn't one for unpivoted QR. On 1.12:

QRCompactWY fieldnames:               (:factors, :T)      # no info field
hasmethod(issuccess, (QRCompactWY,)): false
geqrf! on a matrix with a zero column: info == 0
R diagonal after geqrf!:              [0.0, -1.73, -1.33, -0.882]

geqrf/geqrt only set info < 0 for an illegal argument, never for rank deficiency, which is what geqp3/gelsy are for. LU is the contrasting case and does carry .info (3 on a singular test matrix) plus issuccess, which is what the existing LU fallback keys off.

That is why _notsuccessful(::QRCompactWY) directly above already hand-scans the R diagonal for an exact zero instead of reading an info. This is the same scan with a relative threshold (min(m, n) * eps * max|R[i,i]|, the rcond xGELSY uses) rather than exact equality, which is what catches the near-deficient case. That case is the one that currently returns [4.33e28, -4.33e14, ...] with ReturnCode.Success, so an exact-zero test alone would not close it.

The genuinely simpler option is pivoted QR by default, since geqp3 does give a rank estimate, but that measured 1.5x to 2.8x on typical tall-skinny shapes (1000x100: 2.80 ms vs 7.82 ms), which is why I kept the cheap check on the fast path instead.