# Inplace multiplication by a square matrix

**URL:** <https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702>\
**Category:** General Usage\
**Created:** [January 26, 2017, 2:06pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702 "2017-01-26T14:06:30Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![gdkrmr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdkrmr/32/2791_2.png) [@gdkrmr](https://discourse.julialang.org/u/gdkrmr)\
**Post date:** [January 26, 2017, 2:06pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/1 "2017-01-26T14:06:30Z")

</div>

I have a `d \times d` matrix `R` and a `n \times d` matrix `X`, I want to multiply `R` with `X` and store the result in `X`, is there a way in julia to to do this without making a copy of `X`?

---

<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:** [January 26, 2017, 2:45pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/2 "2017-01-26T14:45:35Z")

</div>

You can operate on a few columns at a time to avoid copying the whole matrix, e.g. to do it one column at a time:

```julia
b = similar(X, size(X,1)) # pre-allocate array for storing results of R * column
for i = 1:size(X, 2)
    X[:,i] = A_mul_B!(R, @view(X[:,i]), b)
end

```

(It would probably be more efficient to multiply a few columns at a time, because that makes better use of BLAS.)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [January 26, 2017, 3:40pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/3 "2017-01-26T15:40:20Z")

</div>

@stevengj, can this convoluted for loop be replaced by the dot syntax you discussed in your post?

> **[More Dots: Syntactic Loop Fusion in Julia](https://julialang.org/blog/2017/01/moredots/)**
>
> More Dots: Syntactic Loop Fusion in Julia | After a lengthy design process (https://github.com/JuliaLang/julia/issues/8450) and preliminary foundations in Julia 0.5 (/blog/2016-10-11-julia-0.5-highlights#vectorized\_function\_calls), Julia 0.6 includes...

---

<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:** [January 26, 2017, 3:44pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/4 "2017-01-26T15:44:39Z")

</div>

> [@juliohm](#):
>
> @stevengj, can this convoluted for loop be replaced by the dot syntax you discussed in your post?

Not without doing something even more convoluted, as far as I can tell.

You could just do

```julia
@views for i = 1:size(X, 2)
    X[:,i] = R * X[:,i]
end

```

in 0.6 with the new `@views` macro, but you need to call `A_mul_B!` if you want to do the `matrix*vector` with a preallocated output vector.

---

<div class="post-metadata">

**Author:** ![wsshin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wsshin/32/360_2.png) [@wsshin](https://discourse.julialang.org/u/wsshin)\
**Post date:** [February 15, 2017, 11:57pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/5 "2017-02-15T23:57:14Z")

</div>

Related to this, I tried a few things and got some questions about the benchmark result.

Below, I have two functions that apply a 2-by-2 matrix `A` to a 2-by-1000 matrix whose each column is a point on the unit circle.

```julia
julia> VERSION
v"0.5.1-pre+31"

julia> function applyA(A)
           θ = linspace(0, 2π, 1001)[1:end-1] # 1000 angles
           X = [cos(θ) sin(θ)]'
           Z = A*X

           return Z
       end
applyA (generic function with 1 method)

julia> function applyAinplace(A)
           n = 1000
           ∆θ = 2π/n
           Z = Array{Float64}(2,n)
           x = Vector{Float64}(2)
           y = Vector{Float64}(2)
           for k = 1:n # 1000 angles
               x[1] = cos(k*∆θ)
               x[2] = sin(k*∆θ)
               BLAS.gemv!(y, 'N', A, x)
               view(Z,:,k) = y
           end

           return Z
       end
applyAinplace (generic function with 1 method)

```

Here is the benchmark result.

```julia
julia> @benchmark applyA($A)
BenchmarkTools.Trial:
  memory estimate: 63.34 KiB
  allocs estimate: 13
  mean time: 57.686 μs (0.78% GC)

julia> @benchmark applyAinplace($A)
BenchmarkTools.Trial:
  memory estimate: 15.94 KiB
  allocs estimate: 3
  mean time: 97.919 μs (1.53% GC)

```

So, even though the in-place version uses less allocations, it is slower. Is this because the in-place version uses level-2 BLAS whereas the other uses level-3 BLAS?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 16, 2017, 12:43am UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/6 "2017-02-16T00:43:27Z")

</div>

> [@stevengj](#):
>
> in 0.6 with the new @views macro, but you need to call A\_mul\_B! if you want to do the matrix\*vector with a preallocated output vector.

Or use this package to make things a little easier:

> **[GitHub - simonbyrne/InplaceOps.jl: Convenient macros for in-place matrix...](https://github.com/simonbyrne/InplaceOps.jl)**
>
> Convenient macros for in-place matrix operations in Julia - GitHub - simonbyrne/InplaceOps.jl: Convenient macros for in-place matrix operations in Julia

---

<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 16, 2017, 1:11pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/7 "2017-02-16T13:11:41Z")

</div>

> [@wsshin](#):
>
> apply a 2-by-2 matrix A to a 2-by-1000 matrix

Just writing out the 2x2 multiply (unrolling the loop for multiplying each column of the 2x1000 matrix — it’s only 16 operations) and sticking `@simd` in front will surely be faster than calling out to BLAS for such a small matrix? And using the StaticArrays package will unroll the loops for you.

Note also that `view(Z,:,k) = y` is unnecessary. `Z[:,k] = y` will write in-place into `Z`. (Slicing on the left-hand side of an assignment operator calls `setindex!` on the slice.) But with only two components it is probably faster to unroll the loop and assign them individually.

i.e. just do

```julia
@simd for k = 1:n
    θ = k*∆θ
    y = A * SVector(cos(θ), sin(θ))
    Z[1,k] = y[1]
    Z[2,k] = y[2]
end

```

where `A` is a 2x2 `SMatrix` (from StaticArrays).

---

<div class="post-metadata">

**Author:** ![wsshin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wsshin/32/360_2.png) [@wsshin](https://discourse.julialang.org/u/wsshin)\
**Post date:** [February 16, 2017, 10:30pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/8 "2017-02-16T22:30:20Z")

</div>

> [@stevengj](#):
>
> Just writing out the 2x2 multiply (unrolling the loop for multiplying each column of the 2x1000 matrix — it’s only 16 operations) and sticking `@simd` in front will surely be faster than calling out to BLAS for such a small matrix? And using the StaticArrays package will unroll the loops for you.
> 
> Note also that `view(Z,:,k) = y` is unnecessary. `Z[:,k] = y` will write in-place into `Z`. (Slicing on the left-hand side of an assignment operator calls `setindex!` on the slice.) But with only two components it is probably faster to unroll the loop and assign them individually.

Didn’t know about `@simd`. Thanks!

A few comments.

1. I put `view(Z,:,k)` because it made a _huge_ difference in the number of allocations. Compare the `allocs estimate` below (≈ 500) with the one shown above (= 3):

```julia
function applyAinplace2(A)
    n = 1000
    ∆θ = 2π/n
    Z = Array{Float64}(2,n)
    x = Vector{Float64}(2)
    y = Vector{Float64}(2)
    for k = 1:n # 1000 angles
        x[1] = cos(k*∆θ)
        x[2] = sin(k*∆θ)
        BLAS.gemv!(y, 'N', A, x)
        Z[:,k] = y
    end

    return Z
end

julia> @benchmark applyAinplace2($A)
BenchmarkTools.Trial:
  memory estimate: 23.58 KiB
  allocs estimate: 492
  mean time: 115.919 μs (2.25% GC)

```

Similar trend in Julia v0.6:

```julia
julia> VERSION
v"0.6.0-dev.2673"

julia> @benchmark applyAinplace(A)
BenchmarkTools.Trial:
  memory estimate: 15.94 KiB
  allocs estimate: 3
  median time: 90.076 μs (0.00% GC)

julia> @benchmark applyAinplace2(A)
BenchmarkTools.Trial:
  memory estimate: 39.20 KiB
  allocs estimate: 1492
  mean time: 140.926 μs (1.55% GC)

```

Should I file an issue on this?

1. `StaticArrays` helps a lot, but `@simd` doesn’t seem to be giving an extra speedup. Am I doing something wrong? I built Julia with `make -j 8`, though not sure if that’s relevant.

```julia
julia> function applyAstatic(A)
           n = 1000
           ∆θ = 2π/n
           Z = Array{Float64}(2,n)
           for k = 1:n
               θ = k*∆θ
               y = A * SVector(cos(θ), sin(θ))
               Z[1,k] = y[1]
               Z[2,k] = y[2]
           end

           return Z
       end
applyAstatic (generic function with 1 method)

julia> function applyAstaticSIMD(A)
           n = 1000
           ∆θ = 2π/n
           Z = Array{Float64}(2,n)
           @simd for k = 1:n
               θ = k*∆θ
               y = A * SVector(cos(θ), sin(θ))
               Z[1,k] = y[1]
               Z[2,k] = y[2]
           end

           return Z
       end
applyAstaticSIMD (generic function with 1 method)

julia> A = rand(2,2); SA = SArray{size(A)}(A);

julia> @benchmark applyAstatic($SA)
BenchmarkTools.Trial:
  memory estimate: 15.75 KiB
  allocs estimate: 1
  mean time: 31.877 μs (0.00% GC)

julia> @benchmark applyAstaticSIMD($SA)
BenchmarkTools.Trial:
  memory estimate: 15.75 KiB
  allocs estimate: 1
  mean time: 30.492 μs (0.00% GC)

```

---

<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 16, 2017, 11:19pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/9 "2017-02-16T23:19:18Z")

</div>

> [@wsshin](#):
>
> StaticArrays helps a lot, but @simd doesn’t seem to be giving an extra speedup

are you passing an SMatrix for A?

---

<div class="post-metadata">

**Author:** ![wsshin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wsshin/32/360_2.png) [@wsshin](https://discourse.julialang.org/u/wsshin)\
**Post date:** [February 16, 2017, 11:31pm UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/10 "2017-02-16T23:31:45Z")

</div>

> [@stevengj](#):
>
> are you passing an SMatrix for A?

Yes. In the code quoted in my previous posting, you can see that I created a static array version of `A` as follows:

```julia
julia> A = rand(2,2); SA = SArray{size(A)}(A);

```

I passed `SA` to the functions.

Technically, `SA` is `SArray` and not `SMatrix`, but I guess `@simd` should work for both cases? I also tried `SA = @SMatrix rand(2,2)`, and the result was no different.

---

<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 17, 2017, 2:15am UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/11 "2017-02-17T02:15:57Z")

</div>

`@simd` just turns on LLVM auto-vectorization optimizations, but compilers are pretty limited in what they can vectorize.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [February 17, 2017, 3:26am UTC](https://discourse.julialang.org/t/inplace-multiplication-by-a-square-matrix/1702/12 "2017-02-17T03:26:59Z")

</div>

Most of the time in your static versions is spent computing sines and cosines, not in matrix products. If you need to compute these every time, try a vector math library (e.g. Yeppp or MKL). SIMD conversions of integers to floats are also not always available - you might need to build a native system image to get that.  
If you want this really fast, precompute the sine/cosine arrays, rearrange to make `k` the first index of `Z`, and go back to ordinary matrices - BLAS libraries have well-tuned code for this case.