# Efficient way for assigning a massive array

**URL:** <https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809>\
**Category:** New to Julia\
**Tags:** array\
**Created:** [January 26, 2020, 2:33pm UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809 "2020-01-26T14:33:15Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![AhmedAlreweny](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@AhmedAlreweny](https://discourse.julialang.org/u/AhmedAlreweny)\
**Post date:** [January 26, 2020, 2:33pm UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809/1 "2020-01-26T14:33:15Z")

</div>

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.

```julia
@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
)

```

```julia
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.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [January 26, 2020, 2:49pm UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809/2 "2020-01-26T14:49:23Z")

</div>

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:

```julia
 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

```

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [January 26, 2020, 10:36pm UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809/3 "2020-01-26T22:36:00Z")

</div>

> [@AhmedAlreweny](#):
>
> massive linear system resulting from a finite difference discretization of an ODE

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

---

<div class="post-metadata">

**Author:** ![AhmedAlreweny](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@AhmedAlreweny](https://discourse.julialang.org/u/AhmedAlreweny)\
**Post date:** [January 27, 2020, 7:54am UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809/4 "2020-01-27T07:54:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![AhmedAlreweny](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@AhmedAlreweny](https://discourse.julialang.org/u/AhmedAlreweny)\
**Post date:** [January 27, 2020, 8:00am UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809/5 "2020-01-27T08:00:07Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [January 27, 2020, 8:31am UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809/6 "2020-01-27T08:31:29Z")

</div>

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.

> [@AhmedAlreweny](#):
>
> The matrix is being filled in blocks, so I don’t really know the `I,J,V` before calling f\_dx

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

```julia
D_matrix[i_range, b_range] = f

```

you do something like

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![AhmedAlreweny](https://avatars.discourse-cdn.com/v4/letter/a/db5fbb/32.png) [@AhmedAlreweny](https://discourse.julialang.org/u/AhmedAlreweny)\
**Post date:** [January 27, 2020, 10:32am UTC](https://discourse.julialang.org/t/efficient-way-for-assigning-a-massive-array/33809/7 "2020-01-27T10:32:33Z")

</div>

I really appreciate your response! Thanks.