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/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 7d56dc7c0..3f2a1052f 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 @@ -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 """ 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")