# Reuse Tridiagonal LU decomposition

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:

``````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.

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> 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> 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.

3 Likes

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.

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})`

2 Likes