why? I am crazy now!!!

[mmquan]

program main
real*8 eps(4)
real*8 a,b,c,d
eps(1) = Z'4037CB020C49BA5E'
eps(2) = Z'403C1E76C8B43958'
eps(3) = Z'BFF943D46B26BF87'
eps(4) = Z'BFF943D46B26BF87'
write(*,'(z17.16)') eps(1)*eps(2) - eps(3)*eps(4)
a = Z'4037CB020C49BA5E'
b = Z'403C1E76C8B43958'
c = Z'BFF943D46B26BF87'
d = Z'BFF943D46B26BF87'
write(*,'(z17.16)') a*b-c*d
end
save this as t.f. I am using Intel fortran compiler for linux 7.1. if I compile this without any optimization like: ifc -O0 -o t t.f then the two outputs get the same. but if I do some optimitation, like: ifc -O2 -o t t.f then it gets different(the last bit). But in Java I always get the second output as I did some optimitation in fortran. And who can tell me which one is more accurate???

[pedro_gonnet]

This sounds like the old 80-bit Intel register story... The double-precision floating point registers on the Pentium chips stores 80 bits intsead of just 64. Therefore, if you do a computation keeping all variables in the registers, you're going to get extra precision in the end. If, however, you write these values back to the sack and read them in again, you truncate them back down to 64 bits and any subsequent operations on those values have "normal" precision. Most Java JIT comilers try to keep all the values in registers and so will the optimized fortran code. The unoptimized code will store each variable after each operation, therefore nukeing any added precision.

This is hell for anybody trying to do numerics up to machine precision :( Is there any compiler flag that tells the FP-unit to knock it off and behave like any normal double-precision unit?

Cheers
Pedro

[Steven_L_Intel1]

-mp

[pedro_gonnet]

Well, no, -mp doesn't do it for me... According to the System User's Guide, -mp does "floating-point user variables declared as floating-point types are not assigned to registers", which is a far cry both conceptually and in performance to simply setting the Floating Point Precision Control Word. Or am I missing something here?

Cheers and thanks
Pedro

[mmquan]

Thanks for your responding!!! But why I do this operation then they get the same results???
program main
real*8 eps(4)
real*8 a,b,c,d
eps(1) = Z'4037CB020C49BA5E'
eps(2) = Z'403C1E76C8B43958'
eps(3) = Z'BFF943D46B26BF87'
eps(4) = Z'BFF943D46B26BF87'
write(*,'(z17.16)') eps(1)*eps(2) - eps(3)*eps(4)
a = eps(1)
b = eps(2)
c = eps(3)
d = eps(4)
write(*,'(z17.16)') a*b-c*d
end
I just assigned the value from array to single variables. Why? Now I want to get the same output in fortran and java, who can tell me how? and the "-mp" does not work, I don't why. after I used -mp the output will be both the unoptimized results, but I need the output like java!!! Or who can tell me how to make java output with the same as fortran? thanks a lot!

quan

[pedro_gonnet]

The documentation on -mp is not quite clear, but apparently it only works on variables and not on temporary values in expressions... This would mean that

write(*,'(z17.16)') eps(1)*eps(2) - eps(3)*eps(4)

is evaluated with 80 bits, whereas with

a = eps(1)
b = eps(2)
c = eps(3)
d = eps(4)
write(*,'(z17.16)') a*b-c*d

the... Aha. I see your point. You're right, the difference in the two printed values is a bug :(

The reason you can't get the same results as with Java? That is what you get for using processors that do not conform to the IEEE double-precision specification :)

[mmquan]

So, you mean if I change another processor, then they should get the same output?

[pedro_gonnet]

Can you run that code on a SPARC-platform or a PowerPC?

[Steven_L_Intel1]

Actually, the IEEE spec does include the 80-bit extended precision format. It's just that X86 processors are unusual in supporting it.

If you compile with -xW or -xN, you'll get SSE2 code that will tend not to use extended-precision.

[mmquan]

Yes, I do have SPARC here. But I have no fortran compiler for SPARC :(

[mmquan]

if I compile with "ifc -xW -o t t.f" it can not run on my linux system. I don't know why? It tells me "** Illegal Instruction **". If I compile with "ifc -xN -o t t.f" then the result is still the same as optimized code!

[Steven_L_Intel1]

You have too many unknowns here - and I also think it is unrealistic to expect absolutely the same floating point computations from Fortran code as you would see in Java. Minor differences in code generation and optimization can cause small differences.

In general, the Fortran computations which use extended precision, or double-precision for single-precision values, give you BETTER results than computations which restrict themselves to single precision.

[mmquan]

Yes, I have said before I am only a starter in programming and know really a little about computer programming. But I don't expect too much, because as you see, throwing Java away, even in fortran the output(if I assign the array item to single variable) is also different. Yes it is a light difference, but if I divided these two slightly different numbers with 0.125 and do 10,000 iteration, then you can see the difference! Now I am doing a simulation program, it's about molecular dynamic in Chemistry. It will run for 3 month before I get the result! So I wonder if I directly use the array elements do multiplication and subtraction, is it more accurate than single variables or not?
This is very strict! Thanks !

quan

[pedro_gonnet]

I found what I was looking for yesterday... With the subroutines GETCONTROLFPQQ and SETCONTROLFPQQ you can set the FP-precision to 53 bits (the size of a 64-bit floating point mantissa).

I have used these subroutines before to set the trapping behaviour, but this makes your fortran code non-portable :(

Steve, you're the Intel-guy :) What was the rationale for not having these values controlled by compiler switches?

Cheers and thanks
Pedro

P.S. the issue of losing/gaining a bit is purely academic. I'm more interested in the implications for numerical algorithms, i.e. termination criteria.

[TimP]

The Windows compilers have initialized the x87 mode to 53 bit precision for several years now. Until recently, the linux compilers did not, because this settingwill break C long double code and produce different results from other linux compiler defaults. Now that SSE2 code has become the usual option, the question is becoming moot. If you don't want 64-bit precision mode, except for C long double, use one of the SSE2 options.

I did see one of the recent ifort compilers setting 53-bit precision at startup, which surprised me. I think this is undesirable in view of the implication that it would break any C long double function which might be called later, unless you knew that GET/SETCONTROLFPQQ are available now in the linux compiler, and should be used in such cases.

Until recently, complex math functions, such as sin(), have been more accurate on linux than in Windows, due to the difference in precision mode setting.