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 [EMAIL PROTECTED]
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'; '[EMAIL PROTECTED]'
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 [EMAIL PROTECTED]
Montana State University
Bozeman, MT 59717-3460
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
