Support Apple Accelerate and improve MKL integration

ext/LinearSolveMKLExt.jl

@@ -16,6 +16,9 @@
     chkstride1(A)
     m, n = size(A)
     lda  = max(1,stride(A, 2))
+    if isempty(ipiv)
+        ipiv = similar(A, BlasInt, min(size(A,1),size(A,2)))
+    end
     ccall((@blasfunc(dgetrf_), MKL_jll.libmkl_rt), Cvoid,
             (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
             Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
@@ -39,7 +42,7 @@
     A = convert(AbstractMatrix, A)
     if cache.isfresh
         cacheval = @get_cacheval(cache, :MKLLUFactorization)
-        fact = LU(getrf!(A)...)
+        fact = LU(getrf!(A; ipiv = cacheval.ipiv)...)
         cache.cacheval = fact
         cache.isfresh = false
     end

src/LinearSolve.jl

@@ -23,6 +23,10 @@ using EnumX
 using Requires
 import InteractiveUtils
 
+using LinearAlgebra: BlasInt, LU
+using LinearAlgebra.LAPACK: require_one_based_indexing, chkfinite, chkstride1, 
+                            @blasfunc, chkargsok
+
 import GPUArraysCore
 import Preferences
 
@@ -87,6 +91,7 @@ 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;
@@ -185,6 +190,7 @@ export CudaOffloadFactorization
 export MKLPardisoFactorize, MKLPardisoIterate
 export PardisoJL
 export MKLLUFactorization
+export AppleAccelerateLUFactorization
 
 export OperatorAssumptions, OperatorCondition
 

src/appleaccelerate.jl

@@ -0,0 +1,65 @@
+# For now, only use BLAS from Accelerate (that is to say, vecLib)
+global const libacc = "/System/Library/Frameworks/Accelerate.framework/Accelerate"
+
+"""
+```julia
+AppleAccelerateLUFactorization()
+```
+
+A wrapper over Apple's Accelerate Library. Direct calls to Acceelrate in a way that pre-allocates workspace
+to avoid allocations and does not require libblastrampoline.
+"""
+struct AppleAccelerateLUFactorization <: AbstractFactorization end
+
+function appleaccelerate_isavailable()
+    libacc_hdl = dlopen_e(libacc)
+    if libacc_hdl == C_NULL
+        return false
+    end
+
+    if dlsym_e(libacc_hdl, "dgemm\$NEWLAPACK\$ILP64") == C_NULL
+        return false
+    end
+    return true
+end
+
+function aa_getrf!(A::AbstractMatrix{<:Float64}; ipiv = similar(A, BlasInt, min(size(A,1),size(A,2))), info = Ref{BlasInt}(), check = false)
+    require_one_based_indexing(A)
+    check && chkfinite(A)
+    chkstride1(A)
+    m, n = size(A)
+    lda  = max(1,stride(A, 2))
+    if isempty(ipiv)
+        ipiv = similar(A, BlasInt, min(size(A,1),size(A,2)))
+    end
+
+    ccall(("dgetrf\$NEWLAPACK\$ILP64", libacc), Cvoid,
+            (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
+            Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
+            m, n, A, lda, ipiv, info)
+    chkargsok(info[])
+    A, ipiv, info[] #Error code is stored in LU factorization type
+end
+
+default_alias_A(::AppleAccelerateLUFactorization, ::Any, ::Any) = false
+default_alias_b(::AppleAccelerateLUFactorization, ::Any, ::Any) = false
+
+function LinearSolve.init_cacheval(alg::AppleAccelerateLUFactorization, A, b, u, Pl, Pr,
+    maxiters::Int, abstol, reltol, verbose::Bool,
+    assumptions::OperatorAssumptions)
+    ArrayInterface.lu_instance(convert(AbstractMatrix, A))
+end
+
+function SciMLBase.solve!(cache::LinearCache, alg::AppleAccelerateLUFactorization;
+    kwargs...)
+    A = cache.A
+    A = convert(AbstractMatrix, A)
+    if cache.isfresh
+        cacheval = @get_cacheval(cache, :AppleAccelerateLUFactorization)
+        fact = LU(aa_getrf!(A; ipiv = cacheval.ipiv)...)
+        cache.cacheval = fact
+        cache.isfresh = false
+    end
+    y = ldiv!(cache.u, @get_cacheval(cache, :AppleAccelerateLUFactorization), cache.b)
+    SciMLBase.build_linear_solution(alg, y, nothing, cache)
+end

test/basictests.jl

@@ -213,6 +213,9 @@ end
                 test_interface(alg, prob1, prob2)
             end
         end
+        if LinearSolve.appleaccelerate_isavailable()
+            test_interface(AppleAccelerateLUFactorization(), prob1, prob2)
+        end
     end
 
     @testset "Generic Factorizations" begin