# Avoiding allocations of small but non-trivial arrays (work array alternative?)

**URL:** <https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084>\
**Category:** Performance\
**Tags:** question\
**Created:** [November 11, 2022, 12:17am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084 "2022-11-11T00:17:02Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 11, 2022, 12:17am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/1 "2022-11-11T00:17:02Z")

</div>

Hi everyone,

I’m in the process of optimizing a portion of my code that has a general problem of too much memory allocation, and consequently poor threading efficiency and overall performance. This is a series of pretty simple functions that are called a gazillion times, and each of which needs to use various small temporary arrays.

I’m able to rely on Static Arrays (which, as I understand, get allocated on the stack) but only when the arrays in question can be initialized in a single line of code. However, in other situations, when I need to perform a small computation to build my array, the only solution I found was defining “work” arrays and passing them around. Now, this can be annoying sometimes, because I have several small arrays that I use.

Consider this MWE:

```julia
function barw_array( xnodes )
    n = length(xnodes)
    barw = Array{Float64}(undef,n)
    for j ∈ 1:n
        barw[j] = 1.0
        for k ∈ 1:n
            if k != j
                barw[j] = barw[j] / (xnodes[j] - xnodes[k])
            end
        end
    end
    return sum(barw)
end

function run_barw_array()
    xnodes = LinRange(0,1,4)

    for i ∈ 1:100000
        u = barw_array( xnodes )
    end
end

@time run_interp()

```

This results in too much memory allocation (that grows with the size of the loop).

```julia
 0.292006 seconds (2.63 M allocations: 109.113 MiB, 7.85% gc time, 15.11% compilation time)

```

The work-array version doesn’t have this problem

```julia
function barw_workvector( barw, xnodes )
    n = length(xnodes)
    for j ∈ 1:n
        barw[j] = 1.0
        for k ∈ 1:n
            if k != j
                barw[j] = barw[j] / (xnodes[j] - xnodes[k])
            end
        end
    end
    return sum(barw)
end

function run_barw_workvector()
    xnodes = LinRange(0,1,4)

    barw = Array{Float64}(undef,length(xnodes))
    for i ∈ 1:100000
        u = barw_workvector( barw, xnodes )
    end

end

@time run_interp_workvector()

```

```julia
0.005938 seconds (2 allocations: 192 bytes)

```

Now, I would love to be able to obtain this functionality with StaticArrays.

I’ve been reading some threads in the forum, like this one:

[Can Julia optimize mutable static arrays to be allocated on the stack?](https://discourse.julialang.org/t/can-julia-optimize-mutable-static-arrays-to-be-allocated-on-the-stack/66360)

but in the proposed solution, they use an `MArray` in conjunction with a work array. If I try to allocate an `MArray` in my original function

```julia
barw = MVector{n,Float64}

```

it seems that these vectors are allocated on the heap, and therefore I have the same poor performance as in the original code.

For a related discussion, this was also an interesting thread from 2018: [Fortran vs Julia stack-allocated arrays](https://discourse.julialang.org/t/fortran-vs-julia-stack-allocated-arrays/18119).

I also tried using StaticArrays with generator constructors that call a function, just to see if that worked

```julia
using StaticArrays

function xdiff(xnodes,i)
    xd = 1.0
    for k ∈ eachindex(xnodes)
        if k != i
            xd = xd / (xnodes[i] - xnodes[k])
        end
    end
    return xd
end

function sarray_constructor( xnodes )
    n = length(xnodes)
    barw = SVector{n,Float64}( ( xdiff(xnodes,i) for i in 1:n )... )
    s = sum( barw )
    return s
end

function run_barw_sarray_constructor()
    xnodes = LinRange(0,1,4)

    for i ∈ 1:100000
        u = barw_sarray_constructor( xnodes )
    end

end

@time run_barw_sarray_constructor()

```

This effectively doesn’t allocate. But regardless of the mess that this would be, there is also a performance penalty – benchmarking with `@btime` I got:

```julia
  7.868 ms (0 allocations: 0 bytes)
  4.956 ms (2 allocations: 192 bytes)

```

Any ideas? Maybe there is a package that enables something along the lines of stack-allocated mutable arrays, or that makes this sort of StaticArray constructors work properly?

Thanks!

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [November 11, 2022, 8:06am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/2 "2022-11-11T08:06:00Z")

</div>

Could you write the timing lines next to all the timing results? It seems obvious for the most part but I’m not sure what the 7.868ms vs 4.956ms `@btime`ings are because the preceding code block only contains 1 `@time` line. Part of my confusion is that you say there is a performance penalty, but the 4.956ms after the `SVector` code is better than the 0.005938s after the work vector code.

I’m also interested to learn the current state of stack-allocating fixed-size mutables, good topic.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 11, 2022, 8:59am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/3 "2022-11-11T08:59:37Z")

</div>

> [@martin.d.maas](#):
>
> it seems that these vectors are allocated on the heap,

Probably because `n` is not known at compile time by the function. One alternative is passing it as a Val parameter, then the MArray version will likely be non allocating as well.

edit:, this:

```julia
julia> function barw_array( xnodes, ::Val{n} ) where n
           barw = zeros(MVector{n,Float64})
           for j ∈ 1:n
               barw[j] = 1.0
               for k ∈ 1:n
                   if k != j
                       barw[j] = barw[j] / (xnodes[j] - xnodes[k])
                   end
               end
           end
           return sum(barw)
       end
barw_array (generic function with 2 methods)

julia> function run_barw_array()
           xnodes = LinRange(0,1,4)
           n = Val(length(xnodes))
           for i ∈ 1:100000
               u = barw_array( xnodes, n)
           end
       end
run_barw_array (generic function with 1 method)

julia> @btime run_barw_array()
  2.395 ms (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 11, 2022, 9:01am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/4 "2022-11-11T09:01:58Z")

</div>

> [@martin.d.maas](#):
>
> ```julia
> barw = SVector{n,Float64}( ( xdiff(xnodes,i) for i in 1:n )... )
> s = sum( barw )
> 
> ```

Maybe off topic, but you can sum directly over the generator, no need for the intermediate array.

---

<div class="post-metadata">

**Author:** ![jmair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmair/32/35117_2.png) [@jmair](https://discourse.julialang.org/u/jmair)\
**Post date:** [November 11, 2022, 9:23am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/5 "2022-11-11T09:23:01Z")

</div>

Is it possible to allocate a single standard array for all the computations and reuse the memory? If you make your functions mutating which reuse the memory, then you can reduce the allocations to 0 inside the hot loop and just allocate once beforehand.

Even if the array needs to change size, if you know the largest size beforehand, then allocate to the largest and use a view to get a reference to a smaller array using the same memory.

If StaticArrays is really what you want, you can do `n = Val(length(xnodes))` as suggested by @lmiq and pass this into a function (so it can know how much memory to allocate on the stack at compile time).

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 11, 2022, 10:28pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/6 "2022-11-11T22:28:01Z")

</div>

> [@lmiq](#):
>
> Probably because `n` is not known at compile time by the function. One alternative is passing it as a Val parameter, then the MArray version will likely be non-allocating as well.

That’s great, this works wonderfully!

I didn’t know about Val parameters, so I will look into them. Here’s the section of the manual for anyone who might be interested: [Value types](https://docs.julialang.org/en/v1/manual/types/#%22Value-types%22)

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 11, 2022, 10:29pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/7 "2022-11-11T22:29:23Z")

</div>

> [@jmair](#):
>
> Is it possible to allocate a single standard array for all the computations and reuse the memory?

Yes, I think this is what I called a “work” array. There are many problems with this approach, for example, in the context of multi-threaded computations these work arrays can give rise to race conditions.

> [@Benny](#):
>
> I’m also interested to learn the current state of stack-allocating fixed-size mutables, good topic.

MArrays and Value-types seem the way to go!

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 11, 2022, 11:14pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/8 "2022-11-11T23:14:37Z")

</div>

> [@martin.d.maas](#):
>
> ```julia
> function barw_array( xnodes )
> n = length(xnodes)
> barw = Array{Float64}(undef,n)
> for j ∈ 1:n
> barw[j] = 1.0
> for k ∈ 1:n
> if k != j
> barw[j] = barw[j] / (xnodes[j] - xnodes[k])
> end
> end
> end
> return sum(barw)
> end
> 
> ```

The array is indeed you problem here, but, at least in this MWE, it looks like it is not needed:

```julia
using LoopVectorization
function barw_array( xnodes )
    n = length(xnodes)
    sum_barw = 0.0
    @turbo for j ∈ 1:n
        barw_j = 1.0
        for k ∈ 1:n
            if k != j
                barw_j *= (xnodes[j] - xnodes[k])
            end
        end
        sum_barw += 1/barw_j
    end
    return sum_barw
end

```

This should be much much faster and should not allocate. I also divided the number of divisions by n, which should improve perf a lot too.

Tell me what you think of the perf

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 11, 2022, 11:33pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/9 "2022-11-11T23:33:32Z")

</div>

Oh, thanks, that looks nice!

The MWE was indeed a simplification of what I need to do (I don’t return the sum of barw, but use the array in conjunction with other data), but it is true that I could try to do something similar to what you’ve done in the real case.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [November 11, 2022, 11:35pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/10 "2022-11-11T23:35:34Z")

</div>

There’s a brand new package that seems like made for you. It’s probably not yet registered (you can still use/try it), nor announced:

> **[GitHub - MasonProtter/Bumper.jl: Bring Your Own Stack](https://github.com/MasonProtter/Bumper.jl)**
>
> Bring Your Own Stack. Contribute to MasonProtter/Bumper.jl development by creating an account on GitHub.

In theory Julia could heap-allocate for you and free before exit of barw\_array, but it would be slower. Or allocate on the stack, but it’s limited so dangerous, you don’t know how large the array is going to be, so I doubt that’s a better policy (or check at runtime the size and only stack-allocate if small).

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 11, 2022, 11:43pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/11 "2022-11-11T23:43:37Z")

</div>

> [@martin.d.maas](#):
>
> ```julia
> function sarray_constructor( xnodes )
> n = length(xnodes)
> barw = SVector{n,Float64}( ( xdiff(xnodes,i) for i in 1:n )... )
> s = sum( barw )
> return s
> end
> 
> ```

There’s no need to pass the length of `xnodes` as a `Val`, just make the input `xnode` a static vector, then the compiler will know `n` statically.

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 11, 2022, 11:55pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/12 "2022-11-11T23:55:06Z")

</div>

> [@lrnv](#):
>
> ```julia
> using LoopVectorization
> function barw_array( xnodes )
> n = length(xnodes)
> sum_barw = 0.0
> @turbo for j ∈ 1:n
> barw_j = 1.0
> for k ∈ 1:n
> if k != j
> barw_j *= (xnodes[j] - xnodes[k])
> end
> end
> sum_barw += 1/barw_j
> end
> return sum_barw
> end
> 
> ```

Were you able to compile this? I get the following error:

```julia
LoadError: Don't know how to handle expression.
if k != j

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 11, 2022, 11:58pm UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/13 "2022-11-11T23:58:12Z")

</div>

Loopvectorization doesn’t handle all Julia code. Try splitting the loop into `1:j-1` and `j+1:n`, perhaps.

Even that might not work, then I’d try moving `@turbo` to the inner loop.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [November 12, 2022, 12:20am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/14 "2022-11-12T00:20:00Z")

</div>

> [@martin.d.maas](#):
>
> in the context of multi-threaded computations these work arrays can give rise to race conditions.

No, unless they are global (only such an array or one you allocate and pass around could have such problems).

The package (@MasonProtter’s) I pointed you to in my last comment does have a global array, and would have problems with threads, except as documented, it also has an array for each thread. Read its docs.

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 12, 2022, 12:41am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/15 "2022-11-12T00:41:12Z")

</div>

Sure, I will take a look at that package (however, I’m always reluctant to add more dependencies).

Having one array for each thread sounds like a sensible choice.

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 12, 2022, 12:46am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/16 "2022-11-12T00:46:12Z")

</div>

As a quick update, here’s actually my complete test (which is not much larger than the original one). It’s the barycentric approximation formula, which I’m hoping to specialize for a fixed degree and arbitrary nodes.

I tested 3 versions, the first one without any array allocations (using `barw_j` instead of `barw[j]`), the second one using MVectors and Val-types, and the third one hard-coding the case `n=4` for a standard StaticArray.

Long story short, MVectors are faster than the version which doesn’t use an array at all (this is surprising), and comes close to the performance of the hard-coded version. In all cases, I’m using multiplications and a single division instead of multiple divisions, which proved faster indeed.

```julia
using LoopVectorization
using StaticArrays
using BenchmarkTools

function Interp1d_no_array( xnodes, fvals, t )
    n = length(xnodes)
    numt, denomt = 0, 0

    @inbounds for j ∈ eachindex(xnodes)
        barw_j = 1.0
        for k ∈ 1:n
            if k != j
            barw_j *= (xnodes[j] - xnodes[k])
            end
        end
        barw_j = 1.0 / barw_j
        tdiff = t - xnodes[j]
        numt = numt + barw_j / tdiff * fvals[j]
        denomt = denomt + barw_j / tdiff
        if ( abs(tdiff) < 1e-14 )
            numt = fvals[j]
            denomt = 1.0
            break
        end
    end
    return numt / denomt
end

function Interp1d_MVector( xnodes, fvals, t, ::Val{n} ) where n 
    
    barw = zeros(MVector{n,Float64})
    for j ∈ 1:n
        barw[j] = 1.0
        for k ∈ 1:n
            if k != j
                barw[j] = barw[j] * (xnodes[j] - xnodes[k])
            end
        end
        barw[j] = 1.0/barw[j]
    end

    numt, denomt = 0, 0
    @inbounds for j ∈ eachindex(xnodes)
        tdiff = t - xnodes[j]
        numt = numt + barw[j] / tdiff * fvals[j]
        denomt = denomt + barw[j] / tdiff
        if ( abs(tdiff) < 1e-14 )
            numt = fvals[j]
            denomt = 1.0
            break
        end
    end
    return numt / denomt
end

function Interp1d_Hardcoded( xn, fvals, t )
    
    barw = SA[1.0/((xn[1]-xn[2])*(xn[1]-xn[3])*(xn[1]-xn[4])), 1.0/((xn[2]-xn[1])*(xn[2]-xn[3])*(xn[2]-xn[4])), 1.0/((xn[3]-xn[1])*(xn[3]-xn[2])*(xn[3]-xn[4])), 1.0/((xn[4]-xn[1])*(xn[4]-xn[2])*(xn[4]-xn[3])) ]

    numt, denomt = 0, 0
    @inbounds for j ∈ eachindex(xn)
        tdiff = t - xn[j]
        numt = numt + barw[j] / tdiff * fvals[j]
        denomt = denomt + barw[j] / tdiff
        if ( abs(tdiff) < 1e-14 )
            numt = fvals[j]
            denomt = 1.0
            break
        end
    end
    return numt / denomt
end

function run_Interp1d_no_array()
    xnodes = LinRange(0,1,4)
    fvals = xnodes.^4
    t = 0.3
    
    for i ∈ 1:100000
        u = Interp1d_no_array( xnodes, fvals, t )
    end
end

function run_Interp1d_MVector()
    xnodes = LinRange(0,1,4)
    fvals = xnodes.^4
    t = 0.3
    
    n = Val(length(xnodes))
    for i ∈ 1:100000
        u = Interp1d_MVector( xnodes, fvals, t, n )
    end
end

function run_Interp1d_Hardcoded()
    xnodes = LinRange(0,1,4)
    fvals = xnodes.^4
    t = 0.3
    
    for i ∈ 1:100000
        u = Interp1d_Hardcoded( xnodes, fvals, t )
    end
end

@btime run_Interp1d_no_array()
4.481 ms (1 allocation: 96 bytes)

@btime run_Interp1d_MVector()
3.745 ms (1 allocation: 96 bytes)

@btime run_Interp1d_Hardcoded()
3.074 ms (1 allocation: 96 bytes)

```

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [November 12, 2022, 2:48am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/17 "2022-11-12T02:48:21Z")

</div>

So if the MArrays are working for you, I’d advise you to just use that. However, in your actual program will `n` always be knowable at compile time, or is it dynamic? Because that can severely change the timings. The `MVector` will also become very slow if `n` gets too big.

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 12, 2022, 4:33am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/18 "2022-11-12T04:33:01Z")

</div>

Yes, the `MArrays` are working fine for me now, in combination with Value-types.

`n` is definitively known at compile time, and is very small (in the range of 4-8). (The reason that performance matters to me is that this routine is called a gazillion times, each time with different data.)

The only thing that I’m leaving for future work is the possibility of using SIMD here. However, I couldn’t get around the `if` statements which seem to preclude its use. I guess I would probably need to find a way to rewrite the expressions.

---

<div class="post-metadata">

**Author:** ![jmair](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmair/32/35117_2.png) [@jmair](https://discourse.julialang.org/u/jmair)\
**Post date:** [November 12, 2022, 6:18am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/19 "2022-11-12T06:18:11Z")

</div>

You can vectorise the loop statement and avoid the branch, the following:

> [@lrnv](#):
>
> ```julia
> for k ∈ 1:n
> if k != j
> barw_j *= (xnodes[j] - xnodes[k])
> end
> end
> 
> ```

Becomes

```julia
for k ∈ 1:n
      barw_j *= (xnodes[j] - xnodes[k]) + (k==j)
 end

```

Which should get rid of the `if` statement. I don’t know how efficient it is but it might be worth checking if it works with `LoopVectorization`

---

<div class="post-metadata">

**Author:** ![martin.d.maas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/martin.d.maas/32/50964_2.png) [@martin.d.maas](https://discourse.julialang.org/u/martin.d.maas)\
**Post date:** [November 12, 2022, 6:41am UTC](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084/20 "2022-11-12T06:41:24Z")

</div>

> [@jmair](#):
>
> Which should get rid of the `if` statement. I don’t know how efficient it is but it might be worth checking if it works with `LoopVectorization`

Good try, but LoopVectorization.jl doesn’t recognize the expression…

```julia
Expression not recognized.
for k = 1:n
    #= interp_turbo.jl:11 =#
    barw[j] = mul_fast(add_fast(sub_fast(xn[j], xn[k]), k == j), barw[j])
    #= interp_turbo.jl:12 =#
end

```

I have another version that gets rid of the if statement and takes advantage of some symmetry (that also won’t work with LoopVectorization).

```julia
    @inbounds for j ∈ 1:n
        for k ∈ j+1:n
            xjk = (xn[j] - xn[k])
            barw[j] *= xjk
            barw[k] *= -xjk
        end
    end

```

but, curiously, this is yet slower than the first “naive” version with the if statement.

So maybe the compiler is, after all, performing some sort of vectorization in the first place.

[Next page](https://discourse.julialang.org/t/avoiding-allocations-of-small-but-non-trivial-arrays-work-array-alternative/90084.md?page=2)