MATLAB and Fortran gives different results (precision)
12 views (last 30 days)
Show older comments
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?
0 Comments
Answers (1)
Walter Roberson
on 26 Dec 2022
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
6 Comments
Walter Roberson
on 27 Dec 2022
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.
See Also
Categories
Find more on Fortran with MATLAB in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!