From renaud at mancala.cbio.uct.ac.za  Fri Sep 10 16:55:29 2010
From: renaud at mancala.cbio.uct.ac.za (Renaud Gaujoux)
Date: Fri, 10 Sep 2010 16:55:29 +0200
Subject: [Rd] Non identical numerical results from R code vs C/C++ code?
In-Reply-To: <D611103AA7EE3B4DAE7F7D49C72B291A02D5AAF7@EXMAIL2.bocad.bank-banque-canada.ca>
References: <4C8A0C86.8030307@cbio.uct.ac.za>
	<D611103AA7EE3B4DAE7F7D49C72B291A02D5AAF7@EXMAIL2.bocad.bank-banque-canada.ca>
Message-ID: <4C8A46E1.1020203@cbio.uct.ac.za>

Thanks Paul for the hints.
After some tests, reducing portion of my code, I found that simply doing 
a naive computation of 'crossprod' in C does NOT give exactly the same 
results as calling the Fortran underlying routine (dsyrk) as used in the 
R source code.
I will try the double 0.0 to see if it makes a difference.

What do you mean by

"You also need to be careful about
conversion from integer to real, to do double conversion."


?
Where are the trap in this type of conversion?
Thanks.

Renaud

On 10/09/2010 16:07, Paul Gilbert wrote:
> With fortran I have always managed to be able to get identical results
> on the same computer with the same OS. You will have trouble if you
> switch OS or hardware, or try the same hardware and OS with different
> math libraries. All the real calculations need to be double, even
> intermediate variables. Also, I've had trouble with arrays not being
> initialized to double 0.0. If you initialize to single 0.0 the
> straggling bits can cause differences. You also need to be careful about
> conversion from integer to real, to do double conversion. I'm not sure
> about C, but I would guess there are some of the same problems.
>
> Paul
>
>    
>> -----Original Message-----
>> From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-
>> project.org] On Behalf Of Renaud Gaujoux
>> Sent: September 10, 2010 6:47 AM
>> To: r-devel at r-project.org
>> Subject: [Rd] Non identical numerical results from R code vs C/C++
>>      
> code?
>    
>> 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.
>>
>> Thank you,
>> Renaud
>>
>> --
>> Renaud Gaujoux
>> Computational Biology - University of Cape Town
>> South Africa
>>
>>
>>
>>
>> ###
>> 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.
>>
>> ###
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>      
> ====================================================================================
>
> La version fran?aise suit le texte anglais.
>
> ------------------------------------------------------------------------------------
>
> This email may contain privileged and/or confidential ...{{dropped:27}}



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