From 818e59d3db8b54ae6d2954bcffad1225174f8412 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 14 Jul 2025 09:06:54 -0400 Subject: [PATCH 1/2] Fix successful return codes for factorization methods Fixes https://github.com/SciML/LinearSolve.jl/issues/630 fixes https://github.com/SciML/LinearSolve.jl/issues/532 --- ext/LinearSolveRecursiveFactorizationExt.jl | 2 +- src/appleaccelerate.jl | 7 ++++++- src/factorization.jl | 18 +++++++++++++++--- src/mkl.jl | 7 ++++++- test/retcodes.jl | 15 +++++++++++---- 5 files changed, 39 insertions(+), 10 deletions(-) diff --git a/ext/LinearSolveRecursiveFactorizationExt.jl b/ext/LinearSolveRecursiveFactorizationExt.jl index c4794c44a..3e8455eb2 100644 --- a/ext/LinearSolveRecursiveFactorizationExt.jl +++ b/ext/LinearSolveRecursiveFactorizationExt.jl @@ -25,7 +25,7 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::RFLUFactorization cache.isfresh = false end y = ldiv!(cache.u, LinearSolve.@get_cacheval(cache, :RFLUFactorization)[1], cache.b) - SciMLBase.build_linear_solution(alg, y, nothing, cache) + SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = ReturnCode.Success) end end diff --git a/src/appleaccelerate.jl b/src/appleaccelerate.jl index 917ad9c4a..92e2f6403 100644 --- a/src/appleaccelerate.jl +++ b/src/appleaccelerate.jl @@ -243,6 +243,11 @@ function SciMLBase.solve!(cache::LinearCache, alg::AppleAccelerateLUFactorizatio res = aa_getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2]) fact = LU(res[1:3]...), res[4] cache.cacheval = fact + + if !LinearAlgebra.issuccess(fact) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Failure) + end cache.isfresh = false end @@ -258,5 +263,5 @@ function SciMLBase.solve!(cache::LinearCache, alg::AppleAccelerateLUFactorizatio aa_getrs!('N', A.factors, A.ipiv, cache.u; info) end - SciMLBase.build_linear_solution(alg, cache.u, nothing, cache) + SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success) end diff --git a/src/factorization.jl b/src/factorization.jl index 8d03a5d2c..d1a496be3 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -150,7 +150,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::LUFactorization; kwargs...) F = @get_cacheval(cache, :LUFactorization) y = _ldiv!(cache.u, F, cache.b) - SciMLBase.build_linear_solution(alg, y, nothing, cache) + SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = ReturnCode.Success) end function do_factorization(alg::LUFactorization, A, b, u) @@ -203,7 +203,7 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::GenericLUFactoriz cache.isfresh = false end y = ldiv!(cache.u, LinearSolve.@get_cacheval(cache, :GenericLUFactorization)[1], cache.b) - SciMLBase.build_linear_solution(alg, y, nothing, cache) + SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = ReturnCode.Success) end function init_cacheval( @@ -953,6 +953,12 @@ A fast factorization which uses a Cholesky factorization on A * A'. Can be much faster than LU factorization, but is not as numerically stable and thus should only be applied to well-conditioned matrices. +!!! warn + `NormalCholeskyFactorization` should only be applied to well-conditioned matrices. As a + method it is not able to easily identify possible numerical issues. As a check it is + recommended that the user checks `A*u-b` is approximately zero, as this may be untrue + even if `sol.retcode === ReturnCode.Success` due to numerical stability issues. + ## Positional Arguments - pivot: Defaults to RowMaximum(), but can be NoPivot() @@ -1010,6 +1016,12 @@ function SciMLBase.solve!(cache::LinearCache, alg::NormalCholeskyFactorization; fact = cholesky(Symmetric((A)' * A), alg.pivot; check = false) end cache.cacheval = fact + + if hasmethod(LinearAlgebra.issuccess, Tuple{typeof(fact)}) && !LinearAlgebra.issuccess(fact) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Failure) + end + cache.isfresh = false end if issparsematrixcsc(A) @@ -1021,7 +1033,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::NormalCholeskyFactorization; else y = ldiv!(cache.u, @get_cacheval(cache, :NormalCholeskyFactorization), A' * cache.b) end - SciMLBase.build_linear_solution(alg, y, nothing, cache) + SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = ReturnCode.Success) end ## NormalBunchKaufmanFactorization diff --git a/src/mkl.jl b/src/mkl.jl index a00882e1d..894306115 100644 --- a/src/mkl.jl +++ b/src/mkl.jl @@ -219,11 +219,16 @@ function SciMLBase.solve!(cache::LinearCache, alg::MKLLUFactorization; res = getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2]) fact = LU(res[1:3]...), res[4] cache.cacheval = fact + + if !LinearAlgebra.issuccess(fact) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Failure) + end cache.isfresh = false end y = ldiv!(cache.u, @get_cacheval(cache, :MKLLUFactorization)[1], cache.b) - SciMLBase.build_linear_solution(alg, y, nothing, cache) + SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = ReturnCode.Success) #= A, info = @get_cacheval(cache, :MKLLUFactorization) diff --git a/test/retcodes.jl b/test/retcodes.jl index af59de6eb..53fd8db4b 100644 --- a/test/retcodes.jl +++ b/test/retcodes.jl @@ -1,4 +1,4 @@ -using LinearSolve, LinearAlgebra, RecursiveFactorization +using LinearSolve, LinearAlgebra, RecursiveFactorization, Test alglist = ( LUFactorization, @@ -14,12 +14,13 @@ alglist = ( @testset "Success" begin for alg in alglist + @show alg A = [2.0 1.0; -1.0 1.0] b = [-1.0, 1.0] prob = LinearProblem(A, b) linsolve = init(prob, alg()) sol = solve!(linsolve) - @test SciMLBase.successful_retcode(sol.retcode) || sol.retcode == ReturnCode.Default # The latter seems off... + @test SciMLBase.successful_retcode(sol.retcode) end end @@ -39,7 +40,13 @@ lualgs = ( prob = LinearProblem(A, b) linsolve = init(prob, alg) sol = solve!(linsolve) - @test !SciMLBase.successful_retcode(sol.retcode) + if alg isa NormalCholeskyFactorization + # This is a known and documented incorrectness in NormalCholeskyFactorization + # due to numerical instability in its method that is fundamental. + @test SciMLBase.successful_retcode(sol.retcode) + else + @test !SciMLBase.successful_retcode(sol.retcode) + end end end @@ -62,4 +69,4 @@ rankdeficientalgs = ( sol = solve!(linsolve) @test SciMLBase.successful_retcode(sol.retcode) end -end \ No newline at end of file +end From f235234cc8445acf7eb6e84f51f02020222be8f9 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 14 Jul 2025 09:19:31 -0400 Subject: [PATCH 2/2] Update test/retcodes.jl --- test/retcodes.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/retcodes.jl b/test/retcodes.jl index 53fd8db4b..42ffd7d7f 100644 --- a/test/retcodes.jl +++ b/test/retcodes.jl @@ -14,7 +14,6 @@ alglist = ( @testset "Success" begin for alg in alglist - @show alg A = [2.0 1.0; -1.0 1.0] b = [-1.0, 1.0] prob = LinearProblem(A, b)