[Rd] gfortran optimization problems

[Dave Roberts]

Bill and Mike,

    Thanks for your help.  I think Mike was right about the 
80-bit/64-bit compare.  As Mike thought, the error occurs at a test for 
equality, i.e.

       if (x .ge. y) then

I know better than to test reals for equality, but this is a closed 
interval test, and it still fails.

       if (x-y .gt. -0.00001)

worked. Bill, thanks for the tip on --ffloat-store.  I tried

gfortran -c alt_duleg.f
gcc -std=gnu99 -ffloat-store -shared -L/usr/local/lib -o alt_duleg.so
alt_duleg.o -lgfortran -lm -lgcc_s

(standard R CMD SHLIB except for the addition of --float-store)

and it seems to have cleared up the problem.  I should have seen this
coming, since I once had a comparison fail on a 4-byte real versus a
double precision that I know are exactly the same in base10.  To fix 
that I canged EVERYTHING to double precision, but I didn't see the 
register problem coning at all.

Now I need to figure out how to implement a makefile for my package,
but that's another story and easily solved I suspect.

Thanks to everyone who took a stab at this one.  I learned a lot.
Sorry to have taken so long to get back to it, but other priorities 
demanded so.

Dave
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
David W. Roberts                                     office 406-994-4548
Professor and Head                                      FAX 406-994-3190
Department of Ecology                         email droberts at montana.edu
Montana State University
Bozeman, MT 59717-3460


William Dunlap wrote:
>  
> 
>> -----Original Message-----
>> From: William Dunlap 
>> Sent: Monday, November 03, 2008 8:34 AM
>> To: 'Mike Prager'; 'droberts at montana.edu'
>> Subject: RE: [Rd] gfortran optimization problems
> ...
>> Another common problem, in C or Fortran, is that optimization 
>> can result in a very small floating point number remaining in 
>> an 80-bit floating-point-unit register instead of being 
>> truncated to the 64-bits available in main memory.  The 
>> difference in roundoff error can result in big differences in 
>> results if your algorithm is ill conditioned or has a test 
>> for an exact 0.0.
>> Adding a debugging WRITE or printf statement generally causes 
>> the variable to be copied to 64-bit main memory.
> 
> A quick way to see if using floating point registers or not
> affects your results is to add the gcc flag (probably a gfortran
> flag as well) '-ffloat-store'.  'man gcc' says:
> 
>        The following options control compiler behavior regarding
> floating
>        point arithmetic.  These options trade off between speed and
> correct-
>        ness.  All must be specifically enabled.
> 
>        -ffloat-store
>            Do not store floating point variables in registers, and
> inhibit
>            other options that might change whether a floating point
> value is
>            taken from a register or memory.
> 
>            This option prevents undesirable excess precision on machines
> such
>            as the 68000 where the floating registers (of the 68881) keep
> more
>            precision than a "double" is supposed to have.  Similarly for
> the
>            x86 architecture.  For most programs, the excess precision
> does
>            only good, but a few programs rely on the precise definition
> of
>            IEEE floating point.  Use -ffloat-store for such programs,
> after
>            modifying them to store all pertinent intermediate
> computations
>            into variables.
> 
> Bill Dunlap
> TIBCO Spotfire Inc
> wdunlap tibco.com 
> 
> 
> 


-- 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
David W. Roberts                                     office 406-994-4548
Professor and Head                                      FAX 406-994-3190
Department of Ecology                         email droberts at montana.edu
Montana State University
Bozeman, MT 59717-3460