# Naive dot product faster in Fortran than in Juila

**URL:** <https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317>\
**Category:** Performance\
**Created:** [June 21, 2021, 5:24pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317 "2021-06-21T17:24:50Z")\
**Posts on this page:** 13\
**Page:** 1

<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:** [June 21, 2021, 5:24pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/1 "2021-06-21T17:24:51Z")

</div>

Does anyone knows why this simple dot product is faster in Fortran then in Julia?

(`LinearAlgebra.dot` is faster and with `LoopVectorization` the Julia code gets also faster. but this is not what I am questioning. I am curious about why without those packages the Fortran compiler is doing a better job than the Julia one).

Here Julia takes about 2x the time of the Fortran execution (compiled with `-O3`).

Julia:

```julia
using BenchmarkTools

function mydot(a, b, N)
  c = 0
  @inbounds @simd for i = 1:N
    c += a[i]*b[i]
  end
  c
end
N = 500_000_000
a = rand(Float32,N) 
b = rand(Float32,N)

@btime mydot($a,$b,$N)

```

Fortran:

```fortran
implicit none

real(real32) :: c
integer, parameter :: N = 500000000
real(real32), allocatable :: a(:), b(:)
real(real32) :: time1, time2

allocate(a(N), b(N))

call random_number(a)
call random_number(b)

call cpu_time(time1)
c = mydot(a, b, N)
call cpu_time(time2)
print *, c
print '("Time: ",f19.17," s")', time2 - time1

contains

function mydot(a, b, N) result(c)
   real(real32), intent(in) :: a(:)
   real(real32), intent(in) :: b(:)
   integer, intent(in) :: N
   real(real32) :: c
   integer :: i
   c = 0.
   do i = 1, N
      c = c + a(i)*b(i)
   enddo
endfunction
endprogram

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 21, 2021, 5:27pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/2 "2021-06-21T17:27:44Z")

</div>

Type stability. You should declare `c=zero(eltype(a))`. The reason this was slowing you down is adding an Int to a Float32 promotes to Float64.

---

<div class="post-metadata">

**Author:** ![jacobadenbaum](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@jacobadenbaum](https://discourse.julialang.org/u/jacobadenbaum)\
**Post date:** [June 21, 2021, 5:34pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/3 "2021-06-21T17:34:44Z")

</div>

Yeah for me that took the execution time down by a full factor of 5

---

<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:** [June 21, 2021, 5:35pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/4 "2021-06-21T17:35:39Z")

</div>

I thought in such a case the compiler would figure that out and do union-splitting automatically. (the problem is clearly shown by `@code_warntype`, which I didn’t because my first feeling was that it this promoting would be harmless).

---

<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:** [June 21, 2021, 5:56pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/5 "2021-06-21T17:56:00Z")

</div>

Those are some long arrays, so this will be completely memory bound.

If you try smaller arrays (once the Julia code is fixed), `gfortran` will be slower unless you dive deep into optimization options, using `-funroll-loops -fvariable-expansion-in-unroller`.

By default, LLVM will unroll the loop and use separate accumulators to split the dependency chain.  
Both these optimizations need to be turned on manually with `gfortran`.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 21, 2021, 6:00pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/6 "2021-06-21T18:00:19Z")

</div>

The problem isn’t the type instability. It’s the promotion type. By using Float64, you lose the ability to use fma, and get halve your vector with.

---

<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:** [June 21, 2021, 6:03pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/7 "2021-06-21T18:03:18Z")

</div>

Ah, I see. That makes sense :-).

---

<div class="post-metadata">

**Author:** ![zdenek\_hurak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zdenek_hurak/32/53118_2.png) [@zdenek\_hurak](https://discourse.julialang.org/u/zdenek_hurak)\
**Post date:** [June 21, 2021, 6:26pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/8 "2021-06-21T18:26:55Z")

</div>

Could you elaborate a bit, please? Only a bit will perhaps suffice.

---

<div class="post-metadata">

**Author:** ![mgkuhn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mgkuhn/32/6276_2.png) [@mgkuhn](https://discourse.julialang.org/u/mgkuhn)\
**Post date:** [June 21, 2021, 6:57pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/9 "2021-06-21T18:57:12Z")

</div>

The [fused multiply–add machine instruction](https://en.wikipedia.org/wiki/Multiply%E2%80%93accumulate_operation) in modern CPUs (with such SIMD extensions) performs several operations of the form `a = b*c + a` in a single clock cycle, but the compiler can only use it if all three variables involved have the same type.

---

<div class="post-metadata">

**Author:** ![mgkuhn](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mgkuhn/32/6276_2.png) [@mgkuhn](https://discourse.julialang.org/u/mgkuhn)\
**Post date:** [June 21, 2021, 7:12pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/10 "2021-06-21T19:12:32Z")

</div>

The following type annotation may prevent you from calling an SIMD-unfriendly method instance of this function:

```julia
function mydot(a::AbstractVector{T},
               b::AbstractVector{T}, N) where T
  c::T = 0
  @inbounds @simd for i = 1:N
    c += a[i]*b[i]
  end
  c
end

```

---

<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:** [July 24, 2021, 9:58pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/11 "2021-07-24T21:58:43Z")

</div>

> [@Oscar\_Smith](#):
>
> The problem isn’t the type instability. It’s the promotion type. By using Float64, you lose the ability to use fma, and get halve your vector with.

While you’re correct that `Float64` will have half throughput of `Float32` I am not sure about your statement regarding the `FMA`. As far as I know, `FMA` for `Float64` is supported on `AVX2` / `AVX512`.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [July 24, 2021, 10:00pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/12 "2021-07-24T22:00:22Z")

</div>

I should have been more clear probably. Float64 fma exists, but there isn’t fma between Float64 and Float32

---

<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:** [July 24, 2021, 11:03pm UTC](https://discourse.julialang.org/t/naive-dot-product-faster-in-fortran-than-in-juila/63317/13 "2021-07-24T23:03:42Z")

</div>

More specifically, in calculating

```julia
c::Float64 = c::Float64 + a::Float32 * b::Float32

```

either

1. promote `a` to `Float64`
2. promote `b` to `Float64`
3. `c = fma(a64, b64, c)`  
or
4. multiply `ab = a*b`
5. promote `ab` to `Float64`
6. add `ab64 + c`.

On most x86 CPUs, the conversion has a reciprocal throughput of 1, while the arithmetic tends to have a r-throughput of 0.5 or 1 (depending on the CPU).  
Thus promoting once tends to be at least as fast as promoting twice.