# Improving performance for matrix assembly

**URL:** <https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252>\
**Category:** General Usage\
**Created:** [February 8, 2023, 4:59am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252 "2023-02-08T04:59:41Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![udax](https://avatars.discourse-cdn.com/v4/letter/u/a587f6/32.png) [@udax](https://discourse.julialang.org/u/udax)\
**Post date:** [February 8, 2023, 4:59am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/1 "2023-02-08T04:59:41Z")

</div>

I have a working code for a large project in scientific computation. Compared to similar existing codes in C, the julia code is around 10 times slower, which is mostly likely due to me being new to Julia. I’ve tried profiling the code but frankly the output is quite overwhelming and I can’t seem to make much traction there.

From what I could distill from the various other posts on this and SO, here is my question with a MWE below: How would you make this code for matrix assembly faster? Please note that this is (obviously) a dummy scaffold, with the real `ker` function being much more complicated.

```julia
#use to set a matrix block
function ker(n)
    return rand(n,n)
end

#main body that assembles a matrix block by block
#each block is of size n x n
function mbody()
    m = 100
    n = 5
    A = zeros(n*m,n*m)
    b = zeros(n*m)
    for i=1:m
        i_idx = 1+(i-1)*n:i*n
        for j=1:m
            j_idx = 1+(j-1)*n:j*n
            A[i_idx,j_idx] = ker(n)
        end
        b[i_idx] = rand(n)
    end
end

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 8, 2023, 5:02am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/2 "2023-02-08T05:02:00Z")

</div>

the slow thing here is that you’re copying a ton of giant arrays and throwing them away. instead, pass a view into ker and fill that.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [February 8, 2023, 5:13am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/3 "2023-02-08T05:13:10Z")

</div>

Concretely,

```julia
julia> function ker!(A)
          for i in eachindex(A)
           A[i] = rand()
          end
       end
ker! (generic function with 2 methods)

julia> function mbody2()
           m = 100
           n = 5
           A = zeros(n*m,n*m)
           b = zeros(n*m)
           for i=1:m
               i_idx = 1+(i-1)*n:i*n
               for j=1:m
                   j_idx = 1+(j-1)*n:j*n
                   ker!(@views A[i_idx,j_idx])
               end
               b[i_idx] = rand(n)
           end
       end
mbody2 (generic function with 1 method)

julia> @btime mbody()
  2.373 ms (10103 allocations: 4.36 MiB)

julia> @btime mbody2()
  1.089 ms (103 allocations: 1.92 MiB)

```

---

<div class="post-metadata">

**Author:** ![udax](https://avatars.discourse-cdn.com/v4/letter/u/a587f6/32.png) [@udax](https://discourse.julialang.org/u/udax)\
**Post date:** [February 8, 2023, 6:09am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/4 "2023-02-08T06:09:11Z")

</div>

Thanks a lot @stillyslalom! The improvement in memory and time is great. When I tried to implement this idea, I ran into some trouble because in my case, the function `ker` accepts multiple argument. Concretely:

(1) If I replace `ker!(@views A[i_idx,j_idx])` by `ker!(@views A[i_idx,j_idx],4)` where `4` is an extra argument (as an example) and modify the function definition to be `ker!(A,z)` I get an error of the sort: `ERROR: MethodError: no method matching setindex!(::Tuple{SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}, Int64}, ::Float64, ::Int64)`.

(2) On the other hand, if I replace `ker!(@views A[i_idx,j_idx])` by `ker!(A,i_idx,j_idx,4)` and within the body of `ker!(A,idx1,idx2,z)` do `A[idx1,idx2] .= rand()+z`, I get no errors and the same improvement in performance as you pointed out in your post.

Can you help me understand why (1) failed, and if (2) is just as good as your original solution? Thanks!

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [February 8, 2023, 6:38am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/5 "2023-02-08T06:38:15Z")

</div>

`@views` seems to be doing something with the other argument.  
You could try to just use `@view` (which should only convert the next expression into a view) or assign the view to a variable before, which usually makes things more readable anyway:

```julia
Aview = @view A[...]
ker!(Aview,4)

```

---

<div class="post-metadata">

**Author:** ![udax](https://avatars.discourse-cdn.com/v4/letter/u/a587f6/32.png) [@udax](https://discourse.julialang.org/u/udax)\
**Post date:** [February 8, 2023, 7:38am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/6 "2023-02-08T07:38:50Z")

</div>

Thanks @Salmon. This does work, but doesn’t seem as optimal. If I follow the approach suggested by me in (2) above, i.e. send in `A` along with the indexing range, then I observe the following (using `@benchmark`):  
a) Memory and allocs estimate is identical  
b) The approach using `Aview` is slower by about 15-20%  
Why would this be happening?

Update: The solution was as simple as to wrap `@view .. ` inside `()` brackets. It has comparable resource utilization as the approach using `Aview`. Thanks to this [link](https://www.juliabloggers.com/the-view-and-views-macros-are-you-sure-you-know-how-they-work/).

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [February 8, 2023, 8:14am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/7 "2023-02-08T08:14:12Z")

</div>

> [@udax](#):
>
> b) The approach using `Aview` is slower by about 15-20%  
> Why would this be happening?

That is weird, afaik the performance should be nearly identical, assuming you are iterating over the correct memory layout in the `ker!` function.

> [@udax](#):
>
> Update: The solution was as simple as to wrap `@view .. ` inside `()` brackets. It has comparable resource utilization as the approach using `Aview`. Thanks to this [link](https://www.juliabloggers.com/the-view-and-views-macros-are-you-sure-you-know-how-they-work/).

Do you mean if you wrap the `@view` in brackets the performance loss is gone?  
Im pretty sure

```julia
ker!( (@view A[...]),4)

```

should be lowered to

```julia
Aview = (@view A[...]
ker!( Aview,4)

```

so Id expect these things to do exactly the same.

---

<div class="post-metadata">

**Author:** ![jmair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmair/32/35117_2.png) [@jmair](https://discourse.julialang.org/u/jmair)\
**Post date:** [February 8, 2023, 8:28am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/8 "2023-02-08T08:28:31Z")

</div>

Building on @stillyslalom’s implementation, I have changed the indexing order (also put everything in the function), but I don’t think it makes much of a difference:

```julia
julia> function mbody3()
                  m = 100
                  n = 5
                  A = zeros(n*m,n*m)
                  b = rand(n*m)
                  @inbounds for j=1:m
                      j_idx = 1+(j-1)*n:j*n
                      for i=1:m
                           i_idx = 1+(i-1)*n:i*n
                           for jj in j_idx
                               for ii in i_idx
                                   A[ii, jj] = rand()
                               end
                           end
                      end
                  end
                  nothing
              end

```

As Julia is column-major, it’s best to get in the habit of iteration over the left-most index first (i.e. swap the i and j from the example). This can make a difference in performance sometimes, but I haven’t run this code on a desktop machine to check.

Also, if you need to set parts of b in the loop, you can get an in-place `rand!` from the `Random` standard library.

---

<div class="post-metadata">

**Author:** ![udax](https://avatars.discourse-cdn.com/v4/letter/u/a587f6/32.png) [@udax](https://discourse.julialang.org/u/udax)\
**Post date:** [February 8, 2023, 8:41am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/9 "2023-02-08T08:41:39Z")

</div>

Sorry, I should have been clearer. What I meant was that `ker!(@views A[i_idx,j_idx],4)` gave an error, but doing `ker!((@views A[i_idx,j_idx]),4)` resolved the error. The performance gap is still there.

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [February 8, 2023, 9:09am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/10 "2023-02-08T09:09:45Z")

</div>

> [@udax](#):
>
> `A[idx1,idx2] .= rand()+ z`

Ah, I think I understand:  
The following two functions do not do the same:

```julia
function ker!(A,z)
    for i in eachindex(A)
        A[i] = rand() + z
    end
end
function ker!(A,z)
    A .= rand()+z
end

```

The first function generates a random number for each index of `A` and writes it to it.  
The second function generates **one** random number and fills `A` with it.  
Since the second function only has to generate one random number, it is likely faster.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 8, 2023, 10:15am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/11 "2023-02-08T10:15:18Z")

</div>

> [@udax](#):
>
> `ker!((@views A[i_idx,j_idx]),4)` resolved the error.

You should use `@view`, not `@views`. You can also use parens:

```julia
ker!(@view(A[i_idx,j_idx]),4)

```

Note that you can use parens the same way as with function calls.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [February 8, 2023, 1:18pm UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/13 "2023-02-08T13:18:31Z")

</div>

Sorry for the confusion–I accidentally copy-pasted that broken definition from the REPL. The version I actually ran for timings includes parentheses around `@views`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 8, 2023, 1:21pm UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/14 "2023-02-08T13:21:12Z")

</div>

> [@DNF](#):
>
> You should use `@view`, not `@views`. You can also use parens:

Or just put `@views` in front of `function` to make every slice in the function use views. Or, if you want to annotate use of views at a finer-grained level, I think it is clearer to write:

```julia
@views ker!(A[i_idx,j_idx],4)

```

(I think a lot of people still don’t grok how the `@views` macro allows you to opt-in to views for arbitrary block of code.)

> [@Salmon](#):
>
> `A .= rand()+z`

Tip: If you change this to

```julia
@. A = rand() + z

```

which is shorthand for `A .= rand.() .+ z`, then it will call `rand()` separately for each element of `A`, just as for your explicit loop.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 8, 2023, 2:14pm UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/15 "2023-02-08T14:14:21Z")

</div>

> [@stevengj](#):
>
> Or just put `@views` in front of `function` to make every slice in the function use views.

Yes, but it is better to use the `@view` macro when applying it locally, than to use `@views` and then awkwardly limiting its scope with `(@views A[i, j])`.

---

<div class="post-metadata">

**Author:** ![uniment](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uniment/32/24532_2.png) [@uniment](https://discourse.julialang.org/u/uniment)\
**Post date:** [February 9, 2023, 2:42am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/16 "2023-02-09T02:42:34Z")

</div>

Is there a practical difference between `@view(A[i, j])` and `@views(A[i, j])`?

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 9, 2023, 6:30am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/17 "2023-02-09T06:30:00Z")

</div>

I don’t think so, it’s about being idiomatic, mostly.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [February 9, 2023, 10:12am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/18 "2023-02-09T10:12:52Z")

</div>

> [@jmair](#):
>
> ```julia
> for jj in j_idx
> for ii in i_idx
> A[ii, jj] = rand()
> end
> end
> 
> ```

The quoted code calls `rand()` for each element. Generating random numbers takes a bit of effort and has specialized code for generating many numbers at once, hence:

```julia
rand(nrows, ncols)

```

is significantly faster for larger matrices. In this case:

```julia
A = rand(n*m,n*m)

```

would accomplish the operation of the whole loop.

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [February 9, 2023, 10:46am UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/19 "2023-02-09T10:46:49Z")

</div>

In this specific case no, but for multi-argument macros, arguments need to be separated by comma instead of space

```julia
(@macroname arg1 arg2 arg3)
@macroname(arg1, arg2, arg3)

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 9, 2023, 12:00pm UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/20 "2023-02-09T12:00:57Z")

</div>

> [@Dan](#):
>
> In this case:
> 
> ```julia
> A = rand(n*m,n*m)
> 
> ```
> 
> would accomplish the operation of the whole loop.

Firstly, the actual `ker` does not return random numbers, so this is slightly beside the point. But, accepting `rand` for the sake of argument:

The above re-assigns the variable `A`, while the loop you quoted works in-place. Of course, you can write `A .= rand(n*m, n*m)`. So this is faster in terms of the `rand` part, but you need to allocate a temporary array, and then copy over the data. It’s not clear to me that this is always faster.

In terms of the sizes in question, `ker` returns a 5x5 matrix, so:

```julia
julia> @btime A[301:305, 301:305] .= rand(5, 5) setup=(A=rand(500,500));
  145.933 ns (1 allocation: 256 bytes)

julia> @btime A[301:305, 301:305] .= rand.() setup=(A=rand(500,500));
  41.616 ns (0 allocations: 0 bytes)

```

Even for larger sizes, this holds up.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [February 9, 2023, 12:10pm UTC](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252/21 "2023-02-09T12:10:05Z")

</div>

> [@DNF](#):
>
> Firstly, the actual `ker` does not return random numbers, so this is slightly beside the point.

The comment was about the general performance different between generating many random values at once vs. one-by-one, which exists. The kernel would obviously not be random in applications.

> [@Dan](#):
>
> `A = rand(n*m,n*m)`

Replaces the whole loop and the allocation of `A` with `zeros`. Therefore there is no need for `.=` syntax for in-place operation. The initial zero-filling is unnecessary as well in this case.

Finally, regarding the `5x5` size, the whole `A` is larger, and my comment specifically says it applies to larger matrices:

```julia
julia> function test()
       A = Matrix{Float64}(undef, 100,100)
       for i in 1:100
           for j in 1:100
               A[j,i] = rand()
           end
       end
       A
       end
test (generic function with 1 method)

julia> @btime test();
  32.231 μs (2 allocations: 78.17 KiB)

julia> @btime rand(100,100);
  11.287 μs (2 allocations: 78.17 KiB)

```

this 3x difference is substantial (even though most of it is due to `@inbounds` missing - though still 50% difference with `@inbounds`).

Finally, this intuition regarding bunching of random generation for efficiency is useful to have. But it may not be pertinent to this case.

[Next page](https://discourse.julialang.org/t/improving-performance-for-matrix-assembly/94252.md?page=2)