# Inplace multiplication of sub-matrices without allocations

**URL:** <https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534>\
**Category:** Performance\
**Tags:** performance, linear-algebra, allocations\
**Created:** [December 11, 2022, 9:23pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534 "2022-12-11T21:23:18Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [December 11, 2022, 9:23pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/1 "2022-12-11T21:23:18Z")

</div>

I am writing a performance-critical part and would like to avoid any memory allocations in that part by creating caches beforehand. I then essentially need to multiply sub-arrays in a loop. Is there any way to do these matrix multiplications of sub-arrays without any allocations?  
Here is a MWE:

```julia
using LinearAlgebra

A = rand(4,10,10)
B = rand(10,10)

C = zeros(4,10,10)

# This is the function I plan to loop over
function foo(A,B,C,dim)
    mul!(C[dim,1:10,1:10], A[dim,1:10,1:10], B)
    nothing
end

foo(A,B,C,1)

@time foo(A,B,C,1)
# 0.000010 seconds (2 allocations: 1.750 KiB)

# This is the function that loops over the first dimension
function wrapper_foo(A,B,C)
    for i in 1:4
        foo(A,B,C,i)
    end
    nothing
end

wrapper_foo(A,B,C)
@time wrapper_foo(A,B,C)
# 0.000037 seconds (12 allocations: 85.500 KiB)

```

In situations where I don’t use sub-arrays, mul!() works without allocating any temporary arrays. Any ideas?

---

<div class="post-metadata">

**Author:** ![pitsianis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitsianis/32/26588_2.png) [@pitsianis](https://discourse.julialang.org/u/pitsianis)\
**Post date:** [December 11, 2022, 9:36pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/2 "2022-12-11T21:36:07Z")

</div>

Use `view` or `@view` with the submatrix expression.

Also, if possible, make `dim` the third dimension in your arrays, to avoid strides.

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [December 11, 2022, 9:53pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/3 "2022-12-11T21:53:56Z")

</div>

Thanks a lot for the fast reply. I already tried to use `@view` but it just increases the number of allocations from 2 to 3 while increasing the amount of allocated memory substantially.  
Thanks for the suggestion about `dim`. My actual problem has a six-dimensional array. With the last two ones not being looped over. Should I change the dimensions such that I loop over dim 3-6 instead?

---

<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:** [December 11, 2022, 9:58pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/4 "2022-12-11T21:58:26Z")

</div>

For multi-dimensional tensor contractions like this, you should probably use `Tullio` with `LoopVectorization`. Decomposing the operation into 2D operations leaves a bunch of performance on the table.

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [December 11, 2022, 10:02pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/5 "2022-12-11T22:02:32Z")

</div>

Here is the code about using `@view`:

```julia
using LinearAlgebra, BenchmarkTools
function foo(A,B,C,dim)
    mul!(C[dim,1:10,1:10], A[dim,1:10,1:10], B)
    nothing
end

function foo2(A,B,C,dim)
    mul!(C[dim,1:10,1:10], @view(A[dim,1:10,1:10]), @view(B[1:10,1:10]))
    nothing
end

foo(A,B,C,1)
foo2(A,B,C,1)

@btime foo(A,B,C,1)
# 504.167 ns (2 allocations: 1.75 KiB)
@btime foo2(A,B,C,1)
# 1.020 μs (3 allocations: 21.38 KiB)  

```

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [December 11, 2022, 10:07pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/6 "2022-12-11T22:07:28Z")

</div>

Thanks for the suggestion. I parallelize the wrapping loop(s) with `Threads.@spawn`. Do you think `LoopVectorization` would be a better fit? Inside the `foo` function happen a few more things than just the matrix multiplication. But most of this results in zero allocations. `Toolio.jl` looks great but I am not sure how I would write the program with it. Would you mind giving me an example with my MVE?  
Thanks a lot!

---

<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:** [December 11, 2022, 10:14pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/7 "2022-12-11T22:14:50Z")

</div>

Tullio works by using [Einstein notation](https://en.wikipedia.org/wiki/Einstein_notation), so this program would turn into

```julia
using LoopVectorization, Tullio
function wrapper_foo(A,B,C)
    @tullio C[i, j, k] = A[i, j, p] * B[p, k]
end

```

---

<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:** [December 11, 2022, 10:15pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/8 "2022-12-11T22:15:25Z")

</div>

View is less performant when taking strides, but more performant when not.

```julia
julia> function foo!(A,B,C,dim)
           mul!(C[dim,1:10,1:10], A[dim,1:10,1:10], B)
           nothing
       end
foo! (generic function with 1 method)

julia> function foo2!(A,B,C,dim)
           mul!(C[dim,1:10,1:10], @view(A[dim,1:10,1:10]), @view(B[1:10,1:10]))
           nothing
       end
foo2! (generic function with 1 method)

julia> function foo3!(A,B,C,dim)
           mul!(C[:,:,dim], @view(A[:,:,dim]), B)
           nothing
       end
foo3! (generic function with 1 method)

julia> function foo4!(A,B,C,dim)
           mul!(@view(C[:,:,dim]), @view(A[:,:,dim]), B)
           nothing
       end
foo4! (generic function with 1 method)

julia> function foo5!(A,B,C,dim)
           mul!(C[:,:,dim], A[:,:,dim], B)
           nothing
       end
foo5! (generic function with 1 method)

julia> @btime foo!(A,B,C,1)
  1.140 μs (2 allocations: 1.75 KiB)

julia> @btime foo2!(A,B,C,1)
  2.022 μs (3 allocations: 21.38 KiB)

julia> @btime foo3!(A,B,C,1)
  450.289 ns (1 allocation: 400 bytes)

julia> @btime foo4!(A,B,C,1)
  324.885 ns (0 allocations: 0 bytes)

julia> @btime foo5!(A,B,C,1)
  526.380 ns (2 allocations: 800 bytes)

```

---

<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:** [December 11, 2022, 10:19pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/9 "2022-12-11T22:19:23Z")

</div>

Out of curiosity, why do we use `view(A, :, :, dim)` instead of `View(A)[:, :, dim]`? It seems like much of the motivation for the `@view` macro is to allow the `begin` and `end` keywords which are valid in `[]` square brackets, but if we simply had an indexable type we get those back.

Playing around, it seems like this should be doable:

```julia
struct View{A<:AbstractArray} a::A end
Base.getindex(v::View, args...) = view(v.a, args...)
Base.axes(v::View, args...) = axes(v.a, args...)
View(A)[:, :, end] # this just works

```

Applying it here:

```julia
julia> function foo6!(A,B,C,dim)
           mul!(View(C)[:,:,dim], View(A)[:,:,dim], B)
           nothing
       end
foo6! (generic function with 1 method)

julia> @btime foo6!(A,B,C,1)
  330.808 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [December 11, 2022, 11:27pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/10 "2022-12-11T23:27:36Z")

</div>

This one can be written as reshape & then one matrix multiplication. Whether that’s faster than Tullio depends on the size (and will depend on your machine & BLAS library).

```julia
function reshape!(A,B,C)
    s1, s2 = size(B)
    mul!(reshape(C, :, s2), reshape(A, :, s1), B)
    C
end

function loop!(A,B,C)
    @views for d in axes(C,1)
        @inbounds mul!(C[d,:,:], A[d,:,:], B)
    end
    C
end

using LoopVectorization, Tullio
function tullio!(A,B,C)
    @tullio C[i, j, k] = A[i, j, p] * B[p, k]
end

using BenchmarkTools, LinearAlgebra
let n = 1 # n=10 changes the ranking?
    A = rand(4n,10n,10n)
    B = rand(10n,10n)
    C1 = zeros(4n,10n,10n)
    C2 = zeros(4n,10n,10n)
    C3 = zeros(4n,10n,10n)
    @btime loop!($A,$B,$C1)
    @btime tullio!($A,$B,$C2)
    @btime reshape!($A,$B,$C3)
    C1 ≈ C2 ≈ C3
end

```

Note that quite a few variants above which do `mul!(C[:,:,dim], ...` etc, without `view`, aren’t changing `C` at all, but just mutating a temporary copy.

> [@CeterisPartybus](#):
>
> My actual problem has a six-dimensional array. With the last two ones not being looped over. Should I change the dimensions such that I loop over dim 3-6 instead?

Can you write the formula for this?

If you must slice then slicing the last dimensions is generally best. But better to avoid it if you can.

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [December 12, 2022, 10:08am UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/11 "2022-12-12T10:08:38Z")

</div>

Awesome! Thanks! But do I understand correctly that I get zero allocations in `foo4!()` only because of switching the dimensions of `A` and `C` and taking the slice with respect to the last and not the first dimension? If so: Why?

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [December 12, 2022, 11:45am UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/12 "2022-12-12T11:45:18Z")

</div>

```julia
mul!(C[dim,1:10,1:10], A[dim,1:10,1:10], B)

```

cannot update the original matrix C. Only (@view C[…]) here can, because C[…] without @view just makes a copy of (part of) C.

---

<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:** [December 12, 2022, 8:44pm UTC](https://discourse.julialang.org/t/inplace-multiplication-of-sub-matrices-without-allocations/91534/13 "2022-12-12T20:44:36Z")

</div>

As @mcabbott and @photor have brought up, taking a slice without `view` and passing it to `map!` causes a new matrix to be allocated, so that `mul!(C[dim,:,:],...` means that the original matrix `C` isn’t even modified (so the code is functionally incorrect). Among the examples I listed, _only_ `foo4!` _and_ `foo6!` _are functional. The rest don’t even work._

You should not accept my comment as an answer to your question; I was just comparing timings. @mcabbott presents some much better solutions. On my machine, for the `n=1` case `tullio!` is fastest, and for the `n=10` case `reshape!` is fastest (even though it causes some small allocations).

> [@CeterisPartybus](#):
>
> If so: Why?

A multi-dimensional array like

```julia
[ 1 4 7
  2 5 8
  3 6 9 ]

```

is stored in memory as a sequence of numbers, and information regarding its “dimensionality” isn’t in the data but in how the code treats it. In Julia, it’s stored in a _column-major_ format, meaning that columns are contiguous like `1 2 3 4 5 6 7 8 9`. Accessing a slice like `A[:,1]` (which gives `[1,2,3]`) is more performant than accessing `A[1,:]` (which gives `[1,4,7]`), because accessing contiguous data from memory is more efficient than taking strides to access non-contiguous data.

As for why exactly a `view` of the non-contiguous data would cause so much allocation, I don’t know the mechanics to answer this.