# Solving a linear system without modifying and allocating arrays

**URL:** <https://discourse.julialang.org/t/solving-a-linear-system-without-modifying-and-allocating-arrays/127315>\
**Category:** Performance\
**Created:** [March 24, 2025, 5:55pm UTC](https://discourse.julialang.org/t/solving-a-linear-system-without-modifying-and-allocating-arrays/127315 "2025-03-24T17:55:36Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![cneverett](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cneverett/32/209565_2.png) [@cneverett](https://discourse.julialang.org/u/cneverett)\
**Post date:** [March 24, 2025, 5:55pm UTC](https://discourse.julialang.org/t/solving-a-linear-system-without-modifying-and-allocating-arrays/127315/1 "2025-03-24T17:55:37Z")

</div>

Hi, I have a setup where I need to repeatedly solve a linear system like A \ B = C every timestep of a system. A is a matrix, B is a vector and both changing each timestep and more importantly A is reused later each timestep for other functions like a matrix multiplication e.g. D = A \* C.

What I assume is a naive function might look like:

```julia
using LinearAlgebra
using BenchmarkTools

A = rand(Float32,10,10)
A_copy = similar(A)
B = rand(Float32,10)
C = similar(B)
D = similar(B)

function lu_div_mul!(A::Array{Float32,2},A_copy::Array{Float32,2},B::Array{Float32,1},C::Array{Float32,1},D::Array{Float32,1})
    copy!(A_copy,A)
    LU = lu!(A_copy)
    ldiv!(C,LU,B)
    mul!(D,A,C)
end

```

```julia-repl
@btime lu_div_mul!($A,$A_copy,$B,$C,$D)
  843.077 ns (2 allocations: 144 bytes)

```

where the two allocations come from the lu factorisation. I would like to avoid these allocations and ideally the unneeded copy of A as in reality the dimensions of A can be rather large. Is there an existing function or method, something like lu!(A\_lu,A), where the factorisation A\_lu can be pre-allocated and reused, avoiding the matrix copy?

---

<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:** [March 24, 2025, 6:08pm UTC](https://discourse.julialang.org/t/solving-a-linear-system-without-modifying-and-allocating-arrays/127315/2 "2025-03-24T18:08:16Z")

</div>

You can probably do that with LinearSolve.jl but it would be nice to have it in `LinearAlgebra`

---

<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:** [March 24, 2025, 6:31pm UTC](https://discourse.julialang.org/t/solving-a-linear-system-without-modifying-and-allocating-arrays/127315/3 "2025-03-24T18:31:12Z")

</div>

> [@cneverett](#):
>
> Is there an existing function or method, something like lu!(A\_lu,A), where the factorisation A\_lu can be pre-allocated and reused, avoiding the matrix copy?

Yes, this feature was recently added: [define lu!(F::LU,A) for general matrices by longemen3000 · Pull Request #1131 · JuliaLang/LinearAlgebra.jl · GitHub](https://github.com/JuliaLang/LinearAlgebra.jl/pull/1131)

I’m not sure how long it will be until this appears in a Julia release, but the code from that PR is pretty short so you could easily copy it.

---

<div class="post-metadata">

**Author:** ![cneverett](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cneverett/32/209565_2.png) [@cneverett](https://discourse.julialang.org/u/cneverett)\
**Post date:** [March 25, 2025, 10:48am UTC](https://discourse.julialang.org/t/solving-a-linear-system-without-modifying-and-allocating-arrays/127315/4 "2025-03-25T10:48:41Z")

</div>

That’s excellent! Based on the docstring looks like this will come with Julia v1.12, until then I shall just copy the code as it is as you say, pretty short.

Thank you!