0

NOTE: Any solutions or work arounds to this issue need to occur in the MATLAB, it is the code where the issues arise. Making the fortran match the MATLAB is counter productive because then neither piece of code will work... I understand that the difference is because of differences in the way the fortran compiler and whatever the hell MATLAB does in place of compiling interpret single versus double precision floats but I am hoping that someone can help me come up with a solution to work around this.

I am working on debugging some code that I have translated from Fortran into MATLAB and have come across something that has me stumped. In both fortran and MATLAB I have the following line

pcnt = 0.9999*(-0.5+z2)

where z2 = 0.51482129528868548. The issue I am having is that there is a difference of 2.4594e-10 in the pcnt calculated in MATLAB and the pcnt calculated in fortran. I have confirmed that z2 is precisely the same (ie z2_matlab-z2_fortran=0) and so I am stumped. both z2 and pcnt are double precision (real*8 for fortran) in both Fortran and MATLAB so as far as I am concerned they should have exactly the same accuracy (since this is being run on the same machine).

Normally I wouldn't care about a difference this small but z2 ends up getting multiplied by a large number, which is then used to calculate an index later on and the small difference ends up being a difference of around 2 for the array index later, leading to huge errors later on in the algorithm (on the order of 1e6 for number that are only at most order 1e7).

Does anyone have any idea why this issue is happening and some way to fix it? I am performing this work on MATLAB R2011a and used the gfortran compiler to compile the fortran all on a 64 bit MacBook pro with an I5 (3rd gen I think) processor.

If someone has any suggestions please let me know because if I can't find a solution to this all of the translation of some 5000 lines of code I performed over the last two weeks would be pretty much worthless otherwise.

Also, any solutions would have to be to the MATLAB code because the Fortran code is the one that currently works.

Thanks in advance,

Andrew

Andrew
  • 693
  • 6
  • 19
  • 1
    It sounds like your real problem may be that the algorithm itself has numerical stability problems. Having said that, what are the actual values that you get for `pcnt` in each case ? – Paul R Jul 18 '14 at 14:11
  • possible duplicate of [Why Are Floating Point Numbers Inaccurate?](http://stackoverflow.com/questions/21895756/why-are-floating-point-numbers-inaccurate) – Alexander Vogt Jul 18 '14 at 14:12
  • Both `0.9999` and `0.5` are declared as *single precision*. When storing a single precision number as double precision, the compiler has to "guess" the last bits. If `z2` is of the order `1e-1`, then a difference of `1e-10` is due to those last bits. Obviously, MATLAB guesses differently the Fortran does. – Alexander Vogt Jul 18 '14 at 14:15
  • @PaulR I agree with you but unfortunately I have no control over that... quite frankly everything about the fortran code is poorly done but again, there's not much I can do about it (its one of the reasons it is being translated into MATLAB. – Andrew Jul 18 '14 at 14:18
  • @AlexanderVogt I disagree that this is a duplicate of that question. It is very similar however I wanted to see if anyone could help me work around this issue. I do agree with your second comment but I am hoping to find a way around that. – Andrew Jul 18 '14 at 14:22
  • A translation from Fortran to MATLAB? I usually go the other way around... is the performance of the code not an issue? – sigma Jul 18 '14 at 15:20
  • @sigma, the readability/maintainability of the code is an issue. Also, I'm just doing what I'm told as an intern so... – Andrew Jul 18 '14 at 19:42
  • Could you provide the answer of Fortran and MatLab you get except only the difference? – EJG89 Jul 21 '14 at 07:55

1 Answers1

5

The Fortran numeric literals are single precision unless the d modifier used, while MATLAB uses double as default numeric literal type. So maybe you should rewrite your pcnt expression like:

pcnt = 0.9999d+0 * (-0.5d+0 + z2)

Conversely, you should convert to single the numeric literals of MATLAB, in order to emulate the Fortran behaviour:

pcnt = single(0.9999) * (single(-0.5) + z2);

Later edit

In extremis, one should not rely on the different compilers' algorithms of interpreting numeric literals; but instead use the native (binary) representation of the said literals:

  • write unformatted from Fortran all the numeric literals that cause you grief (in this formula, 0.9999 is the most probable suspect, because 0.5 can be represented exactly in both FP precision), in a file (see WRITE).
  • load in MATLAB these values stored in the file (see fread()).
  • use the loaded values instead of numeric literals (which should have a suggestive name, as a good programming practice).
  • @AlexanderVogt Yeah, I had underscore typo, sorry. –  Jul 18 '14 at 14:20
  • while this may work to get the fortran to match the MATLAB I need to go the other way around. – Andrew Jul 18 '14 at 14:24
  • Also, I tried going the other way by converting .9999 and 0.5 to single in MATLAB before performing the operation but this had no effect. – Andrew Jul 18 '14 at 14:26
  • Isn't MATLAB always calculating in double precision? – Alexander Vogt Jul 18 '14 at 14:32
  • @Andrew The last resort would be to save the Fortran constants' binary representation in a file; and load them in MATLAB as `single`s. To my knowledge both use the hardware floating point (no FP runtime lib) and the order of operations must be the same. –  Jul 18 '14 at 14:34
  • Yes but you can convert values to single if you want using the `single` command. It may just re upconvert to double when it performs and operation though, I am not sure about that. – Andrew Jul 18 '14 at 14:35
  • @CST-Link how exactly would this work, I have never done something like that? – Andrew Jul 18 '14 at 14:36
  • It is the other way around: `pcnt = double(single(.9999)) * (double(single(-0.5)) + z2)`. This results in the same error of `2.4594e-10`. – sigma Jul 18 '14 at 14:38
  • @AlexanderVogt Yes, it is always converting the expression terms to `double` before evaluating the expression. If one tries, for example, to compare `double(single(0.9999))` with `0.9999` in MATLAB, one will find a difference of approx `1.6594e-08`, which means that the conversion to single then back to double will affect the result. –  Jul 18 '14 at 14:39
  • @sigma Your observation is correct; however, the conversion back to `double` will be done automatically by MATLAB. –  Jul 18 '14 at 14:40
  • 1
    @Andrew The idea would be to write unformatted all the numeric literals that cause you grief (in this formula, `0.9999` is the most probable suspect, because `0.5` can be represented exactly in both FP precision) from Fortran, in a file. See `WRITE`. Then load in MATLAB these values stored in the file, and use them instead of numeric literals (which should have a suggestive name, as a good programming practice). See `fread()`. –  Jul 18 '14 at 14:53
  • Ok, thanks @CST-Link. That what I thought you might mean but I wasn't sure. Also, this solution worked and now gives the correct values so I will accept the above answer. – Andrew Jul 18 '14 at 15:18
  • @Andrew I will add the comment to the answer, if you don't mind me modifying an accepted answer. :-) –  Jul 18 '14 at 15:30