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
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!