MATLAB and Fortran gives different results (precision)

I have 2 small programs in MATLAB and Fortran, respectively, which should be identical

MATLAB version

f=zeros(size(r1,1),3);
for ii=1:size(r1,1)
    for iii=1:size(r2,1)
        if all(r1(ii,:)==r2(iii,:))
            continue
        end
        f(ii,:)=f(ii,:)+ScaFun(r1(ii,:),e1(ii),r2(iii,:),e2(iii));
    end
end

function f=ScaFun(r1,e1,r2,e2)
r21 = r1-r2;
R = sqrt(r21(1)^2+r21(2)^2+r21(3)^2);

f = e1.*e2.*r21./(R^3);
end

Fortran:

function InvSq(r1,e1,r2,e2) result(f)
	
	real*8, allocatable, intent(in) :: r1(:,:), e1(:), r2(:,:), e2(:)
	integer*8 :: n1, n2, ii, jj
	real*8, allocatable :: f(:,:)
	
	n1 = size(r1,1)
	n2 = size(r2,1)
	allocate(f(n1,3),source = 0.0_8)

	do ii = 1,n1
		jj_loop: do jj = 1,n2
			if (all(r1(ii,:)==r2(jj,:))) cycle jj_loop !skip self-to-self interaction (which InvSq_scalar() will return NaN)
			f(ii,:) = f(ii,:)+InvSq_scalar(r1(ii,:),e1(ii),r2(jj,:),e2(jj))
		end do jj_loop
	end do
	
end function InvSq


function InvSq_scalar(r1,e1,r2,e2) result(f21)
	
	real*8, intent(in) :: r1(3), e1, r2(3), e2
	real*8 :: r21(3), R
	real*8 :: f21(3)
	
	r21 = r1-r2
	R = sqrt(r21(1)**2+r21(2)**2+r21(3)**2);
	
	f21 = e1*e2*r21/(R**3)
	
end function InvSq_scalar

However, for randomly generated inputs (e.g., sizein=1e3)

r1=2*rand(sizein,3)-1;
e1=2*rand(sizein,1)-1;
r2=2*rand(sizein,3)-1;
e2=2*rand(sizein,1)-1;

The (absolute) difference between MATLAB and Fortran results are on the order of 1e-13 for most, with a few as high as 1e-12.

Now, I understand that some of the inputs may be badly scaled (i.e., when r1 very close to r2) and this could be problematic for floating point arithmetic. However, they should still use the same ieee doubles and excute the same operations in the same order and should give identical results. Is there something else under the hood makes this happening?

Answer 1

[Walter Roberson]

FORTRAN and MATLAB do not necessarily use the same floating point rounding mode for calculations.

Control over rounding mode is not standardized in FORTRAN. See for example IBM https://www.ibm.com/docs/en/xcafbg/9.0.0?topic=SS3KZ4_9.0.0/com.ibm.xlf111.bg.doc/xlfopg/fpround.html versus Intel https://www.intel.com/content/www/us/en/develop/documentation/fortran-compiler-oneapi-dev-guide-and-reference/top/language-reference/a-to-z-reference/h-to-i/ieee-set-rounding-mode.html

MATLAB has an undocumented function to set rounding behaviour (for code within MATLAB; might not apply to the high performance libraries that are automatically called.) See https://stackoverflow.com/questions/55682034/how-to-change-the-rounding-mode-for-floating-point-operations-in-matlab

[Paul]

At least on the FORTRAN side, might it also be the case that compiler optimizations result in floating point computations executed in a different order than in Matlab, which would also contribute to small differences in results?

[Walter Roberson]

Yes. The default for Fortran compilers should be that the calculations are performed in strict precedence order in the order written, left to right unless the documented precedence says otherwise. However there are typically command line options provided that permit operations to be rearranged. Some of the options possibilities "should" only rearrange if it will not affect rounding, but advanced options often permit reordering that can affect rounding.

MATLAB defines that operations should be executed in precedence order otherwise left to right as written. But... it is known that some code patterns are intercepted and replaced with calls to high performance libraries that might have different internal orders.

[Yi-xiao Liu]

@Walter Roberson

interesting. So this is also the reason when I try to vectorize the code in MATLAB it also produces differences ~1e-12 to 1e-13 versus the loop version?

BTW do you know the default rounding mode (or how to set it) in gfortran? This seems to be the only fortran compiler supported by mex and a quick google dose not yield something similar like the intel one.

[Walter Roberson]

https://stackoverflow.com/questions/53950936/gfortran-linker-argument-for-library-containing-ieee-get-rounding-mode has some code for gfortran

[Walter Roberson]

In MATLAB, sum() over a larger array might not give the same result as looping adding one value at a time. sum() for sufficiently large arrays is passed to the high performance math libraries. The values to be totalled are partitioned, and one core is allocated per partion, and then the partial results are added. This will typically give different rounding results.

[Yi-xiao Liu]

@Walter Roberson

Given that it is unclear how ties are treated when set to IEEE_NEAREST in MATLAB and FORTRAN, I set both to IEEE_TO_ZERO and rerun the tests (complier optimization is also disabled with -O0). However there are still ~1e-12 differences MATLAB vs FORTRAN.

Does MATLAB apply optimizations that could affect the order of summations? I heard a lot was done to improve the for loop performance in recent releases.