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