Migrating Fortran source code to Matlab precision issues

I have a piece of Fortran code that I have recently converted to Matlab. As I am going through and validating the translation, I'm running into an error with precision that is killing the performance of the Matlab code. It is an iterative code that repeats the same operations over many time steps, so the numerical differences that arise from the differences in precision end up building up and lead to wrong pieces of logic.

The issue right now is that in Fortran the single precision float variables XX and YY, are summed and being multiplied by a constant:

       VAR = XX + YY
       T1 = VAR*CONST

So T1 at this point is a float. The next line in Fortran casts the value of T1 to an integer:

       IT1 = T1

In my Matlab scripts I accomplish this with the following:

      var = xx +yy;   
      t1 = var*const;
      it1 = fix(t1);

This works until about the 700th time step and the value of it1 in Matlab diverges from the one in Fortran. The value in Fortran of T1 = 1.0004226, fixing IT1 = 1, and in Matlab t1 = 0.99868, fixing it1 = 0;

My attempts at correcting this problem so far have been to single( ) all the values that go into the calculations of both XX and YY in Matlab. But this gets me to the 800th time step or so, and it falls apart again.

Does anyone have some words of wisdom that can help me with this problem?? I've been wracking my brain about this for the past couple of days and have come up with nothing.

[Geoff]

So the result in t1 is then used in some way for the next calculation? In that case, unless the value is bit-for-bit identical to the one in Fortran, you will surely diverge sooner or later.

Can you post the Fortaran code and test scenario?

Have you tried writing it in C for a second opinion?

[James Tursa]

Ditto. I would like to see the Fortran code and the MATLAB code side by side.

[Chris]

The source code of Fortran is proprietary so I don't have the ability to post either the Fortran or Matlab equivalent on here. The next line of code is something along the lines of

VAL = IT1*CONST2

where CONST2 is a hard coded constant. Val is used in about 5 other different files numerous times in each file.

I was provided with too many files to try converting it all to C. The conversion to Matlab itself took about 2 weeks to finish.

[Geoff]

After having googled myself a 5-second introduction to Fortran, can I ask what type the constant is? If the types are not the same, Fortran does implicit type conversion in much the same manner as C. You need to be careful to mimic this in every single calculation, as MatLab has its own type conversion that behaves a little differently.

[Chris]

Both constants were declared as REAL = single precision floating point

[Geoff]

The disparity between your two programs seems much more significant than just single precision error. I'd be looking for accidental rounding or incorrect typecasting elsewhere in the code.

Answer 1

[John D'Errico]

So it sounds as if in Fortran, you are using SINGLE precision variables, and it is deviation from what you get in MATLAB, which will use double precision by default. The point is, the Fortran computation is probably the one going wrong, as it was done in lower precision.

If you want some clarity, use my HPF tool to verify which values are correct.

http://www.mathworks.com/matlabcentral/fileexchange/36534-hpf-a-big-decimal-class-currently-in-beta-release

[Geoff]

I disagree on principle here =) If an algorithm is designed to use single precision, then it is not necessarily better to throw double precision at it. The way I read this question, it is a curiosity that using single precision in MatLab is giving different results to using single precision in Fortran.

However, I wonder about the integer conversion and whether there is overflow etc. That would be where we tread very carefully in MatLab because it does too much for us behind the scenes when mixing integers and reals.

[Chris]

Well the Fortran is producing the correct answer in the end, whereas the Matlab isn't. So I don't think the Fortran computation is going wrong. I know that each compiler performs operations such as sine and cosine differently.

The way the YY variable above is created is a matrix multiply of a 3x3 and a 3x1 where the elements in the 3x3 matrix are sums or products of sines and cosines of different state variables. So I was thinking maybe that could be an issue.

Similarly, I was wondering if it would be worth pursuing the idea of declaring each of my variables in Fortran as doubles to see if that matches up with Matlab.

[James Tursa]

Well, my 2 cents on this at this point:

- I don't know if the original algorithm was designed for single or not ... maybe it is really old Fortran code that used single for storage or throughput reasons and double precision would have been better from the start. OP doesn't really say one way or the other.

- If I think there is a precision issue then I would generally do what John suggests also, use higher precision variables and compare intermediate output to see what is going on.

- Without giving us any code to look at it is generally hopeless trying to help. The issue could be that the original algorithm is the one going wrong ... could be the MATLAB code is going wrong ... could be they are *both* wrong and need higher precision ... could be the conversion to MATLAB has typos in it or in some other way is not replicating the Fortran algorithm ... etc etc etc.

Without anything concrete to go on, the only thing I can suggest is OP prints a lot of intermediate output on both sides, compare them, and see what is going on and which one is right (if either).

BTW, just because the variables in Fortran are single precision doesn't mean the calculations are done in single precision ... on a PC it is likely that many intermediate calculations are done in 80-bit extended precision and then the result gets converted to single for storage.

[Chris]

Believe me I would like to be able to post up the code to get whatever help I can, but it's just not an option. It is really old Fortran code though.

I am doing a comparison check and printing the state vectors out of Fortran and comparing the results to Fortran is how I have been validating the progress of the simulation so far. There comes a point where the numerical differences in Matlab between the two variables I mentioned above, T1 in Fortran and t1 in Matlab, are

T1 = 1.0001 and t1 = 0.998

Upon fixing these two values to integers I get two different integers. Doing the following:

VAL = IT1*CONST2 !Fortran Code

val = it1*const2; %Matlab Code

produces two different numbers for val. Val is a state variable that is used to update other state variables.

That's interesting though about how intermediate calculations are performed. I'm working on RedHat Linux with running intel processors, and doing my compilation with GNU gfortran, and not a Windows PC ( I'm assuming you were referring to a Windows PC when you said PC)

[John D'Errico]

The problem is, WE can't know what you have done wrong. Have you made a mistake in the re-write into MATLAB? Maybe. Is there a precision problem? Maybe. The only thing you can do (since you want show us the code, probably for a very good reason) is to check each computation. Write it out. Look at the intermediate results to see where the two diverge. My guess is that since this is the result of a long series of iterations when the two diverge, that tells me it is the accumulation of many tiny errors. That is what happens when you may be storing results in single precision, EVEN if they are done in higher precision registers internally. Eventually the loss of precision in the stored intermediate results causes a significant error. Again, only you will be able to verify this, since we cannot inspect what is probably a lengthy code, one where only you knows what it should be doing.