# \[ANN\] FastLapackInterface.jl v1.0.0: Non-allocating LAPACK factorizations

**URL:** <https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354>\
**Category:** Package Announcements\
**Tags:** linearalgebra, lapack, dynare\
**Created:** [June 26, 2022, 11:57am UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354 "2022-06-26T11:57:40Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![louisponet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/louisponet/32/2070_2.png) [@louisponet](https://discourse.julialang.org/u/louisponet)\
**Post date:** [June 26, 2022, 11:57am UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/1 "2022-06-26T11:57:41Z")

</div>

The DynareJulia team and I are excited to announce the first major release of [FastLapackInterface.jl](https://github.com/DynareJulia/FastLapackInterface.jl). The package facilitates the pre-allocation of workspaces to be used with some of the most important LAPACK functions with respect to factorizations. This eliminates all major sources of allocations inside the LAPACK functions. The package for now targets [QR](https://dynarejulia.github.io/FastLapackInterface.jl/dev/LAPACK/#QR), [Schur](https://dynarejulia.github.io/FastLapackInterface.jl/dev/LAPACK/#Schur) and [LU](https://dynarejulia.github.io/FastLapackInterface.jl/dev/LAPACK/#LU) factorizations, but functionality may be added in the future.

We have designed the package to be as consistent with the LAPACK functionality in Base julia, making it highly transparent. The [QRWs](https://dynarejulia.github.io/FastLapackInterface.jl/dev/workspaces/#WorkSpaces) workspace for example can be used with [LinearAlgebra.QR](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.QR) as follows:

```julia
julia> A = [1.2 2.3
            6.2 3.3]
2×2 Matrix{Float64}:
 1.2 2.3
 6.2 3.3

julia> ws = QRWs(A)
QRWs{Float64}
work: 64-element Vector{Float64}
τ: 2-element Vector{Float64}

julia> t = QR(LAPACK.geqrf!(A, ws)...)
QR{Float64, Matrix{Float64}, Vector{Float64}}
Q factor:
2×2 QRPackedQ{Float64, Matrix{Float64}, Vector{Float64}}:
 -0.190022 -0.98178
 -0.98178 0.190022
R factor:
2×2 Matrix{Float64}:
 -6.31506 -3.67692
  0.0 -1.63102

julia> Matrix(t)
2×2 Matrix{Float64}:
 1.2 2.3
 6.2 3.3

```

The package is fully tested and documented, and a suite of benchmarks to compare with Base julia LAPACK functions is included.

All the best,  
Louis

---

<div class="post-metadata">

**Author:** ![Jae-Mo\_Lihm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jae-mo_lihm/32/20145_2.png) [@Jae-Mo\_Lihm](https://discourse.julialang.org/u/Jae-Mo_Lihm)\
**Post date:** [June 29, 2022, 1:35am UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/2 "2022-06-29T01:35:58Z")

</div>

Thanks a lot for this package! It would be very useful in my applications where I need LAPACK operations for many small matrices.

I have two questions:

1. Do you have a plan to integrate more functions? Specifically, I use Hermitian diagonalization (syev / heev) a lot.

2. If I understood the code correctly, currently the non-allocating function throws if the workspace is too small. Would it be possible to have an option to just resize the workspace and run?

---

<div class="post-metadata">

**Author:** ![louisponet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/louisponet/32/2070_2.png) [@louisponet](https://discourse.julialang.org/u/louisponet)\
**Post date:** [June 29, 2022, 8:42am UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/3 "2022-06-29T08:42:11Z")

</div>

1. Yes! In my own work I also mostly run diagonalizations, and I have implemented in the past something similar to a Workspace for it. We are totally open to adding more functionality, and this would be the first one we’d start with.
2. You are correct. It’s not a bad idea to indeed support resizing. If you want, could you open an issue on it to further discuss what would be the best way to do this? I’m not sure if doing that “under the hood” is the best idea since maybe that would lead to some unexpected behaviors for users, i.e. too much magic. Especially considering that for some of the factorizations work buffers are not trivially determinable but rather require a LAPACK call themselves (maybe looking at what LAPACK does to figure out the size and doing it manually in julia is possible though). I think just having a `resize!` function that would take in a new size or matrix would be a good idea.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [June 29, 2022, 8:52am UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/4 "2022-06-29T08:52:00Z")

</div>

`LinearAlgebra` functions often involve two LAPACK calls, one to determine the workspace size and the second to actually perform the calculation after resizing the workspace. Having the resizing step built in might not be too surprising.

On another note, I guess pre-allocating a workspace will not be thread-safe?

---

<div class="post-metadata">

**Author:** ![louisponet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/louisponet/32/2070_2.png) [@louisponet](https://discourse.julialang.org/u/louisponet)\
**Post date:** [June 29, 2022, 1:39pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/5 "2022-06-29T13:39:12Z")

</div>

> [@jishnub](#):
>
> Having the resizing step built in might not be too surprising.

The whole point of the package is to not have allocations occur during the LAPACK calls.

> [@jishnub](#):
>
> On another note, I guess pre-allocating a workspace will not be thread-safe?

Why not? Using the same workspace with multiple threads of course won’t be indeed.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [June 29, 2022, 2:09pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/6 "2022-06-29T14:09:46Z")

</div>

> [@louisponet](#):
>
> The whole point of the package is to not have allocations occur during the LAPACK calls

I see. My (erroneous) impression was that the package allows one to pre-allocate output matrices, and isn’t fastidious about allocations in resizing the workspace. Thanks for clearing this up!

> [@louisponet](#):
>
> Using the same workspace with multiple threads of course won’t be indeed

Yes, I was referring to using the same workspace from multiple threads, which seems unsafe. I guess the solution is to allocate a separate workspace for each thread, and have each task fetch a workspace from a pool.

---

<div class="post-metadata">

**Author:** ![louisponet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/louisponet/32/2070_2.png) [@louisponet](https://discourse.julialang.org/u/louisponet)\
**Post date:** [June 29, 2022, 2:13pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/7 "2022-06-29T14:13:24Z")

</div>

> [@jishnub](#):
>
> Thanks for clearing this up!

No problem!

> [@jishnub](#):
>
> the solution is to allocate a separate workspace for each thread

Exactly!

---

<div class="post-metadata">

**Author:** ![louisponet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/louisponet/32/2070_2.png) [@louisponet](https://discourse.julialang.org/u/louisponet)\
**Post date:** [July 10, 2022, 9:04pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/8 "2022-07-10T21:04:16Z")

</div>

Hi All,

We’re happy to announce an update to `FastLapackInterface` now including Eigen decompositions!

We added 3 new `Workspaces`:

- `EigenWs`
- `HermitianEigenWs`
- `GeneralizedEigenWs`

We also came up with a more uniform API that allows to construct the correct `Workspace` depending on the target `LAPACK` function, for example:

```julia
ws = Workspace(LAPACK.getrf!, A; kwargs...)
factorize!(ws, A)

```

will automatically create a `QRWYWs` and call `getrf!` through `factorize!`.

Let us know what you think!

Cheers

---

<div class="post-metadata">

**Author:** ![louisponet](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/louisponet/32/2070_2.png) [@louisponet](https://discourse.julialang.org/u/louisponet)\
**Post date:** [August 5, 2022, 12:22pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/9 "2022-08-05T12:22:06Z")

</div>

Hi All,

We have another update to announce, `FastLapackInterface` version `1.2.0`.

`Cholesky` and `BunchKaufman` related `Workspaces` were added, and now each `Workspace` can be resized using `resize!`. Automatic resizing functionality while calling the `LAPACK` functions or `factorize!/decompose!` is turned on by default, but can be controlled with the `resize` keyword argument.

i.e.

```julia
A = rand(3,3)
ws = Workspace(LAPACK.getrf!, A; kwargs...)
B = rand(4,4)
factorize!(ws, B)

```

will automatically resize `ws` appropriately. If instead `factorize!(ws, B, resize=false)` were used, an error would be thrown.

If you have any feedback let us know!

Cheers,  
Louis

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 21, 2024, 7:37am UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/10 "2024-01-21T07:37:31Z")

</div>

Is there an example how to use it in the context of solving a linear system.

I have the following code in a loop (`mC` is different per loop):

```julia
ldiv!(bunchkaufman!(mC), vX);

```

What should I do to make this non allocating using `FastLapackInterface.jl`?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 21, 2024, 1:09pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/11 "2024-01-21T13:09:26Z")

</div>

> [@RoyiAvital](#):
>
> I have the following code in a loop (`mC` is different per loop):
> 
> ```julia
> ldiv!(bunchkaufman!(mC), vX);
> 
> ```
> 
> What should I do to make this non allocating using `FastLapackInterface.jl`?

I haven’t used the package, but from glancing [at the manual](https://dynarejulia.github.io/FastLapackInterface.jl/dev), it seems that you first pre-allocate a workspace with `ws = BunchKaufmanWs(A)`, initializing it with a matrix of the right size (assuming all your matrices are the same size), and then you do `ldiv!(factorize!(wc, mC), vX)`.

There is also a lower-level interface [documented here (with an example)](https://dynarejulia.github.io/FastLapackInterface.jl/dev/workspaces/#BunchKaufman-id).

---

<div class="post-metadata">

**Author:** ![MichelJuillard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/micheljuillard/32/10555_2.png) [@MichelJuillard](https://discourse.julialang.org/u/MichelJuillard)\
**Post date:** [January 21, 2024, 2:27pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/12 "2024-01-21T14:27:51Z")

</div>

Unfortunately, the `factorize!()` method is missing for `BunchKaufman`. We will add it soon.

Here is an example showing what is working and what isn’t yet:

```julia
using FastLapackInterface
using LinearAlgebra

# NOT WORKING
function loop_1!(vXs, mCs, ws)
    for mC in mCs    
        # factorization
        F = factorize!(ws, mC)    
        # solving linear systems
        for vX in vXs
            ldiv!(F, vX)
        end 
    end
end

function loop_2!(vXs, mCs, ws)
    for mC in mCs    
        # factorization    
        A, ipiv, info = LAPACK.sytrf!(ws, 'U', mC)        
        F = BunchKaufman(mC, ipiv, 'U', true, false, BLAS.BlasInt(0))
        # solving linear systems
        for vX in vXs
            ldiv!(F, vX)
        end 
    end
end

mCs = []
vXs = []
n = 10
for i = 1:5
    x = randn(n, n)
    mC = (x + x')/2 
    push!(mCs, mC)
    push!(vXs, randn(n))
end

# create workspace
ws = BunchKaufmanWs(mCs[1])

#= NOT WORKING !
mCs_1 = copy(mCs)
vXs_1 = copy(vXs)
loop_1!(vXs_1, mCs_1, ws)
mCs_1 = copy(mCs)
vXs_1 = copy(vXs)
@time loop_1!(vXs_1, mCs_1, ws)
=#

mCs_2 = copy(mCs)
vXs_2 = copy(vXs)
loop_2!(vXs_2, mCs_2, ws)
mCs_2 = copy(mCs)
vXs_2 = copy(vXs)
@time loop_2!(vXs_2, mCs_2, ws)

```

---

<div class="post-metadata">

**Author:** ![MichelJuillard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/micheljuillard/32/10555_2.png) [@MichelJuillard](https://discourse.julialang.org/u/MichelJuillard)\
**Post date:** [January 21, 2024, 2:35pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/13 "2024-01-21T14:35:40Z")

</div>

I opened an issue here: [Add factorize!() method for BunchKaufman · Issue #38 · DynareJulia/FastLapackInterface.jl · GitHub](https://github.com/DynareJulia/FastLapackInterface.jl/issues/38)

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [January 21, 2024, 6:45pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/14 "2024-01-21T18:45:40Z")

</div>

Let me understand this:

```julia
ws = BunchKaufmanWs(mC); #<! Creates a workspace
A, ipiv, info = LAPACK.sytrf!(ws, 'U', mC); #<! Applies the decomposition
F = BunchKaufman(mC, ipiv, 'U', true, false, BLAS.BlasInt(0)); #<! Creates the data structure Julia knows

```

I was aware of the 2 first steps:

1. Allocating the workspace.
2. Applying the decomposition.

I was missing the 3rd step, where you build the data structure as defined in `LinearAlgebra`.

Am I right about the analysis?

---

<div class="post-metadata">

**Author:** ![MichelJuillard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/micheljuillard/32/10555_2.png) [@MichelJuillard](https://discourse.julialang.org/u/MichelJuillard)\
**Post date:** [January 21, 2024, 6:58pm UTC](https://discourse.julialang.org/t/ann-fastlapackinterface-jl-v1-0-0-non-allocating-lapack-factorizations/83354/15 "2024-01-21T18:58:44Z")

</div>

Yes, you are