Efficient way for assigning a massive array

New to Julia
1

Hello,
I am new to Juila (and programming in general). I am working on solving a massive linear system resulting from a finite difference discretization of an ODE. My matrix is almost a bordered matrix (got some elements on the corners because of periodic BCs), so I defined it as a SparseMatrix. To fill up this system, I need multiple calls to one function that computes the Jacobian at each time step. The system size might reach 10 millions. At the beginning, I thought that the linear solver is taking most of the time, but it turns out that the filling up of the system take a MASSIVE amount of time and allocations. This is one example of the resources needed to fill up only one block on the diagonal of the system.

@time (
for j = 1:M
    s=size(f_dx(j), 1)
    e=(1+(j-1)*s)+(N-1)

    D_matrix[1+(j-1)*s:e  ,  1+(j-1)*s:e,
    ] = -f_dx(j)
end
)

57362.172008 seconds (7.79 G allocations: 18.049 TiB, 2.93% gc time)

I need every possible tip to optimize the assigning procedure, please.

2

First, please provide a self-contained example so that the code you post can be run by someone else. This means e.g. giving the definition of f_dx.

How are you initializing D? If it is with something like spzeros then it will indeed take a long time to fill it up like this because you are changing the sparsity pattern all the time. Instead use the I, J, V constructor with sparse. From the manual:

 sparse(I, J, V,[ m, n, combine])

  Create a sparse matrix S of dimensions m x n such that S[I[k], J[k]] = V[k]. The combine function is used to combine duplicates. If m and n are not specified, they are set to maximum(I) and
  maximum(J) respectively. If the combine function is not supplied, combine defaults to + unless the elements of V are Booleans in which case combine defaults to |. All elements of I must satisfy 1 <=
  I[k] <= m, and all elements of J must satisfy 1 <= J[k] <= n. Numerical zeros in (I, J, V) are retained as structural nonzeros; to drop numerical zeros, use dropzeros!.

  For additional documentation and an expert driver, see SparseArrays.sparse!.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> Is = [1; 2; 3];

  julia> Js = [1; 2; 3];

  julia> Vs = [1; 2; 3];

  julia> sparse(Is, Js, Vs)
  3×3 SparseMatrixCSC{Int64,Int64} with 3 stored entries:
    [1, 1]  =  1
    [2, 2]  =  2
    [3, 3]  =  3
6 Likes
3

Would a matrix free iterative method be appropriate in your problem? Setup a matrix-vector product and use krylov methods to solve the system?

4

The next step would be something like that. But I would like to know how far can I push the direct solver to solve this problem. The solver is not the problem, the matrix assignment is the main problem.

1 Like
5

Yes! I am using spzeros to initialize the matrix. The matrix is being filled in blocks, so I don’t really know the I,J,V before calling f_dx (which is by the way returning the Jacobian - NxN sparse Matrix-).

6

The bottom line of this issue is that SparseMatrixCSC is not a good format to incrementally build up a sparse matrix . This is because adding new values to previously non-stored positions is expensive.

You don’t need to know them before calling f_dx. If f_dx is a matrix, instead of doing

D_matrix[i_range, b_range] = f

you do something like

for (a, i) in enumerate(i_range)
    for (b,j) in enumerate(j_range)
        push!(I, i)
        push!(J, j)
        push(V, f[a,b])
    end
end

...

D_matrix = sparse(I, J, V)
7 Likes
7

I really appreciate your response! Thanks.