Expected precision of Fortran FE constant evaluation

[Salvatore Filippone]

That's exactly what I did a long time ago "by hand" (this is a well
know algorithm, you can find it in one of the execises of TAOCP).
However, at the very least the implementation should contain the
algorithm in LAPACK, which does the scaling with just one pass:
---------
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  -- This version written on 25-October-1982.
*>     Modified on 14-October-1993 to inline the call to DLASSQ.
*>     Sven Hammarling, Nag Ltd.
*> \endverbatim
*>
*  =====================================================================
      DOUBLE PRECISION FUNCTION DNRM2(N,X,INCX)
*
*  -- Reference BLAS level1 routine (version 3.8.0) --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2017
*
*     .. Scalar Arguments ..
      INTEGER INCX,N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION X(*)
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION ONE,ZERO
      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ABSXI,NORM,SCALE,SSQ
      INTEGER IX
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC ABS,SQRT
*     ..
      IF (N.LT.1 .OR. INCX.LT.1) THEN
          NORM = ZERO
      ELSE IF (N.EQ.1) THEN
          NORM = ABS(X(1))
      ELSE
          SCALE = ZERO
          SSQ = ONE
*        The following loop is equivalent to this call to the LAPACK
*        auxiliary routine:
*        CALL DLASSQ( N, X, INCX, SCALE, SSQ )
*
          DO 10 IX = 1,1 + (N-1)*INCX,INCX
              IF (X(IX).NE.ZERO) THEN
                  ABSXI = ABS(X(IX))
                  IF (SCALE.LT.ABSXI) THEN
                      SSQ = ONE + SSQ* (SCALE/ABSXI)**2
                      SCALE = ABSXI
                  ELSE
                      SSQ = SSQ + (ABSXI/SCALE)**2
                  END IF
              END IF
   10     CONTINUE
          NORM = SCALE*SQRT(SSQ)
      END IF
*
      DNRM2 = NORM
      RETURN
*
*     End of DNRM2.
*
      END
----------------------

On Mon, Feb 11, 2019 at 7:40 PM Thomas König <tk@tkoenig.net> wrote:
>
> 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...