# Unroll and Vectorize EvalPoly

**URL:** <https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835>\
**Category:** General Usage\
**Created:** [May 21, 2017, 3:09pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835 "2017-05-21T15:09:22Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![goldfita](https://avatars.discourse-cdn.com/v4/letter/g/4491bb/32.png) [@goldfita](https://discourse.julialang.org/u/goldfita)\
**Post date:** [May 21, 2017, 3:09pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/1 "2017-05-21T15:09:22Z")

</div>

I often encounter problems like the following where I want to vectorize a short numerical algorithm but keep the code as high level as possible so I can change the precision, filter length, or other parameters. This function returns a function that takes the dot product with a vector where each element is the evaluation of a polynomial. The loop and polynomial is fully unrolled so there are no unknown lengths when the final function is compiled.

```julia
#https://stackoverflow.com/questions/28077057/julia-evalpoly-macro-with-varargs
function get_farrow_dot(fcoeff::Array{Float64,2})
    fc = fcoeff[:,end:-1:1]
    flen = size(fc,1)
    ex = :0
    for i = 1:flen
        ex = :(muladd(s[$flen-$i+1],@evalpoly(α,$(fc[i,:]...)),$ex))
    end
    eval(:(function (α::Float64,s::Vector{Float64}) @fastmath @inbounds return $ex end))
end

```

I couldn’t get Julia/LLVM to vectorize. There are a few potential issues. If the loop is too small, it won’t be vectorized. Also, the polynomial coefficients are expanded into a vararg tuple of constants; so, I thought maybe LLVM wasn’t seeing the simple matrix structure. However, I came up with the following simplified function which does not vectorize either.

```julia
function test(α::Float64,s::Vector{Float64},fc::Array{Float64,2})
    tot=0.0
   @simd for i=1:200
         @fastmath @inbounds tot+=muladd(s[i],@evalpoly(α,fc[i,1],fc[i,2],fc[i,3],fc[i,4],fc[i,5],fc[i,6],fc[i,7],fc[i,8]),tot)
    end
    tot
end

```

A [straightforward vectorization](https://godbolt.org/g/N3vZvM) without worrying about instruction latency or cache demonstrates this can be vectorized. In fact, GCC does an excellent job without any optimizations from me.

```julia
julia> t=randn(48);
julia> fdot=get_farrow_dot(farrow_coeff);
julia> @btime fdot(-.11,$t)
  166.033 ns (1 allocation: 16 bytes)
1.3992278257126003
julia> @btime fdot_c(-.11,$t,$(farrow_coeff[end:-1:1,:]))
  51.351 ns (0 allocations: 0 bytes)
1.3992278257125998
julia> @btime fdot_copt(-.11,$t,$(farrow_coeff[end:-1:1,:]))
  26.808 ns (0 allocations: 0 bytes)
1.3992278257125998

```

At this point, I’m stuck. I’m not sure if there is a better way to express this so that it yields something like my manually optimized version, or if it’s a limitation of Julia or LLVM. It would be a big deal if I could get this to work. It takes a lot of time even to vectorize something simple like this.

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [May 21, 2017, 3:45pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/2 "2017-05-21T15:45:30Z")

</div>

LLVM vectorization doesn’t seem to handle `muladd` well. Replace it with `*` and `+` and it’ll vectorize. (Also note that in the julia version you are doubling `tot` everytime)

---

<div class="post-metadata">

**Author:** ![goldfita](https://avatars.discourse-cdn.com/v4/letter/g/4491bb/32.png) [@goldfita](https://discourse.julialang.org/u/goldfita)\
**Post date:** [May 21, 2017, 6:44pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/3 "2017-05-21T18:44:41Z")

</div>

The code is nearly identical with muladd replaced by a\*b+c in test. I believe it’s using the bottom 64 bits of the XMM registers.

```julia
        leaq dlatps64_(%r10,%rcx,8), %rdx
        vmovsd (%r10,%rcx,8), %xmm1 # xmm1 = mem[0],zero
        vfmadd213sd (%rax,%rdx), %xmm0, %xmm1
        addq %rax, %rdx
        vfmadd213sd (%rax,%rdx), %xmm0, %xmm1
        addq %rax, %rdx
        vfmadd213sd (%rax,%rdx), %xmm0, %xmm1
        addq %rax, %rdx
        vfmadd213sd (%rax,%rdx), %xmm0, %xmm1
        addq %rax, %rdx
        vfmadd213sd (%rax,%rdx), %xmm0, %xmm1
        addq %rax, %rdx
        vfmadd213sd (%rax,%rdx), %xmm0, %xmm1
        addq %rax, %rdx
        vfmadd213sd (%rax,%rdx), %xmm0, %xmm1
Source line: 4
        vmulsd (%r9,%rcx,8), %xmm1, %xmm1
        vmovq %r8, %xmm2
        vaddsd %xmm2, %xmm1, %xmm1
        vmovq %xmm1, %r8

```

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [May 21, 2017, 7:19pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/4 "2017-05-21T19:19:53Z")

</div>

Because you also need to fix

> [@yuyichao](#):
>
> Also note that in the julia version you are doubling tot everytime

i.e. you are doing `tot += muladd(..., ..., tot)`. This appears to affect the vectorization decision.

---

<div class="post-metadata">

**Author:** ![goldfita](https://avatars.discourse-cdn.com/v4/letter/g/4491bb/32.png) [@goldfita](https://discourse.julialang.org/u/goldfita)\
**Post date:** [May 21, 2017, 7:23pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/5 "2017-05-21T19:23:08Z")

</div>

I had removed the extra +. It still didn’t vectorize. I think the vector is also backwards. I hadn’t meant for test to be functionally correct – just a simplified version to experiment on.

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [May 21, 2017, 8:24pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/6 "2017-05-21T20:24:53Z")

</div>

The following works for me at least…

```julia
julia> function test(α::Float64,s::Vector{Float64},fc::Array{Float64,2})
           tot=0.0
           @simd for i=1:200
                @fastmath @inbounds tot += s[i] * @evalpoly(α,fc[i,1],fc[i,2],fc[i,3],fc[i,4],fc[i,5],fc[i,6],fc[i,7],fc[i,8])
           end
           tot
       end
test (generic function with 1 method)

julia> @code_native test(1.0, Float64[], Matrix{Float64}(0, 0))
        .text
Filename: REPL[1]
        pushq %rbp
        movq %rsp, %rbp
        movq 24(%rsi), %rax
        movq (%rdi), %rcx
Source line: 71
        vbroadcastsd %xmm0, %ymm0
        imulq $56, %rax, %rdx
        addq (%rsi), %rdx
        shlq $3, %rax
        negq %rax
        vxorpd %ymm1, %ymm1, %ymm1
        xorl %esi, %esi
        vxorpd %ymm2, %ymm2, %ymm2
        nopl (%rax,%rax)
Source line: 129
L48:
        leaq (%rdx,%rsi,8), %rdi
        vmovupd (%rdx,%rsi,8), %ymm3
        vmovupd 32(%rdx,%rsi,8), %ymm4
        vfmadd213pd (%rax,%rdi), %ymm0, %ymm3
        vfmadd213pd 32(%rax,%rdi), %ymm0, %ymm4
        leaq (%rdi,%rax), %rdi
        vfmadd213pd (%rax,%rdi), %ymm0, %ymm3
        vfmadd213pd 32(%rax,%rdi), %ymm0, %ymm4
        leaq (%rdi,%rax), %rdi
        vfmadd213pd (%rax,%rdi), %ymm0, %ymm3
        vfmadd213pd 32(%rax,%rdi), %ymm0, %ymm4
        leaq (%rdi,%rax), %rdi
        vfmadd213pd (%rax,%rdi), %ymm0, %ymm3
        vfmadd213pd 32(%rax,%rdi), %ymm0, %ymm4
        leaq (%rdi,%rax), %rdi
        vfmadd213pd (%rax,%rdi), %ymm0, %ymm3
        vfmadd213pd 32(%rax,%rdi), %ymm0, %ymm4
        leaq (%rdi,%rax), %rdi
        vfmadd213pd (%rax,%rdi), %ymm0, %ymm3
        vfmadd213pd 32(%rax,%rdi), %ymm0, %ymm4
        leaq (%rdi,%rax), %rdi
        vfmadd213pd (%rax,%rdi), %ymm0, %ymm3
        vfmadd213pd 32(%rax,%rdi), %ymm0, %ymm4
Source line: 4
        vmulpd (%rcx,%rsi,8), %ymm3, %ymm3
        vmulpd 32(%rcx,%rsi,8), %ymm4, %ymm4
        vaddpd %ymm1, %ymm3, %ymm1
        vaddpd %ymm2, %ymm4, %ymm2
Source line: 74
        addq $8, %rsi
        cmpq $200, %rsi
        jne L48
Source line: 4
        vaddpd %ymm1, %ymm2, %ymm0
        vextractf128 $1, %ymm0, %xmm1
        vaddpd %ymm1, %ymm0, %ymm0
        vhaddpd %ymm0, %ymm0, %ymm0
Source line: 6
        popq %rbp
        vzeroupper
        retq

```

---

<div class="post-metadata">

**Author:** ![goldfita](https://avatars.discourse-cdn.com/v4/letter/g/4491bb/32.png) [@goldfita](https://discourse.julialang.org/u/goldfita)\
**Post date:** [May 21, 2017, 11:27pm UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/7 "2017-05-21T23:27:02Z")

</div>

> [@yuyichao](#):
>
> julia\> @code\_native test(1.0, Float64, Matrix{Float64}(0, 0))  
> .text  
> Filename: REPL[1]  
> pushq %rbp  
> movq %rsp, %rbp  
> movq 24(%rsi), %rax  
> movq (%rdi), %rcx  
> Source line: 71  
> vbroadcastsd %xmm0, %ymm0

Maybe my LLVM is too old?

```julia
Julia Version 0.6.0-rc2.0
Commit 68e911be53 (2017-05-18 02:31 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)

```

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [May 22, 2017, 12:38am UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/8 "2017-05-22T00:38:10Z")

</div>

This is LLVM 4.0 on Linux.

```julia
Julia Version 0.7.0-DEV.176
Commit 50c2a378ba* (2017-05-15 00:13 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E3-1505M v6 @ 3.00GHz
  WORD_SIZE: 64
  BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-4.0.0 (ORCJIT, skylake)

```

---

<div class="post-metadata">

**Author:** ![cstjean](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cstjean/32/1444_2.png) [@cstjean](https://discourse.julialang.org/u/cstjean)\
**Post date:** [May 22, 2017, 3:02am UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/9 "2017-05-22T03:02:09Z")

</div>

I don’t know if it’ll help, but FYI there’s a package for fixed-length arrays with automatically-unrolled operations.

> **[GitHub - JuliaArrays/StaticArrays.jl: Statically sized arrays for Julia](https://github.com/JuliaArrays/StaticArrays.jl)**
>
> Statically sized arrays for Julia. Contribute to JuliaArrays/StaticArrays.jl development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![goldfita](https://avatars.discourse-cdn.com/v4/letter/g/4491bb/32.png) [@goldfita](https://discourse.julialang.org/u/goldfita)\
**Post date:** [May 22, 2017, 3:43am UTC](https://discourse.julialang.org/t/unroll-and-vectorize-evalpoly/3835/10 "2017-05-22T03:43:06Z")

</div>

Thanks for the tip. This will be helpful for other things I’m working on. I can wait for the LLVM version to be bumped up to 4.0 in a future release. It’s also possible to select the version as part of the build, but I’m not sure if that will lead to an unstable Julia.