Make AppleAccelerate the default algorithm when available

The benchmarks show that it's just so much better, and there's no cost to doing this, so we should just always use AppleAccelerate on any machine where it exists.

Author: ChrisRackauckas

src/LinearSolve.jl

@@ -90,6 +90,7 @@ EnumX.@enumx DefaultAlgorithmChoice begin
     SVDFactorization
     CholeskyFactorization
     NormalCholeskyFactorization
+    AppleAccelerateLUFactorization
 end
 
 struct DefaultLinearSolver <: SciMLLinearSolveAlgorithm
@@ -98,6 +99,7 @@ end
 
 include("common.jl")
 include("factorization.jl")
+include("appleaccelerate.jl")
 include("simplelu.jl")
 include("simplegmres.jl")
 include("iterative_wrappers.jl")
@@ -106,7 +108,6 @@ include("solve_function.jl")
 include("default.jl")
 include("init.jl")
 include("extension_algs.jl")
-include("appleaccelerate.jl")
 include("deprecated.jl")
 
 @generated function SciMLBase.solve!(cache::LinearCache, alg::AbstractFactorization;

src/appleaccelerate.jl

@@ -38,7 +38,6 @@ function aa_getrf!(A::AbstractMatrix{<:Float64};
     if isempty(ipiv)
         ipiv = similar(A, Cint, min(size(A, 1), size(A, 2)))
     end
-
     ccall(("dgetrf_", libacc), Cvoid,
         (Ref{Cint}, Ref{Cint}, Ptr{Float64},
             Ref{Cint}, Ptr{Cint}, Ptr{Cint}),
@@ -121,11 +120,16 @@ end
 default_alias_A(::AppleAccelerateLUFactorization, ::Any, ::Any) = false
 default_alias_b(::AppleAccelerateLUFactorization, ::Any, ::Any) = false
 
+const PREALLOCATED_APPLE_LU = begin
+    A = rand(0, 0)
+    luinst = ArrayInterface.lu_instance(A)
+    LU(luinst.factors, similar(A, Cint, 0), luinst.info), Ref{Cint}()
+end
+
 function LinearSolve.init_cacheval(alg::AppleAccelerateLUFactorization, A, b, u, Pl, Pr,
     maxiters::Int, abstol, reltol, verbose::Bool,
     assumptions::OperatorAssumptions)
-    luinst = ArrayInterface.lu_instance(convert(AbstractMatrix, A))
-    LU(luinst.factors, similar(A, Cint, 0), luinst.info), Ref{Cint}()
+    PREALLOCATED_APPLE_LU
 end
 
 function SciMLBase.solve!(cache::LinearCache, alg::AppleAccelerateLUFactorization;

src/default.jl

@@ -1,6 +1,6 @@
 needs_concrete_A(alg::DefaultLinearSolver) = true
 mutable struct DefaultLinearSolverInit{T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12,
-    T13, T14, T15, T16}
+    T13, T14, T15, T16, T17}
     LUFactorization::T1
     QRFactorization::T2
     DiagonalFactorization::T3
@@ -17,6 +17,7 @@ mutable struct DefaultLinearSolverInit{T1, T2, T3, T4, T5, T6, T7, T8, T9, T10,
     SVDFactorization::T14
     CholeskyFactorization::T15
     NormalCholeskyFactorization::T16
+    AppleAccelerateLUFactorization::T17
 end
 
 # Legacy fallback
@@ -157,7 +158,9 @@ function defaultalg(A, b, assump::OperatorAssumptions)
                ArrayInterface.can_setindex(b) &&
                (__conditioning(assump) === OperatorCondition.IllConditioned ||
                 __conditioning(assump) === OperatorCondition.WellConditioned)
-                if length(b) <= 10
+                if appleaccelerate_isavailable()
+                    DefaultAlgorithmChoice.AppleAccelerateLUFactorization
+                elseif length(b) <= 10
                     DefaultAlgorithmChoice.GenericLUFactorization
                 elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500)) &&
                        (A === nothing ? eltype(b) <: Union{Float32, Float64} :
@@ -232,6 +235,8 @@ function algchoice_to_alg(alg::Symbol)
         CholeskyFactorization()
     elseif alg === :NormalCholeskyFactorization
         NormalCholeskyFactorization()
+    elseif alg === :AppleAccelerateLUFactorization
+        AppleAccelerateLUFactorization()
     else
         error("Algorithm choice symbol $alg not allowed in the default")
     end