How to preallocate LU factorization

Specific Domains Numerics
1

OK, I see what you mean.

Suppose, then, that I’m not using StaticArrays and just want to allocate an empty LU with normal Julia arrays—how do I do that?

zeros(LU{Float64, zeros(Matrix{Float64}), zeros(Vector{Int64})})

is an error. (I realize this is a different question.)

Edit: Kind of hacky:

LU(zeros(3,3), [1, 2, 3], 0)
2

The easiest thing is to just do F = lu(rand(n,n)). The only reason to pre-allocate is if you are going to re-use the factorization object many times, in which case doing one extra factorization at the beginning has negligible cost.

3 Likes
3

I don’t think Julia exposes any way to preallocate the permutation vector. Perhaps you already knew this, but the space for the L and U outputs can be “preallocated” by using the lu! function. lu! uses the input array’s storage to store the L and U components. However, the original input matrix is mutilated in the process and, further, its memory is now shared with L and U so can’t be used for other purposes until you’re done with those.
This can be seen in this demonstration

julia> x = reshape(1.0:9,3,3).^2
3×3 Matrix{Float64}:
 1.0  16.0  49.0
 4.0  25.0  64.0
 9.0  36.0  81.0

julia> lu!(x)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
 1.0       0.0   0.0
 0.111111  1.0   0.0
 0.444444  0.75  1.0
U factor:
3×3 Matrix{Float64}:
 9.0  36.0  81.0
 0.0  12.0  40.0
 0.0   0.0  -2.0

julia> x # L in the lower triangle and U in the upper 
3×3 Matrix{Float64}:
 9.0       36.0   81.0
 0.111111  12.0   40.0
 0.444444   0.75  -2.0

There are many ___! functions in LinearAlgebra that follow this same pattern.

4

You may think of the cases you want to build a library based on Julia code and one of the requirements is to be non allocating code.

I can think of many more case the option to have API which exposes all memory allocations is useable.