# When shouldn't we use @simd?

**URL:** <https://discourse.julialang.org/t/when-shouldnt-we-use-simd/18276>\
**Category:** General Usage\
**Created:** [December 4, 2018, 11:28am UTC](https://discourse.julialang.org/t/when-shouldnt-we-use-simd/18276 "2018-12-04T11:28:35Z")\
**Posts on this page:** 1\
**Showing post:** 11

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [December 5, 2018, 4:33pm UTC](https://discourse.julialang.org/t/when-shouldnt-we-use-simd/18276/11 "2018-12-05T16:33:46Z")

</div>

You should not use `@simd` if your algorithm depends on how `+` and `*` and other floating-point operations are associated. Consider the difference between left-to-right summation and left-to-right summation with a `@simd` annotation:

```julia
function mysum(v::Vector{Float64})
    t = 0.0
    for x in v
        t += x
    end
    return t
end

function mysum_simd(v::Vector{Float64})
    t = 0.0
    @simd for x in v
        t += x
    end
    return t
end

```

You already get slightly different results for most random arrays:

```julia
julia> v = rand(10^6);

julia> mysum(v)
499897.00784308364

julia> mysum_simd(v)
499897.00784311077

```

Well, that doesn’t seem terrible. How bad can it get? I addressed this question in this post:

> [@Array ordering and naive summation](https://discourse.julialang.org/t/array-ordering-and-naive-summation/1929):
>
> One of the first things you learn in numerical analysis is that floating-point operations are not associative. A classic example is this: julia\> (0.1 + 0.2) + 0.3 0.6000000000000001 julia\> 0.1 + (0.2 + 0.3) 0.6 I was thinking of ways to make floating-point summation independent of the order of the summands without making the performance much worse (this is a hobby of mine). Julia currently uses a pairwise summation algorithm, which is much better than naive left-to-right reduction while havin…

In short, it can be as bad as possible—you can literally get arbitrarily wrong values. Here’s a minor tweak of the `sumsto` function from that post (changing `realmax` to `floatmax` and reversing the order of the iteration of powers), which really shows how bad this is:

```julia
function sumsto(x::Float64)
    0 <= x < exp2(970) || throw(ArgumentError("sum must be in [0,2^970)"))
    n, p₀ = Base.decompose(x) # integers such that `n*exp2(p₀) == x`
    [floatmax(); [exp2(p) for p in reverse(-1074:969) if iseven(n >> (p-p₀))]
    -floatmax(); [exp2(p) for p in reverse(-1074:969) if isodd(n >> (p-p₀))]]
end

```

```julia
julia> s = rand()*10^6
629436.8129391546

julia> v = sumsto(s)
2046-element Array{Float64,1}:
   1.7976931348623157e308
   4.9896007738368e291
   2.4948003869184e291
   ⋮
   1.862645149230957e-9
   4.656612873077393e-10
   1.1641532182693481e-10

julia> mysum(v)
629436.8129391546

julia> mysum_simd(v)
6.652801031782399e291

julia> abs(mysum_simd(v)-mysum(v))/mysum(v)
1.056945017358796e286

```

Assuming that the result of left-to-right summation is the desired result, that’s a relative error of 10^286! Oops. Julia’s built-in `sum` also gives an entirely different answer by reassociating the `+` operations in a different way:

```julia
julia> sum(v)
0.0

```

Of course, there’s a strong chance that what you really wanted was the _true_ sum, which is 9.9792015476736e291 (rounded to 64-bit precision). None of these algorithms give that, although `mysum_simd` comes closest. Kahan summation does give the correct answer:

```julia
julia> using KahanSummation

julia> sum_kbn(v)
9.9792015476736e291

```

That comes at the cost of being about 2x slower than built-in `sum` or `mysum_simd`:

```julia
julia> @time mysum(v)
  0.000008 seconds (5 allocations: 176 bytes)
629436.8129391546

julia> @time mysum_simd(v)
  0.000005 seconds (5 allocations: 176 bytes)
6.652801031782399e291

julia> @time sum(v)
  0.000005 seconds (5 allocations: 176 bytes)
0.0

julia> @time sum_kbn(v)
  0.000007 seconds (5 allocations: 176 bytes)
9.9792015476736e291

```

Note, however, that Kahan summation is fully correct and still faster than non-SIMD naive summation, so it’s pretty impressive overall.

---

_[View the full topic](https://discourse.julialang.org/t/when-shouldnt-we-use-simd/18276)._