CODEWITHSUNDEEP
mathfloating-point
Alessandro B
Alessandro B

Reputation: 111

Consistent floating point arithmetic in fortran with different compilers on different computers?

How can I get the same results with floating-point arithmetic using different compilers and possibly also different computers?

This is p.f90

program fl
implicit none
real(kind=8) :: a

a=9d0/10d0
write(*,*) a

end program

With gfortran -o p p.f90 I get 0.90000000000000002

With ifort -g -o p p.f90 I get 0.900000000000000

Is there a way to achieve consistency? The small things propagate and become larger.

Best regards Alessandro

Upvotes: 5

Views: 1920

Answers (3)

Thorsten S.
Thorsten S.

Reputation: 4274

As Ali: You cannot.

It is a well known problem of FORTRAN compilers that results cannot be compared with different compilers and architectures. The main reason that it is impossible has to do with a design decision of the FORTRAN community: The computer may change the execution order and use internal computer specialities to optimize speed. As floating-point arithmetic is NOT distributive [ x(y+z) != xy + xz ] and the result DOES depend on the execution order

e.g.

1E308 (MAX) + (1E308 - 2E308) ~= 0E308

(1E308 (MAX) + 1E308) - 2E308 == INFINITY

you are stuck. Some compilers use e.g. the Fused multiply-add for calculating polynominals which give different results than the "normal" and this is quite allowed.

You are stuck. Sorry.

Upvotes: 0

High Performance Mark
High Performance Mark

Reputation: 78324

We'll turn to floating-point arithmetic below, but let me deal first with a misconception you have. Your use of list-directed formatting in the write statement means that the format chosen by the compiler to write out a variable is, well, the compiler's choice; it's not mandated by the language standard. Your use of the second * in write(*,*) tells the compiler to write out the value of a variable as it wishes. So what you have is not evidence that there is something wrong with the arithmetic, but that there are differences between gfortran and ifort. If I modify your write statement to

write(*,'(f21.18)') a

then my Intel Fortran program writes

0.900000000000000022

to the console.

SO is littered with questions arising from unfamiliarity with the details of floating-point arithmetic so I'm not going to write a treatise, just a few observations relevant to your question.

IEEE-754 64-bit floating-point numbers (which is probably what you get with the declaration real(kind=8)) only provide approximately 16 decimal digits of useful information. Actually, since they're binary and there is no easy correspondence between the number of digits in one base and in another, it's actually 15.95 decimal digits and many users round this down and never look at anything after the 15th significant figure in a floating-point number's decimal representation. So both ifort and gfortran are misleading you with their trailing 2s.

IEEE-754 defines not only the formats for f-p numbers, but also some rules for rounding and for arithmetic operations. A carefully-written program which uses only those arithmetic operations (I think square root is also specified) and which takes care about rounding modes and rounding operations should produce the same results on two different processors. Of course, not many useful numerical programs confine themselves to the basic arithmetic operations.

Since the 2003 standard Fortran has included an intrinsic module called ieee_arithmetic which gives the programmer direct access to the underlying IEEE-754 capabilities of the hardware on which it is running -- but note that it does not require that the hardware have any such capabilities. Using ieee_arithmetic and another intrinsic module called ieee_exceptions you should, if your hardware provides the necessary support, be able to write programs which compile under both gfortran and ifort and which produce, when executed, the same results to the last bit of every f-p number.

You'll also need to familiarise yourself with the optimisation options and numerical arithmetic options of the compiler(s) you're using. Most compilers have an option whose meaning is sacrifice IEEE compliance for speed (for ifort I think it's fp-model). In general adherence to IEEE-754 for operations will slow your programs.

Upvotes: 8

Ali
Ali

Reputation: 58471

You cannot.

See under Systems Aspects in What Every Computer Scientist Should Know About Floating-Point Arithmetic.

Upvotes: 1

Related Questions