GESVFactorization + dense LU pivot-buffer reuse (v4.2.0)#1078
Merged
ChrisRackauckas merged 6 commits intoJul 10, 2026
Conversation
…fers Warm refactorizations (cache.A = X then solve!) previously allocated a fresh ipiv vector and LU wrapper via lu!/lu on every A update — and, without alias_A, a fresh O(n^2) factors copy as well. When the cached LU holds size-matched buffers and the matrix is a strided BLAS type with RowMaximum pivoting, factorize through LAPACK.getrf! with the cached ipiv (copying A into the cached factors first when the user's matrix must stay intact). The preallocated-ipiv getrf! method exists on Julia >= 1.11 only; older releases keep the previous allocating path. Best-of-N warm-loop measurements (n = 5/10, Julia 1.12, OpenBLAS single-threaded): alias_A=true refactorization solve drops from 96/144 B to 0 B per call (1.10 -> 1.08 us / 2.72 -> 2.71 us), and alias_A=false from 368/1088 B to 0 B (1.15 -> 1.10 us / 2.89 -> 2.74 us). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
….2.0) New exported dense algorithm mirroring the classic dgesv workflow: on a fresh matrix the right-hand side is copied into u and LAPACK.gesv! factorizes and solves in a single call; the (factors, ipiv) pair is cached so b-only re-solves are one allocation-free LAPACK.getrs! triangular solve. Batched (matrix) right-hand sides work natively through both calls. With alias_A = true the factorization overwrites cache.A directly; otherwise a cacheval-owned buffer is refilled by copyto! and the user's matrix stays pristine. Singular matrices return ReturnCode.Failure instead of leaking the LAPACK throw, matching the existing BandedMatrix-LU catch pattern. Non-strided or non-BlasFloat operators get an informative ArgumentError at init. Warm-cache measurements (n = 5/10, Julia 1.12, OpenBLAS single-threaded) against a raw LAPACK.gesv! call including the matrix copy: refactorize-every-call solves land at parity (1.00 vs 0.93 us, 2.55 vs 2.57 us) and b-only re-solves at 0.60/1.47 us with zero allocations. In the ExponentialUtilities exp! Pade benchmark it is bit-identical to the raw gesv! path and within -1.8% to +2.3% of its runtime for n = 5..250. Tests cover vector/matrix RHS, cache reuse with new A and new b, both alias_A semantics, singular Failure retcodes, allocation bounds, and the unsupported-type error. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
solve! combined getrf!/getrs! instead of gesv! + try/catch: getrf! returns info > 0 for a singular U without raising (gesv! calls chklapackerror and throws), so a singular denominator becomes ReturnCode.Failure directly with no exception on the path. Warm b-only re-solves and the alias_A in-place contract are unchanged. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Follow @ChrisRackauckas's review: instead of catching (or manually info-checking) a throwing LAPACK call, factorize with the idiomatic `lu!(A; check = false)` and branch on `issuccess`, solving with `ldiv!` — the same no-throw pattern the other dense factorizations use. A singular factor becomes ReturnCode.Failure with no exception on the path (and, matching LUFactorization, leaves isfresh set so a re-solve refactorizes). The cached LU makes warm b-only re-solves allocation-free; fresh-A solves stay at raw-gesv! speed (0.62 vs 0.61 µs at n=5, 9.8 vs 10.2 at n=30). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
StableRNGs is not in the test target, so the Downgrade Core job errored
at load ("Package StableRNGs not found"). Seed with MersenneTwister from
the already-available Random stdlib; the rng-threaded rand(rng, ...) calls
and determinism are unchanged.
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Revert the previous stopgap (MersenneTwister) and keep StableRNGs — its cross-version-stable stream is the point — by adding it properly to the [extras], the test target, and [compat = "1"]. Verified the full test target resolves and test/Core/gesv.jl passes (24/24). Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
ChrisRackauckas
marked this pull request as ready for review
July 10, 2026 06:45
This was referenced Jul 10, 2026
ChrisRackauckas
added a commit
that referenced
this pull request
Jul 11, 2026
…1082) The alias_A = true dense LUFactorization refactorization path reuses the cached ipiv through LAPACK.getrf! since #1078, but only for strided BLAS matrices with RowMaximum pivoting; NoPivot/RowNonZero pivoting and non-BLAS element types still allocated a fresh pivot vector and LU wrapper through lu! on every cache.A = X refactorization. Route those through the vendored generic_lufact! with the cached pivot vector, which accepts a preallocated ipiv on every supported Julia version (unlike the getrf! method, which needs >= 1.11). 51x51 Matrix{Float64} warm refactorization (cache.A = X; solve!) with alias_A = true, Julia 1.12.6: NoPivot drops 480 B -> 0 B per refactorization (RowMaximum and the default algorithm were already 0 B via #1078; pre-#1078 all paths allocated 480 B). On Julia 1.10 the NoPivot path drops to 0 B as well, while the BLAS paths keep the previous allocating lu! there by design. Adds test/Core/lu_refactorization.jl covering: zero-allocation warm refactorization cycles (RowMaximum, NoPivot, and the default algorithm), pivot-vector identity across refactorizations, correctness against A \ b across cycles, singular refactorization returning ReturnCode.Failure without throwing and recovering on the next nonsingular matrix, the default algorithm's QR safety fallback surviving a warm singular refactorization, BigFloat (generic eltype) correctness, and pivot reallocation after resize!. Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
ChrisRackauckas
added a commit
that referenced
this pull request
Jul 11, 2026
…ve! (#1081) Fixes the QA group's ExplicitImports no_implicit_imports failure introduced in c75997b (#1078), which added an unqualified issuccess call. Matches the qualified use elsewhere in the file. Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Two commits that were pushed to the #1074 branch after it merged, so they never appeared in that PR — rebased onto main:
GESVFactorization(new exported alg, docs + release notes): fresh A → one combinedLAPACK.gesv!intocache.u(honorsalias_A; non-aliased users keep A pristine via a cacheval buffer); warm b-only re-solves → zero-allocgetrs!on the stored factors; batched matrix RHS native; singular →ReturnCode.Failure, no throw. 24 tests intest/Core/gesv.jl. Profiled: warm b-only re-solve 0.60 µs vs 0.93 µs for a rawgesv!at n=5; in the ExponentialUtilitiesexp!harness it restores raw-LAPACK parity at small n (−1.8% at n=5, ties at n=10–250).LUFactorizationpivot/factors-buffer reuse on refactorization via preallocated-ipivgetrf!— 0 B warm refactorization on Julia ≥ 1.11 (inert on 1.10, which lacks the method; old path kept).Version → 4.2.0 (new public API). First consumer: SciML/ExponentialUtilities.jl#240 (compat
"4.2").🤖 Generated with Claude Code