Handling of getrf! and getrs!

[ChrisRackauckas]

Test case:

using LinearSolve, LinearAlgebra
# using MKL_jll

n = 100
A = rand(n, n)
b1 = rand(n);
b2 = rand(n);

function f(A, b1, b2; alg = LUFactorization())
    prob = LinearProblem(A, b1)

    linsolve = init(prob, alg)
    sol1 = solve!(linsolve)

    s1 = copy(sol1.u)

    linsolve.b = b2
    sol2 = solve!(linsolve)

    s2 = copy(sol2.u)
    norm(s1 + s2)
end

f(A, b1, b2) # Uses BLAS
f(A, b1, b2; alg=RFLUFactorization()) # Uses loops
f(A, b1, b2; alg=MKLLUFactorization()) # Requires `using MKL_jll`

using Enzyme

dA = zero(A)
db1 = zero(b1)
db2 = zero(b2)
Enzyme.autodiff(Reverse, f, Duplicated(A,dA), 
                Duplicated(b1, db1), Duplicated(b2, db2))

The LUFactorization() version uses OpenBLAS, the RFLUFactorization() version uses a pure Julia implementation https://github.com/JuliaLinearAlgebra/RecursiveFactorization.jl/blob/master/src/lu.jl (no BLAS), and MKLLUFactorization swaps to an alternative BLAS with wrappers fully written here: https://github.com/SciML/LinearSolve.jl/blob/main/ext/LinearSolveMKLExt.jl. I assume RFLUFactorization will just work, it's just scalar mutating code, but LUFactorization makes sure base Julia is working and the MKL one is a good test that it's at the right level of libblastrampoline

[wsmoses]

Closing, as the LinearSolve package used above should be fine now