OutOfMemoryError() while solving a linear system

General Usage Performance
1

I am trying to solve a linear system with the direct solver (backslash operator). My matrix is sparse,square and unsymmetric. The system is defined as Ax=b, where A is defined to be a SparseMatrixCSC with dimensions = 803800 x 803800 nnz(A_matrix)=808613995. the b vector is dense and defined to be like that. My A matrix is very close to be a banded matrix, but it is not, due to some elements on the corners, so I believe that the solver is LU factorizing the matrix and then solving the system. I am running my code on a 16 GB RAM laptop. My question is how to know the size limit of my linear system which I can solve on this machine given the sparsity pattern of my matrix?

2

Think of the memory consumption of such a system: you would need

julia> (808613995 * 16) / (2^30)
12.049287483096123

gigabytes to represent the nonzero entries and the row indexes (Float64, and Int).

Maybe someone can suggest something, but I am not sure that doing this in 16GB is feasible.

Consider even with iterative methods with lower memory requirements, eg

but even these would need more RAM I assume.

2 Likes
3

Thanks !

This is true. Unless you have a good preconditioner, IterativeSolver will take ages and also massive memory.

4

Normally, a Krylov method only allocates few n-vectors, where n is the dimension of your matrix. With n = 803800, a n-vector only needs 6.13Mo. So iterative methods could be relevant in your case.

I can recommend you to test DQGMRES method, you can define the number of n-vectors that you want to use during the Arnoldi process. It’s a kind of limited memory variant of GMRES.

(x, stats) = dqgmres(A, b, memory=400, itmax=10000, verbose=true)

You can adjust parameters (set verbose to false and/or reduce/increase memory).

Do you have your linear system available somewhere or a way to generate it that I can give a try too ?

1 Like
5

Thanks for your suggestion! I tried it but it is not actually working with A being sparse. I tried to convert A back to the dense form using Matrix(A) (which I know is not a good idea), and I got an OutOfMemoryError.

6

Have you tried with Pardiso.jl and MKL?

1 Like
7

@AhmedAlreweny, method should work with A being sparse. What is the error that is returned?

using SparseArrays
A = spzeros(100,100)
for i = 1 : 100
    A[i,i] = 100
end
b = ones(100)
dqgmres(A, b)

([0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01  …  0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01], 
Simple stats
  solved: true
  inconsistent: false
  residuals:  [ 1.0e+01  5.4e-16 ]
  Aresiduals: [ ]
  status: solution good enough given atol and rtol
)

Otherwise, you can use MUMPS that has been developed for linear systems of high dimensions.

8

Thanks again! This is the error I got

dqgmres(A_matrix, b_vector)
ERROR: MethodError: no method matching dqgmres(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,2})
Closest candidates are:
dqgmres(::AbstractArray{T<:AbstractFloat,2}, ::SparseVector{T<:AbstractFloat,Ti} where Ti<:Integer, ::Any…; kwargs…) where T<:AbstractFloat

9

Ok, I see your problem.
Your right-hand side b is not a vector but a matrix Array{Float64,2} with one column.
You can convert it to a vector with b = b[:].

1 Like
10

Nope! I will go through the documentation thanks!

11

It seems that a preconditioner is required. It takes very long time without converging to the solution. Thanks for pointing out!

1 Like