Issue correlating exp() from fortran to matlab

Important clarification - I need Matab to match Fortran. I cannot change the Fortran code. I am working in 2018b and the function I am having an issue with is gz = 1.0 / (1+0.0247*exp(0.650*log(0.1089*re))). re is single precision in this case (Fortran). I was able to correlate 0.650*log(0.1089*re) so I know that the issue is with exp. In fortran, exp(0.650*log(0.1089*re)) = 24.09045, (24.0904521942139 if double precision), but in Matlab it is 24.0904541, (24.090454101562500 double precision). How can I modifiy the code to be the same as fortran. I have tried double() and single() everywhere but it didn't work.

The weirdest part: This code is part of a function that is used >500 times. I would say 350/500 times the results match up, but there are random stretches where it is off as described above. There is also a near similar function gx (some of the constants are different) that works 100% of the time. I made sure that I used single() and double() in the same exact spots but still got this issue.

[dpb]

You can't even guarantee the same Fortran code will be completely identical on two systems, what more a comparison between MATLAB and Fortran. The order in which arguments are evaluated can easily cause differences in rounding.

But, 24.09045 and 24.0904541 don't even have same precision being looked at; that's controllable by output formatting.

I don't think you can modify the code to generically guarantee identical results between the two; the only chance I would see would be to convert the Fortran to a dll and call it from MATLAB, but even then the runtime libraries are likely not the same.

[Maxwell]

Basic question (sorry I'm relatively unfamiliar with formatting in matlab). I have format long set for the program and I put

fprintf("single precision gz2 = %f\n",gz2)

in the command line. The output is single precision gz2 = 24.090454.

In fortran I have

print*, 't4 =', t4

and the output is t4 = 24.09045.

Both of these variables are defined in single precision (I have added the single() command to matlab and defined t4 as REAL in fortran). Do both of these print statements show the full value of each variable? If so, why does this single precision value in Matlab have an extra value. Thanks!

[dpb]

No neither show a specific number of digits/precision, either in Fortran or MATLAB.

The asterisk in Fortran PRINT statement is "compiler-dependent" as to how many digits it will print; the "*" is a shorthand for "I don't care about exact precision enough to write a FORMAT statement here, just give me a value".

Similarly, the '%f' format in MATLAB leaves it to the i/o library to determine the number of printed decimal places; I don't recall otomh if the C Standard from which MATLAB i/o is derived defines a specific number of digits or not.

To present the results with identical displayed digits you need something like

write(*,'F15.8') t4

in Fortran

and

fprintf("single precision gz2 = %15.8f\n",gz2)

in MATLAB.

"I have added the single() command to matlab ..."

NOTA BENE: In MATLAB simply defining a variable as single precision does NOT imply that all computations on that variable are in single precision -- only if every element in the operations all the way through are single will operations be performed only in single precision. In your example above of

gz = 1.0 / (1+0.0247*exp(0.650*log(0.1089*re)))

the constants will be converted by the interpreter as double and then the expression will be promoted to double for operations. If put single() around the expression, then you are just casting the result to single precision. If you have previously cast re to single, then it will have been converted from double to single there, but then promoted back to double for the calculation.

In the end, as noted, it's essentially certain you'll never be able to reproduce the two identically and you would find other and probably different results if you moved the Fortran itself to another system or even recompiled it with a new/different compiler.

We had to go through such verification processes at every change of OS and/or major compiler upgrade with the NRC for reactor safety computation codes; even the NRC recognized that one could only produce so much in absolute comparisons; the test suites compared results to be within what were determined a priori to be "computationally equivalent" values for various quantities. Some similar approach here would seem to be the way to proceed.

[Maxwell]

Thank you for your clarification on this. This makes a lot more sense now

[Walter Roberson]

Further to what @dpb was saying: to calculate entirely in single precision, you would use

gz = single(1.0) / (single(1)+single(0.0247)*exp(single(0.650)*log(single(0.1089)*re)))

[dpb]

Indeed, Walter. MATLAB has no equivalent way to set global precision in a code module as does Fortran and other typed language; all numerical quantities are assumed to be class double unless defined otherwise; for quite some time there was no other floating point class at all and actual support for single precision computations is actually a fairly recent innovation.

Where in Fortran 1.0 is single precision (unless one uses some vendor-specific compiler extension) and one must write 1.0D0 to ensure double precision, MATLAB doesn't have an equivalent to the REAL or INTEGER statements in Fortran; that was an on-purpose decision in the beginning as the teaching toolset in that it used the double on 32-bit machines by default and treated everything numeric as double so students didn't have to deal with such programming details but could concentrate directly on solving the problem at hand. Hence, adding the other classes is somewhat of an add-on; you can store variables as single to save the memory, but unless one goes to the extremes as @Walter Roberson shows above, you probably will lose that by creating a double array for working with the array and the extra overhead of the cast to/from besides. The pain in writing as the above is such that one almost never actually codes things fully in anything except the default class in MATLAB.

[dpb]

C1=1; C2=1.0; C3=1E0; C4=1D0;
whos C*

  Name      Size            Bytes  Class     Attributes

  C1        1x1                 8  double              
  C2        1x1                 8  double              
  C3        1x1                 8  double              
  C4        1x1                 8  double              

isequal(C1,C2,C3,C4)

ans = logical
   1

Shows that all the above from an integer to float and explicit single or double exponential form are all converted into doubles by MATLAB.

[Steven Lord]

To address the isequal call in dpb's post, FYI the isequal function doesn't care about data types or complexity.

[isequal(1, single(1)), isequal(1, int32(1)), isequal(1, complex(1, 0))]

ans = 1x3 logical array
   1   1   1

Now if you were using the unit testing framework to test that the variables were equal, those methods are stricter. They do require the data types to match and require complexity to match (unless you specify a tolerance.) [See the Comparison Details entry in the More About section on the documentation page for the constraint used by the verifyEqual method if you're interested in more details.]

test = matlab.unittest.TestCase.forInteractiveUse;
verifyEqual(test, 1, single(1))

Verification failed.
    ---------------------
    Framework Diagnostic:
    ---------------------
    verifyEqual failed.
    --> Classes do not match.
        
        Actual Class:
            double
        Expected Class:
            single
    
    Actual Value:
         1
    Expected Value:
      single
    
         1
    ------------------
    Stack Information:
    ------------------
    In /tmp/Editor_gfvlu/LiveEditorEvaluationHelperEeditorId.m (LiveEditorEvaluationHelperEeditorId) at 3
    In /MATLAB/toolbox/matlab/connector2/interpreter/+connector/+internal/fevalMatlab.p (fevalMatlab) at 0
    In /MATLAB/toolbox/matlab/connector2/interpreter/+connector/+internal/fevalJSON.p (fevalJSON) at 0

[dpb]

" ... the isequal function doesn't care about data types or complexity."

That's why the whos to show that MATLAB created all the Cn as double...there really wasn't too much point in adding the isequal there, granted, it would take a pretty big mistake in the interpreter that they wouldn't all match. :)

But, it is a good point that to prove the type/class needs more strict testing...one could also use

C1=1; C2=1.0; C3=1E0; C4=1D0;
isequal(class(C1),class(C2),class(C3),class(C4)) & isequal(C1,C2,C3,C4)

ans = logical
   1

to mimic it if, like me, one isn't familiar with the unit test framework... :)

isequal(class(C1),class(single(C2)),class(C3),class(C4)) & isequal(C1,C2,C3,C4)

ans = logical
   0

ADDENDUM:

It is, of course, quite cumbersome to have to code each variable individually (which is why I used the whos first since it accepts the asterisk wildcard to get all in one swell foop), so a crude functional form could look something like

function flg=isidentical(varargin)
% returns T if all input arguments match in class and value, F otherwise
  assert(nargin>1,'Must have at least two arguments')
  flg=true;
  for i=2:nargin
    flg=flg & isequal(class(varargin{1}),class(varargin{i})) & isequal(varargin{1},varargin{i});
    if ~flg, return, end        % proven false, no need to continue
  end
end

isidentical(C1,C2,C3,C4)

ans = logical
   1

isidentical(C1,single(C2),C3,C4)

ans = logical
   0

isidentical(C1/2,C2,C3,C4)

ans = logical
   0

[James Tursa]

Another issue is the implementations of exp( ) and log( ). Languages typically defer these calculations to libraries, and the associated language standards do not impose specific alglorithms that must be used for these library functions. These library functions are not held to the same strict results standards as the +, -, *, and / operations are in IEEE. So one cannot expect exact consistency between different languages, or even different versions or compilers of the same language, for these library functions. You might get lucky and have the specific versions of the libraries you are using match up, but it would be a fragile test setup that could easily break if anything got updated in the future. It would be a poor test design to rely on this.

[Walter Roberson]

According to https://www.mathworks.com/help/hdlcoder/ug/ulp-considerations-of-native-floating-point-operators.html both log() and exp() are constrained to be no more than 1 ULP (unit in the last place) different than the true mathematical answer.

[James Tursa]

@Walter Roberson And that leaves a bit of wiggle room in the result. E.g.,

p = randi([-5,5],10,10);
x = sum(2.^p); % some random values that can be represented exactly in double and vpa
X = string(x');
IEEE_hex = string(num2hex(x));
table(X,IEEE_hex) % demonstrate that these are "nice" in binary floating point format

ans = 10x2 table
        X              IEEE_hex     
    __________    __________________

    "120.3125"    "405e140000000000"
    "10.875"      "4025c00000000000"
    "26.34375"    "403a580000000000"
    "34.03125"    "4041040000000000"
    "108.125"     "405b080000000000"
    "56.28125"    "404c240000000000"
    "28.8125"     "403cd00000000000"
    "23.03125"    "4037080000000000"
    "37.25"       "4042a00000000000"
    "31.25"       "403f400000000000"

digits 50
double(exp(vpa(x))) - exp(x)

ans = 1×10

1.0e-03 * 0 0 0 0 0 0 -0.4883 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

It wasn't too hard to find a value x for which exp(x) did not give a result that was correctly rounded "as if the calculation were performed in infinite precision". This is unlike the IEEE requirement for the +, -, *, and / operations which do have this strict requirement.

exp_vpa_x = string(num2hex(double(exp(vpa(x)))));
exp_x = string(num2hex(exp(x)));
table(exp_vpa_x,exp_x)

ans = 10x2 table
        exp_vpa_x               exp_x       
    __________________    __________________

    "4ac7d28912b900dc"    "4ac7d28912b900dc"
    "40e9ccd7d3d55bc5"    "40e9ccd7d3d55bc5"
    "42501110262b98ea"    "42501110262b98ea"
    "43011c005aae8303"    "43011c005aae8303"
    "49afcf51cc64d9fe"    "49afcf51cc64d9fe"
    "4502564116ff8878"    "4502564116ff8878"
    "4287b6b72fba0d83"    "4287b6b72fba0d84"
    "4202ba2f9c4855b1"    "4202ba2f9c4855b1"
    "434abae41f425f94"    "434abae41f425f94"
    "42c0f63a92073371"    "42c0f63a92073371"

The 7th value above differs by 1 ULP. A different library using a different algorithm meeting the same 1 ULP standard as MATLAB could get a different result. Hence the advice to not rely on exact comparisons for library functions.

[dpb]

+1 @James Tursa for elegant exposition...