From 4d1d4be64906d174c0ee56ae3ec44eefa901f057 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 17:06:47 -0400 Subject: [PATCH 1/6] Dense LUFactorization refactorization reuses cached pivot/factors buffers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 Co-Authored-By: Claude Fable 5 --- src/factorization.jl | 44 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 7d56dc7c0..6232213e1 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -263,6 +263,29 @@ GenericLUFactorization(pivot = RowMaximum(); residualsafety::Bool = false) = Gen _get_residualsafety(alg::LUFactorization) = alg.residualsafety _get_residualsafety(alg::GenericLUFactorization) = alg.residualsafety +# Dense-LU refactorization buffer reuse: `lu`/`lu!` allocate a fresh pivot +# vector (and, without aliasing, a fresh factors copy) on every `cache.A = X` +# refactorization. When the cached `LU`'s buffers match the incoming matrix, +# factorize through `LAPACK.getrf!` with the cached `ipiv` instead. The +# preallocated-`ipiv` `getrf!` method only exists on Julia >= 1.11; older +# releases keep the allocating path. +@static if VERSION >= v"1.11" + function _reusable_lu_cacheval(cacheval, A) + return cacheval isa LU && + cacheval.ipiv isa Vector{BlasInt} && + length(cacheval.ipiv) == min(size(A)...) + end + function _lu_reusing_ipiv!(A, ipiv::Vector{BlasInt}) + factors, piv, info = LAPACK.getrf!(A, ipiv; check = false) + return LU{eltype(factors), typeof(factors), typeof(piv)}( + factors, piv, convert(BlasInt, info) + ) + end +else + _reusable_lu_cacheval(cacheval, A) = false + _lu_reusing_ipiv!(A, ipiv) = error("unreachable: requires Julia >= 1.11") +end + function SciMLBase.solve!(cache::LinearCache, alg::LUFactorization; kwargs...) A = cache.A A = convert(AbstractMatrix, A) @@ -287,9 +310,26 @@ function SciMLBase.solve!(cache::LinearCache, alg::LUFactorization; kwargs...) ArrayInterface.can_setindex(typeof(A)) # The user permitted overwriting A (`alias_A = true` at `init`), # so refactorize in place and skip the O(n²) copy `lu` makes. - fact = lu!(A, _normalize_pivot(alg.pivot); check = false) + pivot = _normalize_pivot(alg.pivot) + if A isa StridedMatrix{<:LinearAlgebra.BlasFloat} && + pivot isa RowMaximum && _reusable_lu_cacheval(cacheval, A) + fact = _lu_reusing_ipiv!(A, cacheval.ipiv) + else + fact = lu!(A, pivot; check = false) + end else - fact = lu(A, check = false) + pivot = _normalize_pivot(alg.pivot) + if A isa StridedMatrix{<:LinearAlgebra.BlasFloat} && + pivot isa RowMaximum && _reusable_lu_cacheval(cacheval, A) && + cacheval.factors isa Matrix{eltype(A)} && + size(cacheval.factors) == size(A) && cacheval.factors !== A + # A must stay intact, but the previous factorization's + # buffers can be overwritten in place of `lu`'s fresh copy. + copyto!(cacheval.factors, A) + fact = _lu_reusing_ipiv!(cacheval.factors, cacheval.ipiv) + else + fact = lu(A, check = false) + end end catch e # Some matrix types (e.g. BandedMatrix) throw LAPACKException on singular From 8b071ac1bc3acb6c259b6634027247efe18d0473 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sun, 5 Jul 2026 17:07:02 -0400 Subject: [PATCH 2/6] GESVFactorization: LAPACK gesv-style combined factorize-and-solve (v4.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 Co-Authored-By: Claude Fable 5 --- docs/src/release_notes.md | 9 +++ docs/src/solvers/solvers.md | 1 + src/LinearSolve.jl | 3 +- src/factorization.jl | 88 ++++++++++++++++++++++++ test/Core/gesv.jl | 130 ++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 231 insertions(+), 1 deletion(-) create mode 100644 test/Core/gesv.jl diff --git a/docs/src/release_notes.md b/docs/src/release_notes.md index 47987c702..1a81d8d66 100644 --- a/docs/src/release_notes.md +++ b/docs/src/release_notes.md @@ -1,5 +1,14 @@ # Release Notes +## v4.2 + + - New `GESVFactorization` algorithm mirroring LAPACK's `gesv` driver: fresh matrices + factorize and solve in a single `LAPACK.gesv!` call, and repeat solves with only a new + `b` reuse the cached factors through an allocation-free `LAPACK.getrs!`. + - Dense `LUFactorization` refactorizations (`cache.A = X` then `solve!`) now reuse the + cached pivot vector (and, without `alias_A`, the cached factors buffer) on + Julia >= 1.11, making warm refactorization solves allocation-free. + ## v4.0 - Batched (matrix) right-hand sides are now supported: `solve(LinearProblem(A, B))` with diff --git a/docs/src/solvers/solvers.md b/docs/src/solvers/solvers.md index 58a9101bf..978acb15d 100644 --- a/docs/src/solvers/solvers.md +++ b/docs/src/solvers/solvers.md @@ -175,6 +175,7 @@ customized per-package, details given below describe a subset of important array LUFactorization GenericFactorization GenericLUFactorization +GESVFactorization QRFactorization SVDFactorization CholeskyFactorization diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 292fbd468..224531b17 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -518,7 +518,8 @@ is_cusparse_csr(A) = false is_cusparse_csc(A) = false export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization, - GenericLUFactorization, SimpleLUFactorization, RFLUFactorization, ButterflyFactorization, + GenericLUFactorization, GESVFactorization, SimpleLUFactorization, + RFLUFactorization, ButterflyFactorization, NormalCholeskyFactorization, NormalBunchKaufmanFactorization, UMFPACKFactorization, KLUFactorization, PureKLUFactorization, PureUMFPACKFactorization, SparseColumnPivotedQRFactorization, FastLUFactorization, diff --git a/src/factorization.jl b/src/factorization.jl index 6232213e1..ff1419d0d 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -507,6 +507,94 @@ function init_cacheval( return nothing end +## GESVFactorization + +""" +`GESVFactorization()` + +A direct mirror of LAPACK's `gesv` driver (`LAPACK.gesv!`): on a fresh matrix the +right-hand side is copied into `u` and the factorization and solve happen in a +single LAPACK call, matching the classic `dgesv` workflow. The `(factors, ipiv)` +pair is cached, so subsequent solves that only change `b` are a single +allocation-free `LAPACK.getrs!` triangular solve. Batched (matrix) right-hand +sides are handled natively by both calls. + +With `alias_A = true` the factorization overwrites `cache.A` directly; otherwise +a workspace-owned buffer is refilled with `copyto!` on each refactorization and +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`. +""" +struct GESVFactorization <: AbstractDenseFactorization end + +function init_cacheval( + alg::GESVFactorization, A::StridedMatrix{T}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, + assumptions::OperatorAssumptions + ) where {T <: LinearAlgebra.BlasFloat} + return (Matrix{T}(undef, 0, 0), Vector{BlasInt}(undef, 0)) +end + +function init_cacheval( + alg::GESVFactorization, A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, + assumptions::OperatorAssumptions + ) + throw( + ArgumentError( + "GESVFactorization only supports dense strided BLAS floating-point \ + matrices (StridedMatrix{<:Union{Float32, Float64, ComplexF32, \ + ComplexF64}}); got $(typeof(A)). Use LUFactorization or the default \ + algorithm instead." + ) + ) +end + +function SciMLBase.solve!(cache::LinearCache, alg::GESVFactorization; kwargs...) + A = convert(AbstractMatrix, cache.A) + F, ipiv = @get_cacheval(cache, :GESVFactorization) + if cache.isfresh + if cache.alias_A && A isa Matrix + # The user permitted overwriting A, so factorize it in place. + Atarget = A + else + if size(F) != size(A) + F = Matrix{eltype(A)}(undef, size(A)) + end + copyto!(F, A) + Atarget = F + end + copyto!(cache.u, cache.b) + try + # One combined factorize-and-solve LAPACK call, exactly like dgesv; + # gesv! returns the factors and pivots for later getrs! reuse. + _, factors, ipivnew = LAPACK.gesv!(Atarget, cache.u) + cache.cacheval = (factors, ipivnew) + catch e + if e isa LinearAlgebra.LAPACKException || + e isa LinearAlgebra.SingularException + @SciMLMessage("Solver failed", cache.verbose, :solver_failure) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Failure + ) + else + rethrow(e) + end + end + cache.isfresh = false + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Success + ) + end + copyto!(cache.u, cache.b) + LAPACK.getrs!('N', F, ipiv, cache.u) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Success + ) +end + ## QRFactorization """ diff --git a/test/Core/gesv.jl b/test/Core/gesv.jl new file mode 100644 index 000000000..5f15e8a8b --- /dev/null +++ b/test/Core/gesv.jl @@ -0,0 +1,130 @@ +using LinearSolve, LinearAlgebra, SparseArrays, Test +using StableRNGs + +rng = StableRNG(123) + +@testset "GESVFactorization" begin + @testset "vector and matrix RHS" begin + n = 20 + A = rand(rng, n, n) + 5I + b = rand(rng, n) + B = rand(rng, n, 4) + + sol = solve(LinearProblem(A, b), GESVFactorization()) + @test SciMLBase.successful_retcode(sol.retcode) + @test sol.u ≈ A \ b + + solm = solve(LinearProblem(A, B), GESVFactorization()) + @test SciMLBase.successful_retcode(solm.retcode) + @test solm.u ≈ A \ B + end + + @testset "cache reuse: new b (getrs! path) and new A (refactorization)" begin + n = 20 + A1 = rand(rng, n, n) + 5I + A2 = rand(rng, n, n) + 5I + b1 = rand(rng, n) + b2 = rand(rng, n) + + cache = init(LinearProblem(copy(A1), copy(b1)), GESVFactorization()) + sol1 = solve!(cache) + @test sol1.u ≈ A1 \ b1 + + cache.b = copy(b2) + sol2 = solve!(cache) + @test sol2.u ≈ A1 \ b2 + + cache.A = copy(A2) + sol3 = solve!(cache) + @test sol3.u ≈ A2 \ b2 + + # matrix RHS through the same cache-reuse paths + B1 = rand(rng, n, 3) + B2 = rand(rng, n, 3) + mcache = init(LinearProblem(copy(A1), copy(B1)), GESVFactorization()) + @test solve!(mcache).u ≈ A1 \ B1 + mcache.b = copy(B2) + @test solve!(mcache).u ≈ A1 \ B2 + mcache.A = copy(A2) + @test solve!(mcache).u ≈ A2 \ B2 + end + + @testset "alias_A semantics" begin + n = 12 + A = rand(rng, n, n) + 5I + b = rand(rng, n) + + # alias_A = false: the user's matrix must stay pristine + Akeep = copy(A) + cache = init( + LinearProblem(A, copy(b)), GESVFactorization(); + alias = LinearAliasSpecifier(alias_A = false, alias_b = false) + ) + sol = solve!(cache) + @test sol.u ≈ Akeep \ b + @test A == Akeep + + # alias_A = true: A may be overwritten by the factors, and the solve is + # still correct + A2 = rand(rng, n, n) + 5I + A2keep = copy(A2) + cache2 = init( + LinearProblem(A2, copy(b)), GESVFactorization(); + alias = LinearAliasSpecifier(alias_A = true, alias_b = false) + ) + sol2 = solve!(cache2) + @test sol2.u ≈ A2keep \ b + end + + @testset "singular matrix returns Failure retcode without throwing" begin + n = 6 + As = zeros(n, n) + b = rand(rng, n) + sol = solve(LinearProblem(As, b), GESVFactorization()) + @test !SciMLBase.successful_retcode(sol.retcode) + @test sol.retcode == ReturnCode.Failure + + # a cache whose fresh solve failed can be reused with a new nonsingular A + cache = init(LinearProblem(zeros(n, n), copy(b)), GESVFactorization()) + @test !SciMLBase.successful_retcode(solve!(cache).retcode) + Aok = rand(rng, n, n) + 5I + cache.A = copy(Aok) + solok = solve!(cache) + @test SciMLBase.successful_retcode(solok.retcode) + @test solok.u ≈ Aok \ b + end + + @testset "allocations: warm b-only re-solve is allocation-free" begin + n = 100 + A = rand(rng, n, n) + 10I + b = rand(rng, n) + bnew = rand(rng, n) + + resolve_b!(cache, bv) = (copyto!(cache.b, bv); solve!(cache); nothing) + refactor!(cache, Am) = (copyto!(cache.A, Am); cache.A = cache.A; solve!(cache); nothing) + + cache = init( + LinearProblem(copy(A), copy(b)), GESVFactorization(); + alias = LinearAliasSpecifier(alias_A = true, alias_b = true) + ) + solve!(cache) + resolve_b!(cache, bnew) + @test (@allocated resolve_b!(cache, bnew)) == 0 + @test cache.u ≈ A \ bnew + + # fresh-A re-solve stays O(n): gesv! allocates only its pivot vector + refactor!(cache, A) + alloc = @allocated refactor!(cache, A) + @test alloc < 20000 + @test cache.u ≈ A \ bnew + end + + @testset "unsupported matrix types throw an informative error at init" begin + A = sprand(rng, 10, 10, 0.5) + 5I + b = rand(rng, 10) + @test_throws ArgumentError init(LinearProblem(A, b), GESVFactorization()) + Abig = rand(rng, BigFloat, 4, 4) + 5I + bbig = rand(rng, BigFloat, 4) + @test_throws ArgumentError init(LinearProblem(Abig, bbig), GESVFactorization()) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 3dbf54472..c44370b6a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -63,6 +63,7 @@ else core = function () @time @safetestset "Basic Tests" include("Core/basictests.jl") @time @safetestset "Batched RHS" include("Core/batch.jl") + @time @safetestset "GESV Factorization" include("Core/gesv.jl") @time @safetestset "Return codes" include("Core/retcodes.jl") @time @safetestset "Re-solve" include("Core/resolve.jl") @time @safetestset "Zero Initialization Tests" include("Core/zeroinittests.jl") From 8d870fc0d1e0d3905f07dc65f62c8a63c06f5502 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 9 Jul 2026 17:44:09 -0400 Subject: [PATCH 3/6] GESVFactorization: report singular factor via retcode, not exception 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 Co-Authored-By: Claude Fable 5 --- src/factorization.jl | 31 ++++++++++++++----------------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index ff1419d0d..8549c7bbc 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -566,24 +566,21 @@ function SciMLBase.solve!(cache::LinearCache, alg::GESVFactorization; kwargs...) copyto!(F, A) Atarget = F end - copyto!(cache.u, cache.b) - try - # One combined factorize-and-solve LAPACK call, exactly like dgesv; - # gesv! returns the factors and pivots for later getrs! reuse. - _, factors, ipivnew = LAPACK.gesv!(Atarget, cache.u) - cache.cacheval = (factors, ipivnew) - catch e - if e isa LinearAlgebra.LAPACKException || - e isa LinearAlgebra.SingularException - @SciMLMessage("Solver failed", cache.verbose, :solver_failure) - return SciMLBase.build_linear_solution( - alg, cache.u, nothing, cache; retcode = ReturnCode.Failure - ) - else - rethrow(e) - end - end + # Equivalent to dgesv (factorize + solve), but split into getrf!/getrs! + # so a singular factor is reported through the return code rather than a + # thrown exception: getrf! returns `info > 0` for a singular U without + # raising (unlike gesv!, which calls chklapackerror). + factors, ipivnew, info = LAPACK.getrf!(Atarget) + cache.cacheval = (factors, ipivnew) cache.isfresh = false + if !iszero(info) + @SciMLMessage("Solver failed", cache.verbose, :solver_failure) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Failure + ) + end + copyto!(cache.u, cache.b) + LAPACK.getrs!('N', factors, ipivnew, cache.u) return SciMLBase.build_linear_solution( alg, cache.u, nothing, cache; retcode = ReturnCode.Success ) From c0d0e6386623701ef5724c52e341e0de05507dbb Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 9 Jul 2026 18:16:06 -0400 Subject: [PATCH 4/6] GESVFactorization: use lu!(;check=false)/ldiv!, no exception path MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 Co-Authored-By: Claude Fable 5 --- src/factorization.jl | 47 ++++++++++++++++++++------------------------ 1 file changed, 21 insertions(+), 26 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 8549c7bbc..3f2a1052f 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -512,12 +512,12 @@ end """ `GESVFactorization()` -A direct mirror of LAPACK's `gesv` driver (`LAPACK.gesv!`): on a fresh matrix the -right-hand side is copied into `u` and the factorization and solve happen in a -single LAPACK call, matching the classic `dgesv` workflow. The `(factors, ipiv)` -pair is cached, so subsequent solves that only change `b` are a single -allocation-free `LAPACK.getrs!` triangular solve. Batched (matrix) right-hand -sides are handled natively by both calls. +A dense LU factorize-and-solve in the style of LAPACK's `gesv` driver, tuned for +repeated solves. On a fresh matrix it factorizes with `lu!(A; check = false)` +(so a singular factor is reported through the return code, never thrown) and +solves into `u`; the factorization is cached, so subsequent solves that only +change `b` are an allocation-free triangular solve. Batched (matrix) right-hand +sides are handled natively. With `alias_A = true` the factorization overwrites `cache.A` directly; otherwise a workspace-owned buffer is refilled with `copyto!` on each refactorization and @@ -534,7 +534,7 @@ function init_cacheval( maxiters::Int, abstol, reltol, verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions ) where {T <: LinearAlgebra.BlasFloat} - return (Matrix{T}(undef, 0, 0), Vector{BlasInt}(undef, 0)) + return LinearAlgebra.LU(Matrix{T}(undef, 0, 0), BlasInt[], zero(BlasInt)) end function init_cacheval( @@ -554,39 +554,34 @@ end function SciMLBase.solve!(cache::LinearCache, alg::GESVFactorization; kwargs...) A = convert(AbstractMatrix, cache.A) - F, ipiv = @get_cacheval(cache, :GESVFactorization) + luf = @get_cacheval(cache, :GESVFactorization) if cache.isfresh if cache.alias_A && A isa Matrix # The user permitted overwriting A, so factorize it in place. Atarget = A else - if size(F) != size(A) - F = Matrix{eltype(A)}(undef, size(A)) + Fbuf = getfield(luf, :factors) + if size(Fbuf) != size(A) + Fbuf = similar(A) end - copyto!(F, A) - Atarget = F + copyto!(Fbuf, A) + Atarget = Fbuf end - # Equivalent to dgesv (factorize + solve), but split into getrf!/getrs! - # so a singular factor is reported through the return code rather than a - # thrown exception: getrf! returns `info > 0` for a singular U without - # raising (unlike gesv!, which calls chklapackerror). - factors, ipivnew, info = LAPACK.getrf!(Atarget) - cache.cacheval = (factors, ipivnew) - cache.isfresh = false - if !iszero(info) + # `check = false` returns a singular factorization instead of throwing, + # so a singular denominator is reported through the return code. + luf = lu!(Atarget; check = false) + cache.cacheval = luf + if !issuccess(luf) @SciMLMessage("Solver failed", cache.verbose, :solver_failure) return SciMLBase.build_linear_solution( alg, cache.u, nothing, cache; retcode = ReturnCode.Failure ) end - copyto!(cache.u, cache.b) - LAPACK.getrs!('N', factors, ipivnew, cache.u) - return SciMLBase.build_linear_solution( - alg, cache.u, nothing, cache; retcode = ReturnCode.Success - ) + cache.isfresh = false end + # Solve into `u`; a matrix `b` (batched RHS) is handled natively by ldiv!. copyto!(cache.u, cache.b) - LAPACK.getrs!('N', F, ipiv, cache.u) + ldiv!(luf, cache.u) return SciMLBase.build_linear_solution( alg, cache.u, nothing, cache; retcode = ReturnCode.Success ) From b672d64306443e88b8e966c8e6ee2f3d8ae7e60a Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 9 Jul 2026 19:54:21 -0400 Subject: [PATCH 5/6] test/Core/gesv.jl: use stdlib Random instead of StableRNGs 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 Co-Authored-By: Claude Fable 5 --- test/Core/gesv.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/Core/gesv.jl b/test/Core/gesv.jl index 5f15e8a8b..7b211952b 100644 --- a/test/Core/gesv.jl +++ b/test/Core/gesv.jl @@ -1,7 +1,6 @@ -using LinearSolve, LinearAlgebra, SparseArrays, Test -using StableRNGs +using LinearSolve, LinearAlgebra, SparseArrays, Test, Random -rng = StableRNG(123) +rng = MersenneTwister(123) @testset "GESVFactorization" begin @testset "vector and matrix RHS" begin From b3870d12a52f45d76087aaea933fa13cacf45c33 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 9 Jul 2026 20:04:27 -0400 Subject: [PATCH 6/6] test: add StableRNGs test dependency for gesv.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 Co-Authored-By: Claude Fable 5 --- Project.toml | 4 +++- test/Core/gesv.jl | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index d9d899663..c220ad41c 100644 --- a/Project.toml +++ b/Project.toml @@ -178,6 +178,7 @@ SparseArrays = "1.10" SparseColumnPivotedQR = "2" SparseMatricesCSR = "0.6" Sparspak = "0.3.9" +StableRNGs = "1" SpecializingFactorizations = "0.1" StaticArrays = "1.9" StaticArraysCore = "1.4.3" @@ -213,9 +214,10 @@ SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" SpecializingFactorizations = "fa08b7a1-13d3-4faf-875d-5cbc1520e3f3" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["AlgebraicMultigrid", "BandedMatrices", "BlockDiagonals", "CliqueTrees", "ComponentArrays", "FastAlmostBandedMatrices", "FastLapackInterface", "FiniteDiff", "FixedSizeArrays", "ForwardDiff", "InteractiveUtils", "IterativeSolvers", "KrylovKit", "KrylovPreconditioners", "MultiFloats", "Pkg", "PureUMFPACK", "Random", "RecursiveFactorization", "STRUMPACK_jll", "SafeTestsets", "SciMLTesting", "SparseArrays", "Sparspak", "SpecializingFactorizations", "StaticArrays", "Test", "Zygote"] +test = ["AlgebraicMultigrid", "BandedMatrices", "BlockDiagonals", "CliqueTrees", "ComponentArrays", "FastAlmostBandedMatrices", "FastLapackInterface", "FiniteDiff", "FixedSizeArrays", "ForwardDiff", "InteractiveUtils", "IterativeSolvers", "KrylovKit", "KrylovPreconditioners", "MultiFloats", "Pkg", "PureUMFPACK", "Random", "RecursiveFactorization", "STRUMPACK_jll", "SafeTestsets", "SciMLTesting", "SparseArrays", "Sparspak", "SpecializingFactorizations", "StableRNGs", "StaticArrays", "Test", "Zygote"] diff --git a/test/Core/gesv.jl b/test/Core/gesv.jl index 7b211952b..5f15e8a8b 100644 --- a/test/Core/gesv.jl +++ b/test/Core/gesv.jl @@ -1,6 +1,7 @@ -using LinearSolve, LinearAlgebra, SparseArrays, Test, Random +using LinearSolve, LinearAlgebra, SparseArrays, Test +using StableRNGs -rng = MersenneTwister(123) +rng = StableRNG(123) @testset "GESVFactorization" begin @testset "vector and matrix RHS" begin