WIP: Setup MKL direct factorizations

MWE:

```julia
using LinearSolve, MKL_jll
A = rand(4, 4); b = rand(4); u0 = zeros(4);
lp = LinearProblem(A, b);
truesol = solve(lp, LUFactorization())
mklsol = solve(lp, MKLLUFactorization())
@test truesol ≈ mklsol
```

The segfault can be reproduced just with the triangular solver. MWE without LinearSolve:

```julia
using MKL_jll
using LinearAlgebra: BlasInt, LU
using LinearAlgebra.LAPACK: require_one_based_indexing, chkfinite, chkstride1, chkargsok
const usemkl = MKL_jll.is_available()

function 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))
    ccall((:dgetrf_, MKL_jll.libmkl_rt), 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

A = rand(4,4); b = rand(4)
getrf!(A)
LU(getrf!(A)...) \ b
```

Author: ChrisRackauckas

Project.toml

@@ -29,6 +29,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
+MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
@@ -37,6 +38,7 @@ Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 LinearSolveCUDAExt = "CUDA"
 LinearSolveHYPREExt = "HYPRE"
 LinearSolveIterativeSolversExt = "IterativeSolvers"
+LinearSolveMKLExt = "MKL_jll"
 LinearSolveKrylovKitExt = "KrylovKit"
 LinearSolvePardisoExt = "Pardiso"
 

ext/LinearSolveMKLExt.jl

@@ -0,0 +1,54 @@
+module LinearSolveMKLExt
+
+using MKL_jll
+using LinearAlgebra: BlasInt, LU
+using LinearAlgebra.LAPACK: require_one_based_indexing, chkfinite, chkstride1, chkargsok
+using LinearAlgebra
+const usemkl = MKL_jll.is_available()
+
+using LinearSolve
+using LinearSolve: ArrayInterface, MKLLUFactorization, @get_cacheval, LinearCache, SciMLBase
+
+function 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_, MKL_jll.libmkl_rt), 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(::MKLLUFactorization, ::Any, ::Any) = false
+default_alias_b(::MKLLUFactorization, ::Any, ::Any) = false
+
+function LinearSolve.init_cacheval(alg::MKLLUFactorization, 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::MKLLUFactorization;
+    kwargs...)
+    A = cache.A
+    A = convert(AbstractMatrix, A)
+    if cache.isfresh
+        cacheval = @get_cacheval(cache, :MKLLUFactorization)
+        fact = LU(getrf!(A)...)
+        cache.cacheval = fact
+        cache.isfresh = false
+    end
+    y = ldiv!(cache.u, @get_cacheval(cache, :MKLLUFactorization), cache.b)
+    SciMLBase.build_linear_solution(alg, y, nothing, cache)
+end
+
+end

src/LinearSolve.jl

@@ -181,6 +181,7 @@ export HYPREAlgorithm
 export CudaOffloadFactorization
 export MKLPardisoFactorize, MKLPardisoIterate
 export PardisoJL
+export MKLLUFactorization
 
 export OperatorAssumptions, OperatorCondition
 

src/extension_algs.jl

@@ -337,3 +337,13 @@ A wrapper over the IterativeSolvers.jl MINRES.
 
 """
 function IterativeSolversJL_MINRES end
+
+"""
+```julia
+MKLLUFactorization()
+```
+
+A wrapper over Intel's Math Kernel Library (MKL). Direct calls to MKL in a way that pre-allocates workspace
+to avoid allocations and does not require libblastrampoline.
+"""
+struct MKLLUFactorization <: AbstractFactorization end