Expected precision of Fortran FE constant evaluation

[Thomas König]

Hi Salvatore,

> Right, try this preprint
> https://www.researchgate.net/publication/298896236_Remark_on_Algorithm_539_Carefully_Computing_the_Euclidean_Norm
> Apart from typesetting to an older style it looks essentially the same

Interesting reading.

Regarding the algorithm the authors propose, I have one concern.
The first part of their algorithm is (translating their pseudocode into
Fortran, resulting errors are mine)

   real :: sm, s, t

   sm = 0.
   s = 0.
   do i=1, n
      s = s + x(i)**2
      t = sm
      sm = t + s
      s = s + (t-sm)
   end do
   norm2 = sqrt(sm)

Never translate this with -Ofast (but everybody knew that already ;-)

But this is also a problem for vectorization, because the result of s
is used in the next iteration. So, if we want to have real
high performance on an architecture offering SIMD instructions
with 4 reals in parallel, we will probably have to do something like
(I've left in Fortran's array notation to make the vectors explicit)

   real :: sm(4), s(4), t(4)
   sm = 0.
   s = 0.
   do i=1, n, 4
     s(1:4) = s + x(i:i+3)**2
     t(1:4) = sm(1:4)
     sm(1:4) = t(1:4) + s(1:4)
     s(1:4) = s(1:4) + (t(1:4)-sm(1:4))
   end do
   norm2 = sqrt(sm(1) + sm(2) + sm(3) + sm(4))

and we would have to do all the loop peeling etc by hand.

Hmm...