# Inplace version for LU decomposition in Julia

**URL:** <https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795>\
**Category:** General Usage\
**Tags:** linearalgebra, lapack, linear-algebra\
**Created:** [January 14, 2024, 3:07pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795 "2024-01-14T15:07:08Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![mrVeng](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrveng/32/8836_2.png) [@mrVeng](https://discourse.julialang.org/u/mrVeng)\
**Post date:** [January 14, 2024, 3:07pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/1 "2024-01-14T15:07:08Z")

</div>

Hi there!

I have to compute the lu decomposition of a matrix that is dependent - I was thinking to create a buffer first and then inplace compute the lu decomposition when needed, but after googling and searching through some threads, I still have not found the correct function to use with the `LinearAlgebra.jl` package.

Below is a MWE of my use case. I would like to substitute `LinearAlgebra.lu` with something like `LinearAlgebra.lu!` and find out how exactly to assign the buffer that would work with this.

```julia
using Distributions, LinearAlgebra, Random

T=100
n=3
rng = Xoshiro(1)
d = InverseWishart(10, [1.5 0.8 0.5 ; 0.8 2. .5 ; .5 .5 1.] )

for t in Base.OneTo(T)
    mat = rand(d)
    mat_lu = LinearAlgebra.lu(mat)
    # do some computations with mat_lu
end

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [January 14, 2024, 3:28pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/2 "2024-01-14T15:28:10Z")

</div>

Does the [existing `lu!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.lu!) fit your purpose?

---

<div class="post-metadata">

**Author:** ![mrVeng](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrveng/32/8836_2.png) [@mrVeng](https://discourse.julialang.org/u/mrVeng)\
**Post date:** [January 14, 2024, 3:35pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/3 "2024-01-14T15:35:09Z")

</div>

Thanks! That is a bit embarassing :D.

Does this only work with Sparse Matrices or also out of the box with standard LU decomposed matrices as buffer? I.e., the following would not work (continuation from MWE):

```julia
buffer = LinearAlgebra.lu( rand(d) )
mat = rand(d)
LinearAlgebra.lu!(buffer, mat) # MethodError: no method matching lu!(::LU{Float64, Matrix{Float64}, Vector{Int64}}, ::Matrix{Float64})

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [January 14, 2024, 4:02pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/4 "2024-01-14T16:02:42Z")

</div>

I couldn’t link to it but I think the second docstring for `lu!` on this page is more generic. From what I understand, both the `L` and `U` parts can be stored in the matrix itself

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 14, 2024, 4:53pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/5 "2024-01-14T16:53:52Z")

</div>

I think the 2 argument variant only works for sparse matrices and is intended for factorizing another sparse matrix with the same sparsity structure as the previous matrix.

In your case I think you need to keep the reference to `mat` like so:

```julia
mat = rand(d)
buffer = LinearAlgebra.lu(mat)
mat = rand!(mat)
buffer = LinearAlgebra.lu!(mat)

```

You could maybe wrap this logic into a struct to make it more ergonomic to use.

---

<div class="post-metadata">

**Author:** ![mrVeng](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrveng/32/8836_2.png) [@mrVeng](https://discourse.julialang.org/u/mrVeng)\
**Post date:** [January 15, 2024, 9:02am UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/6 "2024-01-15T09:02:18Z")

</div>

> [@abraemer](#):
>
> ```julia
> mat = rand(d)
> buffer = LinearAlgebra.lu(mat)
> mat = rand!(mat)
> buffer = LinearAlgebra.lu!(mat)
> 
> ```

Thanks a lot! This still seems to be allocating though, i.e.:

```julia
using BenchmarkTools

mat = rand(d)
buffer = LinearAlgebra.lu(mat)
mat = rand!(mat)
buffer = LinearAlgebra.lu!(mat)

@btime LinearAlgebra.lu!($mat) #186.438 ns (1 allocation: 80 bytes)

```

I assume this is because `mat` is a standard matrix, while the `LinearAlgebra.lu!` output is a `LU{Float64, Matrix{Float64}, Vector{Int64}}` output.

Is there a way around this? For a single argument function, I assume I would already need to provide a `LU` buffer, but then I would already need to destructure the time dependent covariance matrix into a LU decomposition, which is what I want to avoid in the first place.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [January 15, 2024, 10:10am UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/7 "2024-01-15T10:10:42Z")

</div>

Judging by the source code of `lu!`, I don’t think there is any way around allocation of `ipiv` without manual hacking.  
Note that strictly speaking, the `!` doesn’t guarantee the absence of allocation: it just warns the user that some inputs might be overwritten.

> <https://github.com/JuliaLang/julia/blob/fc6295df630603bbe4ada01f471316fb457240ae/stdlib/LinearAlgebra/src/lu.jl#L134-L193>

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [January 15, 2024, 2:21pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/8 "2024-01-15T14:21:07Z")

</div>

[`FastLapackInterface.jl`](https://github.com/DynareJulia/FastLapackInterface.jl) provides ways to separate memory allocation from (a subset of) LAPACK computations. It allows to create a “workspace” which can be reused for several calls to the same functions on similar matrices.

The documentation contains an example of this for LU factorization (i.e. re-using the same “workspace” for multiple calls to `getrf!`):  
[https://dynarejulia.github.io/FastLapackInterface.jl/dev/workspaces/#LU-id](https://dynarejulia.github.io/FastLapackInterface.jl/dev/workspaces/#LU-id)

---

<div class="post-metadata">

**Author:** ![j-fu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j-fu/32/11373_2.png) [@j-fu](https://discourse.julialang.org/u/j-fu)\
**Post date:** [January 15, 2024, 2:55pm UTC](https://discourse.julialang.org/t/inplace-version-for-lu-decomposition-in-julia/108795/9 "2024-01-15T14:55:06Z")

</div>

RecursiveFactorizations.jl has a non-allocating `lu!` where you can pass `ipiv`.