Large allocations in `ldiv!` for sparse matrix (with LU factorization)

[superckl]

The method ldiv! has very large allocations happening for sparse LU factorizations. See the example below:

using SparseArrays
using LinearAlgebra

N = 1000

A = sparse(I, N, N) + sprand(N, N, .01)
fact = lu(A)
b = rand(N)
x = zeros(N)

@time ldiv!(x, fact, b)

#Convert to a full matrix and recompute the factorization
A = Matrix(A)
fact = lu(A)

@time ldiv!(x, fact, b);

A characteristic test runs gives

  0.001000 seconds (2 allocations: 46.906 KiB)
  0.000669 seconds (2 allocations: 96 bytes)

The case of the sparse LU decomposition causes a large amount of memory allocation in ldiv! when I would expect it to no different from the case of a dense matrix.

This is on Julia Version 1.6.1.

[ViralBShah]

Moving to SuiteSparse.jl repo.

[superckl]

A quick analysis (with @time) shows the allocations seem to be happening in the call to SuiteSparse here: https://github.com/JuliaLang/SuiteSparse.jl/blob/master/src/umfpack.jl#L414-L416

I see that the method being called here, umfpack_*_solve does indeed allocate two work arrays. The nonallocating version umfpack_*_wsolve requires that one passes those work arrays. I suppose there is no way around this then, since we can't really pass those work arrays through ldiv!. Perhaps they could be stored in the UmfpackLU struct?

[rayegun]

We could allocate them on the Julia side, and if you do a re-solve I suppose they would be reused.

[SobhanMP]

I was considering fixing this.
The umfpack doc mention that there are two arrays that need to be passed, Wi and W. The size of Wi is [n] (same as colptr - 1),
the size of W is a bit more complicated

(page 75, user guide)
i'm not sure what UMFPACK_IRSTEP is. the source says it's control[7] element of the control vector so that's the 8th element in the julia array?

I think given the fact that this changes size given the settings, i think not necessarily a good idea to put in the lu struct. Specially since it's not obvious that using the same decomposition in multiple threads will fail silently. Maybe adding another solve! with 5 argument is a better solution and it's the user's responsibility to allocate the right size arrays?

My main problem with this issue is that it's really bad when the rhs is a matrix not a vector, it typically allocates ~4gb on 10_000x10_000 matrices.

[SobhanMP]

I did a small benchmark with my patch.
old:

new:

i'm not sure why there are 3 allocations now tho :/ will look at it tomorrow.

I think the improvement is important but i still think that any change should not make the lu factorization non thread safe. It would be unexpected. We could have some sort of structure to hold them together. then ldiv! will have 4 arguments. Or a new struct that wraps the actual LU and Wi and W as to not make it thread unsafe by default. This way it's both clear that the lu with Wi and W is not thread safe and we still don't change the number of arguments of ldiv!

Should i implement any of those two alternatives?

code for benchmark:

using Revise, SparseArrays, LinearAlgebra, SuiteSparse, BenchmarkTools, Random

Random.seed!(471);

const A = I + sprandn(1000, 1000, 0.01);
const Af = lu(A);
const b = randn(1000, 1000);
const b1 = randn(1000);

let x = similar(b1)
    @benchmark ldiv!($x, $Af, $b1)
end

let x = similar(b)
    @benchmark ldiv!($x, $Af, $b)
end

[rayegun]

Everything should be held in the UMFPACKFactorization struct. W can be resized if necessary?

The User Guide should give a maximum size. It was my intention to move all of the Control and Info structs into the factorization struct eventually anyway.

[SobhanMP]

Everything should be held in the UMFPACKFactorization struct.

The default solve should not modify the factorization. Think of the poor souls who'll use the same factorization in a thread. nobody expected ldiv! to change the factorization.

[rayegun]

Can you do a PR with your patch? I'm not quite sure what is going on without seeing it. We should be moving the Control and Info structs from global to the UMFPACKFactorization struct, but that's unrelated to this issue for the most part (except thread safety). Outside of that, workspace arrays for solve will indeed have to be passed in or be generated function local (or both: solve!(...; workspace = Vector{T}(...))).

[SobhanMP]

will try