Re: Performance and rounding (?) problem ?

[Christophe Meessen <meessen at cppm dot in2p3 dot fr>]

Andrew Haley a écrit :
> It's possible that
>
>   res[i] += err[i] / k[i];
>
> might give different results on x86 from
>
>   err[i] /= k[i];
>   res[i] += err[i];
>
> due to intermediate 80-bit precision being used for the result.
>   
How could I test this precision issue ? 

I noticed that I didn't mention that k is constant.

One significant difference is that in one form err is modified and in 
the other not. I checked the code and err is not used beyond that 
instruction. In the next iteration err is cleared (memset) before 
filling it with data.

I see a similar problem if the operation is

res[i] *= err[i]/k[i];

I also precomputed 1/k[i] and tried
res[i] *= err[i]*k[i];

Same problem and exactly same intermediate results as with the initial form.


I also tested this and it fails too

const float tmp = err[i] / k[i];
res[i] += tmp;

I guess the compiler optimizes the tmp away.
> However, if you have a 5 Gb matrix I'm assuming you are on a 64-bit system,
> and 64-bit platforms don't use the 80-bit intermediate precision.
>   
Do you mean the hardware or the software ? How is this possible ?
I'm using the 64bit Ubuntu distribution and 64bit g++ and icc.
> What does "uname -a" say?
>   
Linux XXX 2.6.24-19-generic #1 SMP Wed Aug 20 17:53:40 UTC 2008 x86_64 
GNU/Linux