# Fortran benchmark

**URL:** <https://fortran-lang.discourse.group/t/fortran-benchmark/2055>\
**Category:** Meta\
**Created:** [October 12, 2021, 5:58am UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055 "2021-10-12T05:58:43Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![fortran4r](https://avatars.discourse-cdn.com/v4/letter/f/9de0a6/32.png) [@fortran4r](https://fortran-lang.discourse.group/u/fortran4r)\
**Post date:** [October 12, 2021, 5:58am UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/1 "2021-10-12T05:58:43Z")

</div>

In this benchmark, Fortran is slower than C in every program. Is it possible to improve the Fortran code? Since this website is quite popular, if the Fortran code can be improved, that would be a great marketing point.

[https://benchmarksgame-team.pages.debian.net/benchmarksgame/fastest/fortran.html](https://benchmarksgame-team.pages.debian.net/benchmarksgame/fastest/fortran.html)

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 12, 2021, 3:42pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/2 "2021-10-12T15:42:14Z")

</div>

We can add it here: [GitHub - fortran-lang/benchmarks: Fortran benchmarks](https://github.com/fortran-lang/benchmarks) and work on it.

Last time I looked at the nbody benchmark, the bottleneck was 1/sqrt(x). We would have to look at the assembly that compilers generate and see.

Down the road I would like LFortran to generate optimal assembly for things like 1/sqrt(x), those are basic tricks and it’s actually used quite often in any kind of a particle simulation (MD, n-body, etc.).

---

<div class="post-metadata">

**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\
**Post date:** [October 12, 2021, 5:08pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/3 "2021-10-12T17:08:07Z")

</div>

> [@certik](#):
>
> the bottleneck was 1/sqrt(x)

Whoever wants to improve that benchmark should take a look at this then. [Fast inverse square root - Wikipedia](https://en.wikipedia.org/wiki/Fast_inverse_square_root)

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 12, 2021, 5:36pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/4 "2021-10-12T17:36:16Z")

</div>

> [@everythingfunctional](#):
>
> Whoever wants to improve that benchmark should take a look at this then. [Fast inverse square root - Wikipedia](https://en.wikipedia.org/wiki/Fast_inverse_square_root)

That version is actually slow on modern hardware. I assume it was fast on the hardware back then.

Rather, the solution is to emit the same assembly as the C code, which uses optimal assembly. It’s the line ` z = _mm256_rsqrt_pd(z);`.

---

<div class="post-metadata">

**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\
**Post date:** [October 12, 2021, 6:10pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/5 "2021-10-12T18:10:41Z")

</div>

And this is why Knuth’s adage “premature optimization is root of all evil” always seems to hold. Optimizations like this are always premature, because they’ll be out of date when the newer hardware comes out 😜 .

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 12, 2021, 10:16pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/7 "2021-10-12T22:16:46Z")

</div>

> [@kargl](#):
>
> Define “optimal”. Replacing `sqrt(x)` and `1/sqrt(x)` by the numerical literal constant `1` is likely optimal with respect to speed. It is also optimal for `x = 1` . It may not be too optimal for `x /= 1` .

Optimal is getting the same answer as they do to the same number of significant digits, in this case it seems 9 significant digits, as fast as you can:

```auto
PROGRAM OUTPUT:
-0.169075164
-0.169059907

```

So replacing `1/sqrt(x)` by 1 will change the answer, thus one cannot do that.

The program must solve the same equations, i.e. you can change the number of steps and you must get the correct answer out.

We want to keep the Fortran program, it is actually quite readable: [n-body Classic Fortran&nbsp;#6 program (Benchmarks Game)](https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-ifc-6.html), but we want the compiler to generate code that runs as fast as the non-portable C code.

---

<div class="post-metadata">

**Author:** ![themos](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/themos/32/521_2.png) [@themos](https://fortran-lang.discourse.group/u/themos)\
**Post date:** [October 12, 2021, 10:59pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/8 "2021-10-12T22:59:53Z")

</div>

And never uses the COMPLEX data type???

---

<div class="post-metadata">

**Author:** ![themos](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/themos/32/521_2.png) [@themos](https://fortran-lang.discourse.group/u/themos)\
**Post date:** [October 12, 2021, 11:44pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/11 "2021-10-12T23:44:17Z")

</div>

The Fortran mandelbrot program does not use COMPLEX, which I found odd.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 13, 2021, 12:37am UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/13 "2021-10-13T00:37:03Z")

</div>

Here is the relevant part from [https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-ifc-6.html:](https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-ifc-6.html:)

```auto
       d2 = sum(r(:,m)**2)
       distance = 1/sqrt(real(d2))
       distance = distance * (1.5d0 - 0.5d0 * d2 * distance * distance)

```

and it’s the energy printout at the end:

```auto
  e = energy(x, v, mass)
  write (*,'(f12.9)') e

```

That must be correct to 9 decimal places. So what we can do is keep everything else the same and play with the accuracy of sqrt to see what is the minimal accuracy that still gives the correct final answer.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 13, 2021, 2:28am UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/15 "2021-10-13T02:28:36Z")

</div>

Excellent point, I haven’t noticed it. It’s either a mistake, or it is on purpose. Given that it (I assume) returns the correct answer, it means the sqrt can be in single precision (or less) and still recover total energy to 9 digits.

---

<div class="post-metadata">

**Author:** ![vmagnin](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/vmagnin/32/28_2.png) [@vmagnin](https://fortran-lang.discourse.group/u/vmagnin)\
**Post date:** [October 13, 2021, 6:37am UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/16 "2021-10-13T06:37:00Z")

</div>

I do not really understand that part of the algorithm, but if `real(d2)` and `d2` had the same type, the second line would seem useless:

```auto
       distance = 1/sqrt(real(d2))
       distance = distance * (1.5d0 - 0.5d0 * d2 * distance * distance)

```

because it could be analytically simplified and we would obtain (at least mathematically):

```auto
       distance = 1/sqrt(real(d2))
       distance = distance

```

or am I missing something?

---

<div class="post-metadata">

**Author:** ![themos](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/themos/32/521_2.png) [@themos](https://fortran-lang.discourse.group/u/themos)\
**Post date:** [October 13, 2021, 9:45am UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/17 "2021-10-13T09:45:46Z")

</div>

It’s a Newton iteration to improve accuracy.

---

<div class="post-metadata">

**Author:** ![Euler-37](https://avatars.discourse-cdn.com/v4/letter/e/b19c9b/32.png) [@Euler-37](https://fortran-lang.discourse.group/u/Euler-37)\
**Post date:** [October 13, 2021, 11:25am UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/18 "2021-10-13T11:25:54Z")

</div>

It looks strange, because kind of d2 is `dp(kind=8)`, real(d2) will change it to `kind=4 `,so the precision is decrease, and then the iteration is useful maybe. So the Newton iteration looks like for fix the type casting loss of `kind=8` to `kind=4`( 🤔Bug maybe)

Actually,if I change it like this

```fortran
       distance = 1/sqrt(d2)
       !distance = 1/sqrt(real(d2))
       !distance = distance * (1.5d0 - 0.5d0 * d2 * distance * distance)

```

The result is same

* * *

Based on @mecej4 answer , I see, the optimize principle is the single precision sqrt , type casting and 4 multiplication may faster than double precision sqrt. 😮

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [October 13, 2021, 12:09pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/19 "2021-10-13T12:09:25Z")

</div>

It looks like a trick to fool a specific compiler into emitting faster code for a specific FPU. With other compilers/machines, this code may “shoot itself in the foot”.

Here is an interpretation, expanding on the short answer by Themos:

The first of the two lines

```auto
       distance = 1/sqrt(real(d2))
       distance = distance * (1.5d0 - 0.5d0 * d2 * distance * distance)

```

calculates a single-precision approximation, presumably using a machine instruction that is faster than the corresponding double precision instruction. The next line uses Newton iteration, supposedly doubling the number of significant digits that are correct.

If we are using IEEE floats with 23+1 and 52+1 bits for the mantissa, the Newton step takes us from 24 to 48 bits, so the final result may be less precise than the result from a simple use of the double precision square root instruction.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 13, 2021, 3:46pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/20 "2021-10-13T15:46:25Z")

</div>

Indeed, I think that is exactly right. The single precision `1/sqrt(x)` will (presumably) generate the faster single precision instruction, and then the Newton iteration increases the precision, not quite to a full double, but enough to be accurate to 9 digits in the final energy. You can also notice the actual code has a second Newton iteration that is commented out:

```auto
       d2 = sum(r(:,m)**2)
       distance = 1/sqrt(real(d2))
       distance = distance * (1.5d0 - 0.5d0 * d2 * distance * distance)
       !distance = distance * (1.5d0 - 0.5d0 * d2 * distance * distance)
       mag(m) = tstep * distance**3

```

And the reason it is commented out is that we don’t need the extra accuracy. That is a common practice with these things.

I will note that this is such a commonly used technique that compilers like GFortran even have code to do this automatically ([31723 – Use reciprocal and reciprocal square root with -ffast-math](https://gcc.gnu.org/bugzilla/show_bug.cgi?id=31723)) — but the issue is how many iterations the compiler should do (to be safe it really should do enough iterations to recover full double precision accuracy), and since this problem requires lower accuracy, then presumably one has to do it by hand to control the accuracy vs speed better.

Note that the [C version](https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-gcc-9.html) is simply using this routine to do the same:

```auto
// utilize vrsqrtps to compute an approximation of 1/sqrt(x) with float,
// cast back to double and refine using a variation of
// Goldschmidt’s algorithm to get < 1e-9 error
static inline __m256d _mm256_rsqrt_pd(__m256d s) {
    __m128 q = _mm256_cvtpd_ps(s);
    q = _mm_rsqrt_ps(q);
    __m256d x = _mm256_cvtps_pd(q);
    __m256d y = s * x * x;
    __m256d a = _mm256_mul_pd(y, _mm256_set1_pd(0.375));
    a = _mm256_mul_pd(a, y);
    __m256d b = _mm256_mul_pd(y, _mm256_set1_pd(1.25));
    b = _mm256_sub_pd(b, _mm256_set1_pd(1.875));
    y = _mm256_sub_pd(a, b);
    x = _mm256_mul_pd(x, y);
    return x;
}

```

It’s very hard to read, but if you translate the assembly into Fortran line by line, here is what it is doing (except that the assembly is operating on 4 doubles at once, but the Fortran compiler can vectorize in a similar manner):

```fortran
! Computes x = 1/sqrt(s)
real(dp) function rsqrt(s) result(x)
real(dp), intent(in) :: s
real(dp) :: y, a, b
real(sp) :: q
q = real(s, sp)
q = 1/sqrt(q)
x = real(q, dp)
y = s * x * x
a = y * 0.375_dp
a = a * y
b = y * 1.25_dp
b = b - 1.875_dp
y = a - b
x = x * y
end function

```

I tested the last function and it works. It is actually using 2 Newton iterations, an equivalent of:

```fortran
r = 1/sqrt(real(x,sp))
r = r * (1.5_dp - 0.5_dp * x * r**2)
r = r * (1.5_dp - 0.5_dp * x * r**2)

```

Plugging the second equation into the third:

```auto
r = 1/sqrt(real(x,sp))
r = r * (1.5_dp - 0.5_dp * x * r**2) * (1.5_dp - 0.5_dp * x * r**2 * (1.5_dp - 0.5_dp * x * r**2)**2)

```

and expanding we get:

```auto
r = 1/sqrt(real(x,sp))
r = r * (0.0625*x **4*r** 8 - 0.5625*x **3*r** 6 + 1.6875*x **2*r** 4 - 2.4375*x*r**2 + 2.25)

```

And rearranging the `rsqrt` function it is doing:

```auto
r = 1/sqrt(real(x,sp))
r = r * (0.375*x **2*r** 4 - 1.25*x*r**2 + 1.875)

```

Numerically it seems both are getting the same answer, so I am guessing they are somehow related using the r = 1/sqrt(x) relation. [Here](https://gist.github.com/certik/44f34283b69b9e1c9fae7c9e319f696e) is the full code that I used to test it.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [October 13, 2021, 7:48pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/21 "2021-10-13T19:48:19Z")

</div>

Perhaps making the procedure `elemental` and adding

```fortran
!$omp declare simd(rsqrt) [simdlen(n)]

```

would aid in achieving vectorization for the Fortran version. I wonder if setting the optional clause `simdlen(4)` would help the compiler emit a version equivalent to the C version. A few other options are discussed in the article [Explicit Vector Programming in Fortran](https://software.intel.com/content/www/us/en/develop/articles/explicit-vector-programming-in-fortran.html) from Intel.

---

<div class="post-metadata">

**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\
**Post date:** [October 13, 2021, 8:01pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/22 "2021-10-13T20:01:14Z")

</div>

In my experience GFortran correctly vectorizes such functions out of the box most of the time (with `-O3 -march=native -ffast-math -funroll-loops`), sometimes one needs to help it a bit such as by splitting the loop or rewriting the code in various ways. Here is an example where I apply the elemental `rsqrt` function on an array and you can see that GFortran correctly vectorized it: [Compiler Explorer](https://godbolt.org/z/jvsnrx941), from first look the assembly looks very good. I think for simple things like `rsqrt`, it should be possible to make GFortran generate the same assembly as the “C” code (really an assembly code).

---

<div class="post-metadata">

**Author:** ![billlong](https://avatars.discourse-cdn.com/v4/letter/b/71e660/32.png) [@billlong](https://fortran-lang.discourse.group/u/billlong)\
**Post date:** [October 13, 2021, 8:44pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/23 "2021-10-13T20:44:48Z")

</div>

Compilers have been converting 1/sqrt(x) to rsqrt(x) for some time.

---

<div class="post-metadata">

**Author:** ![billlong](https://avatars.discourse-cdn.com/v4/letter/b/71e660/32.png) [@billlong](https://fortran-lang.discourse.group/u/billlong)\
**Post date:** [October 13, 2021, 8:46pm UTC](https://fortran-lang.discourse.group/t/fortran-benchmark/2055/24 "2021-10-13T20:46:05Z")

</div>

Odd indeed. Any Mandelbrot code I ever wrote always used COMPLEX variables.