Add NormalCholeskyFactorization and WellConditioned defaults

Continuing on #289 to close #283

Author: ChrisRackauckas

src/default.jl

@@ -155,17 +155,52 @@ function defaultalg(A, b, assump::OperatorAssumptions{true})
     # whether MKL or OpenBLAS is being used
     if (A === nothing && !(b isa GPUArraysCore.AbstractGPUArray)) || A isa Matrix
         if (A === nothing || eltype(A) <: Union{Float32, Float64, ComplexF32, ComplexF64}) &&
-           ArrayInterface.can_setindex(b) && __conditioning(assump) === OperatorCondition.IllConditioned
+           ArrayInterface.can_setindex(b) && (__conditioning(assump) === OperatorCondition.IllConditioned ||
+           __conditioning(assump) === OperatorCondition.WellConditioned)
+
             if length(b) <= 10
-                alg = GenericLUFactorization()
+                pivot = @static if VERSION < v"1.7beta"
+                    if __conditioning(assump) === OperatorCondition.IllConditioned
+                        Val(true)
+                    else
+                        Val(false)
+                    end
+                    else
+                        if __conditioning(assump) === OperatorCondition.IllConditioned
+                            RowMaximum()
+                        else
+                            RowNonZero()
+                        end
+                    end
+                end
+                alg = GenericLUFactorization(pivot)
             elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500)) &&
                    (A === nothing ? eltype(b) <: Union{Float32, Float64} :
                     eltype(A) <: Union{Float32, Float64})
-                alg = RFLUFactorization()
+                pivot =  if __conditioning(assump) === OperatorCondition.IllConditioned
+                    Val(true)
+                else
+                    Val(false)
+                end
+                alg = RFLUFactorization(;pivot = pivot)
                 #elseif A === nothing || A isa Matrix
                 #    alg = FastLUFactorization()
             else
-                alg = LUFactorization()
+                pivot = @static if VERSION < v"1.7beta"
+                    if __conditioning(assump) === OperatorCondition.IllConditioned
+                        Val(true)
+                    else
+                        Val(false)
+                    end
+                    else
+                        if __conditioning(assump) === OperatorCondition.IllConditioned
+                            RowMaximum()
+                        else
+                            RowNonZero()
+                        end
+                    end
+                end
+                alg = LUFactorization(pivot)
             end
         elseif __conditioning(assump) === OperatorCondition.VeryIllConditioned
             alg = QRFactorization()

src/factorization.jl

@@ -472,6 +472,43 @@ function SciMLBase.solve(cache::LinearCache, alg::RFLUFactorization{P, T};
     SciMLBase.build_linear_solution(alg, y, nothing, cache)
 end
 
+## NormalCholeskyFactorization
+
+struct NormalCholeskyFactorization <: AbstractFactorization
+    pivot::P
+    perm::Bool
+end
+
+function NormalCholeskyFactorization(; pivot = nothing, perm = nothing)
+    if pivot === nothing
+        @static if VERSION < v"1.7beta"
+            Val(true)
+        else
+            RowMaximum()
+        end
+    end
+    NormalCholeskyFactorization(pivot, perm)
+end
+
+function init_cacheval(alg::NormalCholeskyFactorization, A, b, u, Pl, Pr,
+                       maxiters::Int, abstol, reltol, verbose::Bool,
+                       assumptions::OperatorAssumptions)
+    ArrayInterface.cholesky_instance(convert(AbstractMatrix, A), alg.pivot)
+end
+
+function SciMLBase.solve(cache::LinearCache, alg::NormalCholeskyFactorization;
+                         kwargs...)
+    A = cache.A
+    A = convert(AbstractMatrix, A)
+    fact, ipiv = cache.cacheval
+    if cache.isfresh
+        fact = cholesky(Symmetric((A)' * A), alg.pivot)
+        cache = set_cacheval(cache, fact)
+    end
+    y = ldiv!(cache.u, cache.cacheval, A' * cache.b)
+    SciMLBase.build_linear_solution(alg, y, nothing, cache)
+end
+
 ## DiagonalFactorization
 
 struct DiagonalFactorization <: AbstractFactorization end