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..42ffd7d7f 100644
--- a/test/retcodes.jl
+++ b/test/retcodes.jl
@@ -1,4 +1,4 @@
-using LinearSolve, LinearAlgebra, RecursiveFactorization
+using LinearSolve, LinearAlgebra, RecursiveFactorization, Test
 
 alglist = (
     LUFactorization,
@@ -19,7 +19,7 @@ alglist = (
         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 +39,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 +68,4 @@ rankdeficientalgs = (
         sol = solve!(linsolve)
         @test SciMLBase.successful_retcode(sol.retcode)
     end
-end
\ No newline at end of file
+end