OutOfMemoryError on lu decompositions

[alexQueue]

I have a set of sparse matrices of various sizes, and as part of the preconditioner for a larger algorithm (GMRES), I need the LU decomposition of these matrices.

The code is simple.

pc1 = ... # <15654×15654 SparseMatrixCSC{Float64, Int32} with 3327588 stored entries:
pc2 = ... # <31308×31308 SparseMatrixCSC{Float64, Int32} with 13310352 stored entries:

lu(pc1) # works well
lu(pc2) # OutOfMemoryError

Normally this wouldn't be so strange--computers run out of memory operating on large systems. But I attempted this again on a supercomputer with 10x the available memory (190GB), and I still get the same behavior for both matrices.

Additionally, I am able to take the lu decomposition if I first turn it into a dense matrix, so I know that this is a SparseArrays issue rather than a system memory issue.

lu(Matrix(pc2)) # Works, but takes a long time

It seems like there's some internal operation that is running out of memory regardless of the available memory. Is there anything I can do to fix this?

Here is the stacktrace:

ERROR: OutOfMemoryError()
Stacktrace:
 [1] umferror(status::Int32)
   @ SparseArrays.UMFPACK /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:111
 [2] macro expansion
   @ /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:628 [inlined]
 [3] macro expansion
   @ ./lock.jl:267 [inlined]
 [4] umfpack_numeric!(U::SparseArrays.UMFPACK.UmfpackLU{Float64, Int32}; reuse_numeric::Bool, q::Nothing)
   @ SparseArrays.UMFPACK /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:619
 [5] umfpack_numeric!
   @ /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:618 [inlined]
 [6] #lu#7
   @ /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:382 [inlined]
 [7] lu(S::SparseMatrixCSC{Float64, Int32})
   @ SparseArrays.UMFPACK /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:378
 [8] top-level scope
   @ ./timing.jl:273 [inlined]
 [9] top-level scope
   @ ./REPL[53]:0

This is on the 64-bit version of Julia.

[ViralBShah]

Can you provide the matrices?

[SobhanMP]

using SparseArrays, LinearAlgebra
N = 31_308
M = 13_310_352
A = sprandn(N, N, M / N / N);
@time lu(A);

is pretty reasonable?

472.549481 seconds (396.07 k allocations: 35.162 GiB, 0.00% gc time, 0.03% compilation time)

tested on 1.9.0

@Wimmerer don't you agree?

[rayegun]

You're using Int32 indices. It is quite possible that UMFPACK 32-bit is respecting 32-bit memory limits although I will have to verify this by looking through the code.

If indeed that is the issue I can try to update UMFPACK, but it may be a significant assumption made by the library.

Can you retry with 64-bit indices?

[SobhanMP]

It is 64 bits tho 31308×31308 SparseMatrixCSC{Float64, Int64} with 13317263 stored entries:

[alexQueue]

This fixed it, thank you. Converting to 64 bits solved my problem.