early failure on bad factorization

SciML · Apr 22, 2024 · dc7a1bb · dc7a1bb
1 parent 8d6494c
commit dc7a1bb
Show file tree

Hide file tree

Showing 3 changed files with 31 additions and 14 deletions.
diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
@@ -134,26 +134,33 @@ include("extension_algs.jl")
 include("adjoint.jl")
 include("deprecated.jl")
 
+@inline function _issucess(F::LinearAlgebra.QRCompactWY)
+    (m, n) = size(F)
+    U = view(F.factors, 1:min(m, n), 1:n)
+    return any(iszero, Iterators.reverse(@view U[diagind(U)]))
+end
+@inline _issucess(F) = hasmethod(LinearAlgebra.issuccess, (typeof(F),)) ?
+                       LinearAlgebra.issuccess(F) : true
+
 @generated function SciMLBase.solve!(cache::LinearCache, alg::AbstractFactorization;
         kwargs...)
     quote
         if cache.isfresh
             fact = do_factorization(alg, cache.A, cache.b, cache.u)
             cache.cacheval = fact
+
+            # If factorization was not successful, return failure. Don't reset `isfresh`
+            if _issucess(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, $(Meta.quot(defaultalg_symbol(alg)))),
             cache.b)
-
-        #=
-        retcode = if LinearAlgebra.issuccess(fact)
-            SciMLBase.ReturnCode.Success
-        else
-            SciMLBase.ReturnCode.Failure
-        end
-        SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = retcode)
-        =#
-        SciMLBase.build_linear_solution(alg, y, nothing, cache)
+        return SciMLBase.build_linear_solution(alg, y, nothing, cache)
     end
 end
 

diff --git a/src/default.jl b/src/default.jl
@@ -216,9 +216,19 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool})
     elseif assump.condition === OperatorCondition.WellConditioned
         DefaultAlgorithmChoice.NormalCholeskyFactorization
     elseif assump.condition === OperatorCondition.IllConditioned
-        DefaultAlgorithmChoice.QRFactorizationPivoted
+        if is_underdetermined(A)
+            # Underdetermined
+            DefaultAlgorithmChoice.QRFactorizationPivoted
+        else
+            DefaultAlgorithmChoice.QRFactorization
+        end
     elseif assump.condition === OperatorCondition.VeryIllConditioned
-        DefaultAlgorithmChoice.QRFactorizationPivoted
+        if is_underdetermined(A)
+            # Underdetermined
+            DefaultAlgorithmChoice.QRFactorizationPivoted
+        else
+            DefaultAlgorithmChoice.QRFactorization
+        end
     elseif assump.condition === OperatorCondition.SuperIllConditioned
         DefaultAlgorithmChoice.SVDFactorization
     else

diff --git a/test/default_algs.jl b/test/default_algs.jl
@@ -121,7 +121,7 @@ b = [1.0, 0.0, 0.0]
 prob = LinearProblem(A, b)
 sol = solve(prob)
 
-@test sol.u ≈ [0.4, 0.2]
+@test !SciMLBase.successful_retcode(sol.retcode)
 
 ## Show that we cannot select a default alg once by checking the rank, since it might change
 ## later in the cache
@@ -143,4 +143,4 @@ cache.A = [2.0 1.0
 
 sol = solve!(cache)
 
-@test sol.u ≈ [0.4, 0.2]
+@test !SciMLBase.successful_retcode(sol.retcode)