Round-off error due to order of addition using double precision variables.

[desensitized]

Hi. I seem to be having a problem with a Fortran 95 program where I have a very large aaray of double precision floating-point variables. Each element in the aaray are then added together in a random sequence. I am finding that the round-off error's vary significantly (and unacceptably) with multiple executions of the program.

What would be the best work-around for this problem?
- I could convert the effected aarays and variables to quad precision and then convert the anwser back to double precision after computation.
- I could convert my floating-point variable to intergers my multiplying them by a large scalar and converting them to integers. Perform calculation, and then convert back to real (double precision)

I'm sure people have incounterd this problem before, and i'm curious as to what solutions they came up with.

Thank you

[jimdempseyatthecove]

Some explanation of the ranges of the values within the arrays might help.

If the array contains groupings of fairly large numbers and fairly small numbers (one or both with binary fractions requiring significant numbers of bits) then you may experience problems when adding the small numbers to the large numbers (requires more bits of precision if numbers contain a significantly sized binary fraction).

One technique that may show improvement

SumInt = 0.0
SumFraction = 0.0
do i=1,size(A)
DINTpart = DINT(A(i))
FractionPart = A(i) - DINTpart
SumInt = SumInt +DINTpart
SumFraction = SumFraction + FractionPart
enddo
Sum = SumInt + SumFraction

The above loop should be reasonably fast as it has no conditionalbranching

Jim Dempsey

[TimP]

Kahan's algorithm would minimize the rounding errors due to non-optimum order of addition. It is necessary to ensure that optimization doesn't break the algorithm, e.g. ifort -fp:precise.
Kahan's scheme doesn't give as much accuracy as somehow persuading the compiler to sum in x87 extra precision

[jimdempseyatthecove]

Alternate method

SumHigh = 0.0
SumLow= 0.0

Threshold = 10E11_8! you pick threshold
do i=1,size(A)

if(A(i).ge.Threshold) then
SumHigh = SumHigh + A(i)
else
SumLow = SumLow + A(i)
endif
enddo
Sum = SumHigh + SumLow

Picking the threshold will require some knowledge of what is in the array.

Note, you can maintain the individual SumHigh and SumLow too.

Jim Dempsey

[desensitized]

Thanks for the responses. I should mention that the program i am talking about is not my code directly, but my manager's. I am a student at WSU working for a physics lab. With that being said, I do not know the exact magnitude of the numbers in the array other than they are of similar magnitued.

The problem in the code arises when trying to execute it in parallel. If the code is executed in series, the results are repeatable. With parallel execution (on a shared memory system), it is not repeatable.

The program is a simulation of shock loading through sand. The specific algorithm where problem is, occurs when adding the affect (stresses) of all grains on grain 'i'.
In series, all grain contributions 'i - n' , are added to 'i'
In parallel, the contributions are arbitrarly summed and then the value is assigned to grain 'i'.

Error starts off small, but gradually increases with the increasing number of grains and hence grain-stress-contributions.

Thanks,
James

[jimdempseyatthecove]

James,

It may help if you disclose the full loop including your attempt at parallel code statements. There are critical issues relating to when and how best to share variables.

Disclosing this in your first post would have been helpful.

Jim Dempsey