parallel sum different from serial sum

I have implemented a simple program which does the sum of n numbers. I initially implemented a simple serial algorithm. For the parallel version I used the one in the documentation. The code is as follows:

numeri_somma = 1e10;
somma = 0;


tic
for i = 1:numeri_somma
    somma = somma + i;
end
toc

fprintf("%d\n", somma)


x = 0;
tic
parfor i = 1:numeri_somma
    x =  x + i;
end
toc

somma == x

The code gets a speedup however the sum calculated by the serial version and the one calculated by the parallel version does not coincide. In fact, the comparison returns false (logical 0). I wanted to understand why this happens.

Answer 1

[Image Analyst]

In your numerical analysis or linear algebra class did they ever talk about truncation error, where you're adding a much, much smaller number to a huge number? They should have. Might want to review that.

Answer 2

[Bruno Luong]

Consider this example of the sum of three numbers in different order:

(Your parallel change the order of the sum compared to non-parallel)

a = 1e16

a = 1.0000e+16

b = 1

b = 1

c = -a

c = -1.0000e+16

s1 = a + b + c % (a + b) + c

ans = 0

s2 = a + c + b % (a + c) + b

ans = 1

Answer 3

[James Tursa]

Bruno and IA have already given the reasons for this discrepancy. I would simply add that even though you are adding integers and you might have expected the order of adding to not matter in this case, the magnitude of the result still causes rounding errors to creep in which can change the result depending on the order of calculations. E.g., the exact mathematical result is:

n = sym(1e10);
n*(n+1)/2

ans = 

50000000005000000000

You can see that the above number contains 20 significant digits, but double precision numbers can only hold about 15-16 significant digits. So at some point the summing process in double precision will start rounding the results of the intermediate additions because there are values in the 1's digit that you are trying to add. Changing the order of calculations via parallel looping can then easily change the rounding that occurs and get a slightly different result.

For a smaller value of n where the sum doesn't exceed the limits of double precision you will get the same answer regardless of calculation order because all the intermediate integer calculations will be done exactly. But for larger values of n you can get a difference. For your n=1e10 example, the double precision sum doesn't match the exact mathematical result. E.g.,

x = 0;
for i=1:1e10
    x = x + i;
end
fprintf('%f\n',x);

50000000000067862528.000000

In fact, the exact mathematical result can't even be represented in double precision exactly. E.g., here is the closest number to the exact mathematical result that can be represented in double precision:

fprintf('%f\n',50000000005000000000)

50000000005000003584.000000

The double precision spacing between numbers in that range is much greater than 1, so trying to add numbers with values in the 1's digit to numbers in this range simply isn't going to work exactly:

fprintf('%f\n',eps(50000000005000000000))

8192.000000

[Image Analyst]

You might say which is expected to be the more accurate one. I would expect the parallel one to be more accurate since each sum (for each core) won't get as high.