# Inplace axpy! but storing to a third arguement rather than y

**URL:** <https://discourse.julialang.org/t/inplace-axpy-but-storing-to-a-third-arguement-rather-than-y/107188>\
**Category:** Performance\
**Tags:** blas\
**Created:** [December 5, 2023, 11:44pm UTC](https://discourse.julialang.org/t/inplace-axpy-but-storing-to-a-third-arguement-rather-than-y/107188 "2023-12-05T23:44:38Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [December 5, 2023, 11:44pm UTC](https://discourse.julialang.org/t/inplace-axpy-but-storing-to-a-third-arguement-rather-than-y/107188/1 "2023-12-05T23:44:38Z")

</div>

Hello,

from the BLAS package, I know that `axpy!(a,X,Y)` computes Y += a X, where a is a scalar and X and Y are a matrix or a vector.

Now I want instead something like

Z = a X + Y

So, something like `axpyz!(a,X,Y,Z)`.

Is there a BLAS funcion or a Julia function that perform this? If not, how to implement it with the maximum efficiency?

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [December 5, 2023, 11:50pm UTC](https://discourse.julialang.org/t/inplace-axpy-but-storing-to-a-third-arguement-rather-than-y/107188/2 "2023-12-05T23:50:28Z")

</div>

A common pattern for out-of-place versions of otherwise in-place functions is to copy the argument you don’t want overwritten to the output, like `axpy!(a,X,copy!(Z,Y))`.

That said, I imagine that your fastest option here is simply to use broadcasting like `Z .= a .* X .+ Y`. Even slightly better is to allow FMA to be used (for slightly higher speed and accuracy) via `Z .= muladd.(a, X, Y)` (which is almost-certainly what `axpy!` uses). These are better because they don’t have the extra step of `copy!(Z,Y)` like we had to use to make `axpy!` work.

For simple element-wise operations like these, BLAS is unlikely to do notably better than simple broadcast. The broadcast has the additional advantage of keeping everything in Julia, which means that additional types of elements and arrays will “just work” here. BLAS is only compiled for a few common types.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [December 6, 2023, 12:15am UTC](https://discourse.julialang.org/t/inplace-axpy-but-storing-to-a-third-arguement-rather-than-y/107188/3 "2023-12-06T00:15:33Z")

</div>

I tried the following example to test it with `axpy!`

```julia
NN = 400
A = rand(Float64, NN, NN)
B = rand(Float64, NN, NN);
B_cache = deepcopy(B);
B_cache2 = deepcopy(B);

function axpy2!(a,X,Y)
    @. Y += a*X
    Y
end

```

```julia
@benchmark axpy!(0.5, $A, $B_cache)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 62.600 μs … 344.000 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 70.600 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 73.767 μs ± 10.851 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

    ▃ █▅▇▇▄▇▅▄▆▄▄▄▂▁▂▂▁▂▁ ▂ ▂
  █▁█▇████████████████████▇██▇▇██▇▆▇██▇▆▇█▆▆▅▃▃█▆▃▆▄▄▃█▄▁▁▃▄▁▇ █
  62.6 μs Histogram: log(frequency) by time 126 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

```julia
@benchmark axpy2!(0.5, $A, $B_cache2)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 116.100 μs … 541.000 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 132.300 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 139.289 μs ± 29.856 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

     ▆ █ █▄                                                      
  ▅▅▅█▇█▇██▄▃▃▂▂▃▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  116 μs Histogram: frequency by time 267 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

The “home made” version is slightly slower.

Moreover, if I choose random sparse matrices, I have unexpectedly 500 times faster calcs, but I have also some allocations, whose size increase with the dimension of the matrix

```julia
NN = 400
A = sprand(Float64, NN, NN, 0.1/NN)
B = sprand(Float64, NN, NN, 0.1/NN);
B_cache = deepcopy(B);
B_cache2 = deepcopy(B);

```

```julia
@benchmark axpy!(0.5, $A, $B_cache)

BenchmarkTools.Trial: 2252 samples with 1 evaluation.
 Range (min … max): 1.611 ms … 5.790 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 2.186 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 2.206 ms ± 375.599 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

                ▆█                                             
  ▁▃▂▂▂▄▃▆▅▅▅▄▅▄███▆▄▃▃▂▂▂▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  1.61 ms Histogram: frequency by time 3.89 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

```julia
@benchmark axpy2!(0.5, $A, $B_cache2)

BenchmarkTools.Trial: 10000 samples with 8 evaluations.
 Range (min … max): 3.788 μs … 998.751 μs ┊ GC (min … max): 0.00% … 79.00%
 Time (median): 4.162 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 4.876 μs ± 18.765 μs ┊ GC (mean ± σ): 6.49% ± 1.69%

  ▁▅▇██▇▆▅▄▅▄▄▄▃▂▂▁▁▁▁▁ ▂
  ██████████████████████▇█▇▇██▇▇▇▇█▇▇███▇▇███▇▇▇▇███▆▆▆▇▆▆▆▅▄ █
  3.79 μs Histogram: log(frequency) by time 8.12 μs <

 Memory estimate: 4.75 KiB, allocs estimate: 3.

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 6, 2023, 7:57am UTC](https://discourse.julialang.org/t/inplace-axpy-but-storing-to-a-third-arguement-rather-than-y/107188/4 "2023-12-06T07:57:03Z")

</div>

you can use LoopVectorization to speed this up further

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 6, 2023, 10:49am UTC](https://discourse.julialang.org/t/inplace-axpy-but-storing-to-a-third-arguement-rather-than-y/107188/5 "2023-12-06T10:49:32Z")

</div>

In my tests with @albertomercurio 's first example, I get the exact same timings for pretty much everything I try (axpy, manual loop, broadcast, `@turbo`), which is great - it means the base julia is already optimal. I do get better timings (in fact, a superlinear speedup!) when BLAS uses threading: check BLAS.set\_num\_threads. I would have thought this would be too small to benefit from threading, but apparently not…