GESVFactorization + dense LU pivot-buffer reuse (v4.2.0)

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"]

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

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

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,

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
@@ -467,6 +507,86 @@ function init_cacheval(
     return nothing
 end
 
+## GESVFactorization
+
+"""
+`GESVFactorization()`
+
+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
+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 LinearAlgebra.LU(Matrix{T}(undef, 0, 0), BlasInt[], zero(BlasInt))
+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)
+    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
+            Fbuf = getfield(luf, :factors)
+            if size(Fbuf) != size(A)
+                Fbuf = similar(A)
+            end
+            copyto!(Fbuf, A)
+            Atarget = Fbuf
+        end
+        # `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
+        cache.isfresh = false
+    end
+    # Solve into `u`; a matrix `b` (batched RHS) is handled natively by ldiv!.
+    copyto!(cache.u, cache.b)
+    ldiv!(luf, cache.u)
+    return SciMLBase.build_linear_solution(
+        alg, cache.u, nothing, cache; retcode = ReturnCode.Success
+    )
+end
+
 ## QRFactorization
 
 """

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

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")