2
$\begingroup$

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!

$\endgroup$
5
  • $\begingroup$ @njuffaGiven that the result may differ by 30 percent, perhaps I do need to distinguish between them. $\endgroup$
    – H. Zhou
    Mar 4, 2022 at 12:17
  • 2
    $\begingroup$ Firstly real*8 is never the right choice as it is not, and has never been, part of standard Fortran - learn about Fortran kind values which is the correct way to use multiple precisions in Fortran, see e.g. stackoverflow.com/questions/838310/fortran-90-kind-parameter Then as @njuffa says what you choose depends upon the accuracy you require. You might also look at en.wikipedia.org/wiki/Kahan_summation_algorithm which might give you the accuracy you require without having to go to quad precision. $\endgroup$
    – Ian Bush
    Mar 4, 2022 at 12:19
  • $\begingroup$ @IanBush Your comment is very helpful, thanks! Using this algorithm I got same answers for real/real*8 s. I will also have a look at the meaning of real*8! $\endgroup$
    – H. Zhou
    Mar 4, 2022 at 13:18
  • 1
    $\begingroup$ @IanBush Please do not answer in comments. $\endgroup$ Mar 4, 2022 at 15:23
  • $\begingroup$ @OP "which makes them slightly change in base 10" is incorrect. Conversion from binary32 to binary64 is exact. An error of this kind may have been introduced earlier, in conversions from strings or literal constants typed inside your source code. $\endgroup$ Mar 4, 2022 at 15:27

1 Answer 1

10
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thanks a lot for your nice answer! $\endgroup$
    – H. Zhou
    Mar 4, 2022 at 16:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.