[Rd] Non identical numerical results from R code vs C/C++ code?

[Renaud Gaujoux]

Ok.

I quickly tried using LDOUBLE wherever I could, but it did not changed 
the results. I might try harder...
I agree with you Barry, and I will re-double re-check my code.

Thank you both for your help.
Bests,
Renaud

On 10/09/2010 13:24, Barry Rowlingson wrote:
> On Fri, Sep 10, 2010 at 11:46 AM, Renaud Gaujoux
> <renaud at mancala.cbio.uct.ac.za>  wrote:
>    
>> Hi,
>>
>> suppose you have two versions of the same algorithm: one in pure R, the
>> other one in C/C++ called via .Call().
>> Assuming there is no bug in the implementations (i.e. they both do the same
>> thing), is there any well known reason why the C/C++ implementation could
>> return numerical results non identical to the one obtained from the pure R
>> code? (e.g. could it be rounding errors? please explain.)
>> Has anybody had a similar experience?
>>
>> By not identical, I mean very small differences (<  2.4 e-14), but enough to
>> have identical() returning FALSE. Maybe I should not bother, but I want to
>> be sure where the differences come from, at least by mere curiosity.
>>
>> Briefly the R code perform multiple matrix product; the C code is an
>> optimization of those specific products via custom for loops, where entries
>> are not computed in the same order, etc... which improves both memory usage
>> and speed. The result is theoretically the same.
>>      
>   I've had almost exactly this situation recently with an algorithm I
> first implemented in R and then in C. Guess what the problem was? Yes,
> a bug in the C code.
>
>   At first I thought everything was okay because the returned values
> were close-ish, and I thought 'oh, rounding error', but I wasn't happy
> about that. So I dug a bit deeper. This is worth doing even if you are
> sure there no bugs in the C code (remember: there's always one more
> bug).
>
>   First, construct a minimal dataset so you can demonstrate the problem
> with a manageable size of matrix. I went down to 7 data points. Then
> get the C to print out its inputs. Identical, and as expected? Good.
>
>   Now debug intermediate calculations, printing or saving from C and
> checking they are the same as the intermediate calculations from R. If
> possible, try returning intermediate calculations in C and checking
> identical() with R intermediates.
>
>   Eventually you will find out where the first diverge - and this is
> the important bit. It's no use just knowing the final answers come out
> different, it's likely your answer has a sensitive dependence on
> initial conditions. You need to track down when the first bits are
> different.
>
>   Or it could be rounding error...
>
> Barry
>    

 

###
UNIVERSITY OF CAPE TOWN 

This e-mail is subject to the UCT ICT policies and e-mail disclaimer published on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from +27 21 650 4500. This e-mail is intended only for the person(s) to whom it is addressed. If the e-mail has reached you in error, please notify the author. If you are not the intended recipient of the e-mail you may not use, disclose, copy, redirect or print the content. If this e-mail is not related to the business of UCT it is sent by the sender in the sender's individual capacity.

###