From 592dc07fc3ad70fdddf064efa4c643e3b6913926 Mon Sep 17 00:00:00 2001 From: jClugstor Date: Mon, 28 Jul 2025 14:27:21 -0400 Subject: [PATCH 1/6] add test for forwarddiff static arrays --- test/forwarddiff_overloads.jl | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/test/forwarddiff_overloads.jl b/test/forwarddiff_overloads.jl index 3d4a035f3..789e403d5 100644 --- a/test/forwarddiff_overloads.jl +++ b/test/forwarddiff_overloads.jl @@ -193,3 +193,31 @@ overload_x_p = solve(prob, UMFPACKFactorization()) backslash_x_p = A \ b @test ≈(overload_x_p, backslash_x_p, rtol = 1e-9) + + +# Test StaticArrays +# They don't go through the overloads +# But we should test that the overloads don't mess anything up +function static_h(p) + A = @SMatrix [p[1] p[2]+1 p[2]^3; + 3*p[1] p[1]+5 p[2] * p[1]-4; + p[2]^2 9*p[1] p[2]] + + b = SA[p[1] + 1, p[2] * 2, p[1]^2] + + (A, b) +end + +function static_linsolve(p) + A, b = static_h(p) + prob = LinearProblem(A, b) + solve(prob) +end + +function static_backslash(p) + A, b = static_h(p) + A \ b +end + +@test (ForwardDiff.jacobian(static_linsolve, [5.0, 5.0]) ≈ + ForwardDiff.jacobian(static_backslash, [5.0, 5.0])) \ No newline at end of file From 763e60b2f0622514244f901334a68e21c3a4a9a3 Mon Sep 17 00:00:00 2001 From: jClugstor Date: Mon, 28 Jul 2025 14:42:57 -0400 Subject: [PATCH 2/6] add test for cache type --- test/forwarddiff_overloads.jl | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/test/forwarddiff_overloads.jl b/test/forwarddiff_overloads.jl index 789e403d5..d5373c49b 100644 --- a/test/forwarddiff_overloads.jl +++ b/test/forwarddiff_overloads.jl @@ -208,6 +208,16 @@ function static_h(p) (A, b) end +# Check to make sure it's not a DualLinearCache +function static_linsolve_cache_check(p) + A,b = static_h(p) + prob = LinearProblem(A, b) + cache = init(prob) + cache isa LinearCache +end + +@test static_linsolve_cache_check([5.0, 5.0]) + function static_linsolve(p) A, b = static_h(p) prob = LinearProblem(A, b) @@ -220,4 +230,12 @@ function static_backslash(p) end @test (ForwardDiff.jacobian(static_linsolve, [5.0, 5.0]) ≈ - ForwardDiff.jacobian(static_backslash, [5.0, 5.0])) \ No newline at end of file + ForwardDiff.jacobian(static_backslash, [5.0, 5.0])) + +#Test to make sure that the cache is not a DualLinearCache +A, b = static_h([ForwardDiff.Dual(5.0, 1.0, 0.0), ForwardDiff.Dual(5.0, 0.0, 1.0)]) + +static_dual_prob = LinearProblem(A, b) +static_dual_cache = init(static_dual_prob) + +@test static_dual_cache isa LinearSolve.LinearCache \ No newline at end of file From c28a6da887449638a8005bd94ef8422a79511adf Mon Sep 17 00:00:00 2001 From: jClugstor Date: Mon, 28 Jul 2025 15:06:08 -0400 Subject: [PATCH 3/6] import StaticArrays for test --- test/forwarddiff_overloads.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/forwarddiff_overloads.jl b/test/forwarddiff_overloads.jl index d5373c49b..3e8d66297 100644 --- a/test/forwarddiff_overloads.jl +++ b/test/forwarddiff_overloads.jl @@ -2,6 +2,7 @@ using LinearSolve using ForwardDiff using Test using SparseArrays +using StaticArrays function h(p) (A = [p[1] p[2]+1 p[2]^3; From 949e8e8fd404ee3256d0d88a0612b25e00bd3256 Mon Sep 17 00:00:00 2001 From: jClugstor Date: Mon, 28 Jul 2025 15:25:20 -0400 Subject: [PATCH 4/6] remove unnecessary part of test --- test/forwarddiff_overloads.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/test/forwarddiff_overloads.jl b/test/forwarddiff_overloads.jl index 3e8d66297..97fef8d6f 100644 --- a/test/forwarddiff_overloads.jl +++ b/test/forwarddiff_overloads.jl @@ -209,16 +209,6 @@ function static_h(p) (A, b) end -# Check to make sure it's not a DualLinearCache -function static_linsolve_cache_check(p) - A,b = static_h(p) - prob = LinearProblem(A, b) - cache = init(prob) - cache isa LinearCache -end - -@test static_linsolve_cache_check([5.0, 5.0]) - function static_linsolve(p) A, b = static_h(p) prob = LinearProblem(A, b) From 8b4e443ff5c184f07c7f693a8beedd8420a32bc1 Mon Sep 17 00:00:00 2001 From: jClugstor Date: Tue, 29 Jul 2025 09:53:12 -0400 Subject: [PATCH 5/6] add init for StaticLinearProblem --- src/common.jl | 114 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) diff --git a/src/common.jl b/src/common.jl index f447df8e7..0af4821b9 100644 --- a/src/common.jl +++ b/src/common.jl @@ -337,3 +337,117 @@ function SciMLBase.solve(prob::StaticLinearProblem, return SciMLBase.build_linear_solution( alg, u, nothing, prob; retcode = ReturnCode.Success) end + +function SciMLBase.init(prob::StaticLinearProblem, alg::SciMLLinearSolveAlgorithm, + args...; + alias = LinearAliasSpecifier(), + abstol = default_tol(real(eltype(prob.b))), + reltol = default_tol(real(eltype(prob.b))), + maxiters::Int = length(prob.b), + verbose::Bool = false, + Pl = nothing, + Pr = nothing, + assumptions = OperatorAssumptions(issquare(prob.A)), + sensealg = LinearSolveAdjoint(), + kwargs...) + (; A, b, u0, p) = prob + + if haskey(kwargs, :alias_A) || haskey(kwargs, :alias_b) + aliases = LinearAliasSpecifier() + + if haskey(kwargs, :alias_A) + message = "`alias_A` keyword argument is deprecated, to set `alias_A`, + please use an LinearAliasSpecifier, e.g. `solve(prob, alias = LinearAliasSpecifier(alias_A = true))" + Base.depwarn(message, :init) + Base.depwarn(message, :solve) + aliases = LinearAliasSpecifier(alias_A = values(kwargs).alias_A) + end + + if haskey(kwargs, :alias_b) + message = "`alias_b` keyword argument is deprecated, to set `alias_b`, + please use an LinearAliasSpecifier, e.g. `solve(prob, alias = LinearAliasSpecifier(alias_b = true))" + Base.depwarn(message, :init) + Base.depwarn(message, :solve) + aliases = LinearAliasSpecifier( + alias_A = aliases.alias_A, alias_b = values(kwargs).alias_b) + end + else + if alias isa Bool + aliases = LinearAliasSpecifier(alias = alias) + else + aliases = alias + end + end + + if isnothing(aliases.alias_A) + alias_A = default_alias_A(alg, prob.A, prob.b) + else + alias_A = aliases.alias_A + end + + if isnothing(aliases.alias_b) + alias_b = default_alias_b(alg, prob.A, prob.b) + else + alias_b = aliases.alias_b + end + + A = if alias_A || A isa SMatrix + A + elseif A isa Array + copy(A) + elseif issparsematrixcsc(A) + make_SparseMatrixCSC(A) + else + deepcopy(A) + end + + b = if issparsematrix(b) && !(A isa Diagonal) + Array(b) # the solution to a linear solve will always be dense! + elseif alias_b || b isa SVector + b + elseif b isa Array + copy(b) + elseif issparsematrixcsc(b) + # Extension must be loaded if issparsematrixcsc returns true + make_SparseMatrixCSC(b) + else + deepcopy(b) + end + + u0_ = u0 !== nothing ? u0 : __init_u0_from_Ab(A, b) + + # Guard against type mismatch for user-specified reltol/abstol + reltol = real(eltype(prob.b))(reltol) + abstol = real(eltype(prob.b))(abstol) + + precs = if hasproperty(alg, :precs) + isnothing(alg.precs) ? DEFAULT_PRECS : alg.precs + else + DEFAULT_PRECS + end + _Pl, _Pr = precs(A, p) + if isnothing(Pl) + Pl = _Pl + else + # TODO: deprecate once all docs are updated to the new form + #@warn "passing Preconditioners at `init`/`solve` time is deprecated. Instead add a `precs` function to your algorithm." + end + if isnothing(Pr) + Pr = _Pr + else + # TODO: deprecate once all docs are updated to the new form + #@warn "passing Preconditioners at `init`/`solve` time is deprecated. Instead add a `precs` function to your algorithm." + end + cacheval = init_cacheval(alg, A, b, u0_, Pl, Pr, maxiters, abstol, reltol, verbose, + assumptions) + isfresh = true + precsisfresh = false + Tc = typeof(cacheval) + + cache = LinearCache{typeof(A), typeof(b), typeof(u0_), typeof(p), typeof(alg), Tc, + typeof(Pl), typeof(Pr), typeof(reltol), typeof(assumptions.issq), + typeof(sensealg)}( + A, b, u0_, p, alg, cacheval, isfresh, precsisfresh, Pl, Pr, abstol, reltol, + maxiters, verbose, assumptions, sensealg) + return cache +end \ No newline at end of file From ff84c9a28112d8606b9284b4e82a09bf6c9d8519 Mon Sep 17 00:00:00 2001 From: jClugstor Date: Tue, 29 Jul 2025 09:54:39 -0400 Subject: [PATCH 6/6] add explanatory comment --- src/common.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/common.jl b/src/common.jl index 0af4821b9..de0e4d67d 100644 --- a/src/common.jl +++ b/src/common.jl @@ -338,6 +338,7 @@ function SciMLBase.solve(prob::StaticLinearProblem, alg, u, nothing, prob; retcode = ReturnCode.Success) end +# Here to make sure that StaticLinearProblems with Dual elements don't create a Dual linear cache function SciMLBase.init(prob::StaticLinearProblem, alg::SciMLLinearSolveAlgorithm, args...; alias = LinearAliasSpecifier(),