Reuse Tridiagonal LU decomposition

General Usage
1

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.

2

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
3

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.

4

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