# Memory allocation in dot product

**URL:** <https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901>\
**Category:** Performance\
**Tags:** question, memory-allocation\
**Created:** [September 9, 2021, 12:57am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901 "2021-09-09T00:57:00Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Frisus95](https://avatars.discourse-cdn.com/v4/letter/f/77aa72/32.png) [@Frisus95](https://discourse.julialang.org/u/Frisus95)\
**Post date:** [September 9, 2021, 12:57am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/1 "2021-09-09T00:57:00Z")

</div>

Hello everyone.

Is it possible to perform a dot product between multidimensional arrays without allocating memory?

For example:

```julia
function funza(A)
    
return dot(A,A)
    
end

```

```julia
A = rand(5,5,5,5,5);

@btime funza(A)

  428.323 ns (1 allocation: 16 bytes)

1041.8233479785777

```

I suspect that the memory allocation comes from the fact that dot() creates a scalar variable to store the result, but maybe I’m wrong. I cannot find a dot! function in the julia documentation…

In any case, how can I avoid allocating such memory?

Additional and “optional” questions:

- Is there a faster way to do this contraction (maybe using BLAS, unless the latter is used by default)?
- How this can be done on the GPU?

Thanks a lot to everyone!

---

<div class="post-metadata">

**Author:** ![Frisus95](https://avatars.discourse-cdn.com/v4/letter/f/77aa72/32.png) [@Frisus95](https://discourse.julialang.org/u/Frisus95)\
**Post date:** [September 9, 2021, 1:03am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/2 "2021-09-09T01:03:58Z")

</div>

I think I found myself the answer in the documentation of [https://github.com/JuliaSIMD/LoopVectorization.jl#dot-product](https://github.com/JuliaSIMD/LoopVectorization.jl#dot-product)

```julia
function mydotavx(A)
          s = 0.0
          @turbo for i ∈ eachindex(A)
              s += A[i]*A[i]
          end
          s
      end

```

```julia
@btime mydotavx($A)

  105.684 ns (0 allocations: 0 bytes)

1041.8233479785772

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [September 9, 2021, 1:43am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/3 "2021-09-09T01:43:06Z")

</div>

> [@Frisus95](#):
>
> In any case, how can I avoid allocating such memory?

```julia
julia> A = rand(5,5,5,5,5);

julia> @btime mydotavx($A)
  74.698 ns (0 allocations: 0 bytes)
1050.8122689068127

julia> @btime funza($A)
  130.944 ns (0 allocations: 0 bytes)
1050.8122689068127

```

I think it was allocating memory for you because you didn’t interpolate `A`:

> `@btime funza(A)`

Still, LoopVectorization/hand written loops should be much faster than `dot` on most CPUs here because most of the CPU’s execution resources are spent on loading memory. `A[i] * A[i]` needs only a single load instead of two.

If you want multithreading, you’ll need `@tturbo`. For problems too small for multithreading, stick with `@turbo`.

---

<div class="post-metadata">

**Author:** ![jzr](https://avatars.discourse-cdn.com/v4/letter/j/eb9ed0/32.png) [@jzr](https://discourse.julialang.org/u/jzr)\
**Post date:** [September 9, 2021, 3:51am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/4 "2021-09-09T03:51:33Z")

</div>

> [@Elrod](#):
>
> LoopVectorization/hand written loops should be much faster than `dot` on most CPUs here because most of the CPU’s execution resources are spent on loading memory. `A[i] * A[i]` needs only a single load instead of two.

Could LoopVectorization or LinearAlgebra transform `dot(A,A)` into a single-load version by recognizing the identical arguments?

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [September 9, 2021, 4:23am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/5 "2021-09-09T04:23:06Z")

</div>

Currently, LoopVectorization only has access to the information in the scope of the macro, therefore if given code like this:

```julia
function mydotavx_v2(A)
    B = A
    s = 0.0
    @turbo for i ∈ eachindex(A)
        s += A[i]*B[i]
    end
    s
end

```

it has no way of knowing that `A === B`. Maybe this can change once it switches to being a compiler plugin? Not sure what the scope of the analysis capabilities will be there.

In the above example, LoopVectorization’s code will involve two loads.  
LLVM will, however, optimize it into a single load.

```julia
julia> @btime mydotavx_v2($A)
  114.125 ns (0 allocations: 0 bytes)
1050.8122689068127

julia> @code_native syntax=:intel debuginfo=:none mydotavx_v2(A)

```

the loop:

```nohighlight
.LBB0_3: # %L43
                                        # =>This Inner Loop Header: Depth=1
        vmovupd zmm4, zmmword ptr [rax + 8*rdi]
        vmovupd zmm5, zmmword ptr [rax + 8*rdi + 64]
        vmovupd zmm6, zmmword ptr [rax + 8*rdi + 128]
        vmovupd zmm7, zmmword ptr [rax + 8*rdi + 192]
        vfmadd231pd zmm3, zmm4, zmm4 # zmm3 = (zmm4 * zmm4) + zmm3
        vfmadd231pd zmm2, zmm5, zmm5 # zmm2 = (zmm5 * zmm5) + zmm2
        vfmadd231pd zmm0, zmm6, zmm6 # zmm0 = (zmm6 * zmm6) + zmm0
        vfmadd231pd zmm1, zmm7, zmm7 # zmm1 = (zmm7 * zmm7) + zmm1
        add rdi, 32
        cmp rcx, rdi
        jne .LBB0_3

```

4 loads and 4 fmas.

Nonetheless, performance is lost because LoopVectorization’s optimization decisions were made assuming that it has to do two loads per fma. As a result, it assumes less dependency cutting is needed versus the actual single load case that we have here.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 9, 2021, 9:33am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/7 "2021-09-09T09:33:34Z")

</div>

Maybe have some variation of `norm()` with `p = 2` without applying the `sqrt()` operation (Maybe `norm2sqr()`)?  
Then the cases of a dot operation of a linear operator with itself will be efficiently calculated.

---

<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:** [September 9, 2021, 9:45am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/8 "2021-09-09T09:45:59Z")

</div>

> [@RoyiAvital](#):
>
> Maybe `norm2sqr()`

```julia
LinearAlgebra.norm_sqr

```

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 9, 2021, 9:48am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/9 "2021-09-09T09:48:11Z")

</div>

> [@lmiq](#):
>
> `norm_sqr`

Too bad it can’t be found on the [documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) (Where I looked for before posting).

Also:

```julia
help?> LinearAlgebra.norm_sqr
  No documentation found.

```

Moreover, strange result:

```julia
julia> vA = [2; 3]
2-element Vector{Int64}:
 2
 3

julia> LinearAlgebra.norm_sqr(vA)
12.999999999999998

```

There is no `sqrt()` applied, why the result isn’t straight 13?

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [September 9, 2021, 10:24am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/10 "2021-09-09T10:24:17Z")

</div>

Since it is an internal undocumented function, I wouldn’t have too many expectations on exactly how it works.

---

<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:** [September 9, 2021, 11:20am UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/11 "2021-09-09T11:20:04Z")

</div>

> [@Elrod](#):
>
> Maybe this can change once it switches to being a compiler plugin?

That’s a tantalizing remark! For someone who knows next to nothing about ‘compiler plugins’, what would this look like from the user perspective?

---

<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:** [September 9, 2021, 12:35pm UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/12 "2021-09-09T12:35:55Z")

</div>

> [@RoyiAvital](#):
>
> There is no `sqrt()` applied, why the result isn’t straight 13?

This is what the implementation says:

```julia
# faster computation of norm(x)^2, avoiding overflow for integers
norm_sqr(x) = norm(x)^2

```

The overflow part is clear. How that can be faster is beyond my understanding. It certainly isn’t faster than a straightforward implementation, although admittedly the functions give different results when one approaches the overflow scales, and certainly the base function is what is closer to being correct.

```julia
julia> x = rand(1:1000,1000);

julia> @btime LinearAlgebra.norm_sqr($x)
  2.903 μs (0 allocations: 0 bytes)
3.5071326500000006e8

julia> function n2(x)
          n = zero(eltype(x))
          @simd for el in x
             n += el^2
          end
          n
       end
n2 (generic function with 1 method)

julia> @btime n2($x)
  162.360 ns (0 allocations: 0 bytes)
350713265

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [September 9, 2021, 1:54pm UTC](https://discourse.julialang.org/t/memory-allocation-in-dot-product/67901/13 "2021-09-09T13:54:12Z")

</div>

If you have AVX512, converting to floats can be faster:

```julia
julia> function n2f(x)
          n = zero(eltype(x))
          @simd for el in x
             @fastmath n += Float64(el)^2
          end
          n
       end
n2f (generic function with 1 method)

julia> function n2(x)
          n = zero(eltype(x))
          @simd for el in x
             n += el^2
          end
          n
       end
n2 (generic function with 1 method)

julia> x = rand(1:1000,1000);

julia> @btime n2f($x)
  57.579 ns (0 allocations: 0 bytes)
3.36560308e8

julia> @btime n2($x)
  79.610 ns (0 allocations: 0 bytes)
336560308

```

If you don’t, it’ll be slower. `@turbo` will do much better on converting `Int64` to `Float64` than `@simd @fastmath` if you don’t have AVX512, but it’d probably still end up a fair bit slower than the non-converting version.

> [@DNF](#):
>
> That’s a tantalizing remark! For someone who knows next to nothing about ‘compiler plugins’, what would this look like from the user perspective?

I’d also love an answer to that question.  
I’m rewriting LoopVectorization in three parts, in the following order:

1. Loop representation and optimization.
2. Generating code that can be compiled and run from an optimized schedule.
3. Creating a representation from code.

So I’ll wait to start diving into compiler plugins when I actually get to “3.”, which is progressing at a glacial pace so far. Hopefully the work will have advanced and we’ll have some combination of documentation to look at / examples to look at / pandemic dies down so I can more easily ask people (i.e. Keno) questions.