Why are these two calculations which are exactly same giving different results in Fortran using gfortran?

Question

real, dimension(3), parameter :: boxlen = [4.0, 5.0, 7.0]
real, parameter :: mindist = 0.1
integer ::i

  write(*,"(A)") "Operation on array"
  print*, floor(boxlen/mindist)
  write(*,"(/A)") "Operation on individual elements"
  do i=1,3
     print*, i, floor(boxlen(i)/mindist)
  enddo

This is what I get when I run this code.

Operation on array
      40          50          70

Operation on individual elements
       1          39
       2          49
       3          69

Can someone explain why are the two calculations (one using operation on array and another using operation on individual elements) giving different results? I think they should be same.

Which gfortran version are you using? It is working fine on gfortran 5.1! What are your compile options? — Alexander Vogt, Jun 15 '15 at 10:24
With `i=1; print*, floor(boxlen(i)/mindist)-floor(boxlen(1)/mindist)` I get a surprise. That's with gfortran 4.8.1. — francescalus, Jun 15 '15 at 10:25
I am using GNU Fortran (GCC) 4.10.0 20140629 (experimental) [trunk revision 212119]. I compile using gfortran -Wall -std=f2008 -O3 — Yogesh Yadav, Jun 15 '15 at 10:26
The floating point divide by 0.1 is inexact so in principle either outcome is possible. Given obviously the same representation i think we should expect the same result though, — agentp, Jun 15 '15 at 11:36
I cannot reproduce with any of gfortran 4.7.4, 4.9.2 or 5.1.0 running on x86_64-pc-linux-gnu — casey, Jun 15 '15 at 14:09
For completeness, I reproduce with 4.8.1 with mingw-32 but not with 4.6 on x86_64-linux-gnu. — francescalus, Jun 15 '15 at 14:17
I updated my compiler to gfortran 5.1.0 with x86_64-w64-mingw32 built by [Equation Solution](http://www.Equation.com). Now the code runs file and shows expected behaviour. So, was it a bug in gfortran's old version? — Yogesh Yadav, Jun 15 '15 at 14:49
Also, on a side note, does my program's performance depend upon how the compiler was built? Specifically, I want to ask whether the compilers provided by [Equation Solution](http://www.equation.com) are good? Is there any better alternative (for latest gfortran) for Windows? — Yogesh Yadav, Jun 15 '15 at 14:54
@YogeshYadav, the only real "bug" out there is your way of thinking. You're are holding a wrong expectation that 4.0 divided by 0.1 results in 40.0. If you continue to program based on such wrong beliefs you are guaranteed to get in trouble sooner or later. — Wildcat, Jun 15 '15 at 14:54
@Wildcat I understand that floating point calculations are not exact. My question was not that it should be 40.0 but why were the results using operations on array and operations on individual elements different? — Yogesh Yadav, Jun 15 '15 at 14:56
@YogeshYadav, for your question about compilers, from my experience I would better recommend mingw-w64 binaries. — Wildcat, Jun 15 '15 at 14:58
@Wildcat But I can't find a mingw-w64 binary for gfortran 5.1.0. — Yogesh Yadav, Jun 15 '15 at 15:00
@YogeshYadav, I do not answer your question, my response was about your comment about some bug in a particular version gfortran. — Wildcat, Jun 15 '15 at 15:01
@YodeshYadav: I would avoid the equation.com builds. They seem to have plenty of equation-specific bugs, and since they haven't published their modified sources (in violation of the license) noone can help you debug them if they happen to hit you. Mingw-w64 and TDM builds seems fairly popular, if you want an alternative. I don't know whether they have 5.1 builds yet, but unless you're using the more esoteric OOP stuff or coarrays you can probably manage with an older version for the time being. — janneb, Jun 15 '15 at 16:26

score 0 · Answer 1 · answered Jun 15 '15 at 10:52

0

For me with GFortran 4.8.2 on x86_64-linux-gnu I get

Operation on array
          40          50          70

Operation on individual elements
           1          40
           2          50
           3          70

Peering into my crystall ball, if you're using 32-bit x86, there could be a difference in the order of how the computed values are stored to memory and printed.

Also, looking at the generated intermediate code with "-fdump-tree-original", it seems for the array operation, the division is done at compile time, (so should at least be correctly rounded thanks to MPFR).

answered Jun 15 '15 at 10:52

janneb

36,249
2
81
97

I am using 64-bit intel processor. – Yogesh Yadav Jun 15 '15 at 10:56
@YogeshYadav: Good. Are you also using a 64-bit GFortran? I.e. what does "gfortran -v" show for Target:? – janneb Jun 15 '15 at 10:58
For gfortran -v, it gives mingw32 as Target. – Yogesh Yadav Jun 15 '15 at 11:00
@YogeshYadav: Ok, that looks like a 32-bit GFortran then. You might want to try a 64-bit version, see e.g. https://gcc.gnu.org/wiki/GFortranBinaries – janneb Jun 15 '15 at 11:06
1

@janneb, changing a compiler is in my opinion a really bad advice, since the behaviour OP dos not like is totally fine and standard. The floating point arithmetic is inexact, so one should not expect in the first place to get 40.0 dividing 4.0 by 1.0. Consequently, OP actually needs to change his way of thinking about floating point arithmetic and the algorithm, but not the compiler. – Wildcat Jun 15 '15 at 14:50
1

@Wildcat: 4.0/1.0 equaling 40.0 might be standards conforming, but I suspect most users have more stringent requirements on the quality of the implementation. – janneb Jun 15 '15 at 16:16
@Wildcat: More seriously, the point of my suggestion was to narrow down the source of the issue to see whether it was due to x87 quirkyness or due to an actual bug in the compiler. – janneb Jun 15 '15 at 16:19
@Wildcat Your pessimism about floating point is unreasonable. The expressions are equivalent and the compiler should generate the SAME answer for both. Whether or not those answers end in 9 or 0 is debatable, but they must agree. The input data are the same and IEEE 754 floating point is deterministic. – Jeff Hammond Jun 15 '15 at 20:35
@Jeff, please, quote a part of Fortran standard which states that expressions in question are equivalent and that the results must agree. – Wildcat Jun 15 '15 at 21:04
@Jeff, the only thing which I found so far in Fortran 2008 draft is the following. 7.1.4.4 says that: "When an elemental binary operation is applied to a scalar and an array or to two arrays of the same shape, the operation is performed element-by-element on corresponding array elements of the array operands." From reading that I actually would agree that one should indeed expect the same results and that GFortran devs' decision to perform the array division at compile time using MPFR (I trust janneb on that matter) is indeed wrong. – Wildcat Jun 15 '15 at 21:17
Look at IEEE 754 and let me know where it says they won't. – Jeff Hammond Jun 15 '15 at 21:17
@Jeff, I must admit that my criticism of the floating-point arithmetic was wrong. You're right on that matter, IEEE floating-point arithmetic is deterministic. I'm suspicious now about that expressions in question are equivalent. And to the previous quote of the Standard I must add the following. 7.1.5.2.4 states that: "Two expressions of a numeric type are mathematically equivalent if, for all possible values of their primaries, their mathematical values are equal. However, mathematically equivalent expressions of numeric type may produce dierent computational results." – Wildcat Jun 15 '15 at 21:22
Yeah, that caveat is for when e.g. associativity of summation causes differences. I don't think any such transformation is possible here, except perhaps one outside IEEE754, such as a fast approximate reciprocal. – Jeff Hammond Jun 15 '15 at 21:30
So, the thing is, as far as I understand it, that we indeed do have two *mathematically equivalent* expressions out there, but still we are *not* guaranteed to have the same result. For instance, a compiler could replace division by 0.1 with multiplication by 10.0 in the case of the whole array operation and refuse to do so in the case of an explicit do loop. Again, I would not call such decision a smart idea, but it is still would be 100% standard, though odd, behaviour. – Wildcat Jun 15 '15 at 21:30
@Jeff, Fortran standard mentions in NOTE 7.19 that an expression like `A / 5.0` has an allowable alternative form `0.2 * A` that may be used in the evaluation. – Wildcat Jun 15 '15 at 21:35
@Wildcat Fair enough, but it would be rather unsettling for the compiler to transform x/0.1 to 10x inconsistently, particularly when two expressions are trivially equivalent. – Jeff Hammond Jun 15 '15 at 21:38
@Jeff, agree. So, whatever the reason for such inconsistency in GFortran exist, I think it should be eliminated. So it might be a good idea to send a bug report (an enhancement request). – Wildcat Jun 15 '15 at 21:49
@Wildcat: AFAICT the inconsistency that the OP saw is due to an x87 excess precision issue. If you open a bug report about that, it'd just be closed as yet another duplicate of https://gcc.gnu.org/bugzilla/show_bug.cgi?id=323 – janneb Jun 16 '15 at 06:14
The results are indeed same when using mingw64 build of gfortran. But while using 32bit gfortran, the both results should have been wrong ( inexact, actually). My point is that gfortran clearly uses different methods for division (and maybe other operations) for operations on array and operations on individual elements. I think this is a bug which should be reported. – Yogesh Yadav Jun 16 '15 at 07:46

Why are these two calculations which are exactly same giving different results in Fortran using gfortran?

1 Answers1