MATLAB and Fortran gives different results (precision)
7 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Walter Roberson
2022-12-26
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 个评论
Walter Roberson
2022-12-27
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Fortran with MATLAB 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!