# In place GEMM

**URL:** <https://discourse.julialang.org/t/in-place-gemm/56260>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [March 1, 2021, 3:57pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260 "2021-03-01T15:57:21Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Xijiang\_Yu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xijiang_yu/32/14439_2.png) [@Xijiang\_Yu](https://discourse.julialang.org/u/Xijiang_Yu)\
**Post date:** [March 1, 2021, 3:57pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260/1 "2021-03-01T15:57:22Z")

</div>

My LHS is huge and needs to be constructed by parts/sub-matrix. e.g.,

```julia
lhs = [ b11 b12
      b21 b22]

```

Is there a way to specify the results into a sub-matrix, e.g., b22, to save memory? I believe gemm in Julia is also based on blas, which can specify this, though painful.

---

<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:** [March 1, 2021, 4:28pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260/2 "2021-03-01T16:28:29Z")

</div>

> [@Xijiang\_Yu](#):
>
> Is there a way to specify the results into a sub-matrix, e.g., b22, to save memory?

Have you tried just creating a view of b22 (e.g. `b22 = @view lhs[i2:end, j2:end]`) and then doing an in-place operation (`mul!` or `ldiv!`) on that?

---

<div class="post-metadata">

**Author:** ![Pbellive](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pbellive/32/3604_2.png) [@Pbellive](https://discourse.julialang.org/u/Pbellive)\
**Post date:** [March 1, 2021, 5:50pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260/3 "2021-03-01T17:50:38Z")

</div>

Unless I’m missing something or misunderstanding, I don’t think this works. For example:

```julia
A = rand(8,8)
a22 = view(A, 4:8,4:8)
mul!(a22, 10*randn(5,5), a22)

```

errors. This is expected, as the docstring for `mul!` says

```julia
mul!(Y, A, B) -> Y

  Calculates the matrix-matrix or matrix-vector product AB and stores the result in Y, overwriting the existing value of Y. Note that Y must not be aliased with either A or B.

```

To do in-place, overwriting matrix multiplication you need to use `rmul!` or `lmul!`. Unfortunately those functions don’t seems to work with sub-arrays:

```julia
julia> rmul!(a22,4*randn(5,5))
ERROR: MethodError: no method matching rmul!(::SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::Array{Float64,2})
Closest candidates are:
  rmul!(::AbstractArray, ::Number) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/generic.jl:161
  rmul!(::StridedArray{T, 2}, ::LowerTriangular{T,var"#s817"} where var"#s817"<:(StridedArray{T, 2} where T)) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/triangular.jl:749
  rmul!(::StridedArray{T, 2}, ::Transpose{var"#s807",var"#s806"} where var"#s806"<:(LowerTriangular{T,var"#s805"} where var"#s805"<:(StridedArray{T, 2} where T)) where var"#s807") where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/triangular.jl:759

```

Is there a reason why this shouldn’t work for subarrays?

---

<div class="post-metadata">

**Author:** ![Xijiang\_Yu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xijiang_yu/32/14439_2.png) [@Xijiang\_Yu](https://discourse.julialang.org/u/Xijiang_Yu)\
**Post date:** [March 1, 2021, 11:19pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260/4 "2021-03-01T23:19:10Z")

</div>

```julia
t = zeros(1010, 1010)
v = rand(1000, 1000)
c = view(t, 11:1010, 11:1010)
@btime BLAS.gemm!('T', 'N', 1., v, v, 1., c) # 6.135 ms (0 allocations: 0 bytes)

@btime v'v # 5.226 ms (3 allocations: 7.63 MiB)

```

---

<div class="post-metadata">

**Author:** ![Xijiang\_Yu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xijiang_yu/32/14439_2.png) [@Xijiang\_Yu](https://discourse.julialang.org/u/Xijiang_Yu)\
**Post date:** [March 6, 2021, 10:24pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260/5 "2021-03-06T22:24:08Z")

</div>

`Ocatavian.matmul!` is even better. Which

- is faster
- allows mixed precision

---

<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:** [March 6, 2021, 11:30pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260/6 "2021-03-06T23:30:30Z")

</div>

The error message is informative here:

```julia
julia> mul!(a22, 10*randn(5,5), a22)
ERROR: ArgumentError: output matrix must not be aliased with input matrix

```

`mul!` can’t write into the same array it’s reading, as this would over-write data it still needs to read. This is true for any rectangular arrays, whether or not they are views.

> [@Pbellive](#):
>
> To do in-place, overwriting matrix multiplication you need to use `rmul!` or `lmul!` .

And these only work in cases where it is possible to work out the output and over-write the input sequentially, such as various structured matrices. They don’t mind views though:

```julia
julia> lmul!(Diagonal([0,100]), view(rand(4,4), 1:2, 3:4))
2×2 view(::Array{Float64,2}, 1:2, 3:4) with eltype Float64:
  0.0 0.0
 71.1938 90.2883

julia> lmul!(UpperTriangular(rand(2,2)), view(rand(4,4), 1:2, 3:4))
2×2 view(::Array{Float64,2}, 1:2, 3:4) with eltype Float64:
 0.0627902 0.525721
 0.0012771 0.494267

```

---

<div class="post-metadata">

**Author:** ![Pbellive](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pbellive/32/3604_2.png) [@Pbellive](https://discourse.julialang.org/u/Pbellive)\
**Post date:** [March 7, 2021, 11:26pm UTC](https://discourse.julialang.org/t/in-place-gemm/56260/7 "2021-03-07T23:26:04Z")

</div>

Thanks @mcabbott, that helped clarify things.

> [@mcabbott](#):
>
> `mul!` can’t write into the same array it’s reading, as this would over-write data it still needs to read. This is true for any rectangular arrays, whether or not they are views.

For some reason I was thinking BLAS had an overwriting version of GEMM for dense matrices. At least in a naive matrix multiplication algorithm it’s possible to do this by allocating a temporary array that can hold a single row or column of the output. I’m not sure why I thought standard BLAS implementations had such a routine. As far as I can tell they do not, and upon reflection it seems like this would only be useful in rare circumstances.