# Using axpy!

**URL:** <https://discourse.julialang.org/t/using-axpy/7072>\
**Category:** General Usage\
**Tags:** blas, benchmark\
**Created:** [November 14, 2017, 4:55pm UTC](https://discourse.julialang.org/t/using-axpy/7072 "2017-11-14T16:55:59Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 14, 2017, 4:55pm UTC](https://discourse.julialang.org/t/using-axpy/7072/1 "2017-11-14T16:55:59Z")

</div>

The operation of multiplication by a scalar and translation can be done with `Base.LinAlg.axpy!(α, x, y)`. But I can also use `y .+= α*x`.  
I was wondering: is it recommended to use one or the other for writing library code? What does `.+=` actually do? Does it call `axpy!` under the hood?

In microbenchmarks (see below), the first one allocates no memory at all, although their speed is “very very close”.

```julia
using Compat, BenchmarkTools, DataFrames

blas_1(α, x, y) = Base.LinAlg.axpy!(α, x, y)
broadcast_partial(α, x, y) = y .+= α*x
broadcast_full(α, x, y) = @. y += α*x

n = [10^k for k in 1:6]
x, y = [rand(ni) for ni in n], [rand(ni) for ni in n]
z = copy(y); w = copy(y)
α = sqrt(pi)

table = DataFrame(Algorithm=["blas_1", "broadcast_partial", "broadcast_full"])

for (i, ni) in enumerate(n) 
    blas_1_bench = @benchmark blas_1($α, $x[$i], $y[$i])
    broadcast_partial_bench = @benchmark broadcast_partial($α, $x[$i], $z[$i])
    broadcast_full_bench = @benchmark broadcast_full($α, $x[$i], $w[$i])
    times_ms = 1e-6 * [time(blas_1_bench),
                       time(broadcast_partial_bench),
                       time(broadcast_full_bench)]
    table[Symbol("n = $ni (ms)")] = [signif(ti, 3) for ti in times_ms]                                         
end

```

```julia
julia> table
3×7 DataFrames.DataFrame
│ Row │ Algorithm │ n = 10 (ms) │ n = 100 (ms) │ n = 1000 (ms) │ n = 10000 (ms) │ n = 100000 (ms) │ n = 1000000 (ms) │
├─────┼─────────────────────┼─────────────┼──────────────┼───────────────┼────────────────┼─────────────────┼──────────────────┤
│ 1 │ "blas_1" │ 3.0e-5 │ 0.00142 │ 0.00153 │ 0.00516 │ 0.01729 │ 0.31935 │
│ 2 │ "broadcast_partial" │ 6.0e-5 │ 0.00011 │ 0.00078 │ 0.00964 │ 0.11502 │ 2.34353 │
│ 3 │ "broadcast_full" │ 3.0e-5 │ 4.0e-5 │ 0.00027 │ 0.00332 │ 0.05153 │ 0.75011 │

```

_Edit 1:_ added the full loop-fusion version.  
_Edit 2:_ added a benchmark function (use IJulia to see a cute table). Benchmark data on a MacBook Pro (see details below).  
_Edit 3:_ round to three significative digits using `signif(ti, 3)`.

```julia
julia> versioninfo()
Julia Version 0.6.0
Commit 903644385b (2017-06-19 13:05 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

```

---

<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:** [November 14, 2017, 5:16pm UTC](https://discourse.julialang.org/t/using-axpy/7072/2 "2017-11-14T17:16:18Z")

</div>

You want `y .+= α.*x` (note the dot) for the (mostly) non-allocating version, see e.g. [More Dots: Syntactic Loop Fusion in Julia](https://julialang.org/blog/2017/01/moredots). It’s also much faster than `α*x`. axpy! is multithreaded though, so for larger arrays you’re going to see a speed difference.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 14, 2017, 6:24pm UTC](https://discourse.julialang.org/t/using-axpy/7072/3 "2017-11-14T18:24:33Z")

</div>

Got it, thanks! Indeed, the “full dots” version is even faster; i’ve edited the example.

---

<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:** [November 14, 2017, 6:33pm UTC](https://discourse.julialang.org/t/using-axpy/7072/4 "2017-11-14T18:33:24Z")

</div>

Interesting: for me axpy! is 20% faster on 10k elements. For smaller arrays it’s even worse, and axpy is 5x faster for arrays of size 100. For larger arrays the difference levels off, until multithreading kicks off and axpy gets faster again. I guess broadcasting still has some way to go?

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [November 14, 2017, 6:33pm UTC](https://discourse.julialang.org/t/using-axpy/7072/5 "2017-11-14T18:33:57Z")

</div>

You need “$” in front of each variable in the call to the btime macro.

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 14, 2017, 6:42pm UTC](https://discourse.julialang.org/t/using-axpy/7072/6 "2017-11-14T18:42:17Z")

</div>

Like `@btime f($α, $x, $y)`? What I wrote benchmarks something else?

---

<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:** [November 14, 2017, 6:45pm UTC](https://discourse.julialang.org/t/using-axpy/7072/7 "2017-11-14T18:45:58Z")

</div>

I see, thanks! @mforets, see the readme in [GitHub - JuliaCI/BenchmarkTools.jl: A benchmarking framework for the Julia language](https://github.com/JuliaCI/BenchmarkTools.jl). With that the broadcast is slightly faster. For n=100\_000 axpy is faster because of multithreading, for n=1\_000\_000 it still burns all my cores, but actually takes exactly the same time as broadcasting. I officially give up understanding multithreading performance.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [November 14, 2017, 7:59pm UTC](https://discourse.julialang.org/t/using-axpy/7072/8 "2017-11-14T19:59:40Z")

</div>

BLAS level 1 is mostly memory bound and doesn’t benefit that much from multithreading.

@mforets, I don’t see the edited/updated code and benchmarks.

---

<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:** [November 14, 2017, 8:09pm UTC](https://discourse.julialang.org/t/using-axpy/7072/9 "2017-11-14T20:09:44Z")

</div>

I think that statement is architecture-dependent: I’ve been on machines where the memory bandwidth scaled with the number of threads, and machines where it didn’t. But yes that’s probably what’s happening here: for n=100\_000, the arrays are still in cache and multithreading is efficient, and for large sizes it’s memory bandwidth bound anyway so doesn’t speedup with multithreading. Thanks for the explanation!

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 14, 2017, 9:09pm UTC](https://discourse.julialang.org/t/using-axpy/7072/10 "2017-11-14T21:09:14Z")

</div>

There you go ✌

---

<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:** [November 15, 2017, 7:15am UTC](https://discourse.julialang.org/t/using-axpy/7072/11 "2017-11-15T07:15:46Z")

</div>

Just out of curiosity, can you increase n and see if you reproduce what I see, that the benefits of multithreading disappear when you get out of the cache?

---

<div class="post-metadata">

**Author:** ![mforets](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mforets/32/298_2.png) [@mforets](https://discourse.julialang.org/u/mforets)\
**Post date:** [November 15, 2017, 11:42am UTC](https://discourse.julialang.org/t/using-axpy/7072/12 "2017-11-15T11:42:18Z")

</div>

Hi Antoine, i could go up to 10^8, but further than that i get out of swap memory. Here is a plot in log-scale where we can see one crossover:

 ![39](https://global.discourse-cdn.com/julialang/original/3X/7/f/7f9cf419c450360a86d129f24533bf49dcb25dd4.png)

---

<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:** [November 15, 2017, 1:52pm UTC](https://discourse.julialang.org/t/using-axpy/7072/13 "2017-11-15T13:52:49Z")

</div>

In real code (not just benchmarks), be sure you look at `axpy` in the _context_ of any other operations you are doing on the arrays. If you are doing other O(n) operations too, then you’ll probably want to fuse these into a single loop with the `axpy` anyway. If you are also doing Ω(n²) or similarly expensive operations, then aggressively optimizing the O(n) (BLAS-1) stuff is probably not worth it.