Use FastLapackInterface.jl tooling

[ChrisRackauckas]

In Base.LinearAlgebra, the Lapack routines are always called twice, once for computing the necessary workspace and once for executing the function. In FastLapackInterface.jl, I introduce a specialized type for each factorization (LU, Cholesky, QR and Schur) that contains all necessary intermediary matrices and a function that execute the algorithm with preallocated space.

https://git.dynare.org/DynareJulia/FastLapackInterface.jl

Seems like that would be useful to use here since it would match our pre-caching interface.

@MichelJuillard

[MichelJuillard]

Tell me if I can be of any help or if I can improve the current implementation in any way

[MichelJuillard]

LU factorization needs two workspaces:

• lu: a n x n matrix containing L and U triangular matrices
• ipvt: a pivot vector with n elements
  If you want to keep the original matrix, you need the lu workspace, otherwise you only need ipvt
  Does LinearSolve.solve() keep the A matrix intact or is the replaced by the factorization?
  Note also, that in the case of LU factorization, Lapack doesn't need a first pass to compute the size of the workspace as it is self evident here.

[rayegun]

By default it does a deepcopy unless you pass alias_A=true to init in common.jl. Then it does lu!, which (I believe) is destructive to that copy held in the LinearCache.

[ChrisRackauckas]

Yeah for safety it defaults to copying, but the keyword argument alias_A = true turns that off. When that argument is true, we are free to modify the original A how we please.

[MichelJuillard]

I. In factorization.jl:41

• A is the object in which we store the factorization
• we need an additional vector to store the pivots
• then we call xgetrf_ to compute the factorization (small improvement with respect to Base.lu!(), the pivot vector is allocated once for all problems with the same size)

II. If the factorization is stored in A by xgetrf_, we need to call xgetrs_ to actually solve the linear problem instead of Base.ldiv!() that expects a Base.LU object