#8 Should use Kahan's summation in integrator

[Anonymous]

In integrator (RK4EX here), one can see

// a=a+a_k/6
for(long _i1=0; _i1<_main_size*_main_main_ncomponents; _i1++)
_main_main[_i1] += _akfield_main_main[_i1]/6; // HERE

// a_k = a_i + a_k/2
for(long _i1=0; _i1<_main_size*_main_main_ncomponents; _i1++)
_akfield_main_main[_i1] = _aifield_main_main[_i1] + _akfield_main_main[_i1]/2;

t += _step/2; // HERE

I think because we are adding small value to possibly big value many times, we will with high probability yield big errors.

If we are using constant stepsize, things like t = t0 + _i0 * step / 2; will make much better accuracy.

But for adaptative stepsize (is it used in xmds at all?), and for integration components the best way (without changing double to higher and much slower representation, like quad precision (only on ppc), quad-double (via emulation), or gmp (mpf)), is to use Kahan's summation.

Some basic informations can be found on wiki: http://en.wikipedia.org/wiki/Kahan_summation_algorithm

Remember to check if kahan's summation implementation plays well with inlineing and -O3, as many compiler can infere wrongly that it is noop. Using asm or violatile keyword can prevent this.

[Joe Hope]

This roundoff looks like it will be significant for certain simulations. However, compensating for all simulations will be an unnecessary slowdown for most.

Adaptive stepsize algorithms are the 'norm' for deterministic equations, and they in effect maximise the 'small' part of the addition. So the key troublesome simulations will be those that integrate non-trivial equations across extremely large dynamic ranges. Kahan summation may well be important for them.

On the downside, it will be a non-trivial implementation for everything apart from the propagation dimension variable. It will require a lot of new vectors to be defined.

Perhaps we should make this a flag?