Setup MKL direct factorizations

[ChrisRackauckas]

MWE:

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:

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

[ChrisRackauckas]

@staticfloat or @chriselrod maybe you might have an idea of what I'm doing wrong here to cause a segfault?

[codecov[bot]]

Codecov Report

Merging #349 (aaf64d3) into main (98a2292) will decrease coverage by 10.60%. The diff coverage is 88.46%.

@@             Coverage Diff             @@
##             main     #349       +/-   ##
===========================================
- Coverage   75.29%   64.70%   -10.60%     
===========================================
  Files          17       18        +1     
  Lines        1267     1292       +25     
===========================================
- Hits          954      836      -118     
- Misses        313      456      +143     

Files Changed	Coverage Δ
src/LinearSolve.jl	84.09% <50.00%> (-1.63%)	:arrow_down:
ext/LinearSolveMKLExt.jl	91.30% <91.30%> (ø)
src/extension_algs.jl	12.50% <100.00%> (-57.07%)	:arrow_down:

... and 4 files with indirect coverage changes

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

[staticfloat]

You've got MKL in 32-bit mode, but passing it Int64's, so the pivot indices that come back are gibberish, which causes the future LU factorization to have bogus indices, and the \ to segfault:

julia> A, ipiv, info = getrf!(A);

julia> ipiv
4-element Vector{Int64}:
  8589934595
 17179869187
           0
           0

There are a bunch of ways to fix this, but the best is to import @blasfunc (remember that guy we dropped when copy-pasting) from LinearAlgebra.BLAS and use that via @blasfunc(dgetrf_). That's good because you're already polymorphic by using similar(..., BlasInt, ...) when initializing ipiv, so this just makes it match up to the correct version by default.