# Calculating "triple" dot products
**URL:** https://discourse.julialang.org/t/calculating-triple-dot-products/124840
**Category:** Performance
**Created:** [January 16, 2025, 3:52pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840 "2025-01-16T15:52:36Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@kevinli1993](https://discourse.julialang.org/u/kevinli1993)
#### Post date: [January 16, 2025, 3:52pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/1 "2025-01-16T15:52:36Z")
I have Float64 arrays `x`, `y`, `z` and want to compute `sum(x[i] * y[i] * z[i])` over all indices `i`.
I’m curious what is the fastest way to calculate this?
Some ideas:
- a plain for loop
- Einsum.jl
- OMEinsum.jl
- `LinearAlgebra.dot` calls into BLAS.dot, but I’m not sure if there are a 3-vector version of this though.
---
### Author:  [@barucden](https://discourse.julialang.org/u/barucden)
#### Post date: [January 16, 2025, 4:17pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/2 "2025-01-16T16:17:12Z")
Have you already tried some approach that we can compare to?
A straightforward implementation:
```julia
function dot3(x, y, z)
s = zero(promote_type(eltype(x), eltype(y), eltype(z)))
@simd for i in eachindex(x, y, z)
s += x[i] * y[i] * z[i]
end
return s
end
```
I am getting
```julia-repl
julia> x = rand(1000); y = rand(1000); z = rand(1000);
julia> @btime dot3($x, $y, $z);
118.453 ns (0 allocations: 0 bytes)
```
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [January 16, 2025, 5:29pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/3 "2025-01-16T17:29:11Z")
> [@barucden](#):
>
> `eachindex(x, y, z)`
What are the cases in which this is useful?
---
### Author:  [@NonDairyNeutrino](https://discourse.julialang.org/u/NonDairyNeutrino)
#### Post date: [January 16, 2025, 5:43pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/4 "2025-01-16T17:43:32Z")
What about something more Julian like
```
sum(x .* y .* z)
```
```julia
julia> x = rand(1000); y = rand(1000); z = rand(1000);
julia> @btime sum(x .* y .* z);
253.994 ns (6 allocations: 7.95 KiB)
```
---
### Author:  [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [January 16, 2025, 6:02pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/5 "2025-01-16T18:02:25Z")
This `@simd for` version by `barucden` above is probably approximately as fast as you can get. There’s some chance it does some suboptimal vectorization and that something based on `LoopVectorization` (or another) might eek out some more, but I doubt it will make a big difference to your overall program. I suspect that computing `x`, `y`, and `z` in the first place costs significantly more than this sum ever will (for example, even `@btime fill!($x, 1);` costs ~40% of what the sum does).
There just isn’t much optimization that one can do for simple 1-loop calculations.
You can do something simpler like `mapreduce(*, +, x, y, z)` or `sum(prod, zip(x, y, z))`, but those have some open performance issues right now so aren’t competitive with the `@simd for` implementation.
---
### Author:  [@barucden](https://discourse.julialang.org/u/barucden)
#### Post date: [January 16, 2025, 6:02pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/6 "2025-01-16T18:02:46Z")
> [@rafael.guerra](#):
>
> > [@barucden](#):
> >
> > `eachindex(x, y, z)`
>
> What are the cases in which this is useful?
> If you supply more than one `AbstractArray` argument, `eachindex` will create an iterable object that is fast for all arguments
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [January 16, 2025, 6:17pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/7 "2025-01-16T18:17:08Z")
> [@NonDairyNeutrino](#):
>
> What about something more Julian like
>
> ```julia-auto
> sum(x .* y .* z)
>
> ```
I would not call that more Julian. It is very Julian to avoid unnecessary intermediate allocations, so I would avoid this.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [January 16, 2025, 6:47pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/8 "2025-01-16T18:47:01Z")
> [@barucden](#):
>
> If you supply more than one `AbstractArray` argument, `eachindex` will create an iterable object that is fast for all arguments
I asked because for the input example provided with Vector{Float64}, using only `eachindex(x)` gives exactly the same performance.
---
### Author:  [@bertschi](https://discourse.julialang.org/u/bertschi)
#### Post date: [January 16, 2025, 6:57pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/9 "2025-01-16T18:57:55Z")
It’s safer in terms of bound-checking:
```julia-repl
julia> x = rand(3); y = rand(2);
julia> for i in eachindex(x); @show x[i] + y[i] end
x[i] + y[i] = 1.0647502606457395
x[i] + y[i] = 1.0344221061391612
ERROR: BoundsError: attempt to access 2-element Vector{Float64} at index [3]
julia> for i in eachindex(x, y); @show x[i] + y[i] end
ERROR: DimensionMismatch: all inputs to eachindex must have the same indices, got Base.OneTo(3) and Base.OneTo(2)
```
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [January 16, 2025, 9:21pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/10 "2025-01-16T21:21:41Z")
> [@rafael.guerra](#):
>
> `eachindex(x)` gives exactly the same performance.
`x`, `y` and `z` are all similar, so of course there’s no difference. The point is if they are different in some way, especially to catch size differences, but also just to find a workable iterator for all of them.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [January 16, 2025, 9:33pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/11 "2025-01-16T21:33:31Z")
I have read the help docstring but I do not have a good picture of practical examples:
```julia
help?> eachindex
...
If you supply more than one AbstractArray argument, eachindex will create an iterable object that is fast for all
arguments (typically a UnitRange if all inputs have fast linear indexing, a CartesianIndices otherwise).
```
---
### Author:  [@kevinli1993](https://discourse.julialang.org/u/kevinli1993)
#### Post date: [January 22, 2025, 2:32pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/12 "2025-01-22T14:32:05Z")
Thanks all - the responses are great and I learned more about Julia in the process.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [January 22, 2025, 2:43pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/13 "2025-01-22T14:43:39Z")
@rafael.guerra an important aspect is that `eachindex` throws when there’s nothing reasonable to return:
```julia-repl
julia> eachindex(3:5, 3:7)
ERROR: DimensionMismatch: all inputs to eachindex must have the same indices, got Base.OneTo(3) and Base.OneTo(5)
```
Putting a bounds check before a loop helps eliminate bounds checks within the loop body. `eachindex` throwing makes the `eachindex` call a bounds check.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [January 22, 2025, 2:58pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/14 "2025-01-22T14:58:12Z")
@nsajko, thank you. I think @bertschi and others were pointing to that same advantage of the multiple arrays arguments to eachindex.
I was just trying to better understand the help doc string.
Also the bound checks do not seem to impact the performance for the example provided.
> **Benchmark code example**
>
> ```julia
> function dotproduct3a(x, y, z)
> s = zero(eltype(x))
> @simd for i in eachindex(x, y, z) #? note multiple argument eachindex()
> s += x[i]*y[i]*z[i]
> end
> return s
> end
>
> function dotproduct3b(x, y, z)
> s = zero(eltype(x))
> @simd for i in eachindex(x) #? note single argument eachindex()
> s += x[i]*y[i]*z[i]
> end
> return s
> end
>
> n = 100_000
> x = rand(n); y = rand(n); z = rand(n);
>
> using BenchmarkTools
> @btime dotproduct3a($x, $y, $z) # 30.6 μs (0 allocs: 0 bytes)
> @btime dotproduct3b($x, $y, $z) # 30.6 μs (0 allocs: 0 bytes)
>
> ```
> **versioninfo()**
>
> Julia Version 1.11.2
> Commit 5e9a32e7af (2024-12-01 20:02 UTC)
> Build Info:
> Official [https://julialang.org/](https://julialang.org/) release
> Platform Info:
> OS: Windows (x86\_64-w64-mingw32)
> CPU: 20 × 13th Gen Intel(R) Core™ i9-13900H
> WORD\_SIZE: 64
> LLVM: libLLVM-16.0.6 (ORCJIT, goldmont)
> Threads: 1 default, 0 interactive, 1 GC (on 20 virtual cores)
> Environment:
> JULIA\_PKG\_USE\_CLI\_GIT = true
> JULIA\_STACKFRAME\_FUNCTION\_COLOR = blue
> JULIA\_WARN\_COLOR = cyan
> JULIA\_EDITOR = code.cmd -g
> JULIA\_NUM\_THREADS = 8
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [January 22, 2025, 3:29pm UTC](https://discourse.julialang.org/t/calculating-triple-dot-products/124840/15 "2025-01-22T15:29:49Z")
OK, here’s an example where passing all collections to `eachindex` is important for safety. The above example uses `@simd`, which transforms the code of the loop, and I don’t know what transformations exactly does `@simd` do, so in this example I replaced `@simd` with `@fastmath` and `@inbounds`. This keeps the same performance but is easier to understand and discuss, given that there’s no magic code generation:
```julia
function dot3(x, y, z)
s = zero(promote_type(eltype(x), eltype(y), eltype(z)))
@inline @fastmath for i in eachindex(x, y, z)
s += @inbounds x[i] * y[i] * z[i]
end
s
end
```
In this example, it’s important that we have `eachindex(x, y, z)` instead of `eachindex(x)`, because we want an out-of-bounds access to result in a throw, instead of in UB.
With `@simd` it appears that `eachindex(x)` is adequate, at least with the current implementation of `@simd`, but `@simd` is a bit of a spooky black box, not something to base best practices on.