# How to do in-place matrix operation in Julia?

**URL:** <https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [February 14, 2023, 2:57am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601 "2023-02-14T02:57:03Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![J\_I](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_i/32/46880_2.png) [@J\_I](https://discourse.julialang.org/u/J_I)\
**Post date:** [February 14, 2023, 2:57am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/1 "2023-02-14T02:57:03Z")

</div>

Hi, I’m new to Julia, and this is the first post to this forum.

I have problems regarding in-place matrix operations in julia.  
I have to calculate matrix-matrix product in for-loop, and I tried to pre-allocate the memories, before the for-loop, as follows::

```julia
function test_1(N)
    a = rand(Float64, N, N)
    b = rand(Float64, N, N)
    c = Matrix{Float64}(undef, N, N)
    for i in 1:500
        mul!(c, a,b) ## inplace calculation of matrix matrix product
    end
end

```

I confirmed that this mitigates the memory allocation and increases efficiency, compared to the case where allocating memory in each step like below

```julia
function test_2(N)
    a = rand(Float64, N, N)
    b = rand(Float64, N, N)
    for i in 1:500
        c = a*b
    end
end

```

This was an easy task, because of the very simple matrix operation.  
I want to generalize this into more complex matrix operations, like C = A_B - B_A.  
How can I do this?

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [February 14, 2023, 3:02am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/2 "2023-02-14T03:02:43Z")

</div>

You might like [Tullio.jl](https://github.com/mcabbott/Tullio.jl).

---

<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 14, 2023, 3:12am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/3 "2023-02-14T03:12:58Z")

</div>

> [@J\_I](#):
>
> I want to generalize this into more complex matrix operations, like C = A_B - B_A.

Do

```julia
C .= mul!(X1, A, B) .- mul!(X2, B, A)

```

where `X1` and `X2` are additional temporary arrays.

---

<div class="post-metadata">

**Author:** ![J\_I](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_i/32/46880_2.png) [@J\_I](https://discourse.julialang.org/u/J_I)\
**Post date:** [February 14, 2023, 5:10am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/4 "2023-02-14T05:10:23Z")

</div>

Thanks!. I tried Tullio.jl, and it indeed reduces the cost of memory allocation, but the matrix multiplication in this method was much slower than the conventional operation.  
Anyway, I didn’t know this tool, and it might be helpful in other processes.  
Thank you so much!

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [February 14, 2023, 5:17am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/5 "2023-02-14T05:17:01Z")

</div>

> [@stevengj](#):
>
> `C .= mul!(X1, A, B) .- mul!(X2, B, A)`

Or, with even fewer temporary arrays,

```julia
mul!(C, B, A)
mul!(C, A, B, 1.0, -1.0)

```

---

<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 14, 2023, 5:22am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/6 "2023-02-14T05:22:41Z")

</div>

If you want `Tullio` to be faster you should add `LoopVectorization`.

---

<div class="post-metadata">

**Author:** ![J\_I](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_i/32/46880_2.png) [@J\_I](https://discourse.julialang.org/u/J_I)\
**Post date:** [February 14, 2023, 5:25am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/7 "2023-02-14T05:25:06Z")

</div>

Thank you so much. I understand your suggestion. But in that case, we have to prepare several temporary arrays to calculate more complex matrix operations.  
Is there any way to do that more easily and generally while keeping efficiency? (operation AB - BA was just an example, and I have to calculate more complex operations)

---

<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 14, 2023, 5:35am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/8 "2023-02-14T05:35:04Z")

</div>

In general, if you’re doing large matrix multiplication, the allocations won’t matter (since they’re n^2 vs n^3). If you’re doing small matrix multiplication, you might want to check out StaticArrays which will be much faster for smaller than 10x10 or so (also `Vector{SVector}` can be really useful for something like 15000x3.

---

<div class="post-metadata">

**Author:** ![J\_I](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_i/32/46880_2.png) [@J\_I](https://discourse.julialang.org/u/J_I)\
**Post date:** [February 14, 2023, 6:01am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/9 "2023-02-14T06:01:20Z")

</div>

> [@Oscar\_Smith](#):
>
> In general, if you’re doing large matrix multiplication, the allocations won’t matter (since they’re n^2 vs n^3).

yeah, that’s true. I should have given you more information. Basically, I have to calculate the time evolution of the Vector{Matrix} and, in each time step there is matrix multiplication, like below.

```julia
A = Vector{Matrix{ComplexF64}}
B = Vector{Matrix{ComplexF64}}

for t in time_step
   for i in eachindex(A)
        C[i] = A[i] * B[i]
   end
end

```

and the size of each element(Matrix) A is small (~4✕4), while the length of A is about 1000.  
So I thought the cost of allocation might be comparable to that of matrix multiplication itself.

And thank you for suggesting me StaticArrays. I’ve never used it, so let me check it.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [February 14, 2023, 6:12am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/10 "2023-02-14T06:12:55Z")

</div>

Note that static arrays are most useful when the matrix sizes are known at compile time. If this is your use case, the package seems like a perfect fit.

---

<div class="post-metadata">

**Author:** ![J\_I](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_i/32/46880_2.png) [@J\_I](https://discourse.julialang.org/u/J_I)\
**Post date:** [February 14, 2023, 7:43am UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/11 "2023-02-14T07:43:06Z")

</div>

Thank you very much. I confirmed that the matrix-matrix multiplication in SMatrix is much faster than the conventional one. However, my code is rather complicated and a little bit tricky to work with StaticArrays. Let me think about how to incorporate StaticArrays into my code. Thanks again!

---

<div class="post-metadata">

**Author:** ![davidavdav](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/davidavdav/32/1065_2.png) [@davidavdav](https://discourse.julialang.org/u/davidavdav)\
**Post date:** [February 18, 2023, 12:13pm UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/12 "2023-02-18T12:13:22Z")

</div>

A while ago we worked on a package [InplaceLinalg](https://github.com/davidavdav/InplaceLinalg.jl). It lets you write things like

```julia
@inplace C += 2π * A * B'

```

and this will be converted to the correct call `gemm!()`.

I can’t recall if we ever published this package, but you can always install from github.

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [February 18, 2023, 12:49pm UTC](https://discourse.julialang.org/t/how-to-do-in-place-matrix-operation-in-julia/94601/13 "2023-02-18T12:49:29Z")

</div>

There’s also [GitHub - lopezm94/SugarBLAS.jl: Syntactic sugar for BLAS polynomials](https://github.com/lopezm94/SugarBLAS.jl) that does something similar