accuracy of intrinsic SUM function
Steve Kargl
sgk@troutmask.apl.washington.edu
Fri Aug 24 23:02:00 GMT 2012
On Fri, Aug 24, 2012 at 10:24:23PM +0100, N.M. Maclaren wrote:
> On Aug 24 2012, Steve Kargl wrote:
> >
> >The above algorithm probably represents a simple
> >change to gfortran current in-lining of the naive
> >algorithm. Any other algorithm would involve
> >a much more complicated change, and would probably
> >lead to a function call to a library routine.
>
> My brief test indicates that it's a complete waste of effort.
> I have coded up the algorithms, with pairwise needing 8.log_2(N)
> words of scratch space, and here are some results on random
> data (in C++, but so what?):
>
> Please input a size:
> 100000000
> Methods Value Error Time
> Canonical 5.363276169432254392e+03
> Kahan 5.363276169432254392e+03 3.16e-17 0.804
> Extended 5.363276169432254392e+03 1.29e-17 0.213
> Pairwise 5.363276169432252573e+03 3.08e-16 0.555
> Forwards 5.363276169433530413e+03 2.38e-13 0.241
> numeric::accumulate 5.363276169433530413e+03 2.38e-13 0.241
>
> This is on a single AMD system, but the fastest and most accurate
> method is simply to use Intel native 80-bit arithmetic for a
> sequential accumulation. That makes fancy methods futile except
> on non-Intel systems. However, it might well be worth doing if
> it doesn't already (I haven't checked).
>
> I may try with some malice-aforethought data distributions, to see
> if I can make pairwise less accurate on a distribution that is likely
> to occur in practice, but not tonight.
>
> You are welcome to my code, incidentally.
Of course, your 80-bit arithmetic loses at least 49-bits of precision
when you sum the contents of an array declared as 'real(16) x(10000)'.
:-)
--
Steve
More information about the Fortran
mailing list