[Rd] gfortran optimization problems
Dave Roberts
droberts at montana.edu
Thu Nov 13 19:35:51 CET 2008
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
More information about the R-devel
mailing list