# Efficient way for assigning a massive array

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.

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

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

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

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-).

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

I really appreciate your response! Thanks.