[Numpy-discussion] numpy and roundoff(?)

[Damian Eads]

Lisandro Dalcin wrote:
> On 3/1/08, Charles R Harris <charlesr.harris at gmail.com> wrote:
>> So they differ in the least significant bit. Not surprising, I expect the
>> Fortran compiler might well perform operations in different order,
>> accumulate in different places, etc. It might also accumulate in higher
>> precision registers or round differently depending on hardware and various
>> flags.
> 
> Of course, but a completely unrelated but equivalent C implementation
> of this problem, as you can check in line 313 at this link
> 
> http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/snes/examples/tutorials/ex5.c.html
> 
> behaves almost the same that my F90 implemented residual. Perhaps
> Fortran compiler (gfortran) will generate the same code as the C one,
> but I'm not sure, Fortran compilers can be smarter that C compilers
> for this kind of looping.
> 
>> The exp functions in Fortran and C might also return slightly
>> different results.
> 
> I believe this is not the source of the problem, I've tried commenting
> that term, and differences are still there.
> 
>> I don't think the differences are significant, but if you
>> really want to compare results you will need a higher precision solution to
>> compare against.
> 
> I agree, the differences are not significant, but they end up having a
> noticeable impact. I'm still surprised!.
> 
> Let's stop all this now. I'll be back as soon as I can produce some
> self-contained code to show and reproducing the problem.

At work we noticed a significant difference in results occurring in two 
versions of our code, an earlier version written in C++, and a later 
version written in Python/numpy. The algorithms were structured about 
the same. My colleague found the cause of the discrepancy; it turned out 
to be a difference in the way numpy and the C++ program were compiled. 
One used -mfpmath=sse, and the other, -mfpmath=387. Keeping them both 
the same cleared the discrepancy.

Damian