# Reuse Tridiagonal LU decomposition

**URL:** <https://discourse.julialang.org/t/reuse-tridiagonal-lu-decomposition/115149>\
**Category:** General Usage\
**Created:** [June 4, 2024, 3:24am UTC](https://discourse.julialang.org/t/reuse-tridiagonal-lu-decomposition/115149 "2024-06-04T03:24:38Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![bclyons12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bclyons12/32/33734_2.png) [@bclyons12](https://discourse.julialang.org/u/bclyons12)\
**Post date:** [June 4, 2024, 3:24am UTC](https://discourse.julialang.org/t/reuse-tridiagonal-lu-decomposition/115149/1 "2024-06-04T03:24:38Z")

</div>

I need to perform a loop of linear solves with different matrices and vectors, but the same Tridiagonal structure. I’d like to do this without allocations in the loop. Is there anyway to preallocate the Tridiagonal LU data structure and then reuse those arrays inside the loop. Here’s a MWE showing what I would like to work, but how it must work now:

```julia
using LinearAlgebra

N = 10
Dl, D, Du = rand(N-1), rand(N), rand(N-1)
A = Tridiagonal(Dl, D, Du)
L = lu(A)
x = zeros(N)
b = rand(N)
for k in 1:3
    D .= rand(N) # this updates the diagonal of A
    b .= rand(N)
    L = lu(A) # this allocates; ** would like lu!(L, A)
    ldiv!(x, L, b)
end

```

Perhaps this is better done with LinearSolve.jl, but I haven’t found a way to make that completely allocation free either.

---

<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 4, 2024, 6:32am UTC](https://discourse.julialang.org/t/reuse-tridiagonal-lu-decomposition/115149/2 "2024-06-04T06:32:06Z")

</div>

Perhaps you could try `lu!(A)`, which operates in-place on `A`? This isn’t allocation-free either, but it allocates less. Also, you may initialize the tridiagonal matrix with a second off-diagonal, which is used in the `lu!` computation:

```julia
julia> A = Tridiagonal(Dl, D, Du, similar(Du, length(Du)-1));

julia> @btime lu!($A);
  350.357 ns (1 allocation: 144 bytes)

```

Unrelated, but if allocations are your concern, you may also use `rand!(D)` instead of `D .= rand(N)`, as the vector won’t be allocated in the former. Similarly, for `b`.

Incidentally, julia v1.11 will include an `lu!` that accepts the factorization and updates it in-place, so this will be (almost) allocation-free.

```julia
julia> A = Tridiagonal(Dl, D, Du);

julia> L = lu(A);

julia> @btime lu!($L, $A);
  146.500 ns (1 allocation: 48 bytes)

```

I’m not quite sure why there’s the one allocation, but it’s almost negligible.

---

<div class="post-metadata">

**Author:** ![bclyons12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bclyons12/32/33734_2.png) [@bclyons12](https://discourse.julialang.org/u/bclyons12)\
**Post date:** [June 7, 2024, 7:17pm UTC](https://discourse.julialang.org/t/reuse-tridiagonal-lu-decomposition/115149/3 "2024-06-07T19:17:10Z")

</div>

Thanks @jishnub . The julia v1.11 version looks like what I’d like. Just doing `lu!(A)` won’t work because I want to update A in a loop, so I’d need to reallocate the `Tridiagonal` each iteration.

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [June 8, 2024, 6:20am UTC](https://discourse.julialang.org/t/reuse-tridiagonal-lu-decomposition/115149/4 "2024-06-08T06:20:59Z")

</div>

> [@jishnub](#):
>
> I’m not quite sure why there’s the one allocation, but it’s almost negligible.

the allocation is one `Vector{Int}`(for the pivot indices). at the moment, there is no julia function that allows to do `lu!(A,piv,ipiv::Vector{Int})`