Merge pull request #386 from SciML/appleaccelerate_default

Make AppleAccelerate the default algorithm when available

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,18 @@ end
 default_alias_A(::AppleAccelerateLUFactorization, ::Any, ::Any) = false
 default_alias_b(::AppleAccelerateLUFactorization, ::Any, ::Any) = false
 
+const PREALLOCATED_APPLE_LU = @static if VERSION >= v"1.8" 
+    A = rand(0, 0)
+    luinst = ArrayInterface.lu_instance(A)
+    LU(luinst.factors, similar(A, Cint, 0), luinst.info), Ref{Cint}()
+else
+    nothing
+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
@@ -159,6 +160,8 @@ function defaultalg(A, b, assump::OperatorAssumptions)
                 __conditioning(assump) === OperatorCondition.WellConditioned)
                 if length(b) <= 10
                     DefaultAlgorithmChoice.GenericLUFactorization
+                elseif VERSION >= v"1.8" && appleaccelerate_isavailable()
+                    DefaultAlgorithmChoice.AppleAccelerateLUFactorization
                 elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500)) &&
                        (A === nothing ? eltype(b) <: Union{Float32, Float64} :
                         eltype(A) <: 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