Accuracy loss in single-precision Euclidean norm computation

Question

I do hydrodynamics simulations with Fortran and recently I met with this issue: I have a single-precision array b of length n=512**3, and I wish to compute the root-mean-squre of its elements. All I do is

s=0
do i=1,n
  s=s+b(i)**2
enddo

and then compute sqrt(s/n).

I found that it makes about $30\%$ difference in the result depending on whether I have declared s as single-precision (real) or double-precision (real*8). My question is which should be the correct way to do?

My naive understanding is that, when I use real*8 s, the code will first convert b(i) into double precision and then do the computation. In that case, single-precision values in b will be padded with zeros in base 2, which makes them slightly change in base 10. Even so, I thought it was still accurate in base 2 and there should be no harm of promoting them to double precision before the computation. On the other hand, if I stick to single precision then each time I take the square b(i)**2, there is some tiny error introduced. Based on these arguments I though real*8 should be the right choice.

If my thoughts above were correct, does it mean one need to go for quadruple precision when dealing with double-precision data?

I'm pretty new to numerical computations and Fortran, so any answers/comments will be helpful!

Answer 1

The squares are harmless (as long as they don't overflow/underflow), because the relative perturbation they introduce is of the order of the machine precision $u\approx 10^{-8}$. Your troubles here come from the $n\approx 10^8$ sums you are performing in sequence, which introduce an error of the order of $nu\approx 1$. You can confirm this by replacing each machine-sum $a\oplus b$ with $(a+b)(1+\delta_i)$ and bounding the error with $|\delta_i|\leq u$, or by checking the error analysis of scalar products and summation in Chapters 3 and 4 of Higham's accuracy and stability of numerical algorithms.

In any case, your possibilities are three:

• conversion to higher precision, as you have already tried.
• compensated summation (Kahan's algorithm), which has been suggested in the comments.
• pairwise summation, in which you sum numbers two by two recursively according by a binary tree. This reduces the $nu$ factor to $u\log_2 n$. It is the only strategy that does not require to increase the total number of floating-point operations, but the book-keeping involved may be more expensive than switching to the other methods.

This is explained also in Chapter 4 of Higham's book, which I suggest you to read in detail if you are interested in the details.